jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7852 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2012-02-28 17:44:54 +00:00
parent c2777b94c2
commit 77377e36da
4 changed files with 224 additions and 222 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/KSPACE/pair_lj_cut_coul_long.cpp
View File

@ -732,7 +732,9 @@ double PairLJCutCoulLong::init_one(int i, int j)
cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
} }
double cut = MAX(cut_lj[i][j],cut_coul); // include TIP4P qdist in full cutoff, qdist = 0.0 if not TIP4P
double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);

1
src/KSPACE/pair_lj_cut_coul_long.h
View File

@ -52,6 +52,7 @@ class PairLJCutCoulLong : public Pair {
double **epsilon,**sigma; double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset; double **lj1,**lj2,**lj3,**lj4,**offset;
double *cut_respa; double *cut_respa;
double qdist; // TIP4P distance from O site to negative charge
double g_ewald; double g_ewald;
double tabinnersq; double tabinnersq;

440
src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
View File

@ -123,238 +123,236 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j]; jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
// LJ interation based on true rsq
if (rsq < cut_ljsq[itype][jtype]) {
r2inv = 1.0/rsq; r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj *= factor_lj * r2inv;
f[i][0] += delx*forcelj;
f[i][1] += dely*forcelj;
f[i][2] += delz*forcelj;
f[j][0] -= delx*forcelj;
f[j][1] -= dely*forcelj;
f[j][2] -= delz*forcelj;
if (eflag) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,forcelj,delx,dely,delz);
}
if (rsq < cut_ljsq[itype][jtype]) { // adjust rsq and delxyz for off-site O charge(s) if necessary
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj *= factor_lj * r2inv;
f[i][0] += delx*forcelj; if (itype == typeO || jtype == typeO) {
f[i][1] += dely*forcelj; if (jtype == typeO) {
f[i][2] += delz*forcelj; find_M(j,jH1,jH2,xjM);
f[j][0] -= delx*forcelj; x2 = xjM;
f[j][1] -= dely*forcelj; } else x2 = x[j];
f[j][2] -= delz*forcelj; delx = x1[0] - x2[0];
dely = x1[1] - x2[1];
delz = x1[2] - x2[2];
rsq = delx*delx + dely*dely + delz*delz;
}
if (eflag) { // Coulombic interaction based on modified rsq
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
if (evflag) ev_tally(i,j,nlocal,newton_pair, if (rsq < cut_coulsq) {
evdwl,0.0,forcelj,delx,dely,delz); r2inv = 1 / rsq;
} if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
// adjust rsq and delxyz for off-site O charge(s) grij = g_ewald * r;
expm2 = exp(-grij*grij);
if (itype == typeO || jtype == typeO) { t = 1.0 / (1.0 + EWALD_P*grij);
if (jtype == typeO) { erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
find_M(j,jH1,jH2,xjM); prefactor = qqrd2e * qtmp*q[j]/r;
x2 = xjM; forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
} else x2 = x[j]; if (factor_coul < 1.0) {
delx = x1[0] - x2[0]; forcecoul -= (1.0-factor_coul)*prefactor;
dely = x1[1] - x2[1];
delz = x1[2] - x2[2];
rsq = delx*delx + dely*dely + delz*delz;
}
// test current rsq against cutoff and compute Coulombic force
if (rsq < cut_coulsq) {
r2inv = 1 / rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) {
forcecoul -= (1.0-factor_coul)*prefactor;
}
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
} }
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
}
cforce = forcecoul * r2inv;
cforce = forcecoul * r2inv; // if i,j are not O atoms, force is applied directly
// if i or j are O atoms, force is on fictitious atom & partitioned
// force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999)
// f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f
// preserves total force and torque on water molecule
// virial = sum(r x F) where each water's atoms are near xi and xj
// vlist stores 2,4,6 atoms whose forces contribute to virial
n = 0;
if (itype != typeO) {
f[i][0] += delx * cforce;
f[i][1] += dely * cforce;
f[i][2] += delz * cforce;
if (vflag) {
v[0] = x[i][0] * delx * cforce;
v[1] = x[i][1] * dely * cforce;
v[2] = x[i][2] * delz * cforce;
v[3] = x[i][0] * dely * cforce;
v[4] = x[i][0] * delz * cforce;
v[5] = x[i][1] * delz * cforce;
vlist[n++] = i;
}
} else {
// if i,j are not O atoms, force is applied directly fd[0] = delx*cforce;
// if i or j are O atoms, force is on fictitious atom & partitioned fd[1] = dely*cforce;
// force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999) fd[2] = delz*cforce;
// f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f
// preserves total force and torque on water molecule
// virial = sum(r x F) where each water's atoms are near xi and xj
// vlist stores 2,4,6 atoms whose forces contribute to virial
n = 0; delxOM = x[i][0] - x1[0];
delyOM = x[i][1] - x1[1];
delzOM = x[i][2] - x1[2];
if (itype != typeO) { ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
f[i][0] += delx * cforce;
f[i][1] += dely * cforce;
f[i][2] += delz * cforce;
if (vflag) {
v[0] = x[i][0] * delx * cforce;
v[1] = x[i][1] * dely * cforce;
v[2] = x[i][2] * delz * cforce;
v[3] = x[i][0] * dely * cforce;
v[4] = x[i][0] * delz * cforce;
v[5] = x[i][1] * delz * cforce;
vlist[n++] = i;
}
} else {
fd[0] = delx*cforce;
fd[1] = dely*cforce;
fd[2] = delz*cforce;
delxOM = x[i][0] - x1[0];
delyOM = x[i][1] - x1[1];
delzOM = x[i][2] - x1[2];
ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
(qdist*qdist); (qdist*qdist);
f1[0] = ddotf * delxOM; f1[0] = ddotf * delxOM;
f1[1] = ddotf * delyOM; f1[1] = ddotf * delyOM;
f1[2] = ddotf * delzOM; f1[2] = ddotf * delzOM;
fO[0] = fd[0] - alpha * (fd[0] - f1[0]); fO[0] = fd[0] - alpha * (fd[0] - f1[0]);
fO[1] = fd[1] - alpha * (fd[1] - f1[1]); fO[1] = fd[1] - alpha * (fd[1] - f1[1]);
fO[2] = fd[2] - alpha * (fd[2] - f1[2]); fO[2] = fd[2] - alpha * (fd[2] - f1[2]);
fH[0] = 0.5 * alpha * (fd[0] - f1[0]); fH[0] = 0.5 * alpha * (fd[0] - f1[0]);
fH[1] = 0.5 * alpha * (fd[1] - f1[1]); fH[1] = 0.5 * alpha * (fd[1] - f1[1]);
fH[2] = 0.5 * alpha * (fd[2] - f1[2]); fH[2] = 0.5 * alpha * (fd[2] - f1[2]);
f[i][0] += fO[0]; f[i][0] += fO[0];
f[i][1] += fO[1]; f[i][1] += fO[1];
f[i][2] += fO[2]; f[i][2] += fO[2];
f[iH1][0] += fH[0]; f[iH1][0] += fH[0];
f[iH1][1] += fH[1]; f[iH1][1] += fH[1];
f[iH1][2] += fH[2]; f[iH1][2] += fH[2];
f[iH2][0] += fH[0]; f[iH2][0] += fH[0];
f[iH2][1] += fH[1]; f[iH2][1] += fH[1];
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[i],x[iH1],xH1); domain->closest_image(x[i],x[iH1],xH1);
domain->closest_image(x[i],x[iH2],xH2); domain->closest_image(x[i],x[iH2],xH2);
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
vlist[n++] = i; vlist[n++] = i;
vlist[n++] = iH1; vlist[n++] = iH1;
vlist[n++] = iH2; vlist[n++] = iH2;
}
} }
if (jtype != typeO) {
f[j][0] -= delx * cforce;
f[j][1] -= dely * cforce;
f[j][2] -= delz * cforce;
if (vflag) {
v[0] -= x[j][0] * delx * cforce;
v[1] -= x[j][1] * dely * cforce;
v[2] -= x[j][2] * delz * cforce;
v[3] -= x[j][0] * dely * cforce;
v[4] -= x[j][0] * delz * cforce;
v[5] -= x[j][1] * delz * cforce;
vlist[n++] = j;
}
} else {
fd[0] = -delx*cforce;
fd[1] = -dely*cforce;
fd[2] = -delz*cforce;
delxOM = x[j][0] - x2[0];
delyOM = x[j][1] - x2[1];
delzOM = x[j][2] - x2[2];
ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
(qdist*qdist);
f1[0] = ddotf * delxOM;
f1[1] = ddotf * delyOM;
f1[2] = ddotf * delzOM;
fO[0] = fd[0] - alpha * (fd[0] - f1[0]);
fO[1] = fd[1] - alpha * (fd[1] - f1[1]);
fO[2] = fd[2] - alpha * (fd[2] - f1[2]);
fH[0] = 0.5 * alpha * (fd[0] - f1[0]);
fH[1] = 0.5 * alpha * (fd[1] - f1[1]);
fH[2] = 0.5 * alpha * (fd[2] - f1[2]);
f[j][0] += fO[0];
f[j][1] += fO[1];
f[j][2] += fO[2];
f[jH1][0] += fH[0];
f[jH1][1] += fH[1];
f[jH1][2] += fH[2];
f[jH2][0] += fH[0];
f[jH2][1] += fH[1];
f[jH2][2] += fH[2];
if (vflag) {
domain->closest_image(x[j],x[jH1],xH1);
domain->closest_image(x[j],x[jH2],xH2);
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
vlist[n++] = j;
vlist[n++] = jH1;
vlist[n++] = jH2;
}
}
if (eflag) {
if (!ncoultablebits || rsq <= tabinnersq)
ecoul = prefactor*erfc;
else {
table = etable[itable] + fraction*detable[itable];
ecoul = qtmp*q[j] * table;
}
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
if (evflag) ev_tally_list(n,vlist,ecoul,v);
} }
if (jtype != typeO) {
f[j][0] -= delx * cforce;
f[j][1] -= dely * cforce;
f[j][2] -= delz * cforce;
if (vflag) {
v[0] -= x[j][0] * delx * cforce;
v[1] -= x[j][1] * dely * cforce;
v[2] -= x[j][2] * delz * cforce;
v[3] -= x[j][0] * dely * cforce;
v[4] -= x[j][0] * delz * cforce;
v[5] -= x[j][1] * delz * cforce;
vlist[n++] = j;
}
} else {
fd[0] = -delx*cforce;
fd[1] = -dely*cforce;
fd[2] = -delz*cforce;
delxOM = x[j][0] - x2[0];
delyOM = x[j][1] - x2[1];
delzOM = x[j][2] - x2[2];
ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
(qdist*qdist);
f1[0] = ddotf * delxOM;
f1[1] = ddotf * delyOM;
f1[2] = ddotf * delzOM;
fO[0] = fd[0] - alpha * (fd[0] - f1[0]);
fO[1] = fd[1] - alpha * (fd[1] - f1[1]);
fO[2] = fd[2] - alpha * (fd[2] - f1[2]);
fH[0] = 0.5 * alpha * (fd[0] - f1[0]);
fH[1] = 0.5 * alpha * (fd[1] - f1[1]);
fH[2] = 0.5 * alpha * (fd[2] - f1[2]);
f[j][0] += fO[0];
f[j][1] += fO[1];
f[j][2] += fO[2];
f[jH1][0] += fH[0];
f[jH1][1] += fH[1];
f[jH1][2] += fH[2];
f[jH2][0] += fH[0];
f[jH2][1] += fH[1];
f[jH2][2] += fH[2];
if (vflag) {
domain->closest_image(x[j],x[jH1],xH1);
domain->closest_image(x[j],x[jH2],xH2);
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
vlist[n++] = j;
vlist[n++] = jH1;
vlist[n++] = jH2;
}
}
if (eflag) {
if (!ncoultablebits || rsq <= tabinnersq)
ecoul = prefactor*erfc;
else {
table = etable[itable] + fraction*detable[itable];
ecoul = qtmp*q[j] * table;
}
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
if (evflag) ev_tally_list(n,vlist,ecoul,v);
} }
} }
} }
@ -397,9 +395,11 @@ void PairLJCutCoulLongTIP4P::init_style()
if (atom->tag_enable == 0) if (atom->tag_enable == 0)
error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires atom IDs"); error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires atom IDs");
if (!force->newton_pair) if (!force->newton_pair)
error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires newton pair on"); error->all(FLERR,
"Pair style lj/cut/coul/long/tip4p requires newton pair on");
if (!atom->q_flag) if (!atom->q_flag)
error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires atom attribute q"); error->all(FLERR,
"Pair style lj/cut/coul/long/tip4p requires atom attribute q");
if ( (strcmp(force->kspace_style,"pppm/tip4p") != 0) && if ( (strcmp(force->kspace_style,"pppm/tip4p") != 0) &&
(strcmp(force->kspace_style,"pppm/tip4p/proxy") != 0) ) (strcmp(force->kspace_style,"pppm/tip4p/proxy") != 0) )
error->all(FLERR,"Pair style is incompatible with KSpace style"); error->all(FLERR,"Pair style is incompatible with KSpace style");

1
src/KSPACE/pair_lj_cut_coul_long_tip4p.h
View File

@ -37,7 +37,6 @@ class PairLJCutCoulLongTIP4P : public PairLJCutCoulLong {
protected: protected:
int typeH,typeO; // atom types of TIP4P water H and O atoms int typeH,typeO; // atom types of TIP4P water H and O atoms
int typeA,typeB; // angle and bond types of TIP4P water int typeA,typeB; // angle and bond types of TIP4P water
double qdist; // distance from O site to negative charge
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P
void find_M(int, int &, int &, double *); void find_M(int, int &, int &, double *);