Homec4science
Diffusion lammps (efficient_neuronet) Browse

lammps/src/MANYBODY/pair_comb3.cpp
f0cd77adefcdefficient_neuronet

Last Change
Raw File

pair_comb3.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Ray Shan (Sandia, tnshan@sandia.gov)
Updates and debug: Tao Liang (U Florida, liang75@ufl.edu)
Dundar Yilmaz (dundar.yilmaz@zirve.edu.tr)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_comb3.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "group.h"
#include "update.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
#define MAXNEIGH 24
/* ---------------------------------------------------------------------- */
PairComb3::PairComb3(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
ghostneigh = 1;
nmax = 0;
NCo = NULL;
bbij = NULL;
map = NULL;
esm = NULL;
nelements = 0;
elements = NULL;
nparams = 0;
maxparam = 0;
params = NULL;
elem2param = NULL;
intype = NULL;
afb = NULL;
dafb = NULL;
fafb = NULL;
dfafb = NULL;
ddfafb = NULL;
phin = NULL;
dphin = NULL;
erpaw = NULL;
vvdw = NULL;
vdvdw = NULL;
dpl = NULL;
xcctmp = NULL;
xchtmp = NULL;
xcotmp = NULL;
sht_num = NULL;
sht_first = NULL;
ipage = NULL;
pgsize = oneatom = 0;
cflag = 0;
// set comm size needed by this Pair
comm_forward = 1;
comm_reverse = 1;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairComb3::~PairComb3()
{
memory->destroy(NCo);
if (elements)
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
memory->sfree(params);
memory->destroy(elem2param);
memory->destroy(afb);
memory->destroy(dpl);
memory->destroy(dafb);
memory->destroy(fafb);
memory->destroy(phin);
memory->destroy(bbij);
memory->destroy(vvdw);
memory->destroy(vdvdw);
memory->destroy(dphin);
memory->destroy(erpaw);
memory->destroy(dfafb);
memory->destroy(ddfafb);
memory->destroy(xcctmp);
memory->destroy(xchtmp);
memory->destroy(xcotmp);
memory->destroy(intype);
memory->destroy(sht_num);
memory->sfree(sht_first);
delete [] ipage;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cutghost);
delete [] map;
delete [] esm;
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutghost,n+1,n+1,"pair:cutghost");
map = new int[n+1];
esm = new double[n];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairComb3::settings(int narg, char **arg)
{
if (strcmp(arg[0],"polar_on") == 0) {
pol_flag = 1;
if (comm->me == 0) fprintf(screen,
" PairComb3: polarization is on \n");
} else if (strcmp(arg[0],"polar_off") == 0) {
pol_flag = 0;
if (comm->me == 0) fprintf(screen,
" PairComb3: polarization is off \n");
} else {
error->all(FLERR,"Illegal pair_style command");
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairComb3::coeff(int narg, char **arg)
{
int i,j,n;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
// nelements = # of unique elements
// elements = list of element names
if (elements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
}
elements = new char*[atom->ntypes];
for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
nelements = 0;
for (i = 3; i < narg; i++) {
if ((strcmp(arg[i],"C") == 0) && (cflag == 0)) {
if( comm->me == 0) fprintf(screen,
" PairComb3: Found C: reading additional library file\n");
read_lib();
cflag = 1;
}
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
map[i-2] = j;
if (j == nelements) {
n = strlen(arg[i]) + 1;
elements[j] = new char[n];
strcpy(elements[j],arg[i]);
nelements++;
}
}
// read potential file and initialize potential parameters
read_file(arg[2]);
setup();
n = atom->ntypes;
// generate Wolf 1/r energy and van der Waals look-up tables
tables();
// clear setflag since coeff() called once with I,J = * *
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairComb3::init_style()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Pair style COMB3 requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style COMB3 requires newton pair on");
if (!atom->q_flag)
error->all(FLERR,"Pair style COMB3 requires atom attribute q");
// need a full neighbor list
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->ghost = 1;
// local Comb neighbor list
// create pages if first time or if neighbor pgsize/oneatom has changed
int create = 0;
if (ipage == NULL) create = 1;
if (pgsize != neighbor->pgsize) create = 1;
if (oneatom != neighbor->oneatom) create = 1;
if (create) {
delete [] ipage;
pgsize = neighbor->pgsize;
oneatom = neighbor->oneatom;
int nmypage = comm->nthreads;
ipage = new MyPage<int>[nmypage];
for (int i = 0; i < nmypage; i++)
ipage[i].init(oneatom,pgsize);
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairComb3::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
cutghost[j][i] = cutghost[i][j] = cutmax;
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairComb3::read_lib()
{
unsigned int maxlib = 1024;
int i,j,k,l,nwords,m;
int ii,jj,kk,ll,mm,iii;
char s[maxlib];
char **words = new char*[80];
MPI_Comm_rank(world,&comm->me);
// open libraray file on proc 0
if (comm->me == 0) {
FILE *fp = force->open_potential("lib.comb3");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open COMB3 lib.comb3 file");
error->one(FLERR,str);
}
// read and store at the same time
fgets(s,maxlib,fp);
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ccutoff[0] = atof(words[0]);
ccutoff[1] = atof(words[1]);
ccutoff[2] = atof(words[2]);
ccutoff[3] = atof(words[3]);
ccutoff[4] = atof(words[4]);
ccutoff[5] = atof(words[5]);
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ch_a[0] = atof(words[0]);
ch_a[1] = atof(words[1]);
ch_a[2] = atof(words[2]);
ch_a[3] = atof(words[3]);
ch_a[4] = atof(words[4]);
ch_a[5] = atof(words[5]);
ch_a[6] = atof(words[6]);
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
nsplpcn = atoi(words[0]);
nsplrad = atoi(words[1]);
nspltor = atoi(words[2]);
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
maxx = atoi(words[0]);
maxy = atoi(words[1]);
maxz = atoi(words[2]);
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
maxxc = atoi(words[0]);
maxyc = atoi(words[1]);
maxconj = atoi(words[2]);
for (l=0; l<nsplpcn; l++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
maxxcn[l] = atoi(words[1]);
vmaxxcn[l] = atof(words[2]);
dvmaxxcn[l] = atof(words[3]);
}
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ntab = atoi(words[0]);
for (i=0; i<ntab+1; i++){
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
pang[i] = atof(words[1]);
dpang[i] = atof(words[2]);
ddpang[i] = atof(words[3]);
}
for (l=0; l<nsplpcn; l++)
for (i=0; i<maxx+1; i++)
for (j=0; j<maxy+1; j++)
for (k=0; k<maxz+1; k++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ll = atoi(words[0])-1;
ii = atoi(words[1]);
jj = atoi(words[2]);
kk = atoi(words[3]);
pcn_grid[ll][ii][jj][kk] = atof(words[4]);
pcn_gridx[ll][ii][jj][kk] = atof(words[5]);
pcn_gridy[ll][ii][jj][kk] = atof(words[6]);
pcn_gridz[ll][ii][jj][kk] = atof(words[7]);
}
for (l=0; l<nsplpcn; l++)
for (i=0; i<maxx; i++)
for (j=0; j<maxy; j++)
for (k=0; k<maxz; k++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ll = atoi(words[0])-1;
ii = atoi(words[1]);
jj = atoi(words[2]);
kk = atoi(words[3]);
for(iii=0; iii<2; iii++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
for(m=0; m<32 ; m++) {
mm=iii*32+m;
pcn_cubs[ll][ii][jj][kk][mm] = atof(words[m]);
}
}
}
for (l=0; l<nsplrad; l++)
for (i=0; i<maxxc+1; i++)
for (j=0; j<maxyc+1; j++)
for (k=0; k<maxconj; k++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ll = atoi(words[0])-1;
ii = atoi(words[1]);
jj = atoi(words[2]);
kk = atoi(words[3])-1;
rad_grid[ll][ii][jj][kk] = atof(words[4]);
rad_gridx[ll][ii][jj][kk] = atof(words[5]);
rad_gridy[ll][ii][jj][kk] = atof(words[6]);
rad_gridz[ll][ii][jj][kk] = atof(words[7]);
}
for (l=0; l<nsplrad; l++)
for (i=0; i<maxxc; i++)
for (j=0; j<maxyc; j++)
for (k=0; k<maxconj-1; k++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ll = atoi(words[0])-1;
ii = atoi(words[1]);
jj = atoi(words[2]);
kk = atoi(words[3])-1;
for (iii=0; iii<2; iii++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
for(m=0; m<32 ; m++){
mm=iii*32+m;
rad_spl[ll][ii][jj][kk][mm] = atof(words[m]);
}
}
}
for (l=0; l<nspltor; l++)
for (i=0; i<maxxc+1; i++)
for (j=0; j<maxyc+1; j++)
for (k=0; k<maxconj; k++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ll = atoi(words[0])-1;
ii = atoi(words[1]);
jj = atoi(words[2]);
kk = atoi(words[3])-1;
tor_grid[ll][ii][jj][kk] = atof(words[4]);
tor_gridx[ll][ii][jj][kk] = atof(words[5]);
tor_gridy[ll][ii][jj][kk] = atof(words[6]);
tor_gridz[ll][ii][jj][kk] = atof(words[7]);
}
for (l=0; l<nspltor; l++)
for (i=0; i<maxxc; i++)
for (j=0; j<maxyc; j++)
for (k=0; k<maxconj-1; k++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
ll = atoi(words[0])-1;
ii = atoi(words[1]);
jj = atoi(words[2]);
kk = atoi(words[3])-1;
for(iii=0; iii<2; iii++) {
fgets(s,maxlib,fp);
nwords = 0;
words[nwords++] = strtok(s," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
for (m=0; m<32 ; m++){
mm=iii*32+m;
tor_spl[ll][ii][jj][kk][mm] = atof(words[m]);
}
}
}
fclose(fp);
}
k = 0;
for (i=0; i<4; i++)
for (j=0; j<4; j++) {
iin2[k][0] = i;
iin2[k][1] = j;
k ++;
}
l = 0;
for (i=0; i<4; i++)
for (j=0; j<4; j++)
for (k=0; k<4; k++) {
iin3[l][0] = i;
iin3[l][1] = j;
iin3[l][2] = k;
l ++;
}
MPI_Bcast(&ccutoff[0],6,MPI_DOUBLE,0,world);
MPI_Bcast(&ch_a[0],7,MPI_DOUBLE,0,world);
MPI_Bcast(&nsplpcn,1,MPI_INT,0,world);
MPI_Bcast(&nsplrad,1,MPI_INT,0,world);
MPI_Bcast(&nspltor,1,MPI_INT,0,world);
MPI_Bcast(&maxx,1,MPI_INT,0,world);
MPI_Bcast(&maxy,1,MPI_INT,0,world);
MPI_Bcast(&maxz,1,MPI_INT,0,world);
MPI_Bcast(&maxxc,1,MPI_INT,0,world);
MPI_Bcast(&maxyc,1,MPI_INT,0,world);
MPI_Bcast(&maxconj,1,MPI_INT,0,world);
MPI_Bcast(&maxxcn,4,MPI_INT,0,world);
MPI_Bcast(&vmaxxcn,4,MPI_DOUBLE,0,world);
MPI_Bcast(&dvmaxxcn,4,MPI_DOUBLE,0,world);
MPI_Bcast(&ntab,1,MPI_INT,0,world);
MPI_Bcast(&pang[0],20001,MPI_DOUBLE,0,world);
MPI_Bcast(&dpang[0],20001,MPI_DOUBLE,0,world);
MPI_Bcast(&ddpang[0],20001,MPI_DOUBLE,0,world);
MPI_Bcast(&pcn_grid[0][0][0][0],500,MPI_DOUBLE,0,world);
MPI_Bcast(&pcn_gridx[0][0][0][0],500,MPI_DOUBLE,0,world);
MPI_Bcast(&pcn_gridy[0][0][0][0],500,MPI_DOUBLE,0,world);
MPI_Bcast(&pcn_gridz[0][0][0][0],500,MPI_DOUBLE,0,world);
MPI_Bcast(&pcn_cubs[0][0][0][0][0],16384,MPI_DOUBLE,0,world);
MPI_Bcast(&rad_grid[0][0][0][0],825,MPI_DOUBLE,0,world);
MPI_Bcast(&rad_gridx[0][0][0][0],825,MPI_DOUBLE,0,world);
MPI_Bcast(&rad_gridy[0][0][0][0],825,MPI_DOUBLE,0,world);
MPI_Bcast(&rad_gridz[0][0][0][0],825,MPI_DOUBLE,0,world);
MPI_Bcast(&rad_spl[0][0][0][0][0],30720,MPI_DOUBLE,0,world);
MPI_Bcast(&tor_grid[0][0][0][0],275,MPI_DOUBLE,0,world);
MPI_Bcast(&tor_gridx[0][0][0][0],275,MPI_DOUBLE,0,world);
MPI_Bcast(&tor_gridy[0][0][0][0],275,MPI_DOUBLE,0,world);
MPI_Bcast(&tor_gridz[0][0][0][0],275,MPI_DOUBLE,0,world);
MPI_Bcast(&tor_spl[0][0][0][0][0],10240,MPI_DOUBLE,0,world);
MPI_Bcast(&iin2[0][0],32,MPI_INT,0,world);
MPI_Bcast(&iin3[0][0],192,MPI_INT,0,world);
delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairComb3::read_file(char *file)
{
int params_per_line = 74;
char **words = new char*[params_per_line+1];
if (params) delete [] params;
params = NULL;
nparams = 0;
// open file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = fopen(file,"r");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open COMB3 potential file %s",file);
error->one(FLERR,str);
}
}
// read each line out of file, skipping blank lines or leading '#'
// store line of params if all 3 element tags are in element list
int n,nwords,ielement,jelement,kelement;
char line[MAXLINE],*ptr;
int eof = 0;
nwords=0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != params_per_line){
error->all(FLERR,"Incorrect format in COMB3 potential file");
}
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((nwords <= params_per_line)
&& (words[nwords++] = strtok(NULL," \t\n\r\f"))) {
continue;
}
// ielement,jelement,kelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next line
for (ielement = 0; ielement < nelements; ielement++)
if (strcmp(words[0],elements[ielement]) == 0) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (strcmp(words[1],elements[jelement]) == 0) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (strcmp(words[2],elements[kelement]) == 0) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].ielementgp = atoi(words[3]);
params[nparams].jelementgp = atoi(words[4]);
params[nparams].kelementgp = atoi(words[5]);
params[nparams].ang_flag = atoi(words[6]);
params[nparams].pcn_flag = atoi(words[7]);
params[nparams].rad_flag = atoi(words[8]);
params[nparams].tor_flag = atoi(words[9]);
params[nparams].vdwflag = atof(words[10]);
params[nparams].powerm = atof(words[11]);
params[nparams].veps = atof(words[12]);
params[nparams].vsig = atof(words[13]);
params[nparams].paaa = atof(words[14]);
params[nparams].pbbb = atof(words[15]);
params[nparams].lami = atof(words[16]);
params[nparams].alfi = atof(words[17]);
params[nparams].powern = atof(words[18]);
params[nparams].QL = atof(words[19]);
params[nparams].QU = atof(words[20]);
params[nparams].DL = atof(words[21]);
params[nparams].DU = atof(words[22]);
params[nparams].qmin = atof(words[23]);
params[nparams].qmax = atof(words[24]);
params[nparams].chi = atof(words[25]);
params[nparams].dj = atof(words[26]);
params[nparams].dk = atof(words[27]);
params[nparams].dl = atof(words[28]);
params[nparams].esm = atof(words[29]);
params[nparams].cmn1 = atof(words[30]);
params[nparams].cmn2 = atof(words[31]);
params[nparams].pcmn1 = atof(words[32]);
params[nparams].pcmn2 = atof(words[33]);
params[nparams].coulcut = atof(words[34]);
params[nparams].polz = atof(words[35]);
params[nparams].curl = atof(words[36]);
params[nparams].curlcut1 = atof(words[37]);
params[nparams].curlcut2 = atof(words[38]);
params[nparams].curl0 = atof(words[39]);
params[nparams].alpha1 = atof(words[40]);
params[nparams].bigB1 = atof(words[41]);
params[nparams].alpha2 = atof(words[42]);
params[nparams].bigB2 = atof(words[43]);
params[nparams].alpha3 = atof(words[44]);
params[nparams].bigB3 = atof(words[45]);
params[nparams].lambda = atof(words[46]);
params[nparams].bigA = atof(words[47]);
params[nparams].beta = atof(words[48]);
params[nparams].bigr = atof(words[49]);
params[nparams].bigd = atof(words[50]);
params[nparams].pcos6 = atof(words[51]);
params[nparams].pcos5 = atof(words[52]);
params[nparams].pcos4 = atof(words[53]);
params[nparams].pcos3 = atof(words[54]);
params[nparams].pcos2 = atof(words[55]);
params[nparams].pcos1 = atof(words[56]);
params[nparams].pcos0 = atof(words[57]);
params[nparams].pcna = atof(words[58]);
params[nparams].pcnb = atof(words[59]);
params[nparams].pcnc = atof(words[60]);
params[nparams].pcnd = atof(words[61]);
params[nparams].p6p0 = atof(words[62]);
params[nparams].p6p1 = atof(words[63]);
params[nparams].p6p2 = atof(words[64]);
params[nparams].p6p3 = atof(words[65]);
params[nparams].p6p4 = atof(words[66]);
params[nparams].p6p5 = atof(words[67]);
params[nparams].p6p6 = atof(words[68]);
params[nparams].ptork1=atof(words[69]);
params[nparams].ptork2=atof(words[70]);
params[nparams].addrepr=atof(words[71]);
params[nparams].addrep=atof(words[72]);
params[nparams].pcross = atof(words[73]);
params[nparams].powermint = int(params[nparams].powerm);
// parameter sanity checks
if (params[nparams].lambda < 0.0 || params[nparams].powern < 0.0 ||
params[nparams].beta < 0.0 || params[nparams].alpha1 < 0.0 ||
params[nparams].bigB1< 0.0 || params[nparams].bigA< 0.0 ||
params[nparams].bigB2< 0.0 || params[nparams].alpha2 <0.0 ||
params[nparams].bigB3< 0.0 || params[nparams].alpha3 <0.0 ||
params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
params[nparams].bigd > params[nparams].bigr ||
params[nparams].powerm - params[nparams].powermint != 0.0 ||
params[nparams].addrepr < 0.0 || params[nparams].powermint < 1.0 ||
params[nparams].QL > 0.0 || params[nparams].QU < 0.0 ||
params[nparams].DL < 0.0 || params[nparams].DU > 0.0 ||
params[nparams].pcross < 0.0 ||
params[nparams].esm < 0.0 || params[nparams].veps < 0.0 ||
params[nparams].vsig < 0.0 || params[nparams].vdwflag < 0.0
)
error->all(FLERR,"Illegal COMB3 parameter");
nparams++;
}
delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairComb3::setup()
{
int i,j,k,m,n;
// set elem2param for all element triplet combinations
// must be a single exact match to lines read from file
// do not allow for ACB in place of ABC
memory->destroy(elem2param);
memory->create(elem2param,nelements,nelements,nelements,"pair:elem2param");
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
n = -1;
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement &&
k == params[m].kelement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
elem2param[i][j][k] = n;
}
// compute parameter values derived from inputs
for (m = 0; m < nparams; m++) {
params[m].cut = params[m].bigr + params[m].bigd;
params[m].cutsq = params[m].cut*params[m].cut;
params[m].c1 = pow(2.0*params[m].powern*1.0e-16,-1.0/params[m].powern);
params[m].c2 = pow(2.0*params[m].powern*1.00e-8,-1.0/params[m].powern);
params[m].c3 = 1.0/params[m].c2;
params[m].c4 = 1.0/params[m].c1;
params[m].Qo = (params[m].QU+params[m].QL)/2.0; // (A22)
params[m].dQ = (params[m].QU-params[m].QL)/2.0; // (A21)
params[m].aB = 1.0 /
(1.0-pow(fabs(params[m].Qo/params[m].dQ),10)); // (A20)
params[m].bB = pow(fabs(params[m].aB),0.1)/params[m].dQ; // (A19)
params[m].nD = log(params[m].DU/(params[m].DU-params[m].DL))/
log(params[m].QU/(params[m].QU-params[m].QL));
params[m].bD = (pow((params[m].DL-params[m].DU),(1.0/params[m].nD)))/
(params[m].QU-params[m].QL);
params[m].lcut = params[m].coulcut;
params[m].lcutsq = params[m].lcut*params[m].lcut;
}
// set cutmax to max of all params
cutmin = cutmax = 0.0;
polar = 0;
for (m = 0; m < nparams; m++) {
if (params[m].cutsq > cutmin) cutmin = params[m].cutsq + 2.0;
if (params[m].lcut > cutmax) cutmax = params[m].lcut;
}
chicut1 = 7.0;
chicut2 = cutmax;
}
/* ---------------------------------------------------------------------- */
void PairComb3::Short_neigh()
{
int nj,*neighptrj,icontrol;
int iparam_ij,*ilist,*jlist,*numneigh,**firstneigh;
int inum,jnum,i,j,ii,jj,itype,jtype;
double rr1,rsq1,delrj[3];
double **x = atom->x;
int *type = atom->type;
if (atom->nmax > nmax) {
memory->sfree(sht_first);
nmax = atom->nmax;
sht_first = (int **) memory->smalloc(nmax*sizeof(int *),
"pair:sht_first");
memory->grow(dpl,nmax,3,"pair:dpl");
memory->grow(xcctmp,nmax,"pair:xcctmp");
memory->grow(xchtmp,nmax,"pair:xchtmp");
memory->grow(xcotmp,nmax,"pair:xcotmp");
memory->grow(NCo,nmax,"pair:NCo");
memory->grow(sht_num,nmax,"pair:sht_num");
memory->grow(bbij,nmax,MAXNEIGH,"pair:bbij");
}
inum = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// create COMB neighbor list
ipage->reset();
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
dpl[i][0] = dpl[i][1] = dpl[i][2] = 0.0;
nj = 0;
neighptrj = ipage->vget();
itype = map[type[i]];
jlist = firstneigh[i];
jnum = numneigh[i];
NCo[i] = 0.0;
xcctmp[i] = 0.0;
xchtmp[i] = 0.0;
xcotmp[i] = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj] & NEIGHMASK;
delrj[0] = x[i][0] - x[j][0];
delrj[1] = x[i][1] - x[j][1];
delrj[2] = x[i][2] - x[j][2];
rsq1 = vec3_dot(delrj,delrj);
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype][jtype];
if (rsq1 > cutmin) continue;
neighptrj[nj++] = j;
rr1 = sqrt(rsq1);
NCo[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;
icontrol = params[iparam_ij].jelementgp;
if( icontrol == 1)
xcctmp[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;
if (icontrol == 2)
xchtmp[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;
if (icontrol == 3)
xcotmp[i] += comb_fc(rr1,&params[iparam_ij]) * params[iparam_ij].pcross;
}
sht_first[i] = neighptrj;
sht_num[i] = nj;
ipage->vgot(nj);
if (ipage->status())
error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
// communicating coordination number to all nodes
pack_flag = 2;
comm->forward_comm_pair(this);
}
/* ---------------------------------------------------------------------- */
void PairComb3::compute(int eflag, int vflag)
{
int i,ii,k,kk,j,jj,im,inum,jnum,itype,jtype,ktype;
int iparam_i,iparam_ij,iparam_ji;
int iparam_ijk,iparam_jik,iparam_ikj,iparam_jli,iparam_ikl;
int sht_jnum,*sht_jlist,sht_lnum,*sht_llist;
int sht_mnum,*sht_mlist,sht_pnum,*sht_plist;
int *ilist,*jlist,*numneigh,**firstneigh,mr1,mr2,mr3,inty,nj;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rsq1,rsq2,rsq3,iq,jq,yaself;
double eng_tmp,vionij,fvionij,sr1,sr2,sr3;
double zeta_ij,prefac_ij1,prefac_ij2,prefac_ij3,prefac_ij4,prefac_ij5;
double zeta_ji,prefac_ji1,prefac_ji2,prefac_ji3,prefac_ji4,prefac_ji5;
double delrj[3],delrk[3],fi[3],fj[3],fk[3],fl[3];
double ep6p_ij,ep6p_ji,fip6p[3],fjp6p[3],fkp6p[3],flp6p[3];
double potal,fac11,fac11e;
tagint itag, jtag;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
tagint *tag = atom->tag;
int *type = atom->type;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
// coordination terms
double xcn, ycn;
double kcn, lcn;
int torindx;
// torsion and radical variables
int l, ll, ltype, m, mm, mtype, p, pp, ptype;
int iparam_jil, iparam_ijl, iparam_ki, iparam_lj;
int iparam_jl, iparam_ik, iparam_km, iparam_lp;
double kconjug, lconjug, kradtot, lradtot;
double delrl[3], delrm[3], delrp[3], ddprx[3], srmu;
double zet_addi,zet_addj;
evdwl = eng_tmp = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = vflag_atom = 0;
// Build short range neighbor list
Short_neigh();
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
yaself = vionij = fvionij = fpair = 0.0;
// self energy correction term: potal
potal_calc(potal,fac11,fac11e);
// generate initial dipole tensor
if (pol_flag )
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
iq = q[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj] & NEIGHMASK;
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < x[i][2]) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
jtype = map[type[j]];
jq = q[j];
delrj[0] = x[j][0] - xtmp;
delrj[1] = x[j][1] - ytmp;
delrj[2] = x[j][2] - ztmp;
rsq = vec3_dot(delrj,delrj);
iparam_ij = elem2param[itype][jtype][jtype];
iparam_ji = elem2param[jtype][itype][itype];
if (rsq > params[iparam_ij].lcutsq) continue;
tri_point(rsq, mr1, mr2, mr3, sr1, sr2, sr3);
dipole_init(&params[iparam_ij],&params[iparam_ji],fac11,delrj,
rsq,mr1,mr2,mr3,sr1,sr2,sr3,iq,jq,i,j);
}
}
// loop over full neighbor list of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
iq = q[i];
nj = 0;
iparam_i = elem2param[itype][itype][itype];
// self energy, only on i atom
yaself = self(&params[iparam_i],iq);
// dipole self energy
if (pol_flag)
yaself += dipole_self(&params[iparam_i],i);
if (evflag) ev_tally(i,i,nlocal,0,0.0,yaself,0.0,0.0,0.0,0.0);
// two-body interactions (long:R + A, short: only R)
jlist = firstneigh[i];
jnum = numneigh[i];
sht_jlist = sht_first[i];
sht_jnum = sht_num[i];
// long range interactions
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj] & NEIGHMASK;
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < x[i][2]) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
jtype = map[type[j]];
jq = q[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
iparam_ij = elem2param[itype][jtype][jtype];
iparam_ji = elem2param[jtype][itype][itype];
if (rsq > params[iparam_ij].lcutsq) continue;
inty = intype[itype][jtype];
// three-point interpolation
tri_point(rsq, mr1, mr2, mr3, sr1, sr2, sr3);
// Q-indenpendent: van der Waals
vdwaals(inty,mr1,mr2,mr3,rsq,sr1,sr2,sr3,eng_tmp,fpair);
evdwl = eng_tmp;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
if (evflag)
ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
// Q-dependent: Coulombic, field, polarization
// 1/r energy and forces
direct(&params[iparam_ij], &params[iparam_ji],
mr1, mr2, mr3, rsq, sr1, sr2, sr3, iq, jq,
fac11, fac11e, eng_tmp, fvionij, i, j);
vionij = eng_tmp;
// field correction to self energy
field(&params[iparam_ij], &params[iparam_ji],rsq,iq,jq,
eng_tmp,fvionij);
vionij += eng_tmp;
// sums up long range Q-dependent forces (excluding dipole)
f[i][0] += delx*fvionij;
f[i][1] += dely*fvionij;
f[i][2] += delz*fvionij;
f[j][0] -= delx*fvionij;
f[j][1] -= dely*fvionij;
f[j][2] -= delz*fvionij;
// sums up long range Q-dependent energies (excluding dipole)
if (evflag)
ev_tally(i,j,nlocal,newton_pair,0.0,vionij,fvionij,delx,dely,delz);
// polarization field
if (pol_flag) {
dipole_calc(&params[iparam_ij], &params[iparam_ji],fac11,
delx,dely,delz,rsq,mr1,mr2,mr3,
sr1,sr2,sr3,iq,jq,i,j,eng_tmp,fvionij,ddprx);
vionij = eng_tmp;
// sums up dipole energies
if (evflag)
ev_tally(i,j,nlocal,newton_pair,0.0,vionij,fvionij,delx,dely,delz);
// sums up dipole forces
f[i][0] += (ddprx[0] + delx*fvionij);
f[i][1] += (ddprx[1] + dely*fvionij);
f[i][2] += (ddprx[2] + delz*fvionij);
f[j][0] -= (ddprx[0] + delx*fvionij);
f[j][1] -= (ddprx[1] + dely*fvionij);
f[j][2] -= (ddprx[2] + delz*fvionij);
}
if (rsq > params[iparam_ij].cutsq) continue;
repulsive(&params[iparam_ij], &params[iparam_ji], rsq,
fpair, eflag, eng_tmp, iq, jq);
evdwl = eng_tmp;
// repulsion is pure two-body, sums up pair repulsive forces
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
if (evflag)
ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
}
// many-body interactions: start of short-range
xcn = NCo[i];
for (jj = 0; jj < sht_jnum; jj++) {
j = sht_jlist[jj];
sht_llist = sht_first[j];
sht_lnum = sht_num[j];
jtag = tag[j];
if ( jtag <= itag ) continue ;
ycn = NCo[j];
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype][jtype];
iparam_ji = elem2param[jtype][itype][itype];
delrj[0] = x[j][0] - xtmp;
delrj[1] = x[j][1] - ytmp;
delrj[2] = x[j][2] - ztmp;
rsq1 = vec3_dot(delrj,delrj);
if (rsq1 > params[iparam_ij].cutsq) continue;
nj ++;
// this Qj for q-dependent BSi
jq = q[j];
// accumulate bondorder zeta for each i-j interaction via k and l loops
zeta_ij = 0.0;
bbtor = 0.0;
kconjug = 0.0;
for (kk = 0; kk < sht_jnum; kk++) { // kk is neighbor of ii
k = sht_jlist[kk];
if (j == k) continue;
ktype = map[type[k]];
iparam_ijk = elem2param[itype][jtype][ktype];
iparam_ikj = elem2param[itype][ktype][jtype];
iparam_jik = elem2param[jtype][itype][ktype];
iparam_ik = elem2param[itype][ktype][ktype];
delrk[0] = x[k][0] - xtmp;
delrk[1] = x[k][1] - ytmp;
delrk[2] = x[k][2] - ztmp;
rsq2 = vec3_dot(delrk,delrk);
if (rsq2 > params[iparam_ik].cutsq) continue;
// 3-body zeta in bond order
zeta_ij += zeta(&params[iparam_ijk], &params[iparam_ik],
rsq1, rsq2, delrj, delrk, i, xcn);
// radical initialization: apply only to CC,CO,OC bonds
if (params[iparam_ij].rad_flag > 0 &&
params[iparam_ik].ielementgp == 1 &&
params[iparam_ik].jelementgp == 1) {
iparam_ki = elem2param[ktype][itype][itype];
kcn=NCo[k];
kconjug += rad_init(rsq2,&params[iparam_ki],i,kradtot,kcn);
}
// torsion: i-j-k-l: apply to all C-C bonds
if( params[iparam_ij].tor_flag != 0 ) {
srmu = vec3_dot(delrj,delrk)/(sqrt(rsq1*rsq2));
srmu = sqrt(1.0-srmu*srmu);
if(srmu > 0.1) {
for (ll = 0; ll < sht_lnum; ll++) { // ll is neighbor of jj
l = sht_llist[ll];
if(l==i || l==j || l==k) continue;
ltype = map[type[l]];
delrl[0] = x[l][0] - x[j][0];
delrl[1] = x[l][1] - x[j][1];
delrl[2] = x[l][2] - x[j][2];
rsq3 = vec3_dot(delrl,delrl);
iparam_jl = elem2param[jtype][ltype][ltype];
if (rsq3 > params[iparam_jl].cutsq) continue;
iparam_ikl = elem2param[itype][ktype][ltype];
torindx = params[iparam_ij].tor_flag;
bbtor += bbtor1(torindx, &params[iparam_ikl],&params[iparam_jl],
rsq1,rsq2,rsq3,delrj,delrk,delrl,srmu);
}
}
}
}
zeta_ji = 0.0;
lconjug = 0.0;
for (ll = 0; ll < sht_lnum; ll++) {
l = sht_llist[ll];
if (l == i) continue;
ltype = map[type[l]];
iparam_jil = elem2param[jtype][itype][ltype];
iparam_ijl = elem2param[itype][jtype][ltype];
iparam_jl = elem2param[jtype][ltype][ltype];
iparam_lj = elem2param[ltype][jtype][jtype];
delrk[0] = x[l][0] - x[j][0];
delrk[1] = x[l][1] - x[j][1];
delrk[2] = x[l][2] - x[j][2];
rsq2 = vec3_dot(delrk,delrk);
delrl[0] = x[l][0] - x[j][0];
delrl[1] = x[l][1] - x[j][1];
delrl[2] = x[l][2] - x[j][2];
rsq2 = vec3_dot(delrl,delrl);
if (rsq2 > params[iparam_jl].cutsq) continue;
vec3_scale(-1,delrj,delrl); // ji_hat is -(ij_hat)
zeta_ji += zeta(&params[iparam_jil], &params[iparam_jl]
, rsq1, rsq2, delrl, delrk, j, ycn);
// radical initialization: apply only to CC,CO,OC bonds
if(params[iparam_ji].rad_flag > 0
&& params[iparam_jl].ielementgp == 1
&& params[iparam_jl].jelementgp == 1) {
iparam_lj = elem2param[ltype][jtype][jtype];
lcn=NCo[l];
lconjug += rad_init(rsq2,&params[iparam_lj],j,lradtot,lcn);
}
}
force_zeta(&params[iparam_ij], &params[iparam_ji],
rsq1, xcn, ycn, zeta_ij, zeta_ji, fpair,
prefac_ij1, prefac_ij2, prefac_ij3, prefac_ij4, prefac_ij5,
prefac_ji1, prefac_ji2, prefac_ji3, prefac_ji4, prefac_ji5,
eflag, eng_tmp, iq, jq, i, j, nj, bbtor, kconjug, lconjug);
evdwl = eng_tmp;
selfp6p(&params[iparam_ij],&params[iparam_ji],rsq1,eng_tmp,fpair);
evdwl += eng_tmp;
f[i][0] += delrj[0]*fpair;
f[i][1] += delrj[1]*fpair;
f[i][2] += delrj[2]*fpair;
f[j][0] -= delrj[0]*fpair;
f[j][1] -= delrj[1]*fpair;
f[j][2] -= delrj[2]*fpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,-fpair,-delrj[0],-delrj[1],-delrj[2]);
// attractive term via loop over k (3-body forces: i-j-k)
zet_addi=0;
zet_addj=0;
for (kk = 0; kk < sht_jnum; kk++) {
k = sht_jlist[kk];
if (j == k) continue;
sht_mlist = sht_first[k];
sht_mnum = sht_num[k];
ktype = map[type[k]];
iparam_ijk = elem2param[itype][jtype][ktype];
iparam_ikj = elem2param[itype][ktype][jtype];
iparam_jik = elem2param[jtype][itype][ktype];
iparam_ik = elem2param[itype][ktype][ktype];
delrk[0] = x[k][0] - xtmp;
delrk[1] = x[k][1] - ytmp;
delrk[2] = x[k][2] - ztmp;
rsq2 = vec3_dot(delrk,delrk);
if (rsq2 > params[iparam_ik].cutsq) continue;
// BO-dependent 3-body E & F
attractive(&params[iparam_ijk], &params[iparam_jik],&params[iparam_ikj],
prefac_ij1, prefac_ij2, prefac_ij3, prefac_ij4, prefac_ij5,
rsq1,rsq2,delrj,delrk,fi,fj,fk,i,xcn);
ep6p_ij = ep6p(&params[iparam_ijk],&params[iparam_ikj],rsq1,rsq2,delrj,delrk,zet_addi);
fp6p(&params[iparam_ijk],&params[iparam_ikj],rsq1,rsq2,delrj,delrk,fip6p,fjp6p,fkp6p);
// Sums up i-j-k forces: LP contribution
for (im = 0; im < 3; im++) {
fi[im] += fip6p[im];
fj[im] += fjp6p[im];
fk[im] += fkp6p[im];
}
// Sums up i-j-k forces: Tallies into global force vector
for (im = 0; im < 3; im++) {
f[i][im] += fi[im];
f[j][im] += fj[im];
f[k][im] += fk[im];
}
// torsion and radical: apply to all C-C bonds
if( params[iparam_ijk].tor_flag != 0 && fabs(ptorr)>1.0e-8) {
srmu = vec3_dot(delrj,delrk)/(sqrt(rsq1*rsq2));
srmu = sqrt(1.0-srmu*srmu);
if(srmu > 0.1) {
for (ll = 0; ll < sht_lnum; ll++) { // ll is neighbor of jj
l = sht_llist[ll];
if (l==i||l==j||l==k) continue;
ltype = map[type[l]];
delrl[0] = x[l][0] - x[j][0];
delrl[1] = x[l][1] - x[j][1];
delrl[2] = x[l][2] - x[j][2];
rsq3 = vec3_dot(delrl,delrl);
iparam_jl = elem2param[jtype][ltype][ltype];
if (rsq3 > params[iparam_jl].cutsq) continue;
iparam_ikl = elem2param[itype][ktype][ltype];
torindx = params[iparam_ij].tor_flag;
tor_force(torindx, &params[iparam_ikl], &params[iparam_jl],srmu,
rsq1,rsq2,rsq3,delrj,delrk,delrl);
for (im = 0; im < 3; im++) {
f[i][im] += fi_tor[im];
f[j][im] += fj_tor[im];
f[k][im] += fk_tor[im];
f[l][im] += fl_tor[im];
}
}
}
}
if( params[iparam_ijk].rad_flag>=1 &&
params[iparam_ijk].ielementgp==1 &&
params[iparam_ijk].kelementgp==1) {
iparam_ki = elem2param[ktype][itype][itype];
kcn=NCo[k];
double rik=sqrt(rsq2);
kradtot = -comb_fc(rik,&params[iparam_ki])*params[iparam_ki].pcross+kcn;
rad_forceik(&params[iparam_ki],rsq2,delrk,kconjug,kradtot);
for (im = 0; im < 3; im++) {
f[i][im] += fi_rad[im];
f[k][im] += fk_rad[im];
}
if (fabs(radtmp) > 1.0e-12) {
for (mm = 0; mm < sht_mnum; mm++) { // mm is neighbor of kk
m = sht_mlist[mm];
if (m == k) continue;
mtype = map[type[m]];
delrm[0] = x[m][0] - x[k][0];
delrm[1] = x[m][1] - x[k][1];
delrm[2] = x[m][2] - x[k][2];
rsq3 = vec3_dot(delrm,delrm);
iparam_km = elem2param[ktype][mtype][mtype];
iparam_ki = elem2param[ktype][itype][itype];
if (rsq3 > params[iparam_km].cutsq) continue;
rad_force(&params[iparam_km],rsq3,delrm,radtmp);
for (im = 0; im < 3; im++) {
f[k][im] += fj_rad[im];
f[m][im] += fk_rad[im];
}
}
}
}
if (evflag)
ev_tally(i,j,nlocal,newton_pair,ep6p_ij,0.0,0.0,0.0,0.0,0.0);
if (vflag_atom)
v_tally3(i,j,k,fj,fk,delrj,delrk);
} // k-loop
// attractive term via loop over l (3-body forces: j-i-l)
for (ll = 0; ll < sht_lnum; ll++) {
l = sht_llist[ll];
if (l == i) continue;
sht_plist = sht_first[l];
sht_pnum = sht_num[l];
ltype = map[type[l]];
iparam_jil = elem2param[jtype][itype][ltype];
iparam_jli = elem2param[jtype][ltype][itype];
iparam_ijl = elem2param[itype][jtype][ltype];
iparam_jl = elem2param[jtype][ltype][ltype];
delrk[0] = x[l][0] - x[j][0];
delrk[1] = x[l][1] - x[j][1];
delrk[2] = x[l][2] - x[j][2];
rsq2 = vec3_dot(delrk,delrk);
if (rsq2 > params[iparam_jl].cutsq) continue;
vec3_scale(-1,delrj,delrl);
attractive(&params[iparam_jil],&params[iparam_ijl],&params[iparam_jli],
prefac_ji1,prefac_ji2,prefac_ji3,prefac_ji4,prefac_ji5,
rsq1,rsq2,delrl,delrk,fj,fi,fl,j,ycn);
// BO-independent 3-body j-i-l LP and BB correction and forces
ep6p_ji = ep6p(&params[iparam_jil],&params[iparam_jli],rsq1,rsq2,delrl,delrk,zet_addj);
fp6p(&params[iparam_jil],&params[iparam_jli],rsq1,rsq2,delrl,delrk,fjp6p,fip6p,flp6p);
if (evflag)
ev_tally(j,i,nlocal,newton_pair,ep6p_ji,0.0,0.0,0.0,0.0,0.0);
// BO-dependent 3-body E & F
for (im = 0; im < 3; im++) {
fj[im] += fjp6p[im];
fi[im] += fip6p[im];
fl[im] += flp6p[im];
}
// Sums up j-i-l forces: Tallies into global force vector
for (im = 0; im < 3; im++) {
f[j][im] += fj[im];
f[i][im] += fi[im];
f[l][im] += fl[im];
}
// radical i-j-l-p: apply to all CC,CO,OC bonds
if( params[iparam_jil].rad_flag >= 1 &&
params[iparam_jil].ielementgp == 1 &&
params[iparam_jil].kelementgp == 1 ) {
iparam_lj = elem2param[ltype][jtype][jtype];
lcn=NCo[l];
double rjl=sqrt(rsq2);
lradtot=-comb_fc(rjl,&params[iparam_lj])*params[iparam_lj].pcross +lcn;
rad_forceik(&params[iparam_lj],rsq2,delrk,lconjug,lradtot);
for (im = 0; im < 3; im++) {
f[j][im] += fi_rad[im];
f[l][im] += fk_rad[im];
}
if (fabs(radtmp)>1.0e-12) {
for (pp = 0; pp < sht_pnum; pp++) { // pp is neighbor of ll
p = sht_plist[pp];
if (p == l) continue;
ptype = map[type[p]];
delrp[0] = x[p][0] - x[l][0];
delrp[1] = x[p][1] - x[l][1];
delrp[2] = x[p][2] - x[l][2];
rsq3 = vec3_dot(delrp,delrp);
iparam_lp = elem2param[ltype][ptype][ptype];
if (rsq3 > params[iparam_lp].cutsq) continue;
vec3_scale(-1,delrj,delrj);
rad_force(&params[iparam_lp],rsq3,delrp,radtmp);
vec3_scale(-1,delrj,delrj);
for (im = 0; im < 3; im++) {
f[l][im] += fj_rad[im];
f[p][im] += fk_rad[im];
}
}
}
}
if (vflag_atom)
v_tally3(j,i,l,fi,fl,delrl,delrk);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairComb3::repulsive(Param *parami, Param *paramj, double rsq,
double &fforce,int eflag, double &eng, double iq, double jq)
{
double r,tmp_fc,tmp_fc_d,Di,Dj;
double caj,vrcs,fvrcs;
double LamDiLamDj,fcdA,rlm1,bigA;
double romi = parami->addrep;
double rrcs = parami->bigr + parami->bigd;
double addr = parami->addrepr;
r = sqrt(rsq);
if (r > rrcs) return ;
tmp_fc = comb_fc(r,parami);
tmp_fc_d = comb_fc_d(r,parami);
Di = parami->DU + pow(fabs(parami->bD*(parami->QU-iq)),parami->nD);
Dj = paramj->DU + pow(fabs(paramj->bD*(paramj->QU-jq)),paramj->nD);
bigA = parami->bigA;
rlm1 = parami->lambda;
fcdA = tmp_fc_d - tmp_fc * rlm1;
LamDiLamDj = exp(0.5*(parami->lami*Di+paramj->lami*Dj)-rlm1*r);
caj = bigA * LamDiLamDj;
fforce = -caj * fcdA;
// additional repulsion
vrcs = 1.0; fvrcs = 0.0;
if (romi != 0.0 && r < addr) {
vrcs += romi * pow((1.0-r/addr),2.0);
fvrcs = romi * 2.0 * (r/addr-1.0)/addr;
fforce = fforce*vrcs - caj * tmp_fc * vrcs * fvrcs;
}
fforce /= r;
// eng = repulsive energy
eng = caj * tmp_fc * vrcs;
}
/* ---------------------------------------------------------------------- */
double PairComb3::zeta(Param *parami, Param *paramj, double rsqij,
double rsqik, double *delrij, double *delrik, int i, double xcn)
{
double rij,rik,costheta,arg,ex_delr,rlm3;
rij = sqrt(rsqij);
if (rij > parami->bigr+parami->bigd) return 0.0;
rik = sqrt(rsqik);
costheta = vec3_dot(delrij,delrik) / (rij*rik);
rlm3 = parami->beta;
arg = pow(rlm3*(rij-rik),int(parami->powermint));
if (arg > 69.0776) ex_delr = 1.e30;
else if (arg < -69.0776) ex_delr = 0.0;
else ex_delr = exp(arg);
return comb_fc(rik,paramj) * comb_gijk(costheta,parami,xcn) * ex_delr;
}
/* ---------------------------------------------------------------------- */
void PairComb3::selfp6p(Param *parami, Param *paramj, double rsq,
double &eng, double &force)
{
double r,comtti,comttj,fcj,fcj_d;
r=sqrt(rsq);
fcj=comb_fc(r,parami);
fcj_d=comb_fc_d(r,parami);
comtti = comttj = 0.0;
double pilp0 = parami->p6p0;
double pilp1 = parami->p6p1, pilp2 = parami->p6p2, pilp3 = parami->p6p3;
double pilp4 = parami->p6p4, pilp5 = parami->p6p5, pilp6 = parami->p6p6;
comtti = pilp0 + pilp1 + pilp2 + pilp3 + pilp4 + pilp5 + pilp6;
double pjlp0 = paramj->p6p0;
double pjlp1 = paramj->p6p1, pjlp2 = paramj->p6p2, pjlp3 = paramj->p6p3;
double pjlp4 = paramj->p6p4, pjlp5 = paramj->p6p5, pjlp6 = paramj->p6p6;
comttj = pjlp0 + pjlp1 + pjlp2 + pjlp3 + pjlp4 + pjlp5 + pjlp6;
eng = 0.5 * fcj * (comtti + comttj);
force += 0.5 * fcj_d * (comtti + comttj)/r;
}
/* ---------------------------------------------------------------------- */
double PairComb3::ep6p(Param *paramj, Param *paramk, double rsqij, double rsqik,
double *delrij, double *delrik , double &zet_add)
{
double comtt;
double pplp0 = paramj->p6p0;
double pplp1 = paramj->p6p1, pplp2 = paramj->p6p2, pplp3 = paramj->p6p3;
double pplp4 = paramj->p6p4, pplp5 = paramj->p6p5, pplp6 = paramj->p6p6;
double rij,rik,costheta,lp0,lp1,lp2,lp3,lp4,lp5,lp6;
double rmu,rmu2,rmu3,rmu4,rmu5,rmu6,fcj,fck;
comtt=0.0;
rij = sqrt(rsqij);
rik = sqrt(rsqik);
costheta = vec3_dot(delrij,delrik) / (rij*rik);
fcj = comb_fc(rij,paramj);
fck = comb_fc(rik,paramk);
rmu = costheta;
rmu2 = rmu *rmu; rmu3 = rmu2*rmu; rmu4 = rmu3*rmu;
rmu5 = rmu4*rmu; rmu6 = rmu5*rmu;
lp0 = pplp0;
lp1 = pplp1*rmu;
lp2 = pplp2*rmu2;
lp3 = pplp3*rmu3;
lp4 = pplp4*rmu4;
lp5 = pplp5*rmu5;
lp6 = pplp6*rmu6;
comtt = lp0 + lp1 + lp2 + lp3 + lp4 + lp5 + lp6;
return 0.5 * fck * comtt *fcj;
}
/*---------------------------------------------------------------------- */
void PairComb3::fp6p(Param *paramij,Param *paramik, double rsqij, double rsqik,
double *delrij, double *delrik, double *drilp,
double *drjlp, double *drklp)
{
double pplp0 = paramij->p6p0;
double pplp1 = paramij->p6p1, pplp2 = paramij->p6p2, pplp3 = paramij->p6p3;
double pplp4 = paramij->p6p4, pplp5 = paramij->p6p5, pplp6 = paramij->p6p6;
double ffj1,ffj2,ffk1,ffk2;
double rij,rik,costheta;
double rmu,comtt,comtt_d,com4k,com5,com5k,fcj,fcj_d,fck,fck_d;
double lp0,lp1,lp2,lp3,lp4,lp5,lp6;
double lp1_d,lp2_d,lp3_d,lp4_d,lp5_d,lp6_d;
double rmu2, rmu3, rmu4, rmu5, rmu6;
ffj1 = 0.0, ffj2 = 0.0;
ffk1 = 0.0, ffk2 = 0.0;
rij = sqrt(rsqij); rik = sqrt(rsqik);
costheta = vec3_dot(delrij,delrik) / (rij*rik);
fcj = comb_fc(rij,paramij);
fck = comb_fc(rik,paramik);
fcj_d = comb_fc_d(rij,paramij);
fck_d = comb_fc_d(rik,paramik);
rmu = costheta;
rmu2 = rmu *rmu; rmu3 = rmu2*rmu;
rmu4 = rmu3*rmu; rmu5 = rmu4*rmu; rmu6 = rmu5*rmu;
lp0 = pplp0;
lp1 = pplp1*rmu;
lp2 = pplp2*rmu2;
lp3 = pplp3*rmu3;
lp4 = pplp4*rmu4;
lp5 = pplp5*rmu5;
lp6 = pplp6*rmu6;
lp1_d = pplp1;
lp2_d = pplp2*2.0*rmu;
lp3_d = pplp3*3.0*rmu2;
lp4_d = pplp4*4.0*rmu3;
lp5_d = pplp5*5.0*rmu4;
lp6_d = pplp6*6.0*rmu5;
comtt = lp0 + lp1 + lp2 + lp3 + lp4 + lp5 + lp6;
comtt_d = lp1_d + lp2_d + lp3_d + lp4_d + lp5_d + lp6_d;
com4k = fcj * fck_d * comtt;
com5 = fcj * fck * comtt_d;
com5k = fck * comtt * fcj_d;
ffj1 = 0.5*(-com5/(rij*rik));
ffj2 = 0.5*(com5*rmu/rsqij-com5k/rij);
ffk1 = ffj1;
ffk2 = 0.5*(-com4k/rik+com5*rmu/rsqik);
// j-atom
vec3_scale(ffj1,delrik,drjlp);
vec3_scaleadd(ffj2,delrij,drjlp,drjlp);
// k-atom
vec3_scale(ffk1,delrij,drklp);
vec3_scaleadd(ffk2,delrik,drklp,drklp);
// i-atom
vec3_add(drjlp,drklp,drilp);
vec3_scale(-1.0,drilp,drilp);
}
/* ---------------------------------------------------------------------- */
void PairComb3::force_zeta(Param *parami, Param *paramj, double rsq,
double xcn, double ycn, double &zeta_ij, double &zeta_ji, double &fforce,
double &prefac_ij1, double &prefac_ij2, double &prefac_ij3,
double &prefac_ij4, double &prefac_ij5,
double &prefac_ji1, double &prefac_ji2, double &prefac_ji3,
double &prefac_ji4, double &prefac_ji5,
int eflag, double &eng, double iq, double jq,
int i, int j, int nj, double bbtor, double kconjug, double lconjug)
{
double r,att_eng,att_force,bij; // att_eng is -cbj
double boij, dbij1, dbij2, dbij3, dbij4, dbij5;
double boji, dbji1, dbji2, dbji3, dbji4, dbji5;
double pradx, prady;
r = sqrt(rsq);
if (r > parami->bigr + parami->bigd) return;
comb_fa(r, parami, paramj, iq, jq, att_eng, att_force);
comb_bij_d(zeta_ij,parami,r,i,boij,dbij1,dbij2,dbij3,dbij4,dbij5,xcn);
comb_bij_d(zeta_ji,paramj,r,j,boji,dbji1,dbji2,dbji3,dbji4,dbji5,ycn);
bij = 0.5*(boij + boji);
// radical energy
if ( parami->rad_flag>0 ) {
rad_calc( r, parami, paramj, kconjug, lconjug, i, j, xcn, ycn);
bij += brad[0];
pradx = brad[1]*att_eng;
prady = brad[2]*att_eng;
brad[3] = 1.0 * brad[3]*att_eng;
}
// torsion energy
if ( parami->tor_flag!=0) {
tor_calc( r, parami, paramj, kconjug, lconjug, i, j, xcn, ycn);
bij += btor[0] * bbtor;
ptorr = att_eng * btor[0];
pradx += 1.0 * btor[1] * bbtor * att_eng;
prady += 1.0 * btor[2] * bbtor * att_eng;
brad[3]+= 1.0 * btor[3] * bbtor * att_eng;
}
fforce = 1.0*bij*att_force/r; // divide by r will done compute
bbij[i][nj] = bij;
prefac_ij1 = -0.5*att_eng*dbij1; // prefac_ij1 = -pfij
prefac_ij2 = -0.5*att_eng*dbij2; // prefac_ij2 = -pfij1
prefac_ij3 = -0.5*att_eng*dbij3; // prefac_ij3 = -pfij2
prefac_ij4 = -0.5*att_eng*dbij4; // prefac_ij4 = -pfij3
prefac_ij5 = -0.5*att_eng*dbij5; // prefac_ij5 = -pfij4
prefac_ji1 = -0.5*att_eng*dbji1; // prefac_ji1 = -pfji
prefac_ji2 = -0.5*att_eng*dbji2; // prefac_ji2 = -pfji1
prefac_ji3 = -0.5*att_eng*dbji3; // prefac_ji3 = -pfji2
prefac_ji4 = -0.5*att_eng*dbji4; // prefac_ji4 = -pfji3
prefac_ji5 = -0.5*att_eng*dbji5; // prefac_ji5 = -pfji4
// combines com6 & com7 below
if ( parami->rad_flag>0 || parami->tor_flag!=0 ) {
prefac_ij2-=pradx;
prefac_ji2-=prady;
}
// eng = attraction energy
if (eflag) eng = 1.0*bij*att_eng;
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fc(double r, Param *param)
{
double r_inn = param->bigr - param->bigd;
double r_out = param->bigr + param->bigd;
if (r <= r_inn) return 1.0;
if (r >= r_out) return 0.0;
return 0.5*(1.0 + cos(MY_PI*(r-r_inn)/(r_out-r_inn)));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fc_d(double r, Param *param)
{
double r_inn = param->bigr - param->bigd;
double r_out = param->bigr + param->bigd;
if (r <= r_inn) return 0.0;
if (r >= r_out) return 0.0;
return -MY_PI2/(r_out-r_inn)*sin(MY_PI*(r-r_inn)/(r_out-r_inn));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fccc(double xcn)
{
double cut1 = ccutoff[0];
double cut2 = ccutoff[1];
if (xcn <= cut1) return 1.0;
if (xcn >= cut2) return 0.0;
return 0.5*(1.0 + cos(MY_PI*(xcn-cut1)/(cut2-cut1)));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fccc_d(double xcn)
{
double cut1 = ccutoff[0];
double cut2 = ccutoff[1];
if (xcn <= cut1) return 0.0;
if (xcn >= cut2) return 0.0;
return -MY_PI2/(cut2-cut1)*sin(MY_PI*(xcn-cut1)/(cut2-cut1));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fcch(double xcn)
{
double cut1 = ccutoff[2];
double cut2 = ccutoff[3];
if (xcn <= cut1) return 1.0;
if (xcn >= cut2) return 0.0;
return 0.5*(1.0 + cos(MY_PI*(xcn-cut1)/(cut2-cut1)));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fcch_d(double xcn)
{
double cut1 = ccutoff[2];
double cut2 = ccutoff[3];
if (xcn <= cut1) return 0.0;
if (xcn >= cut2) return 0.0;
return -MY_PI2/(cut2-cut1)*sin(MY_PI*(xcn-cut1)/(cut2-cut1));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fccch(double xcn)
{
double cut1 = ccutoff[4];
double cut2 = ccutoff[5];
if (xcn <= cut1) return 1.0;
if (xcn >= cut2) return 0.0;
return 0.5*(1.0 + cos(MY_PI*(xcn-cut1)/(cut2-cut1)));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fccch_d(double xcn)
{
double cut1 = ccutoff[4];
double cut2 = ccutoff[5];
if (xcn <= cut1) return 0.0;
if (xcn >= cut2) return 0.0;
return -MY_PI2/(cut2-cut1)*sin(MY_PI*(xcn-cut1)/(cut2-cut1));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fcsw(double rsq)
{
double r = sqrt(rsq);
if (r <= chicut1) return 1.0;
if (r >= chicut2) return 0.0;
return 0.5*(1.0 + cos(MY_PI*(r-chicut1)/(chicut2-chicut1)));
}
/* ---------------------------------------------------------------------- */
double PairComb3::self(Param *param, double qi)
{
double self_tmp, cmin, cmax, qmin, qmax;
double s1=param->chi, s2=param->dj, s3=param->dk, s4=param->dl;
self_tmp = 0.0;
qmin = param->qmin;
qmax = param->qmax;
cmin = cmax = 100.0;
self_tmp = qi*(s1+qi*(s2+qi*(s3+qi*s4)));
if (qi < qmin) self_tmp += cmin * pow((qi-qmin),4);
if (qi > qmax) self_tmp += cmax * pow((qi-qmax),4);
return self_tmp;
}
/* ---------------------------------------------------------------------- */
void PairComb3::comb_fa(double r, Param *parami, Param *paramj, double iq,
double jq, double &att_eng, double &att_force)
{
double Bsi;
double qi,qj,Di,Dj;
double AlfDiAlfDj, YYBn, YYBj;
double alfij1= parami->alpha1;
double alfij2= parami->alpha2;
double alfij3= parami->alpha3;
double pbij1= parami->bigB1;
double pbij2= parami->bigB2;
double pbij3= parami->bigB3;
if (r > parami->bigr + parami->bigd) Bsi = 0.0;
qi = iq; qj = jq;
Di = Dj = Bsi = 0.0;
Di = parami->DU + pow(fabs(parami->bD*(parami->QU-qi)),parami->nD);
Dj = paramj->DU + pow(fabs(paramj->bD*(paramj->QU-qj)),paramj->nD);
YYBn = (parami->aB-fabs(pow(parami->bB*(qi-parami->Qo),10)));
YYBj = (paramj->aB-fabs(pow(paramj->bB*(qj-paramj->Qo),10)));
if (YYBn*YYBj > 0.0 ) {
AlfDiAlfDj = exp(0.5*(parami->alfi*Di+paramj->alfi*Dj));
Bsi = (pbij1*exp(-alfij1*r)+pbij2*exp(-alfij2*r)+pbij3*exp(-alfij3*r))*
sqrt(YYBn*YYBj)*AlfDiAlfDj; // Bsi is cbj
att_eng = -Bsi * comb_fc(r,parami);
att_force = -(Bsi*comb_fc_d(r,parami)-comb_fc(r,parami)*sqrt(YYBn*YYBj)*
AlfDiAlfDj*(alfij1*pbij1*exp(-alfij1*r)+
alfij2*pbij2*exp(-alfij2*r)+alfij3*pbij3*exp(-alfij3*r)));
} else {
att_eng = 0.0;
att_force = 0.0;
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::comb_bij_d(double zet, Param *param, double r, int i,
double &tbij, double &tbij1, double &tbij2,
double &tbij3, double &tbij4, double &tbij5, double xcn)
{
double pcorn,dpcorn,dxccij,dxchij,dxcoij;
double zeta = zet;
double zetang,tmp_tbij, pow_n;
pcorn = dpcorn = dxccij = dxchij = dxcoij = 0.0;
coord(param,r,i,pcorn,dpcorn,dxccij,dxchij,dxcoij,xcn); // coordination term
zetang=zeta;
pow_n=param->powern;
zeta = pow(zetang,pow_n)+pcorn;
tmp_tbij=pow_n*pow(zetang,(pow_n-1.0));
if ((1.0 + zeta) < 0.1 ){
zeta=0.1-1.0;
tbij = pow(1.0 + zeta, -0.5/pow_n);
tbij1=0.0;
}
else if (zeta > param->c1) {
tbij = pow(zeta,-0.5/pow_n);
tbij1 = -0.5/pow_n*pow(zeta,(-0.5/pow_n-1.0));
} else if (zeta > param->c2) {
tbij = pow(zeta,-0.5/pow_n)-0.5/pow_n*pow(zeta,(-0.5/pow_n-1.0));
tbij1 = -0.5/pow_n/zeta;
} else if (fabs(zeta) < param->c4) {
tbij = 1.0;
tbij1 = 0.0;
} else if (fabs(zeta) < param->c3) {
tbij = 1.0 - zeta/(2.0*pow_n);
tbij1 = -1/(2.0*pow_n);
} else {
tbij = pow(1.0 + zeta, -0.5/pow_n);
tbij1 = -0.5/pow_n * pow(1.0 + zeta,(-1.0-0.5/pow_n));
}
tbij2 = tbij1 * dpcorn;
tbij3 = tbij1 * dxccij;
tbij4 = tbij1 * dxchij;
tbij5 = tbij1 * dxcoij;
tbij1 = tbij1 * tmp_tbij;
}
/* ---------------------------------------------------------------------- */
void PairComb3::coord(Param *param, double r, int i,
double &pcorn, double &dpcorn, double &dxccij,
double &dxchij, double &dxcoij, double xcn)
{
int ixmin,iymin,izmin;
double xcntot,xcccn,xchcn,xcocn;
int tri_flag= param-> pcn_flag;
int jele_gp= param->jelementgp;
double pan = param->pcna;
double pbn = param->pcnb;
double pcn = param->pcnc;
double pdn = param->pcnd;
xcccn = xchcn = xcocn = 0.0;
xcccn = xcctmp[i];
xchcn = xchtmp[i];
xcocn = xcotmp[i];
xcntot = -comb_fc(r,param)*param->pcross + xcn;
pcorn = dpcorn = dxccij = dxchij = dxcoij = 0.0;
pcorn = 0.0; dpcorn = 0.0;
if(xcntot < 0.0) xcntot = 0.0;
if (tri_flag>0) {
if(jele_gp==1) xcccn = xcccn-comb_fc(r,param)*param->pcross;
if(jele_gp==2) xchcn = xchcn-comb_fc(r,param)*param->pcross;
if(jele_gp==3) xcocn = xcocn-comb_fc(r,param)*param->pcross;
if(xcccn < 0.0) xcccn = 0.0;
if(xchcn < 0.0) xchcn = 0.0;
if(xcocn < 0.0) xcocn = 0.0;
if(xcccn > maxx) xcccn = maxx;
if(xchcn > maxy) xchcn = maxy;
if(xcocn > maxz) xcocn = maxz;
double xcntritot=xcccn+xchcn+xcocn;
if(xcntritot > maxxcn[tri_flag-1]) {
pcorn = vmaxxcn[tri_flag-1]+(xcntot-maxxcn[tri_flag-1])*dvmaxxcn[tri_flag-1];
dxccij = dxchij = dxcoij = dvmaxxcn[tri_flag-1];
}
else {
ixmin=int(xcccn+1.0e-12);
iymin=int(xchcn+1.0e-12);
izmin=int(xcocn+1.0e-12);
if (fabs(float(ixmin)-xcccn)>1.0e-8 ||
fabs(float(iymin)-xchcn)>1.0e-8 ||
fabs(float(izmin)-xcocn)>1.0e-8) {
cntri_int(tri_flag,xcccn,xchcn,xcocn,ixmin,iymin,izmin,
pcorn,dxccij,dxchij,dxcoij,param);
}
else {
pcorn = pcn_grid[tri_flag-1][ixmin][iymin][izmin];
dxccij = pcn_gridx[tri_flag-1][ixmin][iymin][izmin];
dxchij = pcn_gridy[tri_flag-1][ixmin][iymin][izmin];
dxcoij = pcn_gridz[tri_flag-1][ixmin][iymin][izmin];
}
}
} else {
pcorn = pan*xcntot+pbn*exp(pcn*xcntot)+pdn;
dpcorn = pan+pbn*pcn*exp(pcn*xcntot);
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::cntri_int(int tri_flag, double xval, double yval,
double zval, int ixmin, int iymin, int izmin, double &vval,
double &dvalx, double &dvaly, double &dvalz, Param *param)
{
double x;
vval = 0.0; dvalx = 0.0; dvaly = 0.0; dvalz = 0.0;
if(ixmin >= maxx-1) { ixmin=maxx-1; }
if(iymin >= maxy-1) { iymin=maxy-1; }
if(izmin >= maxz-1) { izmin=maxz-1; }
for (int j=0; j<64; j++) {
x = pcn_cubs[tri_flag-1][ixmin][iymin][izmin][j]
*pow(xval,iin3[j][0])*pow(yval,iin3[j][1])
*pow(zval,iin3[j][2]);
vval += x;
if(xval>1.0e-8) {dvalx += x*iin3[j][0]/xval;}
if(yval>1.0e-8) {dvaly += x*iin3[j][1]/yval;}
if(zval>1.0e-8) {dvalz += x*iin3[j][2]/zval;}
}
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_gijk(double costheta, Param *param, double nco_tmp)
{
double rmu1 = costheta;
double rmu2 = rmu1*rmu1;
double rmu3 = rmu2*rmu1;
double rmu4 = rmu3*rmu1;
double rmu5 = rmu4*rmu1;
double rmu6 = rmu5*rmu1;
double co6 = param->pcos6*rmu6;
double co5 = param->pcos5*rmu5;
double co4 = param->pcos4*rmu4;
double co3 = param->pcos3*rmu3;
double co2 = param->pcos2*rmu2;
double co1 = param->pcos1*rmu1;
double co0 = param->pcos0;
double pcross = param->pcross;
double gmu;
if (param->ang_flag==1) {
double qtheta, gmu1, gmu2, rrmu, astep;
int k;
qtheta = comb_fccc(nco_tmp);
astep = 2.0/ntab;
rrmu = (rmu1+1.0)/astep;
k = int(rrmu);
gmu1 = co6+co5+co4+co3+co2+co1+co0;
gmu2 = pang[k]+(pang[k+1]-pang[k])*(rrmu-k);
gmu = gmu2+qtheta*(gmu1-gmu2);
return gmu*pcross;
} else if (param->ang_flag==2){
double qtheta, gmu1, gmu2;
double ch6 = ch_a[6]*rmu6;
double ch5 = ch_a[5]*rmu5;
double ch4 = ch_a[4]*rmu4;
double ch3 = ch_a[3]*rmu3;
double ch2 = ch_a[2]*rmu2;
double ch1 = ch_a[1]*rmu1;
double ch0 = ch_a[0];
qtheta = comb_fccch(nco_tmp);
gmu1 = co6+co5+co4+co3+co2+co1+co0;
gmu2 = ch6+ch5+ch4+ch3+ch2+ch1+ch0;
gmu = gmu2+qtheta*(gmu1-gmu2);
return gmu*pcross;
} else {
gmu = co6+co5+co4+co3+co2+co1+co0;
return gmu*pcross;
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::comb_gijk_d(double costheta, Param *param, double nco_tmp,
double &gijk_d, double &com3jk)
{
double rmu1 = costheta;
double rmu2 = rmu1*rmu1;
double rmu3 = rmu2*rmu1;
double rmu4 = rmu3*rmu1;
double rmu5 = rmu4*rmu1;
double rmu6 = rmu5*rmu1;
double co6 = param->pcos6; //*rmu5*6.0;
double co5 = param->pcos5; //*rmu4*5.0;
double co4 = param->pcos4; //*rmu3*4.0;
double co3 = param->pcos3; //*rmu2*3.0;
double co2 = param->pcos2; //*rmu1*2.0;
double co1 = param->pcos1;
double co0 = param->pcos0;
double pcross = param->pcross;
gijk_d = com3jk = 0.0;
if (param->ang_flag==1) {
double qtheta, dqtheta, gmu1, gmu2, dgmu1,dgmu2, rrmu, astep;
int k;
qtheta = comb_fccc(nco_tmp);
dqtheta = comb_fccc_d(nco_tmp);
astep = 2.0/ntab;
rrmu = (rmu1+1.0)/astep;
k = int(rrmu);
gmu1 =rmu6*co6+rmu5*co5+rmu4*co4
+rmu3*co3+rmu2*co2+rmu1*co1+co0;
dgmu1 =6.0*rmu5*co6+5.0*rmu4*co5+4.0*rmu3*co4
+3.0*rmu2*co3+2.0*rmu1*co2+co1;
gmu2 = pang[k]+(pang[k+1]-pang[k])*(rrmu-k);
dgmu2 = dpang[k]+(dpang[k+1]-dpang[k])*(rrmu-k);
gijk_d = pcross*(dgmu2+qtheta*(dgmu1-dgmu2));
com3jk = dqtheta * (gmu1-gmu2);
} else if(param->ang_flag==2) {
double qtheta, dqtheta, gmu1, gmu2, dgmu1,dgmu2;
double ch6 = ch_a[6];
double ch5 = ch_a[5];
double ch4 = ch_a[4];
double ch3 = ch_a[3];
double ch2 = ch_a[2];
double ch1 = ch_a[1];
double ch0 = ch_a[0];
qtheta = comb_fccch(nco_tmp);
dqtheta = comb_fccch_d(nco_tmp);
gmu1 =rmu6*co6+rmu5*co5+rmu4*co4
+rmu3*co3+rmu2*co2+rmu1*co1+co0;
dgmu1 =6.0*rmu5*co6+5.0*rmu4*co5+4.0*rmu3*co4
+3.0*rmu2*co3+2.0*rmu1*co2+co1;
gmu2 =rmu6*ch6+rmu5*ch5+rmu4*ch4
+rmu3*ch3+rmu2*ch2+rmu1*ch1+ch0;
dgmu2 =6.0*rmu5*ch6+5.0*rmu4*ch5+4.0*rmu3*ch4
+3.0*rmu2*ch3+2.0*rmu1*ch2+ch1;
gijk_d = pcross*(dgmu2+qtheta*(dgmu1-dgmu2));
com3jk = dqtheta * (gmu1-gmu2);
} else {
gijk_d = pcross*(6.0*rmu5*co6+5.0*rmu4*co5+4.0*rmu3*co4
+3.0*rmu2*co3+2.0*rmu1*co2+co1);
com3jk = 0.0;
}
}
/*------------------------------------------------------------------------- */
void PairComb3::attractive(Param *parami, Param *paramj , Param *paramk, double prefac_ij1,
double prefac_ij2, double prefac_ij3, double prefac_ij4,
double prefac_ij5, double rsqij, double rsqik, double *delrij,
double *delrik, double *fi, double *fj,double *fk, int i, double xcn)
{
double rij_hat[3],rik_hat[3];
double rij,rijinv,rik,rikinv;
rij = sqrt(rsqij);
rijinv = 1.0/rij;
vec3_scale(rijinv,delrij,rij_hat);
rik = sqrt(rsqik);
rikinv = 1.0/rik;
vec3_scale(rikinv,delrik,rik_hat);
comb_zetaterm_d(prefac_ij1, prefac_ij2, prefac_ij3, prefac_ij4, prefac_ij5,
rij_hat, rij,rik_hat, rik, fi, fj, fk, parami, paramj, paramk,xcn);
}
/* ---------------------------------------------------------------------- */
void PairComb3::comb_zetaterm_d(double prefac_ij1, double prefac_ij2,
double prefac_ij3, double prefac_ij4, double prefac_ij5,
double *rij_hat, double rij, double *rik_hat, double rik, double *dri,
double *drj, double *drk, Param *parami, Param *paramj, Param *paramk, double xcn)
{
double gijk,gijk_d,ex_delr,ex_delr_d,fc_k,cos_theta,tmp,rlm3;
double dcosdri[3],dcosdrj[3],dcosdrk[3],dfc_i,dfc_k;
double com6, com3j, com3k, com3jk;
int mint = int(parami->powermint);
double pcrossi = parami->pcross;
double pcrossj = paramj->pcross;
double pcrossk = paramk->pcross;
int icontrol = parami->pcn_flag;
dfc_i = comb_fc_d(rij,parami);
fc_k = comb_fc(rik,paramk);
dfc_k = comb_fc_d(rik,paramk);
rlm3 = parami->beta;
tmp = pow(rlm3*(rij-rik),mint);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
ex_delr *= pcrossi;
cos_theta = vec3_dot(rij_hat,rik_hat);
gijk = comb_gijk(cos_theta,parami,xcn);
comb_gijk_d(cos_theta,parami,xcn,gijk_d,com3jk);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// com6 & com7
if(icontrol > 0){
if(parami->kelementgp==1) {com6 = prefac_ij3*pcrossk*dfc_k;}
if(parami->kelementgp==2) {com6 = prefac_ij4*pcrossk*dfc_k;}
if(parami->kelementgp==3) {com6 = prefac_ij5*pcrossk*dfc_k;}
if(parami->rad_flag>=1 || parami->tor_flag!=0)
{com6+=prefac_ij2*pcrossk*dfc_k;}
} else {
com6 = prefac_ij2*pcrossi*dfc_k;
}
if (parami->ang_flag==1 || parami->ang_flag==2) {
com3j = com3jk*ex_delr*pcrossk*pcrossj*fc_k*dfc_i;
com3k = com3jk*ex_delr*pcrossk*pcrossk*fc_k*dfc_k;
} else {
com3j = 0.0;
com3k = 0.0;
}
ex_delr_d = mint*pow(rlm3,mint)*pow((rij-rik),(mint-1))*ex_delr; // com3
vec3_scale(-dfc_k*gijk*ex_delr,rik_hat,dri); // com1
vec3_scaleadd(fc_k*gijk_d*ex_delr,dcosdri,dri,dri); // com2
vec3_scaleadd(fc_k*gijk*ex_delr_d,rik_hat,dri,dri); // com3 cont'd
vec3_scaleadd(-fc_k*gijk*ex_delr_d,rij_hat,dri,dri); // com3 sums j
vec3_scaleadd(-com3k,rik_hat,dri,dri); // com3k
vec3_scaleadd(-com3j,rij_hat,dri,dri); // com3j
vec3_scale(prefac_ij1,dri,dri);
vec3_scaleadd(-com6,rik_hat,dri,dri); // com6
vec3_scale(fc_k*gijk_d*ex_delr,dcosdrj,drj); // com2
vec3_scaleadd(fc_k*gijk*ex_delr_d,rij_hat,drj,drj); // com3 cont'd
vec3_scaleadd(com3j,rij_hat,drj,drj); // com3j
vec3_scale(prefac_ij1,drj,drj);
vec3_scale(dfc_k*gijk*ex_delr,rik_hat,drk); // com1
vec3_scaleadd(fc_k*gijk_d*ex_delr,dcosdrk,drk,drk); // com2
vec3_scaleadd(-fc_k*gijk*ex_delr_d,rik_hat,drk,drk); // com3 cont'd
vec3_scaleadd(com3k,rik_hat,drk,drk); // com3k
vec3_scale(prefac_ij1,drk,drk);
vec3_scaleadd(com6,rik_hat,drk,drk); // com6
}
/* ---------------------------------------------------------------------- */
void PairComb3::costheta_d(double *rij_hat, double rij, double *rik_hat,
double rik, double *dri, double *drj, double *drk)
{
double cos_theta = vec3_dot(rij_hat,rik_hat);
vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj);
vec3_scale(1.0/rij,drj,drj);
vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk);
vec3_scale(1.0/rik,drk,drk);
vec3_add(drj,drk,dri);
vec3_scale(-1.0,dri,dri);
}
/* ---------------------------------------------------------------------- */
void PairComb3::tables()
{
int i,j,k,m, nntypes, ncoul,nnbuf, ncoul_lim, inty, itype, jtype;
int iparam_i, iparam_ij, iparam_ji;
double r,dra,drin,drbuf,rc,z,zr,zrc,ea,eb,ea3,eb3,alf;
double exp2er,exp2ersh,fafash,dfafash,F1,dF1,ddF1,E1,E2,E3,E4;
double exp2ear,exp2ebr,exp2earsh,exp2ebrsh,fafbsh,dfafbsh;
double afbshift, dafbshift, exp2ershift;
int n = nelements;
dra = 0.001;
drin = 0.100;
drbuf = 0.100;
nnbuf = int(drbuf/dra) +1;
rc = cutmax;
alf = 0.20;
nmax = atom->nmax;
nntypes = int((n+1)*n/2.0)+1;
ncoul = int((rc-drin)/dra)+ nnbuf;
ncoul_lim = int(ncoul * 1.20);
// allocate arrays
memory->create(intype,n,n,"pair:intype");
memory->create(erpaw,ncoul_lim,3,"pair:erpaw");
memory->create(fafb,ncoul_lim,nntypes,"pair:fafb");
memory->create(dfafb,ncoul_lim,nntypes,"pair:dfafb");
memory->create(ddfafb,ncoul_lim,nntypes,"pair:ddfafb");
memory->create(phin,ncoul_lim,nntypes,"pair:phin");
memory->create(dphin,ncoul_lim,nntypes,"pair:dphin");
memory->create(afb,ncoul_lim,nntypes,"pair:afb");
memory->create(dafb,ncoul_lim,nntypes,"pair:dafb");
memory->create(vvdw,ncoul,nntypes,"pair:vvdw");
memory->create(vdvdw,ncoul,nntypes,"pair:vdvdw");
memory->create(dpl,nmax,3,"pair:dpl");
memory->create(bbij,nmax,MAXNEIGH,"pair:bbij");
memory->create(xcctmp,nmax,"pair:xcctmp");
memory->create(xchtmp,nmax,"pair:xchtmp");
memory->create(xcotmp,nmax,"pair:xcotmp");
memory->create(NCo,nmax,"pair:NCo");
memory->create(sht_num,nmax,"pair:sht_num");
sht_first = (int **) memory->smalloc(nmax*sizeof(int *),
"pair:sht_first");
// set interaction number: 0-0=0, 1-1=1, 0-1=1-0=2
m = 0; k = n;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
if (j == i) {
intype[i][j] = m;
m += 1;
} else if (j != i && j > i) {
intype[i][j] = k;
k += 1;
} else if (j != i && j < i) {
intype[i][j] = intype[j][i];
}
}
}
// default arrays to zero
for (i = 0; i < ncoul; i ++) {
for (j = 0; j < nntypes; j ++) {
fafb[i][j] = 0.0;
dfafb[i][j] = 0.0;
ddfafb[i][j] = 0.0;
phin[i][j] = 0.0;
dphin[i][j] = 0.0;
afb[i][j] = 0.0;
dafb[i][j] = 0.0;
}
}
// direct 1/r energy with Slater 1S orbital overlap
for (i = 0; i < n; i++) {
if (map[i+1] < 0) continue;
r = drin - dra;
itype = map[i+1];
iparam_i = elem2param[itype][itype][itype];
z = params[iparam_i].esm;
exp2ershift = exp(-2.0*z*rc);
afbshift = -exp2ershift*(z+1.0/rc);
dafbshift = exp2ershift*(2.0*z*z+2.0*z/rc+1.0/(rc*rc));
if (comm->me == 0 && screen)
fprintf(screen," element[%d] = %-2s, z = %g\n",i+1,elements[map[i+1]],z);
for (j = 0; j < ncoul; j++) {
exp2er = exp(-2.0 * z * r);
phin[j][i] = 1.0 - exp2er * (1.0 + 2.0 * z * r * (1.0 + z * r));
dphin[j][i] = (4.0 * exp2er * z * z * z * r * r);
afb[j][i] = -exp2er*(z+1.0/r)-afbshift-(r-rc)*dafbshift;
dafb[j][i] = -(exp2er*(2.0*z*z+2.0*z/r+1.0/(r*r))-dafbshift);
r += dra;
}
}
for (i = 0; i < n; i ++) {
if (map[i+1] < 0) continue;
for (j = 0; j < n; j ++) {
if (map[j+1] < 0) continue;
r = drin - dra;
if (j == i) {
itype = map[i+1];
inty = intype[itype][itype];
iparam_i = elem2param[itype][itype][itype];
z = params[iparam_i].esm;
zrc = z * rc;
exp2ersh = exp(-2.0 * zrc);
fafash = -exp2ersh * (1.0 / rc +
z * (11.0/8.0 + 3.0/4.0*zrc + zrc*zrc/6.0));
dfafash = exp2ersh * (1.0/(rc*rc) + 2.0*z/rc +
z*z*(2.0 + 7.0/6.0*zrc + zrc*zrc/3.0));
for (k = 0; k < ncoul; k ++) {
zr = z * r;
exp2er = exp(-2.0*zr);
F1 = -exp2er * (1.0 / r +
z * (11.0/8.0 + 3.0/4.0*zr + zr*zr/6.0));
dF1 = exp2er * (1.0/(r*r) + 2.0*z/r +
z*z*(2.0 + 7.0/6.0*zr + zr*zr/3.0));
ddF1 = -exp2er * (2.0/(r*r*r) + 4.0*z/(r*r) + 4.0*z*z/r +
z*z*z/3.0*(17.0/2.0 + 5.0*zr + 2.0*zr*zr));
fafb[k][inty] = F1-fafash-(r-rc)*dfafash;
dfafb[k][inty] = -(dF1 - dfafash);
ddfafb[k][inty] = ddF1;
r += dra;
}
} else if (j != i) {
itype = map[i+1];
jtype = map[j+1];
inty = intype[itype][jtype];
iparam_ij = elem2param[itype][jtype][jtype];
ea = params[iparam_ij].esm;
ea3 = ea*ea*ea;
iparam_ji = elem2param[jtype][itype][itype];
eb = params[iparam_ji].esm;
eb3 = eb*eb*eb;
E1 = ea*eb3*eb/((ea+eb)*(ea+eb)*(ea-eb)*(ea-eb));
E2 = eb*ea3*ea/((ea+eb)*(ea+eb)*(eb-ea)*(eb-ea));
E3 = (3.0*ea*ea*eb3*eb-eb3*eb3) /
((ea+eb)*(ea+eb)*(ea+eb)*(ea-eb)*(ea-eb)*(ea-eb));
E4 = (3.0*eb*eb*ea3*ea-ea3*ea3) /
((ea+eb)*(ea+eb)*(ea+eb)*(eb-ea)*(eb-ea)*(eb-ea));
exp2earsh = exp(-2.0*ea*rc);
exp2ebrsh = exp(-2.0*eb*rc);
fafbsh = -exp2earsh*(E1 + E3/rc)-exp2ebrsh*(E2 + E4/rc);
dfafbsh =
exp2earsh*(2.0*ea*(E1+E3/rc)+E3/(rc*rc)) +
exp2ebrsh*(2.0*eb*(E2+E4/rc)+E4/(rc*rc));
for (k = 0; k < ncoul; k ++) {
exp2ear = exp(-2.0*ea*r);
exp2ebr = exp(-2.0*eb*r);
fafb[k][inty] =
- exp2ear*(E1+E3/r) - exp2ebr*(E2+E4/r)
- fafbsh - (r-rc) * dfafbsh;
dfafb[k][inty] = -(exp2ear*(2.0*ea*(E1+E3/r) + E3/(r*r))
+ exp2ebr*(2.0*eb*(E2+E4/r) + E4/(r*r))- dfafbsh);
ddfafb[k][inty] = -exp2ear*(4.0*ea*ea*(E1+E3/r)+4.0*ea*E3/(r*r)
+2.0*E3/(r*r*r))
-exp2ebr*(4.0*eb*eb*(E2+E4/r)+4.0*eb*E4/(r*r)
+2.0*E4/(r*r*r));
r += dra;
}
}
}
}
for (i = 0; i < ncoul_lim; i ++) {
r = dra * (i-1) + drin;
erpaw[i][0] = erfc(r*alf);
erpaw[i][1] = exp(-r*r*alf*alf);
}
// end wolf summation
// van der Waals
int ii,jj;
double **rvdw, *cc2, *cc3, *vrc, *rrc;
double r6, r7, r12, r13, rf6, rf12, drf7, drf13;
double drcc, temp6, temp7, temp12, temp13;
double vsigt, vepst, vdwt, dvdwt;
vrc = new double[13];
rrc = new double[13];
cc2 = new double[nntypes];
cc3 = new double[nntypes];
memory->create(rvdw,2,nntypes,"pair:rvdw");
vrc[0] = rc;
for (i=1; i<13; i++) {
vrc[i] = vrc[i-1] * vrc[0];
}
// generate spline coefficients for CC, CH, HH vdw
for (ii = 0; ii < n; ii ++) {
for (jj = ii; jj < n; jj ++) {
itype = ii;
jtype = jj;
inty = intype[itype][jtype];
iparam_ij = elem2param[itype][jtype][jtype];
// parameter check: eps > 0
if(params[iparam_ij].vdwflag > 0) {
if(params[iparam_ij].vdwflag==1){
rvdw[0][inty] = params[iparam_ij].bigr + params[iparam_ij].bigd;
}
else {
rvdw[0][inty] = params[iparam_ij].bigr - params[iparam_ij].bigd;
}
rvdw[1][inty] = params[iparam_ij].vsig * 0.950;
// radius check: outer radius vs. sigma
if( rvdw[0][inty] > rvdw[1][inty] )
error->all(FLERR,"Error in vdw spline: inner radius > outer radius");
rrc[0] = rvdw[1][inty];
for (i=1; i<13; i++)
rrc[i] = rrc[i-1] * rrc[0];
drcc = rrc[0] - rvdw[0][inty];
temp6 = 1.0/rrc[5]-1.0/vrc[5]+6.0*(rrc[0]-vrc[0])/vrc[6];
temp7 = 6.0*(1.0/vrc[6]-1.0/rrc[6]);
temp12 = 1.0/rrc[11]-1.0/vrc[11]+(rrc[0]-vrc[0])*12.0/vrc[12];
temp13 = 12.0*(1.0/vrc[12]-1.0/rrc[12]);
vsigt = params[iparam_ij].vsig;
vepst = params[iparam_ij].veps;
vsigt = vsigt*vsigt*vsigt*vsigt*vsigt*vsigt;
vdwt = vepst*(vsigt*vsigt*temp12-vsigt*temp6);
dvdwt = vepst*(vsigt*vsigt*temp13-vsigt*temp7);
cc2[inty] = (3.0/drcc*vdwt-dvdwt)/drcc;
cc3[inty] = (vdwt/(drcc*drcc)-cc2[inty] )/drcc;
}
}
}
// generate vdw look-up table
for (ii = 0; ii < n; ii ++) {
for (jj = ii; jj < n; jj ++) {
itype = ii;
jtype = jj;
inty = intype[itype][jtype];
iparam_ij = elem2param[itype][jtype][jtype];
r = drin;
for (k = 0; k < ncoul; k ++) {
r6 = r*r*r*r*r*r;
r7 = r6 * r;
rf6 = 1.0/r6-1.0/vrc[5]+(r-vrc[0])*6.0/vrc[6];
drf7 = 6.0*(1.0/vrc[6]-1.0/r7);
vsigt = params[iparam_ij].vsig;
vepst = params[iparam_ij].veps;
vsigt = vsigt*vsigt*vsigt*vsigt*vsigt*vsigt;
if(params[iparam_ij].vdwflag>0) {
if(r <= rvdw[0][inty]) {
vvdw[k][inty] = 0.0;
vdvdw[k][inty] = 0.0;
}
else if ( r > rvdw[0][inty] && r <= rvdw[1][inty]) {
drcc = r-rvdw[0][inty];
vvdw[k][inty] = drcc*drcc*(drcc*cc3[inty]+cc2[inty]);
vdvdw[k][inty] = drcc*(3.0*drcc*cc3[inty]+2.0*cc2[inty]);
} else {
r12 = r6*r6;
r13 = r6*r7;
rf12 = 1.0/r12-1.0/vrc[11]+(r-vrc[0])*12.0/vrc[12];
drf13= 12.0*(1.0/vrc[12]-1.0/r13);
vvdw[k][inty] = vepst*(vsigt*vsigt*rf12-vsigt*rf6);
vdvdw[k][inty] = vepst*(vsigt*vsigt*drf13-vsigt*drf7);
}
} else {
vvdw[k][inty]=0.0;
vdvdw[k][inty]=0.0;
}
r += dra;
}
}
}
delete [] vrc;
delete [] rrc;
delete [] cc2;
delete [] cc3;
memory->destroy(rvdw);
}
/* ---------------------------------------------------------------------- */
void PairComb3::potal_calc(double &calc1, double &calc2, double &calc3)
{
double alf,rcoul,esucon;
int m;
rcoul = 0.0;
for (m = 0; m < nparams; m++)
if (params[m].lcut > rcoul) rcoul = params[m].lcut;
alf = 0.20;
esucon = force->qqr2e;
calc2 = (erfc(rcoul*alf)/rcoul/rcoul+2.0*alf/MY_PIS*
exp(-alf*alf*rcoul*rcoul)/rcoul)*esucon/rcoul;
calc3 = (erfc(rcoul*alf)/rcoul)*esucon;
calc1 = -(alf/MY_PIS*esucon+calc3*0.5);
}
/* ---------------------------------------------------------------------- */
void PairComb3::tri_point(double rsq, int &mr1, int &mr2,
int &mr3, double &sr1, double &sr2, double &sr3)
{
double r, rin, dr, dd, rr1, rridr, rridr2;
rin = 0.1000; dr = 0.0010;
r = sqrt(rsq);
if (r < rin + 2.0*dr) r = rin + 2.0*dr;
if (r > cutmax - 2.0*dr) r = cutmax - 2.0*dr;
rridr = (r-rin)/dr;
mr1 = int(rridr) ;
dd = rridr - float(mr1);
if (dd > 0.5) mr1 += 1;
rr1 = float(mr1)*dr;
rridr = (r - rin - rr1)/dr;
rridr2 = rridr * rridr;
sr1 = (rridr2 - rridr) * 0.50;
sr2 = 1.0 - rridr2;
sr3 = (rridr2 + rridr) * 0.50;
mr2 = mr1 + 1;
mr3 = mr1 + 2;
}
/* ---------------------------------------------------------------------- */
void PairComb3::vdwaals(int inty, int mr1, int mr2, int mr3, double rsq,
double sr1, double sr2, double sr3,
double &eng, double &fforce)
{
double r = sqrt(rsq);
eng = 1.0*(sr1*vvdw[mr1-1][inty]+sr2*vvdw[mr2-1][inty]+sr3*vvdw[mr3-1][inty]);
fforce = -1.0/r*(sr1*vdvdw[mr1-1][inty]+sr2*vdvdw[mr2-1][inty]+sr3*vdvdw[mr3-1][inty]);
}
/* ---------------------------------------------------------------------- */
void PairComb3::direct(Param *parami, Param *paramj, int mr1,
int mr2, int mr3, double rsq, double sr1, double sr2, double sr3,
double iq, double jq, double fac11, double fac11e,
double &pot_tmp, double &for_tmp, int i, int j)
{
double r,erfcc,fafbnl,potij,esucon;
double r3,erfcd,dfafbnl,smf2,dvdrr,alf,alfdpi;
double afbn,afbj,sme1n,sme1j,sme1,sme2,dafbn, dafbj,smf1n,smf1j;
double curli = parami->curl;
double curlj = paramj->curl;
int inti = parami->ielement;
int intj = paramj->ielement;
int inty = intype[inti][intj];
double curlij0 = parami->curl0;
double curlji0 = paramj->curl0;
double curlij1,curlji1,dcurlij,dcurlji;
double fcp1j,xcoij,xcoji;
int icurl, jcurl;
int ielegp = parami->ielementgp;
int jelegp = paramj->ielementgp;
r = sqrt(rsq);
r3 = r * rsq;
alf = 0.20;
alfdpi = 2.0*alf/MY_PIS;
esucon = force->qqr2e;
pot_tmp = for_tmp = 0.0;
icurl=jcurl=0;
if(ielegp==2 && curli>curlij0) {
icurl=1;
curlij1=curli;
}
if(jelegp==2 && curlj>curlji0) {
jcurl=1;
curlji1=curlj;
}
if(icurl==1 || jcurl ==1) {
xcoij = xcotmp[i];
xcoji = xcotmp[j];
fcp1j = comb_fc_d(r,parami);
if(icurl==1) {
curli=curlij1+(curlij0-curlij1)*comb_fc_curl(xcoij,parami);
dcurlij=fcp1j*(curlij0-curlij1)*comb_fc_curl_d(xcoij,parami);
}
if(jcurl==1) {
curlj=curlji1+(curlji0-curlji1)*comb_fc_curl(xcoji,paramj);
dcurlji=fcp1j*(curlji0-curlji1)*comb_fc_curl_d(xcoji,paramj);
}
}
erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
afbn = sr1*afb[mr1][inti] + sr2*afb[mr2][inti] + sr3*afb[mr3][inti];
afbj = sr1*afb[mr1][intj] + sr2*afb[mr2][intj] + sr3*afb[mr3][intj];
fafbnl= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
potij = (erfcc/r * esucon - fac11e);
sme1n = iq*curlj*(afbn-fafbnl)*esucon;
sme1j = jq*curli*(afbj-fafbnl)*esucon;
sme1 = sme1n + sme1j;
sme2 = (potij + fafbnl * esucon) * iq * jq;
pot_tmp = 1.0 * (sme1+sme2);
// 1/r force (wrt r)
erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
dafbn = sr1*dafb[mr1][inti] + sr2*dafb[mr2][inti] + sr3*dafb[mr3][inti];
dafbj = sr1*dafb[mr1][intj] + sr2*dafb[mr2][intj] + sr3*dafb[mr3][intj];
dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
smf1n = iq * curlj * (dafbn-dfafbnl)*esucon/r;
smf1j = jq * curli * (dafbj-dfafbnl)*esucon/r;
if(jcurl==1 && ielegp == 3 && dcurlji != 0.0){
smf1n += dcurlji*iq*(afbn-fafbnl)*esucon/r;
}
if(icurl==1 && jelegp == 3 && dcurlij != 0.0){
smf1j += dcurlij*jq*(afbj-fafbnl)*esucon/r;
}
smf2 = dvdrr + dfafbnl * esucon/r;
for_tmp = 1.0 * iq * jq * smf2 + smf1n + smf1j;
}
/* ---------------------------------------------------------------------- */
void PairComb3::field(Param *parami, Param *paramj, double rsq, double iq,
double jq, double &eng_tmp,double &for_tmp)
{
double r,r3,r4,r5,rc,rc2,rc3,rc4,rc5;
double cmi1,cmi2,cmj1,cmj2,pcmi1,pcmi2;
double rf3i,rcf3i,rf5i,rcf5i;
double drf3i,drcf3i,drf5i,drcf5i;
double rf3,rf5,drf4,drf6;
double smpn,smpl,rfx1,rfx2;
r = sqrt(rsq);
r3 = r * r * r;
r4 = r3 * r;
r5 = r4 * r;
rc = parami->lcut;
rc2 = rc * rc;
rc3 = rc*rc*rc;
rc4 = rc3 * rc;
rc5 = rc4 * rc;
cmi1 = parami->cmn1;
cmi2 = parami->cmn2;
cmj1 = paramj->cmn1;
cmj2 = paramj->cmn2;
pcmi1 = parami->pcmn1;
pcmi2 = parami->pcmn2;
rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));
drf3i = 3/r*rf3i-6*rsq*rf3i*rf3i;
drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
drf5i = 5/r*rf5i-10*r4*rf5i*rf5i;
drcf5i = 5/rc*rcf5i-10*rc4*rcf5i*rcf5i;
rf3 = rf3i-rcf3i-(r-rc)*drcf3i;
rf5 = rf5i-rcf5i-(r-rc)*drcf5i;
drf4 = drf3i - drcf3i;
drf6 = drf5i - drcf5i;
// field correction energy
smpn = jq*(cmi1*rf3+jq*cmi2*rf5);
smpl = iq*(cmj1*rf3+iq*cmj2*rf5);
eng_tmp = 1.0 * (smpn + smpl);
// field correction force
rfx1 = jq*(cmi1*drf4+jq*cmi2*drf6)/r;
rfx2 = iq*(cmj1*drf4+iq*cmj2*drf6)/r;
for_tmp -= 1.0 * (rfx1 + rfx2);
}
/* ---------------------------------------------------------------------- */
double PairComb3::rad_init(double rsq2,Param *param,int i,
double &radtot, double cnconj)
{
double r, fc1k, radcut;
r = sqrt(rsq2);
fc1k = comb_fc(r,param);
radtot = -fc1k * param->pcross + cnconj;
radcut = comb_fcch(radtot);
return fc1k * param->pcross * radcut;
}
/* ---------------------------------------------------------------------- */
void PairComb3::rad_calc(double r, Param *parami, Param *paramj,
double kconjug, double lconjug, int i, int j, double xcn, double ycn)
{
int ixmin, iymin, izmin;
int radindx;
double xrad, yrad, zcon, vrad, pradx, prady, pradz;
vrad = pradx = prady = pradz = 0.0;
xrad = -comb_fc(r,parami)*parami->pcross + xcn;
yrad = -comb_fc(r,paramj)*paramj->pcross + ycn;
zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);
if(xrad < 0.0) xrad = 0.0;
if(yrad < 0.0) yrad = 0.0;
if(zcon < 1.0) zcon = 1.0;
if(xrad > maxxc) xrad = maxxc;
if(yrad > maxyc) yrad = maxyc;
if(zcon > maxconj) zcon = maxconj;
ixmin = int(xrad+1.0e-12);
iymin = int(yrad+1.0e-12);
izmin = int(zcon+1.0e-12);
radindx=parami->rad_flag-1;
if (fabs(float(ixmin)-xrad)>1.0e-8 ||
fabs(float(iymin)-yrad)>1.0e-8 ||
fabs(float(izmin)-zcon)>1.0e-8) {
rad_int(radindx,xrad,yrad,zcon,ixmin,iymin,izmin,
vrad,pradx,prady,pradz);
} else {
vrad = rad_grid[radindx][ixmin][iymin][izmin-1];
pradx = rad_gridx[radindx][ixmin][iymin][izmin-1];
prady = rad_gridy[radindx][ixmin][iymin][izmin-1];
pradz = rad_gridz[radindx][ixmin][iymin][izmin-1];
}
brad[0] = vrad;
brad[1] = pradx;
brad[2] = prady;
brad[3] = pradz;
}
/* ---------------------------------------------------------------------- */
void PairComb3::rad_int(int radindx,double xrad, double yrad, double zcon, int l,
int m, int n, double &vrad, double &pradx, double &prady,
double &pradz)
{
int j;
double x;
vrad = pradx = prady = pradz = 0.0;
if(l >= maxxc-1) { l=maxxc-1;}
if(m >= maxyc-1) { m=maxyc-1; }
if(n >= maxconj-1) { n=maxconj-1;}
for (j=0; j<64; j++) {
x = rad_spl[radindx][l][m][n-1][j] * pow(xrad,iin3[j][0])
* pow(yrad,iin3[j][1]) * pow(zcon,iin3[j][2]);
vrad += x;
if(xrad > 1.0e-8) pradx += x*iin3[j][0]/xrad;
if(yrad > 1.0e-8) prady += x*iin3[j][1]/yrad;
if(zcon > 1.0e-8) pradz += x*iin3[j][2]/zcon;
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::rad_forceik(Param *paramk, double rsq2, double *delrk,
double conjug, double radtot)
{
int nm;
double rik, fc1k, fcp1k;
double pradk, ffkk2, fktmp[3];
double radcut = comb_fcch(radtot);
double dradcut = comb_fcch_d(radtot);
for (nm=0; nm<3; nm++) {
fi_rad[nm] = fk_rad[nm] = 0.0;
}
radtmp =0.0;
rik = sqrt(rsq2);
fc1k = comb_fc(rik, paramk);
fcp1k = comb_fc_d(rik,paramk);
pradk = brad[3]*fcp1k*radcut*paramk->pcross*2.0*conjug;
radtmp= brad[3]*fc1k*dradcut*paramk->pcross*2.0*conjug;
ffkk2 = -pradk/rik;
for (nm=0; nm<3; nm++) {
fktmp[nm] = - ffkk2 * delrk[nm];
}
for (nm=0; nm<3; nm++) {
fi_rad[nm] = fktmp[nm];
fk_rad[nm] = -fktmp[nm];
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::rad_force(Param *paramm, double rsq3,
double *delrm, double dpradk)
{
int nm;
double rkm, fcp1m;
double comkm, ffmm2, fkm[3];
for (nm=0; nm<3; nm++) {
fj_rad[nm] = fk_rad[nm] = 0.0;
fkm[nm]=0.0;
}
rkm = sqrt(rsq3);
fcp1m = comb_fc_d(rkm, paramm);
comkm = dpradk * fcp1m * paramm->pcross;
ffmm2 = -comkm/rkm;
for (nm=0; nm<3; nm++) {
fkm[nm] = -ffmm2 * delrm[nm];
}
for (nm=0; nm<3; nm++) {
fj_rad[nm] = fkm[nm];
fk_rad[nm] = -fkm[nm];
}
}
/* ---------------------------------------------------------------------- */
double PairComb3::bbtor1(int torindx, Param *paramk, Param *paraml,
double rsq1, double rsq2, double rsq3, double *delrj,
double *delrk, double *delrl, double srmu)
{
double rmul, rij, rik, rjl;
rij = sqrt(rsq1);
rik = sqrt(rsq2);
rjl = sqrt(rsq3);
vec3_scale(-1.0,delrl,delrl);
rmul = vec3_dot(delrj,delrl)/(rij*rjl);
vec3_scale(-1.0,delrl,delrl);
rmul = sqrt(1.0-rmul*rmul);
if(rmul > 0.1 ) {
double fc1k, fc1l, TT1, TT2, rmut, btt, tork[3], torl[3];
fc1k = comb_fc(rik,paramk);
fc1l = comb_fc(rjl,paraml);
TT1 = rik*rjl*rij*rij*srmu*rmul;
tork[0] = delrk[1]*delrj[2] - delrk[2]*delrj[1];
torl[0] = delrj[1]*delrl[2] - delrj[2]*delrl[1];
tork[1] = delrk[2]*delrj[0] - delrk[0]*delrj[2];
torl[1] = delrj[2]*delrl[0] - delrj[0]*delrl[2];
tork[2] = delrk[0]*delrj[1] - delrk[1]*delrj[0];
torl[2] = delrj[0]*delrl[1] - delrj[1]*delrl[0];
TT2 = vec3_dot(tork,torl);
rmut = pow((TT2/TT1),2);
if(torindx>=1) {
btt = 1.0 - rmut;
return btt * fc1k * fc1l;
}
else {
btt=paramk->ptork1-TT2/TT1;
btt=paramk->ptork2*pow(btt,2);
return btt * fc1k * fc1l;
}
} else {
return 0.0;
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::tor_calc(double r, Param *parami, Param *paramj,
double kconjug, double lconjug, int i, int j, double xcn, double ycn)
{
int ixmin, iymin, izmin;
double vtor, dtorx, dtory, dtorz;
double xtor, ytor, zcon;
int torindx;
vtor = dtorx = dtory = dtorz = 0.0;
torindx=parami->tor_flag;
if(torindx<0){
vtor=1.0;
dtorx=0.0;
dtory=0.0;
dtorz=0.0;
} else {
xtor = -comb_fc(r, parami) * parami->pcross + xcn;
ytor = -comb_fc(r, paramj) * paramj->pcross + ycn;
zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);
if (xtor < 0.0) xtor = 0.0;
if (ytor < 0.0) ytor = 0.0;
if (zcon < 1.0) zcon = 1.0;
if (xtor > maxxc) xtor = maxxc;
if (ytor > maxyc) ytor = maxyc;
if (zcon > maxconj) zcon = maxconj;
ixmin = int(xtor+1.0e-12);
iymin = int(ytor+1.0e-12);
izmin = int(zcon+1.0e-12);
torindx=torindx-1;
if (fabs(float(ixmin)-xtor)>1.0e-8 ||
fabs(float(iymin)-ytor)>1.0e-8 ||
fabs(float(izmin)-zcon)>1.0e-8) {
tor_int(torindx,xtor,ytor,zcon,ixmin,iymin,izmin,
vtor,dtorx,dtory,dtorz);
} else {
vtor = tor_grid[torindx][ixmin][iymin][izmin-1];
dtorx = tor_gridx[torindx][ixmin][iymin][izmin-1];
dtory = tor_gridy[torindx][ixmin][iymin][izmin-1];
dtorz = tor_gridz[torindx][ixmin][iymin][izmin-1];
}
}
btor[0] = vtor;
btor[1] = dtorx;
btor[2] = dtory;
btor[3] = dtorz;
}
/* ---------------------------------------------------------------------- */
void PairComb3::tor_int(int torindx,double xtor, double ytor, double zcon, int l,
int m, int n, double &vtor, double &dtorx, double &dtory, double &dtorz)
{
int j;
double x;
vtor = dtorx = dtory = dtorz = 0.0;
if(l >= maxxc-1) { l=maxxc-1; } //boundary condition changed
if(m >= maxyc-1) { m=maxyc-1; }
if(n >= maxconj-1) { n=maxconj-1; }
for (j=0; j<64; j++) {
x = tor_spl[torindx][l][m][n-1][j] * pow(xtor,iin3[j][0])
* pow(ytor,iin3[j][1]) * pow(zcon,iin3[j][2]);
vtor += x;
if(xtor > 1.0e-8 ) dtorx += x*iin3[j][0]/xtor;
if(ytor > 1.0e-8 ) dtory += x*iin3[j][1]/ytor;
if(zcon > 1.0e-8 ) dtorz += x*iin3[j][2]/zcon;
}
}
/* ---------------------------------------------------------------------- */
void PairComb3::tor_force(int torindx, Param *paramk, Param *paraml,
double srmu, double rsq1,double rsq2, double rsq3,
double *delrj, double *delrk, double *delrl)
{
int nm;
double rmu, rmul, srmul, rij, rik, rjl;
for (nm=0; nm<3; nm++) {
fi_tor[nm] = fj_tor[nm] = fk_tor[nm] = fl_tor[nm] = 0.0;
}
rij = sqrt(rsq1);
rik = sqrt(rsq2);
rjl = sqrt(rsq3);
rmu = vec3_dot(delrj,delrk)/(rij*rik);
vec3_scale(-1.0,delrl,delrl);
rmul = vec3_dot(delrj,delrl)/(rij*rjl);
vec3_scale(-1.0,delrl,delrl);
srmul = sqrt(1.0-rmul*rmul);
if(acos(rmul) > MY_PI) srmul = -srmul;
if(srmul > 0.1 ) {
double fc1k, fcp1k, fc1l, fcp1l, srmul2, dt1dik, dt1djl;
double TT1, TT2, rmut, btt, tork[3], torl[3];
double dt2dik[3], dt2djl[3], dt2dij[3], AA, AA2;
double tfij[4], tfik[2], tfjl[2], tjx[3], tjy[3], tjz[3];
double tkx[2], tky[2], tkz[2], tlx[2], tly[2], tlz[2];
fc1k = comb_fc(rik,paramk);
fcp1k = comb_fc_d(rik,paramk);
fc1l = comb_fc(rjl,paraml);
fcp1l = comb_fc_d(rjl,paraml);
srmul2 = pow(srmul,2);
TT1 = rik*rjl*rij*rij*srmu*srmul;
dt1dik = -rmu/pow(srmu,2);
dt1djl = -rmul/srmul2;
tork[0] = delrk[1]*delrj[2] - delrk[2]*delrj[1];
torl[0] = delrj[1]*delrl[2] - delrj[2]*delrl[1];
tork[1] = delrk[2]*delrj[0] - delrk[0]*delrj[2];
torl[1] = delrj[2]*delrl[0] - delrj[0]*delrl[2];
tork[2] = delrk[0]*delrj[1] - delrk[1]*delrj[0];
torl[2] = delrj[0]*delrl[1] - delrj[1]*delrl[0];
TT2 = vec3_dot(tork,torl);
dt2dik[0] = -delrj[1]*torl[2] + delrj[2]*torl[1];
dt2dik[1] = -delrj[2]*torl[0] + delrj[0]*torl[2];
dt2dik[2] = -delrj[0]*torl[1] + delrj[1]*torl[0];
dt2djl[0] = delrj[1]*tork[2] - delrj[2]*tork[1];
dt2djl[1] = delrj[2]*tork[0] - delrj[0]*tork[2];
dt2djl[2] = delrj[0]*tork[1] - delrj[1]*tork[0];
dt2dij[0] = -delrk[2]*torl[1] + delrl[2]*tork[1]
+ delrk[1]*torl[2] - delrl[1]*tork[2];
dt2dij[1] = -delrk[0]*torl[2] + delrl[0]*tork[2]
+ delrk[2]*torl[0] - delrl[2]*tork[0];
dt2dij[2] = -delrk[1]*torl[0] + delrl[1]*tork[0]
+ delrk[0]*torl[1] - delrl[0]*tork[1];
rmut = TT2/TT1;
if(torindx>=1) {
btt = 1.0 - pow(rmut,2);
AA = -2.0 * ptorr * rmut * fc1k * fc1l / TT1;
}
else {
btt=paramk->ptork1-rmut;
btt=paramk->ptork2*pow(btt,2);
AA = -2.0 * ptorr * paramk->ptork2 *
(paramk->ptork1-rmut) * fc1k * fc1l /TT1;
}
AA2 = AA * TT2;
tfij[0] = -(dt1dik*AA2)/rij/rik;
tfij[1] = AA2/rij/rij - dt1dik*AA2*rmu/rij/rij;
tfij[2] = -dt1djl*AA2/rij/rjl;
tfij[3] = AA2/rij/rij - dt1djl*AA2*rmul/rij/rij;
tfik[0] = tfij[0];
tfik[1] = (AA2/rik - btt*ptorr*fc1l*fcp1k)/rik -
dt1dik*AA2*rmu/rik/rik;
tfjl[0] = tfij[2];
tfjl[1] = (AA2/rjl - btt*ptorr*fc1k*fcp1l)/rjl -
dt1djl*AA2*rmul/rjl/rjl;
tjx[0] = tfij[0]*delrk[0] - tfij[1]*delrj[0];
tjy[0] = tfij[0]*delrk[1] - tfij[1]*delrj[1];
tjz[0] = tfij[0]*delrk[2] - tfij[1]*delrj[2];
tjx[1] = -tfij[2]*delrl[0] - tfij[3]*delrj[0];
tjy[1] = -tfij[2]*delrl[1] - tfij[3]*delrj[1];
tjz[1] = -tfij[2]*delrl[2] - tfij[3]*delrj[2];
tjx[2] = -dt2dij[0] * AA;
tjy[2] = -dt2dij[1] * AA;
tjz[2] = -dt2dij[2] * AA;
tkx[0] = tfik[0]*delrj[0] - tfik[1]*delrk[0];
tky[0] = tfik[0]*delrj[1] - tfik[1]*delrk[1];
tkz[0] = tfik[0]*delrj[2] - tfik[1]*delrk[2];
tkx[1] = -dt2dik[0] * AA;
tky[1] = -dt2dik[1] * AA;
tkz[1] = -dt2dik[2] * AA;
tlx[0] = -tfjl[0]*delrj[0] - tfjl[1]*delrl[0];
tly[0] = -tfjl[0]*delrj[1] - tfjl[1]*delrl[1];
tlz[0] = -tfjl[0]*delrj[2] - tfjl[1]*delrl[2];
tlx[1] = -dt2djl[0] * AA;
tly[1] = -dt2djl[1] * AA;
tlz[1] = -dt2djl[2] * AA;
fi_tor[0] = tjx[0]+tjx[1]+tjx[2]+tkx[0]+tkx[1];
fi_tor[1] = tjy[0]+tjy[1]+tjy[2]+tky[0]+tky[1];
fi_tor[2] = tjz[0]+tjz[1]+tjz[2]+tkz[0]+tkz[1];
fj_tor[0] = -tjx[0]-tjx[1]-tjx[2]+tlx[0]+tlx[1];
fj_tor[1] = -tjy[0]-tjy[1]-tjy[2]+tly[0]+tly[1];
fj_tor[2] = -tjz[0]-tjz[1]-tjz[2]+tlz[0]+tlz[1];
fk_tor[0] = -tkx[0]-tkx[1];
fk_tor[1] = -tky[0]-tky[1];
fk_tor[2] = -tkz[0]-tkz[1];
fl_tor[0] = -tlx[0]-tlx[1];
fl_tor[1] = -tly[0]-tly[1];
fl_tor[2] = -tlz[0]-tlz[1];
}
}
/* ---------------------------------------------------------------------- */
double PairComb3::combqeq(double *qf_fix, int &igroup)
{
int i,j,ii, jj,itype,jtype,jnum;
int iparam_i,iparam_ji,iparam_ij;
int *ilist,*jlist,*numneigh,**firstneigh;
int mr1,mr2,mr3,inty,nj;
double xtmp,ytmp,ztmp,rsq1,delrj[3];
double iq,jq,fqi,fqij,fqji,sr1,sr2,sr3;
double potal,fac11,fac11e;
int sht_jnum,*sht_jlist;
tagint itag, jtag;
double **x = atom->x;
double *q = atom->q;
tagint *tag = atom->tag;
int *type = atom->type;
int inum = list->inum;
int *mask = atom->mask;
int groupbit = group->bitmask[igroup];
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
qf = qf_fix;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
qf[i] = 0.0;
dpl[i][0] = dpl[i][1] = dpl[i][2] = 0.0;
}
}
// communicating charge force to all nodes, first forward then reverse
pack_flag = 1;
comm->forward_comm_pair(this);
// self energy correction term: potal
potal_calc(potal,fac11,fac11e);
// loop over full neighbor list of my atoms
fqi = fqij = fqji = 0.0;
for (ii = 0; ii < inum; ii ++) {
i = ilist[ii];
itag = tag[i];
nj = 0;
if (mask[i] & groupbit) {
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
iq = q[i];
iparam_i = elem2param[itype][itype][itype];
// charge force from self energy
fqi = qfo_self(&params[iparam_i],iq);
jlist = firstneigh[i];
jnum = numneigh[i];
sht_jlist = sht_first[i];
sht_jnum = sht_num[i];
// two-body interactions
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj] & NEIGHMASK;
jtag = tag[j];
if (itag >= jtag) continue;
jtype = map[type[j]];
inty = intype[itype][jtype];
jq = q[j];
delrj[0] = xtmp - x[j][0];
delrj[1] = ytmp - x[j][1];
delrj[2] = ztmp - x[j][2];
rsq1 = vec3_dot(delrj,delrj);
iparam_ij = elem2param[itype][jtype][jtype];
iparam_ji = elem2param[jtype][itype][itype];
// long range q-dependent
if (rsq1 > params[iparam_ij].lcutsq) continue;
// polynomial three-point interpolation
tri_point(rsq1,mr1,mr2,mr3,sr1,sr2,sr3);
// 1/r charge forces
qfo_direct(&params[iparam_ij],&params[iparam_ji],
mr1,mr2,mr3,rsq1,sr1,sr2,sr3,fac11e,fqij,fqji,
iq,jq,i,j);
fqi += fqij; qf[j] += fqji;
// field correction to self energy and charge force
qfo_field(&params[iparam_ij],&params[iparam_ji],rsq1,
iq,jq,fqij,fqji);
fqi += fqij; qf[j] += fqji;
// polarization field charge force
if (pol_flag) {
qfo_dipole(fac11,mr1,mr2,mr3,inty,rsq1,delrj,sr1,sr2,sr3,
fqij,fqji,i,j);
fqi += fqij; qf[j] += fqji;
}
}
for (jj = 0; jj < sht_jnum; jj++) {
j = sht_jlist[jj];
jtag = tag[j];
if (itag >= jtag) continue;
jtype = map[type[j]];
inty = intype[itype][jtype];
jq = q[j];
delrj[0] = xtmp - x[j][0];
delrj[1] = ytmp - x[j][1];
delrj[2] = ztmp - x[j][2];
rsq1 = vec3_dot(delrj,delrj);
iparam_ij = elem2param[itype][jtype][jtype];
iparam_ji = elem2param[jtype][itype][itype];
if (rsq1 >= params[iparam_ij].cutsq) continue;
nj ++;
// charge force in Aij and Bij
qfo_short(&params[iparam_ij],&params[iparam_ji],
rsq1,iq,jq,fqij,fqji,i,j,nj);
fqi += fqij; qf[j] += fqji;
}
qf[i] += fqi;
}
}
comm->reverse_comm_pair(this);
// sum charge force on each node and return it
double eneg = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit){
eneg += qf[i];
itag=tag[i];
}
}
MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Bcast(&enegtot,1,MPI_DOUBLE,0,world);
return enegtot;
}
/* ---------------------------------------------------------------------- */
double PairComb3::qfo_self(Param *param, double qi)
{
double self_d,cmin,cmax,qmin,qmax;
double s1 = param->chi;
double s2 = param->dj;
double s3 = param->dk;
double s4 = param->dl;
self_d = 0.0;
qmin = param->qmin;
qmax = param->qmax;
cmin = cmax = 100.0;
self_d = s1+qi*(2.0*s2+qi*(3.0*s3+qi*4.0*s4));
if (qi < qmin) self_d += 4.0 * cmin * pow((qi-qmin),3);
if (qi > qmax) self_d += 4.0 * cmax * pow((qi-qmax),3);
return self_d;
}
/* ---------------------------------------------------------------------- */
void PairComb3::qfo_direct(Param *parami, Param *paramj, int mr1,
int mr2, int mr3, double rsq, double sr1, double sr2,
double sr3, double fac11e, double &fqij, double &fqji,
double iq, double jq, int i, int j)
{
double r, erfcc, fafbnl, vm, vmfafb, esucon;
double afbn, afbj, sme1n, sme1j;
double curli = parami->curl;
double curlj = paramj->curl;
int inti = parami->ielement;
int intj = paramj->ielement;
int inty = intype[inti][intj];
double curlij0 = parami->curl0;
double curlji0 = paramj->curl0;
double curlij1,curlji1;
int icurl, jcurl;
int ielegp = parami->ielementgp;
int jelegp = paramj->ielementgp;
r = sqrt(rsq);
esucon=force->qqr2e;
icurl = jcurl = 0;
if(ielegp==2 && curli>curlij0) {
icurl=1;
curlij1=curli;
}
if(jelegp==2 && curlj>curlji0) {
jcurl=1;
curlji1=curlj;
}
if(icurl==1 || jcurl ==1) {
double xcoij= xcotmp[i];
double xcoji= xcotmp[j];
if(icurl==1) {
curli=curlij1+(curlij0-curlij1)*comb_fc_curl(xcoij,parami);
}
if(jcurl==1) {
curlj=curlji1+(curlji0-curlji1)*comb_fc_curl(xcoji,paramj);
}
}
// 1/r force (wrt q)
erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
fafbnl= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
afbn = sr1*afb[mr1][inti] + sr2*afb[mr2][inti] + sr3*afb[mr3][inti];
afbj = sr1*afb[mr1][intj] + sr2*afb[mr2][intj] + sr3*afb[mr3][intj];
vm = (erfcc/r * esucon - fac11e);
vmfafb = vm + esucon * fafbnl;
sme1n = curlj * (afbn - fafbnl) * esucon;
sme1j = curli * (afbj - fafbnl) * esucon;
fqij = 1.0 * (jq * vmfafb + sme1n);
fqji = 1.0 * (iq * vmfafb + sme1j);
}
/* ---------------------------------------------------------------------- */
void PairComb3::qfo_field(Param *parami, Param *paramj, double rsq,
double iq,double jq, double &fqij, double &fqji)
{
double r,r3,r5,rc,rc2,rc3,rc4,rc5;
double cmi1,cmi2,cmj1,cmj2,pcmi1,pcmi2;
double rf3i,rcf3i,rf5i,rcf5i;
double drcf3i,drcf5i,rf3,rf5;
r = sqrt(rsq);
r3 = r * rsq;
r5 = r3 * rsq;
rc = parami->lcut;
rc2= rc*rc;
rc3 = rc*rc*rc;
rc4 = rc3 * rc;
rc5 = rc4 * rc;
cmi1 = parami->cmn1;
cmi2 = parami->cmn2;
cmj1 = paramj->cmn1;
cmj2 = paramj->cmn2;
pcmi1 = parami->pcmn1;
pcmi2 = parami->pcmn2;
rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));
drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
drcf5i = 5/rc*rcf5i-10*rc4*rcf5i*rcf5i;
rf3 = rf3i-rcf3i-(r-rc)*drcf3i;
rf5 = rf5i-rcf5i-(r-rc)*drcf5i;
// field correction charge force
fqij = 1.0 * cmj1*rf3+2.0*iq*cmj2*rf5;
fqji = 1.0 * cmi1*rf3+2.0*jq*cmi2*rf5;
}
/* ---------------------------------------------------------------------- */
void PairComb3::qfo_dipole(double fac11, int mr1, int mr2, int mr3,
int inty, double rsq, double *delrj, double sr1, double sr2,
double sr3, double &fqij, double &fqji, int i, int j)
{
double erfcc, erfcd, dvdrr, dfafbnl, smf2;
double r, r3, alfdpi, esucon;
r = sqrt(rsq);
r3 = r * rsq;
alfdpi = 0.4/MY_PIS;
esucon = force->qqr2e;
erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
smf2 = (dvdrr + dfafbnl*esucon)/r;
fqij = dpl[i][0]*delrj[0] + dpl[i][1]*delrj[1] +dpl[i][2]*delrj[2];
fqji = dpl[j][0]*delrj[0] + dpl[j][1]*delrj[1] +dpl[j][2]*delrj[2];
fqij *= smf2;
fqji *= smf2;
}
/* ---------------------------------------------------------------------- */
void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
double iq, double jq, double &fqij, double &fqji,
int i, int j, int nj)
{
double r, tmp_fc;
double Di, Dj, dDi, dDj, Bsi, Bsj, dBsi, dBsj;
double QUchi, QOchi, QUchj, QOchj;
double bij, caj, cbj, caqpn, caqpj, cbqpn, cbqpj;
double LamDiLamDj, AlfDiAlfDj;
double rlm1 = parami->lambda;
double alfij1= parami->alpha1;
double alfij2= parami->alpha2;
double alfij3= parami->alpha3;
double pbij1= parami->bigB1;
double pbij2= parami->bigB2;
double pbij3= parami->bigB3;
caj = cbj = caqpn = caqpj = cbqpn = cbqpj = 0.0;
r = sqrt(rsq);
tmp_fc = comb_fc(r,parami);
bij = bbij[i][nj];
QUchi = (parami->QU - iq) * parami->bD;
QUchj = (paramj->QU - jq) * paramj->bD;
QOchi = (iq - parami->Qo) * parami->bB;
QOchj = (jq - paramj->Qo) * paramj->bB;
if (iq < parami->QL-0.2) {
iq = parami->QL-0.2;
Di = parami->DL;
dDi = Bsi = dBsi = 0.0;
} else if (iq > parami->QU+0.2) {
iq = parami->QU+0.2;
Di = parami->DU;
dDi = Bsi = dBsi = 0.0;
} else {
Di = parami->DU + pow(QUchi,parami->nD); // YYDin
dDi = -parami->nD * parami->bD * pow(QUchi,(parami->nD-1.0)); // YYDiqp
Bsi = parami->aB - pow(QOchi,10); // YYBsin
dBsi = -parami->bB * 10.0 * pow(QOchi,9.0); // YYBsiqp
}
if (jq < paramj->QL-0.2) {
jq = paramj->QL-0.2;
Dj = paramj->DL;
dDj = Bsj = dBsj = 0.0;
} else if (jq > paramj->QU+0.2) {
jq = paramj->QU+0.2;
Dj = paramj->DU;
dDj = Bsj = dBsj = 0.0;
} else {
Dj = paramj->DU + pow(QUchj,paramj->nD); // YYDij
dDj = -paramj->nD * paramj->bD * pow(QUchj,(paramj->nD-1.0)); // YYDiqpj
Bsj = paramj->aB - pow(QOchj,10); // YYBsij
dBsj = -paramj->bB * 10.0 * pow(QOchj,9.0); // YYBsiqpj
}
LamDiLamDj = exp(0.5*(parami->lami*Di+paramj->lami*Dj)-rlm1*r);
caj = 0.5 * tmp_fc * parami->bigA * LamDiLamDj;
if (Bsi*Bsj > 0.0) {
AlfDiAlfDj = exp(0.5*(parami->alfi*Di+paramj->alfi*Dj));
cbj=-0.5*tmp_fc*bij*sqrt(Bsi*Bsj)*AlfDiAlfDj*
(pbij1*exp(-alfij1*r)+pbij2*exp(-alfij2*r)+pbij3*exp(-alfij3*r));
cbqpn = cbj * (parami->alfi * dDi + dBsi/Bsi);
cbqpj = cbj * (paramj->alfi * dDj + dBsj/Bsj);
} else {
cbj = cbqpn = cbqpj = 0.0;
}
caqpn = caj * parami->lami * dDi;
caqpj = caj * paramj->lami * dDj;
fqij = 1.0 * (caqpn + cbqpn);
fqji = 1.0 * (caqpj + cbqpj);
}
/* ---------------------------------------------------------------------- */
void PairComb3::dipole_init(Param *parami, Param *paramj, double fac11,
double *delrj, double rsq, int mr1, int mr2, int mr3, double sr1,
double sr2, double sr3, double iq, double jq, int i, int j)
{
double erfcc, erfcd, dvdrr, dfafbnl, smf2, phinn, phinj, efn, efj;
double r, r3, alfdpi, esucon;
double rcd, rct, tmurn, tmurj, poln[3], polj[3], Qext[3];
int nm;
int inti = parami->ielement;
int intj = paramj->ielement;
int inty = intype[inti][intj];
for(nm=0; nm<3; nm++) Qext[nm] = 0.0;
r = sqrt(rsq);
r3 = r * rsq;
rcd = 1.0/(r3);
rct = 3.0*rcd/rsq;
alfdpi = 0.4/MY_PIS;
esucon = force->qqr2e;
erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
smf2 = dvdrr/esucon + dfafbnl/r;
phinn = sr1*phin[mr1][inti] + sr2*phin[mr2][inti] + sr3*phin[mr3][inti];
phinj = sr1*phin[mr1][intj] + sr2*phin[mr2][intj] + sr3*phin[mr3][intj];
efn = jq * smf2;
efj = iq * smf2;
tmurn = dpl[i][0]*delrj[0] + dpl[i][1]*delrj[1] + dpl[i][2]*delrj[2];
tmurj = dpl[j][0]*delrj[0] + dpl[j][1]*delrj[1] + dpl[j][2]*delrj[2];
for (nm=0; nm<3; nm++) {
poln[nm] = (tmurj*delrj[nm]*rct - dpl[j][nm]*rcd)*phinj;
polj[nm] = (tmurn*delrj[nm]*rct - dpl[i][nm]*rcd)*phinn;
}
for (nm=0; nm<3; nm++) {
dpl[i][nm] += (Qext[nm]/esucon + delrj[nm]*efn + poln[nm])*parami->polz*0.50;
dpl[j][nm] += (Qext[nm]/esucon - delrj[nm]*efj + polj[nm])*paramj->polz*0.50;
}
}
/* ---------------------------------------------------------------------- */
double PairComb3::dipole_self(Param *parami, int i)
{
double esucon = force->qqr2e;
double apn = parami->polz;
double selfdpV = 0.0;
if (apn != 0.0) {
selfdpV= (dpl[i][0]*dpl[i][0]+dpl[i][1]*dpl[i][1]+dpl[i][2]*dpl[i][2])
*esucon/(2.0*apn); }
return selfdpV;
}
/* ---------------------------------------------------------------------- */
void PairComb3::dipole_calc(Param *parami, Param *paramj, double fac11,
double delx, double dely, double delz, double rsq,
int mr1, int mr2, int mr3, double sr1, double sr2, double sr3,
double iq, double jq, int i, int j, double &vionij,
double &fvionij, double *ddprx)
{
double erfcc, erfcd, dvdrr, dfafbnl, ef, phinn, phinj, efn, efj;
double r, r3, alf, alfdpi, esucon, dphinn, dphinj, ddfafbnl;
double def, defn, defj, tmun, tmuj, emuTmu, edqn, edqj, ddvdrr;
double rcd, rct, tmurn, tmurj, tmumu, poln[3], polj[3], delr1[3];
double demuTmu, ddpr, dcoef;
int nm;
int inti = parami->ielement;
int intj = paramj->ielement;
int inty = intype[inti][intj];
r = sqrt(rsq);
r3 = r * rsq;
esucon = force->qqr2e;
rcd = esucon/r3;
rct = 3.0*rcd/rsq;
alf = 0.2;
alfdpi = 2.0*alf/MY_PIS;
delr1[0] = delx;
delr1[1] = dely;
delr1[2] = delz;
// generate energy & force information from tables
erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
ddvdrr = (2.0*erfcc/r3 + 2.0*alfdpi*erfcd*(1.0/rsq+alf*alf))*esucon;
dfafbnl= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
phinn = sr1*phin[mr1][inti] + sr2*phin[mr2][inti] + sr3*phin[mr3][inti];
phinj = sr1*phin[mr1][intj] + sr2*phin[mr2][intj] + sr3*phin[mr3][intj];
dphinn = sr1*dphin[mr1][inti] + sr2*dphin[mr2][inti] + sr3*dphin[mr3][inti];
dphinj = sr1*dphin[mr1][intj] + sr2*dphin[mr2][intj] + sr3*dphin[mr3][intj];
ddfafbnl= sr1*ddfafb[mr1][inty] + sr2*ddfafb[mr2][inty] + sr3*ddfafb[mr3][inty];
ef = (dvdrr + dfafbnl * esucon)/r;
efn = jq * ef;
efj = -iq * ef;
def = (ddvdrr + ddfafbnl * esucon)/r;
defn = jq * def;
defj = -iq * def;
// dipole - dipole field tensor (Tij)
tmurn = dpl[i][0]*delr1[0] + dpl[i][1]*delr1[1] + dpl[i][2]*delr1[2];
tmurj = dpl[j][0]*delr1[0] + dpl[j][1]*delr1[1] + dpl[j][2]*delr1[2];
tmumu = dpl[i][0]*dpl[j][0] + dpl[i][1]*dpl[j][1] + dpl[i][2]*dpl[j][2];
for (nm=0; nm<3; nm++) {
poln[nm] = (tmurj*delr1[nm]*rct - dpl[j][nm]*rcd);
polj[nm] = (tmurn*delr1[nm]*rct - dpl[i][nm]*rcd);
}
tmun = dpl[j][0]*polj[0] + dpl[j][1]*polj[1] + dpl[j][2]*polj[2];
tmuj = dpl[i][0]*poln[0] + dpl[i][1]*poln[1] + dpl[i][2]*poln[2];
// dipole - dipole energy
emuTmu = -0.5*(tmun*phinn+tmuj*phinj);
// dipole - charge energy
edqn = -0.5 * (tmurn * efn);
edqj = -0.5 * (tmurj * efj);
// overall dipole energy
vionij = emuTmu + edqn + edqj;
// dipole - dipole force
demuTmu = (tmun*dphinn + tmuj*dphinj)/r;
ddpr = 5.0*tmurn*tmurj/rsq - tmumu;
dcoef = rct * (phinn+phinj);
for (nm = 0; nm < 3; nm ++) {
ddprx[nm] = dcoef * (ddpr*delr1[nm] - tmurn*dpl[j][nm] - tmurj*dpl[i][nm])
+ demuTmu * delr1[nm];
}
// dipole - charge force
fvionij = -tmurn*defn - tmurj*defj;
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fc_curl(double rocn, Param *param)
{
double r_inn = param->curlcut1;
double r_out = param->curlcut2;
if (rocn <= r_inn) return 1.0;
if (rocn >= r_out) return 0.0;
return 0.5*(1.0 + cos(MY_PI*(rocn-r_inn)/(r_out-r_inn)));
}
/* ---------------------------------------------------------------------- */
double PairComb3::comb_fc_curl_d(double rocn, Param *param)
{
double r_inn = param->curlcut1;
double r_out = param->curlcut2;
if (rocn <= r_inn) return 0.0;
if (rocn >= r_out) return 0.0;
return -MY_PI2/(r_out-r_inn)*sin(MY_PI*(rocn-r_inn)/(r_out-r_inn));
}
/* ---------------------------------------------------------------------- */
int PairComb3::heaviside(double rr)
{
if (rr <= 0.0) return 0;
else return 1;
}
/* ---------------------------------------------------------------------- */
double PairComb3::switching(double rr)
{
if (rr <= 0.0) return 1.0;
else if (rr >= 1.0) return 0.0;
else return heaviside(-rr)+heaviside(rr)*heaviside(1.0-rr)
* (1.0-(3.0-2.0*rr)*rr*rr);
}
/* ---------------------------------------------------------------------- */
double PairComb3::switching_d(double rr)
{
if (rr <= 0.0) return 0.0;
else if (rr >= 1.0) return 0.0;
else return heaviside(rr)*heaviside(1.0-rr)
* 6.0*rr*(rr-1.0);
}
/* ---------------------------------------------------------------------- */
int PairComb3::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
if (pack_flag == 1) {
for (i = 0; i < n; i ++) {
j = list[i];
buf[m++] = qf[j];
}
} else if (pack_flag == 2) {
for (i = 0; i < n; i ++) {
j = list[i];
buf[m++] = NCo[j];
}
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairComb3::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n ;
if (pack_flag == 1) {
for (i = first; i < last; i++)
qf[i] = buf[m++];
} else if (pack_flag == 2) {
for (i = first; i < last; i++)
NCo[i] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int PairComb3::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (pack_flag == 1) {
for (i = first; i < last; i++)
buf[m++] = qf[i];
} else if (pack_flag == 2) {
for (i = first; i < last; i++)
buf[m++] = NCo[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairComb3::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
if (pack_flag == 1) {
for (i = 0; i < n; i++) {
j = list[i];
qf[j] += buf[m++];
}
} else if (pack_flag == 2) {
for (i = 0; i < n; i++) {
j = list[i];
NCo[j] += buf[m++];
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairComb3::memory_usage()
{
double bytes = maxeatom * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
bytes += nmax * sizeof(int);
bytes += nmax * 8.0 * sizeof(double);
bytes += 25000*2*sizeof(double);
for (int i = 0; i < comm->nthreads; i++)
bytes += ipage[i].size();
return bytes;
}
c4science · Help