forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11440 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
8031ade7de
commit
7bf0ad314c
|
|
@ -323,7 +323,7 @@ void PairComb3::read_lib()
|
|||
FILE *fp = fopen("lib.comb3","r");
|
||||
if (fp == NULL) {
|
||||
char str[128];
|
||||
sprintf(str,"Cannot open COMB3 C library file");
|
||||
sprintf(str,"Cannot open COMB3 C library file \n");
|
||||
error->one(FLERR,str);
|
||||
}
|
||||
|
||||
|
|
@ -750,17 +750,17 @@ void PairComb3::read_file(char *file)
|
|||
params[nparams].pcnb = atof(words[59]);
|
||||
params[nparams].pcnc = atof(words[60]);
|
||||
params[nparams].pcnd = atof(words[61]);
|
||||
params[nparams].plp1 = atof(words[62]);
|
||||
params[nparams].plp2 = atof(words[63]);
|
||||
params[nparams].plp3 = atof(words[64]);
|
||||
params[nparams].plp4 = atof(words[65]);
|
||||
params[nparams].plp5 = atof(words[66]);
|
||||
params[nparams].plp6 = atof(words[67]);
|
||||
params[nparams].addrep = atof(words[68]);
|
||||
params[nparams].pbb1=atof(words[69]);
|
||||
params[nparams].pbb2=atof(words[70]);
|
||||
params[nparams].ptork1=atof(words[71]);
|
||||
params[nparams].ptork2=atof(words[72]);
|
||||
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);
|
||||
|
||||
|
|
@ -774,7 +774,7 @@ void PairComb3::read_file(char *file)
|
|||
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].addrep < 0.0 || params[nparams].powermint < 1.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 ||
|
||||
|
|
@ -860,11 +860,12 @@ void PairComb3::Short_neigh()
|
|||
int n,nj,*neighptrj,icontrol;
|
||||
int iparam_ij,iparam_ji,iparam_jk,*ilist,*jlist,*numneigh,**firstneigh;
|
||||
int inum,jnum,i,j,k,l,ii,jj,kk,ll,itype,jtype,ktype,ltype;
|
||||
tagint itag,jtag;
|
||||
tagint itag, jtag;
|
||||
double rr1,rr2,rsq1,rsq2,delrj[3],delrk[3];
|
||||
int sht_knum,*sht_klist;
|
||||
|
||||
double **x = atom->x;
|
||||
tagint *tag = atom->tag;
|
||||
int *type = atom->type;
|
||||
int ntype = atom->ntypes;
|
||||
int nlocal = atom->nlocal;
|
||||
|
|
@ -895,6 +896,7 @@ void PairComb3::Short_neigh()
|
|||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
itag = tag[i];
|
||||
dpl[i][0] = dpl[i][1] = dpl[i][2] = 0.0;
|
||||
|
||||
nj = 0;
|
||||
|
|
@ -911,6 +913,7 @@ void PairComb3::Short_neigh()
|
|||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
jtag = tag[j];
|
||||
|
||||
delrj[0] = x[i][0] - x[j][0];
|
||||
delrj[1] = x[i][1] - x[j][1];
|
||||
|
|
@ -960,15 +963,15 @@ void PairComb3::compute(int eflag, int vflag)
|
|||
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;
|
||||
tagint itag,jtag;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double rr,rsq,rsq1,rsq2,rsq3,iq,jq,yaself;
|
||||
double fqi,fqij,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],delri[3],fi[3],fj[3],fk[3],fl[3];
|
||||
double elp_ij,elp_ji,filp[3],fjlp[3],fklp[3],fllp[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;
|
||||
|
|
@ -1332,7 +1335,7 @@ void PairComb3::compute(int eflag, int vflag)
|
|||
eflag, eng_tmp, iq, jq, i, j, nj, bbtor, kconjug, lconjug);
|
||||
|
||||
evdwl = eng_tmp;
|
||||
selflp(¶ms[iparam_ij],¶ms[iparam_ji],rsq1,eng_tmp,fpair);
|
||||
selfp6p(¶ms[iparam_ij],¶ms[iparam_ji],rsq1,eng_tmp,fpair);
|
||||
|
||||
evdwl += eng_tmp;
|
||||
f[i][0] += delrj[0]*fpair;
|
||||
|
|
@ -1371,14 +1374,14 @@ void PairComb3::compute(int eflag, int vflag)
|
|||
prefac_ij1, prefac_ij2, prefac_ij3, prefac_ij4, prefac_ij5,
|
||||
rsq1,rsq2,delrj,delrk,fi,fj,fk,i,xcn);
|
||||
|
||||
elp_ij = elp(¶ms[iparam_ijk],¶ms[iparam_ikj],rsq1,rsq2,delrj,delrk,zet_addi);
|
||||
flp(¶ms[iparam_ijk],¶ms[iparam_ikj],rsq1,rsq2,delrj,delrk,filp,fjlp,fklp);
|
||||
ep6p_ij = ep6p(¶ms[iparam_ijk],¶ms[iparam_ikj],rsq1,rsq2,delrj,delrk,zet_addi);
|
||||
fp6p(¶ms[iparam_ijk],¶ms[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] += filp[im];
|
||||
fj[im] += fjlp[im];
|
||||
fk[im] += fklp[im];
|
||||
fi[im] += fip6p[im];
|
||||
fj[im] += fjp6p[im];
|
||||
fk[im] += fkp6p[im];
|
||||
}
|
||||
|
||||
// Sums up i-j-k forces: Tallies into global force vector
|
||||
|
|
@ -1465,7 +1468,7 @@ void PairComb3::compute(int eflag, int vflag)
|
|||
}
|
||||
|
||||
if (evflag)
|
||||
ev_tally(i,j,nlocal,newton_pair,elp_ij,0.0,0.0,0.0,0.0,0.0);
|
||||
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);
|
||||
|
||||
|
|
@ -1500,17 +1503,17 @@ void PairComb3::compute(int eflag, int vflag)
|
|||
rsq1,rsq2,delrl,delrk,fj,fi,fl,j,ycn);
|
||||
|
||||
// BO-independent 3-body j-i-l LP and BB correction and forces
|
||||
elp_ji = elp(¶ms[iparam_jil],¶ms[iparam_jli],rsq1,rsq2,delrl,delrk,zet_addj);
|
||||
flp(¶ms[iparam_jil],¶ms[iparam_jli],rsq1,rsq2,delrl,delrk,fjlp,filp,fllp);
|
||||
ep6p_ji = ep6p(¶ms[iparam_jil],¶ms[iparam_jli],rsq1,rsq2,delrl,delrk,zet_addj);
|
||||
fp6p(¶ms[iparam_jil],¶ms[iparam_jli],rsq1,rsq2,delrl,delrk,fjp6p,fip6p,flp6p);
|
||||
|
||||
if (evflag)
|
||||
ev_tally(j,i,nlocal,newton_pair,elp_ji,0.0,0.0,0.0,0.0,0.0);
|
||||
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] += fjlp[im];
|
||||
fi[im] += filp[im];
|
||||
fl[im] += fllp[im];
|
||||
fj[im] += fjp6p[im];
|
||||
fi[im] += fip6p[im];
|
||||
fl[im] += flp6p[im];
|
||||
}
|
||||
|
||||
// Sums up j-i-l forces: Tallies into global force vector
|
||||
|
|
@ -1583,6 +1586,7 @@ void PairComb3::repulsive(Param *parami, Param *paramj, double rsq,
|
|||
|
||||
double romi = parami->addrep;
|
||||
double rrcs = parami->bigr + parami->bigd;
|
||||
double addr = parami->addrepr;
|
||||
|
||||
r = sqrt(rsq);
|
||||
if (r > rrcs) return ;
|
||||
|
|
@ -1605,9 +1609,9 @@ void PairComb3::repulsive(Param *parami, Param *paramj, double rsq,
|
|||
// additional repulsion
|
||||
|
||||
vrcs = 1.0; fvrcs = 0.0;
|
||||
if (romi > 0.0) {
|
||||
vrcs += romi * pow((1.0-r/rrcs),2.0);
|
||||
fvrcs = romi * 2.0 * (r/rrcs-1.0)/rrcs;
|
||||
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;
|
||||
|
|
@ -1639,7 +1643,7 @@ double PairComb3::zeta(Param *parami, Param *paramj, double rsqij,
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairComb3::selflp(Param *parami, Param *paramj, double rsq,
|
||||
void PairComb3::selfp6p(Param *parami, Param *paramj, double rsq,
|
||||
double &eng, double &force)
|
||||
{
|
||||
double r,comtti,comttj,fcj,fcj_d;
|
||||
|
|
@ -1650,21 +1654,15 @@ void PairComb3::selflp(Param *parami, Param *paramj, double rsq,
|
|||
fcj_d=comb_fc_d(r,parami);
|
||||
comtti = comttj = 0.0;
|
||||
|
||||
if (fabs(parami->plp1) > 1.0e-6 || fabs(parami->plp2) > 1.0e-6
|
||||
|| fabs(parami->plp3) > 1.0e-6 || fabs(parami->plp4) > 1.0e-6
|
||||
|| fabs(parami->plp5) > 1.0e-6 || fabs(parami->plp6) > 1.0e-6 ) {
|
||||
double pilp1 = parami->plp1, pilp2 = parami->plp2, pilp3 = parami->plp3;
|
||||
double pilp4 = parami->plp4, pilp5 = parami->plp5, pilp6 = parami->plp6;
|
||||
comtti = pilp1 + pilp2 + pilp3 + pilp4 + pilp5 + pilp6;
|
||||
}
|
||||
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;
|
||||
|
||||
if (fabs(paramj->plp1) > 1.0e-6 || fabs(paramj->plp2) > 1.0e-6
|
||||
|| fabs(paramj->plp3) > 1.0e-6 || fabs(paramj->plp4) > 1.0e-6
|
||||
|| fabs(paramj->plp5) > 1.0e-6 || fabs(paramj->plp6) > 1.0e-6) {
|
||||
double pjlp1 = paramj->plp1, pjlp2 = paramj->plp2, pjlp3 = paramj->plp3;
|
||||
double pjlp4 = paramj->plp4, pjlp5 = paramj->plp5, pjlp6 = paramj->plp6;
|
||||
comttj = pjlp1 + pjlp2 + pjlp3 + pjlp4 + pjlp5 + pjlp6;
|
||||
}
|
||||
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;
|
||||
|
|
@ -1672,79 +1670,56 @@ void PairComb3::selflp(Param *parami, Param *paramj, double rsq,
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairComb3::elp(Param *paramj, Param *paramk, double rsqij, double rsqik,
|
||||
double PairComb3::ep6p(Param *paramj, Param *paramk, double rsqij, double rsqik,
|
||||
double *delrij, double *delrik , double &zet_add)
|
||||
{
|
||||
int con_flag, plp_flag;
|
||||
con_flag = plp_flag = 0;
|
||||
if (fabs(paramj->pbb2) > 1.0e-6) con_flag = 1;
|
||||
if (fabs(paramj->plp1) > 1.0e-6 || fabs(paramj->plp2) > 1.0e-6
|
||||
|| fabs(paramj->plp3) > 1.0e-6 || fabs(paramj->plp4) > 1.0e-6
|
||||
|| fabs(paramj->plp5) > 1.0e-6 || fabs(paramj->plp6) > 1.0e-6 )
|
||||
plp_flag = 1;
|
||||
|
||||
if(con_flag || plp_flag) {
|
||||
double rij,rik,costheta,lp1,lp2,lp3,lp4,lp5,lp6;
|
||||
double rmu,rmu2,rmu3,rmu4,rmu5,rmu6,comtt,fcj,fck,fcj_d;
|
||||
double pplp1 = paramj->plp1, pplp2 = paramj->plp2, pplp3 = paramj->plp3;
|
||||
double pplp4 = paramj->plp4, pplp5 = paramj->plp5, pplp6 = paramj->plp6;
|
||||
|
||||
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,fcj_d;
|
||||
comtt=0.0;
|
||||
rij = sqrt(rsqij);
|
||||
rik = sqrt(rsqik);
|
||||
costheta = vec3_dot(delrij,delrik) / (rij*rik);
|
||||
fcj = comb_fc(rij,paramj);
|
||||
fcj_d = comb_fc_d(rij,paramj);
|
||||
fck = comb_fc(rik,paramk);
|
||||
|
||||
rmu = costheta;
|
||||
|
||||
// Legendre Polynomial functions
|
||||
if (plp_flag) {
|
||||
rmu2 = rmu *rmu; rmu3 = rmu2*rmu; rmu4 = rmu3*rmu;
|
||||
rmu5 = rmu4*rmu; rmu6 = rmu5*rmu;
|
||||
lp1 = pplp1*rmu;
|
||||
lp2 = pplp2*0.5*(3.0*rmu2-1.0);
|
||||
lp3 = pplp3*0.5*(5.0*rmu3-3.0*rmu);
|
||||
lp4 = pplp4*(35.0*rmu4-30.0*rmu2+3.0)/8.0;
|
||||
lp5 = pplp5*(63.0*rmu5-70.0*rmu3+15.0*rmu)/8.0;
|
||||
lp6 = pplp6*(231.0*rmu6-315.0*rmu4+105.0*rmu2-5.0)/16.0;
|
||||
comtt = lp1 + lp2 + lp3 + lp4 + lp5 + lp6;
|
||||
}
|
||||
else {comtt = 0.0;}
|
||||
|
||||
// bond-bending terms
|
||||
|
||||
if (con_flag) {
|
||||
double c123 = paramj->pbb1;
|
||||
comtt += paramj->pbb2 *(rmu-c123)*(rmu-c123);
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
/*---------------------------------------------------------------------- */
|
||||
|
||||
void PairComb3::flp(Param *paramij,Param *paramik, double rsqij, double rsqik,
|
||||
void PairComb3::fp6p(Param *paramij,Param *paramik, double rsqij, double rsqik,
|
||||
double *delrij, double *delrik, double *drilp,
|
||||
double *drjlp, double *drklp)
|
||||
{
|
||||
int con_flag, plp_flag;
|
||||
double ffj1,ffj2,ffk1,ffk2,com5k;
|
||||
|
||||
con_flag = plp_flag = 0;
|
||||
ffj1 = 0.0; ffj2 = 0.0; ffk1 = 0.0; ffk2 = 0.0;
|
||||
if (fabs(paramij->pbb2) > 1.0e-6) con_flag = 1;
|
||||
if (fabs(paramij->plp1) > 1.0e-6 || fabs(paramij->plp2) > 1.0e-6
|
||||
|| fabs(paramij->plp3) > 1.0e-6 || fabs(paramij->plp4) > 1.0e-6
|
||||
|| fabs(paramij->plp5) > 1.0e-6 || fabs(paramij->plp6) > 1.0e-6 )
|
||||
|
||||
plp_flag=1;
|
||||
|
||||
if(con_flag || plp_flag) {
|
||||
double rij,rik,costheta;
|
||||
double rmu,comtt,comtt_d,com4k,com5,fcj,fcj_d,fck,fck_d;
|
||||
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);
|
||||
|
|
@ -1752,57 +1727,34 @@ void PairComb3::flp(Param *paramij,Param *paramik, double rsqij, double rsqik,
|
|||
fcj_d = comb_fc_d(rij,paramij);
|
||||
fck_d = comb_fc_d(rik,paramik);
|
||||
rmu = costheta;
|
||||
|
||||
// Legendre Polynomial functions and derivatives
|
||||
|
||||
if (plp_flag) {
|
||||
double 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;
|
||||
double pplp1 = paramij->plp1, pplp2 = paramij->plp2, pplp3 = paramij->plp3;
|
||||
double pplp4 = paramij->plp4, pplp5 = paramij->plp5, pplp6 = paramij->plp6;
|
||||
|
||||
rmu2 = rmu *rmu; rmu3 = rmu2*rmu;
|
||||
rmu4 = rmu3*rmu; rmu5 = rmu4*rmu; rmu6 = rmu5*rmu;
|
||||
lp0 = pplp0;
|
||||
lp1 = pplp1*rmu;
|
||||
lp2 = pplp2*0.5*(3.0*rmu2-1.0);
|
||||
lp3 = pplp3*0.5*(5.0*rmu3-3.0*rmu);
|
||||
lp4 = pplp4*(35.0*rmu4-30.0*rmu2+3.0)/8.0;
|
||||
lp5 = pplp5*(63.0*rmu5-70.0*rmu3+15.0*rmu)/8.0;
|
||||
lp6 = pplp6*(231.0*rmu6-315.0*rmu4+105.0*rmu2-5.0)/16.0;
|
||||
lp2 = pplp2*rmu2;
|
||||
lp3 = pplp3*rmu3;
|
||||
lp4 = pplp4*rmu4;
|
||||
lp5 = pplp5*rmu5;
|
||||
lp6 = pplp6*rmu6;
|
||||
lp1_d = pplp1;
|
||||
lp2_d = pplp2*3.0*rmu;
|
||||
lp3_d = pplp3*(7.5*rmu2-1.5);
|
||||
lp4_d = pplp4*(140.0*rmu3-60.0*rmu)/8.0;
|
||||
lp5_d = pplp5*(315.0*rmu4-210.0*rmu2+15.0)/8.0;
|
||||
lp6_d = pplp6*(1386.0*rmu5-1260.0*rmu3+210.0)/16.0;
|
||||
comtt = lp1 + lp2 + lp3 + lp4 + lp5 + lp6;
|
||||
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;
|
||||
} else {
|
||||
comtt = 0.0;
|
||||
comtt_d = 0.0;
|
||||
}
|
||||
|
||||
// bond-bending terms derivatives
|
||||
|
||||
if (con_flag) {
|
||||
double c123 = paramij->pbb1;
|
||||
comtt += paramij->pbb2 *(rmu-c123)*(rmu-c123);
|
||||
comtt_d += 2.0*paramij->pbb2*(rmu-c123);
|
||||
}
|
||||
|
||||
com4k = fcj * fck_d * comtt;
|
||||
com5 = fcj * fck * comtt_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);
|
||||
|
||||
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);
|
||||
} else {
|
||||
ffj1 = 0.0; ffj2 = 0.0;
|
||||
ffk1 = 0.0; ffk2 = 0.0;
|
||||
}
|
||||
|
||||
// j-atom
|
||||
vec3_scale(ffj1,delrik,drjlp);
|
||||
vec3_scaleadd(ffj2,delrij,drjlp,drjlp);
|
||||
|
|
@ -1831,12 +1783,8 @@ void PairComb3::force_zeta(Param *parami, Param *paramj, double rsq,
|
|||
double boij, dbij1, dbij2, dbij3, dbij4, dbij5;
|
||||
double boji, dbji1, dbji2, dbji3, dbji4, dbji5;
|
||||
double pradx, prady;
|
||||
|
||||
int inti=parami->ielement;
|
||||
int intj=paramj->ielement;
|
||||
tagint *tag=atom->tag;
|
||||
tagint itag=tag[i];
|
||||
tagint jtag=tag[j];
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (r > parami->bigr + parami->bigd) return;
|
||||
|
|
@ -2648,9 +2596,9 @@ void PairComb3::tables()
|
|||
|
||||
rvdw[1][inty] = params[iparam_ij].vsig * 0.950;
|
||||
|
||||
// radius check: outer radius vs. sigma
|
||||
// radius check: outter radius vs. sigma
|
||||
if( rvdw[0][inty] > rvdw[1][inty] )
|
||||
error->all(FLERR,"Error in vdw spline: inner radius > outer radius");
|
||||
error->all(FLERR,"Error in vdw spline: inner radius > outter radius");
|
||||
|
||||
rrc[0] = rvdw[1][inty];
|
||||
|
||||
|
|
@ -3347,16 +3295,16 @@ void PairComb3::tor_force(int torindx, Param *paramk, Param *paraml,
|
|||
|
||||
double PairComb3::combqeq(double *qf_fix, int &igroup)
|
||||
{
|
||||
int i,j,ii,jj,itype,jtype,jnum;
|
||||
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;
|
||||
tagint itag,jtag;
|
||||
double xtmp,ytmp,ztmp,rr,rsq,rsq1,rsq2,delrj[3],zeta_ij;
|
||||
double iq,jq,fqi,fqj,fqij,fqji,yaself,yaself_d,sr1,sr2,sr3;
|
||||
double rr_sw,ij_sw,ji_sw,fq_swi,fq_swj;
|
||||
double potal,fac11,fac11e;
|
||||
int sht_jnum,*sht_jlist;
|
||||
tagint itag, jtag;
|
||||
|
||||
double **x = atom->x;
|
||||
double *q = atom->q;
|
||||
|
|
@ -3507,7 +3455,7 @@ double PairComb3::combqeq(double *qf_fix, int &igroup)
|
|||
i = ilist[ii];
|
||||
if (mask[i] & groupbit){
|
||||
eneg += qf[i];
|
||||
itag=tag[i];
|
||||
itag=tag[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -3694,6 +3642,7 @@ void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
|
|||
double QUchi, QOchi, QUchj, QOchj;
|
||||
double bij, caj, cbj, caqpn, caqpj, cbqpn, cbqpj;
|
||||
double LamDiLamDj, AlfDiAlfDj;
|
||||
double addr = parami->addrepr;
|
||||
double romi = parami->addrep;
|
||||
double rrcs = parami->bigr + parami->bigd;
|
||||
double rlm1 = parami->lambda;
|
||||
|
|
@ -3710,7 +3659,7 @@ void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
|
|||
bij = bbij[i][nj];
|
||||
|
||||
// additional repulsion
|
||||
if (romi > 0.0) vrcs = romi * pow((1.0-r/rrcs),2.0);
|
||||
if (romi != 0.0 && r < addr) vrcs = romi * pow((1.0-r/addr),2.0);
|
||||
|
||||
QUchi = (parami->QU - iq) * parami->bD;
|
||||
QUchj = (paramj->QU - jq) * paramj->bD;
|
||||
|
|
|
|||
|
|
@ -47,9 +47,9 @@ class PairComb3 : public Pair {
|
|||
double pcos6,pcos5,pcos4,pcos3,pcos2,pcos1,pcos0;
|
||||
double gamma,powerm,powern,bigA,bigB1,bigB2,bigB3;
|
||||
double bigd,bigr,cut,cutsq,c1,c2,c3,c4;
|
||||
double plp1,plp2,plp3,plp4,plp5,plp6;
|
||||
double pbb1,pbb2,ptork1,ptork2;
|
||||
double addrep, vdwflag;
|
||||
double p6p0,p6p1,p6p2,p6p3,p6p4,p6p5,p6p6;
|
||||
double ptork1,ptork2;
|
||||
double addrepr,addrep, vdwflag;
|
||||
double QU,QL,DU,DL,Qo,dQ,aB,bB,nD,bD,qmin,qmax;
|
||||
double chi,dj,dk,dl,dm,esm,cmn1,cmn2,pcmn1,pcmn2;
|
||||
double coulcut, lcut, lcutsq;
|
||||
|
|
@ -182,9 +182,9 @@ class PairComb3 : public Pair {
|
|||
double &, double &, double &, double &, Param *);
|
||||
|
||||
// Legendre polynomials
|
||||
void selflp(Param *, Param *, double, double &, double &);
|
||||
double elp(Param *, Param *, double, double, double *, double * ,double &);
|
||||
void flp(Param *, Param *, double, double, double *, double *, double *,
|
||||
void selfp6p(Param *, Param *, double, double &, double &);
|
||||
double ep6p(Param *, Param *, double, double, double *, double * ,double &);
|
||||
void fp6p(Param *, Param *, double, double, double *, double *, double *,
|
||||
double *, double *);
|
||||
|
||||
// long range q-dependent terms
|
||||
|
|
@ -280,39 +280,34 @@ E: Incorrect args for pair coefficients
|
|||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style COMB3 requires atom IDs
|
||||
E: Pair style COMB requires atom IDs
|
||||
|
||||
This is a requirement to use this potential.
|
||||
This is a requirement to use the AIREBO potential.
|
||||
|
||||
E: Pair style COMB3 requires newton pair on
|
||||
E: Pair style COMB requires newton pair on
|
||||
|
||||
See the newton command. This is a restriction to use the COMB3
|
||||
See the newton command. This is a restriction to use the COMB
|
||||
potential.
|
||||
|
||||
E: Pair style COMB3 requires atom attribute q
|
||||
E: Pair style COMB requires atom attribute q
|
||||
|
||||
The atom style defined does not have this attribute.
|
||||
Self-explanatory.
|
||||
|
||||
E: All pair coeffs are not set
|
||||
|
||||
All pair coefficients must be set in the data file or by the
|
||||
pair_coeff command before running a simulation.
|
||||
|
||||
E: Cannot open COMB3 C library file
|
||||
|
||||
The extra lib.comb3 file for carbon cannot be opened. Check that it
|
||||
exists.
|
||||
|
||||
E: Cannot open COMB3 potential file %s
|
||||
E: Cannot open COMB potential file %s
|
||||
|
||||
The specified COMB potential file cannot be opened. Check that the
|
||||
path and name are correct.
|
||||
|
||||
E: Incorrect format in COMB3 potential file
|
||||
E: Incorrect format in COMB potential file
|
||||
|
||||
Incorrect number of words per line in the potential file.
|
||||
|
||||
E: Illegal COMB3 parameter
|
||||
E: Illegal COMB parameter
|
||||
|
||||
One or more of the coefficients defined in the potential file is
|
||||
invalid.
|
||||
|
|
@ -327,14 +322,12 @@ E: Potential file is missing an entry
|
|||
The potential file for a SW or Tersoff potential does not have a
|
||||
needed entry.
|
||||
|
||||
E: Neighbor list overflow, boost neigh_modify one
|
||||
W: Pair COMB charge %.10f with force %.10f hit min barrier
|
||||
|
||||
There are too many neighbors of a single atom. Use the neigh_modify
|
||||
command to increase the max number of neighbors allowed for one atom.
|
||||
You may also want to boost the page size.
|
||||
Something is possibly wrong with your model.
|
||||
|
||||
E: Error in vdw spline: inner radius > outer radius
|
||||
W: Pair COMB charge %.10f with force %.10f hit max barrier
|
||||
|
||||
Self-explanatory.
|
||||
Something is possibly wrong with your model.
|
||||
|
||||
*/
|
||||
|
|
|
|||
Loading…
Reference in New Issue