ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

pair TIP4P bug fix for cutoffs >> box size

Browse Source
This commit is contained in:
Steve Plimpton
2017-01-06 09:57:27 -07:00
parent 62dea1bb63
commit f6a819580c
45 changed files with 663 additions and 416 deletions
Download Patch File Download Diff File Expand all files Collapse all files

70
src/KOKKOS/modify_kokkos.cpp
View File

@ -81,7 +81,7 @@ void ModifyKokkos::setup_pre_exchange()
atomKK->sync(fix[list_min_pre_exchange[i]]->execution_space, atomKK->sync(fix[list_min_pre_exchange[i]]->execution_space,
fix[list_min_pre_exchange[i]]->datamask_read); fix[list_min_pre_exchange[i]]->datamask_read);
if (!fix[list_min_pre_exchange[i]]->kokkosable) lmp->kokkos->auto_sync = 1; if (!fix[list_min_pre_exchange[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
fix[list_min_pre_exchange[i]]->min_setup_pre_exchange(); fix[list_min_pre_exchange[i]]->setup_pre_exchange();
lmp->kokkos->auto_sync = 0; lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_exchange[i]]->execution_space, atomKK->modified(fix[list_min_pre_exchange[i]]->execution_space,
fix[list_min_pre_exchange[i]]->datamask_modify); fix[list_min_pre_exchange[i]]->datamask_modify);
@ -110,7 +110,7 @@ void ModifyKokkos::setup_pre_neighbor()
atomKK->sync(fix[list_min_pre_neighbor[i]]->execution_space, atomKK->sync(fix[list_min_pre_neighbor[i]]->execution_space,
fix[list_min_pre_neighbor[i]]->datamask_read); fix[list_min_pre_neighbor[i]]->datamask_read);
if (!fix[list_min_pre_neighbor[i]]->kokkosable) lmp->kokkos->auto_sync = 1; if (!fix[list_min_pre_neighbor[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
fix[list_min_pre_neighbor[i]]->min_setup_pre_neighbor(); fix[list_min_pre_neighbor[i]]->setup_pre_neighbor();
lmp->kokkos->auto_sync = 0; lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_neighbor[i]]->execution_space, atomKK->modified(fix[list_min_pre_neighbor[i]]->execution_space,
fix[list_min_pre_neighbor[i]]->datamask_modify); fix[list_min_pre_neighbor[i]]->datamask_modify);
@ -139,13 +139,42 @@ void ModifyKokkos::setup_pre_force(int vflag)
atomKK->sync(fix[list_min_pre_force[i]]->execution_space, atomKK->sync(fix[list_min_pre_force[i]]->execution_space,
fix[list_min_pre_force[i]]->datamask_read); fix[list_min_pre_force[i]]->datamask_read);
if (!fix[list_min_pre_force[i]]->kokkosable) lmp->kokkos->auto_sync = 1; if (!fix[list_min_pre_force[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
fix[list_min_pre_force[i]]->min_setup_pre_force(vflag); fix[list_min_pre_force[i]]->setup_pre_force(vflag);
lmp->kokkos->auto_sync = 0; lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_force[i]]->execution_space, atomKK->modified(fix[list_min_pre_force[i]]->execution_space,
fix[list_min_pre_force[i]]->datamask_modify); fix[list_min_pre_force[i]]->datamask_modify);
} }
} }
/* ----------------------------------------------------------------------
setup pre_reverse call, only for fixes that define pre_reverse
called from Verlet, RESPA, Min
------------------------------------------------------------------------- */
void ModifyKokkos::setup_pre_reverse(int eflag, int vflag)
{
if (update->whichflag == 1)
for (int i = 0; i < n_pre_reverse; i++) {
atomKK->sync(fix[list_pre_reverse[i]]->execution_space,
fix[list_pre_reverse[i]]->datamask_read);
if (!fix[list_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
fix[list_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_pre_reverse[i]]->execution_space,
fix[list_pre_reverse[i]]->datamask_modify);
}
else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_reverse; i++) {
atomKK->sync(fix[list_min_pre_reverse[i]]->execution_space,
fix[list_min_pre_reverse[i]]->datamask_read);
if (!fix[list_min_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
fix[list_min_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_reverse[i]]->execution_space,
fix[list_min_pre_reverse[i]]->datamask_modify);
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
1st half of integrate call, only for relevant fixes 1st half of integrate call, only for relevant fixes
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -231,6 +260,23 @@ void ModifyKokkos::pre_force(int vflag)
} }
} }
/* ----------------------------------------------------------------------
pre_reverse call, only for relevant fixes
------------------------------------------------------------------------- */
void ModifyKokkos::pre_reverse(int eflag, int vflag)
{
for (int i = 0; i < n_pre_reverse; i++) {
atomKK->sync(fix[list_pre_reverse[i]]->execution_space,
fix[list_pre_reverse[i]]->datamask_read);
if (!fix[list_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
fix[list_pre_reverse[i]]->pre_reverse(eflag,vflag);
lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_pre_reverse[i]]->execution_space,
fix[list_pre_reverse[i]]->datamask_modify);
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
post_force call, only for relevant fixes post_force call, only for relevant fixes
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -476,6 +522,23 @@ void ModifyKokkos::min_pre_force(int vflag)
} }
} }
/* ----------------------------------------------------------------------
minimizer pre-reverse call, only for relevant fixes
------------------------------------------------------------------------- */
void ModifyKokkos::min_pre_reverse(int eflag, int vflag)
{
for (int i = 0; i < n_min_pre_reverse; i++) {
atomKK->sync(fix[list_min_pre_reverse[i]]->execution_space,
fix[list_min_pre_reverse[i]]->datamask_read);
if (!fix[list_min_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
fix[list_min_pre_reverse[i]]->min_pre_reverse(eflag,vflag);
lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_reverse[i]]->execution_space,
fix[list_min_pre_reverse[i]]->datamask_modify);
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
minimizer force adjustment call, only for relevant fixes minimizer force adjustment call, only for relevant fixes
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -658,4 +721,3 @@ int ModifyKokkos::min_reset_ref()
} }
return itmpall; return itmpall;
} }

3
src/KOKKOS/modify_kokkos.h
View File

@ -26,12 +26,14 @@ class ModifyKokkos : public Modify {
void setup_pre_exchange(); void setup_pre_exchange();
void setup_pre_neighbor(); void setup_pre_neighbor();
void setup_pre_force(int); void setup_pre_force(int);
void setup_pre_reverse(int, int);
void initial_integrate(int); void initial_integrate(int);
void post_integrate(); void post_integrate();
void pre_decide(); void pre_decide();
void pre_exchange(); void pre_exchange();
void pre_neighbor(); void pre_neighbor();
void pre_force(int); void pre_force(int);
void pre_reverse(int,int);
void post_force(int); void post_force(int);
void final_integrate(); void final_integrate();
void end_of_step(); void end_of_step();
@ -48,6 +50,7 @@ class ModifyKokkos : public Modify {
void min_pre_exchange(); void min_pre_exchange();
void min_pre_neighbor(); void min_pre_neighbor();
void min_pre_force(int); void min_pre_force(int);
void min_pre_reverse(int,int);
void min_post_force(int); void min_post_force(int);
double min_energy(double *); double min_energy(double *);

40
src/KSPACE/pair_lj_cut_tip4p_long.cpp
View File

@ -88,8 +88,8 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
double r,r2inv,r6inv,forcecoul,forcelj,cforce; double r,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj; double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; double fO[3],fH[3],fd[3],v[6];
double *x1,*x2; double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
double rsq; double rsq;
@ -146,14 +146,20 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(tag[i] + 1); iH1 = atom->map(tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(tag[i] + 2); iH2 = atom->map(tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -216,14 +222,20 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(tag[j] + 1); jH1 = atom->map(tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(tag[j] + 2); jH2 = atom->map(tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -327,9 +339,8 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[i],x[iH1],xH1); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
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];
@ -385,9 +396,8 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; 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[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[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@ -565,12 +575,10 @@ void PairLJCutTIP4PLong::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);

40
src/KSPACE/pair_lj_long_tip4p_long.cpp
View File

@ -86,8 +86,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
double r,r2inv,forcecoul,forcelj,cforce; double r,r2inv,forcecoul,forcelj,cforce;
double factor_coul; double factor_coul;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];// f1[3]; double fO[3],fH[3],fd[3],v[6];
double *x1,*x2; double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
double rsq; double rsq;
@ -145,14 +145,20 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
itype = type[i]; itype = type[i];
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(tag[i] + 1); iH1 = atom->map(tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(tag[i] + 2); iH2 = atom->map(tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -253,14 +259,20 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
if (itype == typeO || jtype == typeO) { if (itype == typeO || jtype == typeO) {
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(tag[j] + 1); jH1 = atom->map(tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(tag[j] + 2); jH2 = atom->map(tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -365,9 +377,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[i],x[iH1],xH1); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
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];
@ -423,9 +434,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; 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[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[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@ -1548,12 +1558,10 @@ void PairLJLongTIP4PLong::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);

40
src/KSPACE/pair_tip4p_long.cpp
View File

@ -84,8 +84,8 @@ void PairTIP4PLong::compute(int eflag, int vflag)
double r,r2inv,forcecoul,cforce; double r,r2inv,forcecoul,cforce;
double factor_coul; double factor_coul;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; double fO[3],fH[3],fd[3],v[6];
double *x1,*x2; double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
double rsq; double rsq;
@ -140,14 +140,20 @@ void PairTIP4PLong::compute(int eflag, int vflag)
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(tag[i] + 1); iH1 = atom->map(tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(tag[i] + 2); iH2 = atom->map(tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -184,14 +190,20 @@ void PairTIP4PLong::compute(int eflag, int vflag)
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(tag[j] + 1); jH1 = atom->map(tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(tag[j] + 2); jH2 = atom->map(tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -295,9 +307,8 @@ void PairTIP4PLong::compute(int eflag, int vflag)
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[i],x[iH1],xH1); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
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];
@ -353,9 +364,8 @@ void PairTIP4PLong::compute(int eflag, int vflag)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; 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[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[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@ -488,12 +498,10 @@ void PairTIP4PLong::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);

7
src/KSPACE/pppm_disp_tip4p.cpp
View File

@ -502,17 +502,20 @@ void PPPMDispTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
double **x = atom->x; double **x = atom->x;
double delx1 = x[iH1][0] - x[i][0]; double delx1 = x[iH1][0] - x[i][0];
double dely1 = x[iH1][1] - x[i][1]; double dely1 = x[iH1][1] - x[i][1];
double delz1 = x[iH1][2] - x[i][2]; double delz1 = x[iH1][2] - x[i][2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2][0] - x[i][0]; double delx2 = x[iH2][0] - x[i][0];
double dely2 = x[iH2][1] - x[i][1]; double dely2 = x[iH2][1] - x[i][1];
double delz2 = x[iH2][2] - x[i][2]; double delz2 = x[iH2][2] - x[i][2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2); xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2); xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);

7
src/KSPACE/pppm_tip4p.cpp
View File

@ -493,17 +493,20 @@ void PPPMTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
double **x = atom->x; double **x = atom->x;
double delx1 = x[iH1][0] - x[i][0]; double delx1 = x[iH1][0] - x[i][0];
double dely1 = x[iH1][1] - x[i][1]; double dely1 = x[iH1][1] - x[i][1];
double delz1 = x[iH1][2] - x[i][2]; double delz1 = x[iH1][2] - x[i][2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2][0] - x[i][0]; double delx2 = x[iH2][0] - x[i][0];
double dely2 = x[iH2][1] - x[i][1]; double dely2 = x[iH2][1] - x[i][1];
double delz2 = x[iH2][2] - x[i][2]; double delz2 = x[iH2][2] - x[i][2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2); xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2); xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);

18
src/MOLECULE/fix_cmap.cpp
View File

@ -62,10 +62,13 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
crosstermlist(NULL), num_crossterm(NULL), crossterm_type(NULL), crossterm_atom1(NULL), Fix(lmp, narg, arg),
crossterm_atom2(NULL), crossterm_atom3(NULL), crossterm_atom4(NULL), crossterm_atom5(NULL), crosstermlist(NULL), num_crossterm(NULL), crossterm_type(NULL),
g_axis(NULL), cmapgrid(NULL), d1cmapgrid(NULL), d2cmapgrid(NULL), d12cmapgrid(NULL) crossterm_atom1(NULL), crossterm_atom2(NULL), crossterm_atom3(NULL),
crossterm_atom4(NULL), crossterm_atom5(NULL),
g_axis(NULL), cmapgrid(NULL), d1cmapgrid(NULL), d2cmapgrid(NULL),
d12cmapgrid(NULL)
{ {
if (narg != 4) error->all(FLERR,"Illegal fix cmap command"); if (narg != 4) error->all(FLERR,"Illegal fix cmap command");
@ -208,6 +211,13 @@ void FixCMAP::setup_pre_neighbor()
/* --------------------------------------------------------------------- */ /* --------------------------------------------------------------------- */
void FixCMAP::setup_pre_reverse(int eflag, int vflag)
{
pre_reverse(eflag,vflag);
}
/* --------------------------------------------------------------------- */
void FixCMAP::min_setup(int vflag) void FixCMAP::min_setup(int vflag)
{ {
pre_neighbor(); pre_neighbor();

1
src/MOLECULE/fix_cmap.h
View File

@ -31,6 +31,7 @@ class FixCMAP : public Fix {
void init(); void init();
void setup(int); void setup(int);
void setup_pre_neighbor(); void setup_pre_neighbor();
void setup_pre_reverse(int, int);
void min_setup(int); void min_setup(int);
void pre_neighbor(); void pre_neighbor();
void pre_reverse(int, int); void pre_reverse(int, int);

42
src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
View File

@ -87,8 +87,8 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
int n,vlist[6]; int n,vlist[6];
int iH1,iH2,jH1,jH2; int iH1,iH2,jH1,jH2;
double cforce; double cforce;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; double fO[3],fH[3],fd[3],v[6];
double *x1,*x2; double *x1,*x2,*xH1,*xH2;
evdwl = ecoul = 0.0; evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag); if (eflag || vflag) ev_setup(eflag,vflag);
@ -139,14 +139,20 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(tag[i] + 1); iH1 = atom->map(tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(tag[i] + 2); iH2 = atom->map(tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -209,14 +215,20 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(tag[j] + 1); jH1 = atom->map(tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(tag[j] + 2); jH2 = atom->map(tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -294,10 +306,9 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
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); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
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];
@ -353,9 +364,8 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; 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[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[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@ -711,12 +721,10 @@ void PairLJCutTIP4PCut::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);

40
src/MOLECULE/pair_tip4p_cut.cpp
View File

@ -75,8 +75,8 @@ void PairTIP4PCut::compute(int eflag, int vflag)
int n,vlist[6]; int n,vlist[6];
int iH1,iH2,jH1,jH2; int iH1,iH2,jH1,jH2;
double cforce; double cforce;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; double fO[3],fH[3],fd[3],v[6];
double *x1,*x2; double *x1,*x2,*xH1,*xH2;
ecoul = 0.0; ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag); if (eflag || vflag) ev_setup(eflag,vflag);
@ -125,14 +125,20 @@ void PairTIP4PCut::compute(int eflag, int vflag)
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(tag[i] + 1); iH1 = atom->map(tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(tag[i] + 2); iH2 = atom->map(tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -169,14 +175,20 @@ void PairTIP4PCut::compute(int eflag, int vflag)
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(tag[j] + 1); jH1 = atom->map(tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(tag[j] + 2); jH2 = atom->map(tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -255,9 +267,8 @@ void PairTIP4PCut::compute(int eflag, int vflag)
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if(vflag) { if(vflag) {
domain->closest_image(x[i],x[iH1],xH1); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
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];
@ -313,9 +324,8 @@ void PairTIP4PCut::compute(int eflag, int vflag)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; 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[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[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@ -526,12 +536,10 @@ void PairTIP4PCut::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);

38
src/OPT/pair_lj_cut_tip4p_long_opt.cpp
View File

@ -111,9 +111,9 @@ void PairLJCutTIP4PLongOpt::eval()
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce; double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj; double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double v[6],xH1[3],xH2[3]; double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz; double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
const double *x1,*x2; const double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,inum,jnum,itype,jtype,itable,key; int i,j,ii,jj,inum,jnum,itype,jtype,itable,key;
@ -155,14 +155,21 @@ void PairLJCutTIP4PLongOpt::eval()
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(tag[i] + 1); iH1 = atom->map(tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(tag[i] + 2); iH2 = atom->map(tag[i] + 2);
hneigh[i][2] = 1; hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_opt(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite_opt(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -226,14 +233,21 @@ void PairLJCutTIP4PLongOpt::eval()
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(tag[j] + 1); jH1 = atom->map(tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(tag[j] + 2); jH2 = atom->map(tag[j] + 2);
hneigh[j][2] = 1; hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_opt(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite_opt(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -339,9 +353,8 @@ void PairLJCutTIP4PLongOpt::eval()
f[iH2][2] += fHz; f[iH2][2] += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(x[i],x[iH1],xH1); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
v[0] = x[i][0]*fOx + xH1[0]*fHx + xH2[0]*fHx; v[0] = x[i][0]*fOx + xH1[0]*fHx + xH2[0]*fHx;
v[1] = x[i][1]*fOy + xH1[1]*fHy + xH2[1]*fHy; v[1] = x[i][1]*fOy + xH1[1]*fHy + xH2[1]*fHy;
v[2] = x[i][2]*fOz + xH1[2]*fHz + xH2[2]*fHz; v[2] = x[i][2]*fOz + xH1[2]*fHz + xH2[2]*fHz;
@ -399,9 +412,8 @@ void PairLJCutTIP4PLongOpt::eval()
f[jH2][2] += fHz; f[jH2][2] += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fOx + xH1[0]*fHx + xH2[0]*fHx; v[0] += x[j][0]*fOx + xH1[0]*fHx + xH2[0]*fHx;
v[1] += x[j][1]*fOy + xH1[1]*fHy + xH2[1]*fHy; v[1] += x[j][1]*fOy + xH1[1]*fHy + xH2[1]*fHy;
v[2] += x[j][2]*fOz + xH1[2]*fHz + xH2[2]*fHz; v[2] += x[j][2]*fOz + xH1[2]*fHz + xH2[2]*fHz;
@ -448,12 +460,10 @@ void PairLJCutTIP4PLongOpt::compute_newsite_opt(const double * xO,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM[0] = xO[0] + prefac * (delx1 + delx2); xM[0] = xO[0] + prefac * (delx1 + delx2);

7
src/QEQ/fix_qeq.cpp
View File

@ -310,13 +310,6 @@ void FixQEq::setup_pre_force_respa(int vflag, int ilevel)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixQEq::min_setup_pre_force(int vflag)
{
setup_pre_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixQEq::init_storage() void FixQEq::init_storage()
{ {
nlocal = atom->nlocal; nlocal = atom->nlocal;

1
src/QEQ/fix_qeq.h
View File

@ -33,7 +33,6 @@ class FixQEq : public Fix {
void setup_pre_force(int); void setup_pre_force(int);
void setup_pre_force_respa(int, int); void setup_pre_force_respa(int, int);
void pre_force_respa(int, int, int); void pre_force_respa(int, int, int);
void min_setup_pre_force(int);
void min_pre_force(int); void min_pre_force(int);
// derived child classes must provide these functions // derived child classes must provide these functions

5
src/REPLICA/compute_event_displace.cpp
View File

@ -48,7 +48,7 @@ ComputeEventDisplace::ComputeEventDisplace(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Distance must be > 0 for compute event/displace"); error->all(FLERR,"Distance must be > 0 for compute event/displace");
displace_distsq = displace_dist * displace_dist; displace_distsq = displace_dist * displace_dist;
// fix event ID will be set later by PRD // fix event ID will be set later by accelerated dynamics method
id_event = NULL; id_event = NULL;
} }
@ -75,7 +75,8 @@ void ComputeEventDisplace::init()
fix_event = (FixEvent*) modify->fix[ifix]; fix_event = (FixEvent*) modify->fix[ifix];
if (strcmp(fix_event->style,"EVENT/PRD") != 0 && if (strcmp(fix_event->style,"EVENT/PRD") != 0 &&
strcmp(fix_event->style,"EVENT/TAD") != 0) strcmp(fix_event->style,"EVENT/TAD") != 0 &&
strcmp(fix_event->style,"EVENT/HYPER") != 0)
error->all(FLERR,"Compute event/displace has invalid fix event assigned"); error->all(FLERR,"Compute event/displace has invalid fix event assigned");
} }

2
src/REPLICA/fix_neb.cpp
View File

@ -395,7 +395,7 @@ void FixNEB::inter_replica_comm()
// ----------------------------------------------------- // -----------------------------------------------------
// single proc per replica // single proc per replica
// all atoms are NEB atoms and no atom sorting // all atoms are NEB atoms and no atom sorting is enabled
// direct comm of x -> xprev and x -> xnext // direct comm of x -> xprev and x -> xnext
if (cmode == SINGLE_PROC_DIRECT) { if (cmode == SINGLE_PROC_DIRECT) {

23
src/USER-ATC/fix_atc.cpp
View File

@ -632,16 +632,6 @@ void FixATC::min_pre_exchange()
throw; throw;
} }
} }
void FixATC::min_setup_pre_exchange()
{
try {
atc_->setup_pre_exchange();
}
catch (ATC::ATC_Error& atcError) {
ATC::LammpsInterface::instance()->print_msg(atcError.error_description());
throw;
}
}
double FixATC::memory_usage() double FixATC::memory_usage()
{ {
@ -896,19 +886,6 @@ void FixATC::setup_pre_neighbor()
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixATC::min_setup_pre_neighbor()
{
if (atc_->is_initialized()) {
try {
atc_->pre_neighbor();
}
catch (ATC::ATC_Error& atcError) {
ATC::LammpsInterface::instance()->print_msg(atcError.error_description());
throw;
}
}
}
/* ---------------------------------------------------------------------- */
void FixATC::min_pre_force(int vflag) void FixATC::min_pre_force(int vflag)
{ {
try { try {

2
src/USER-ATC/fix_atc.h
View File

@ -52,7 +52,6 @@ namespace LAMMPS_NS {
and is called before domain->pbc() and comm->exchange(). */ and is called before domain->pbc() and comm->exchange(). */
void setup_pre_exchange(); void setup_pre_exchange();
void pre_exchange(); void pre_exchange();
void min_setup_pre_exchange();
void min_pre_exchange(); void min_pre_exchange();
double memory_usage(); double memory_usage();
@ -84,7 +83,6 @@ namespace LAMMPS_NS {
neighbor->build(). */ neighbor->build(). */
void pre_neighbor(); void pre_neighbor();
void setup_pre_neighbor(); void setup_pre_neighbor();
void min_setup_pre_neighbor();
/** pre/post_force is used to modify fix-specific data /** pre/post_force is used to modify fix-specific data
and is before/after the various force computations. */ and is before/after the various force computations. */

7
src/USER-DPD/fix_shardlow.cpp
View File

@ -178,13 +178,6 @@ void FixShardlow::min_pre_exchange()
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixShardlow::min_setup_pre_exchange()
{
memset(atom->ssaAIR, 0, sizeof(int)*atom->nlocal);
}
/* ---------------------------------------------------------------------- */
void FixShardlow::setup(int vflag) void FixShardlow::setup(int vflag)
{ {
bool fixShardlow = false; bool fixShardlow = false;

1
src/USER-DPD/fix_shardlow.h
View File

@ -37,7 +37,6 @@ class FixShardlow : public Fix {
virtual void initial_integrate(int); virtual void initial_integrate(int);
void setup_pre_exchange(); void setup_pre_exchange();
void pre_exchange(); void pre_exchange();
void min_setup_pre_exchange();
void min_pre_exchange(); void min_pre_exchange();
void grow_arrays(int); void grow_arrays(int);

40
src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp
View File

@ -86,8 +86,8 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
double factor_coul,factor_lj; double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double denc, denlj, r4sig6; double denc, denlj, r4sig6;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; double fO[3],fH[3],fd[3],v[6];
double *x1,*x2; double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
double rsq; double rsq;
@ -143,14 +143,20 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1); iH1 = atom->map(atom->tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2); iH2 = atom->map(atom->tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -216,14 +222,20 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1); jH1 = atom->map(atom->tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2); jH2 = atom->map(atom->tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -314,9 +326,8 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[i],x[iH1],xH1); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
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];
@ -372,9 +383,8 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; 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[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[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@ -549,12 +559,10 @@ void PairLJCutTIP4PLongSoft::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);

40
src/USER-FEP/pair_tip4p_long_soft.cpp
View File

@ -85,8 +85,8 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
double factor_coul; double factor_coul;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double denc; double denc;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; double fO[3],fH[3],fd[3],v[6];
double *x1,*x2; double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
double rsq; double rsq;
@ -140,14 +140,20 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
if (itype == typeO) { if (itype == typeO) {
if (hneigh[i][0] < 0) { if (hneigh[i][0] < 0) {
hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1); iH1 = atom->map(atom->tag[i] + 1);
hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2); iH2 = atom->map(atom->tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
hneigh[i][0] = iH1;
hneigh[i][1] = iH2;
hneigh[i][2] = 1;
} else { } else {
iH1 = hneigh[i][0]; iH1 = hneigh[i][0];
iH2 = hneigh[i][1]; iH2 = hneigh[i][1];
@ -184,14 +190,20 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh[j][0] < 0) { if (hneigh[j][0] < 0) {
hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1); jH1 = atom->map(atom->tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2); jH2 = atom->map(atom->tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
hneigh[j][0] = jH1;
hneigh[j][1] = jH2;
hneigh[j][2] = 1;
} else { } else {
jH1 = hneigh[j][0]; jH1 = hneigh[j][0];
jH2 = hneigh[j][1]; jH2 = hneigh[j][1];
@ -282,9 +294,8 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[i],x[iH1],xH1); xH1 = x[iH1];
domain->closest_image(x[i],x[iH2],xH2); xH2 = x[iH2];
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];
@ -340,9 +351,8 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (vflag) { if (vflag) {
domain->closest_image(x[j],x[jH1],xH1); xH1 = x[jH1];
domain->closest_image(x[j],x[jH2],xH2); xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; 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[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[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@ -473,12 +483,10 @@ void PairTIP4PLongSoft::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0]; double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1]; double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2]; double delz1 = xH1[2] - xO[2];
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0]; double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1]; double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2]; double delz2 = xH2[2] - xO[2];
domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);

7
src/USER-INTEL/fix_intel.cpp
View File

@ -356,6 +356,13 @@ void FixIntel::setup(int vflag)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixIntel::setup_pre_reverse(int eflag, int vflag)
{
pre_reverse(eflag,vflag);
}
/* ---------------------------------------------------------------------- */
void FixIntel::pair_init_check(const bool cdmessage) void FixIntel::pair_init_check(const bool cdmessage)
{ {
#ifdef INTEL_VMASK #ifdef INTEL_VMASK

2
src/USER-INTEL/fix_intel.h
View File

@ -43,6 +43,8 @@ class FixIntel : public Fix {
virtual int setmask(); virtual int setmask();
virtual void init(); virtual void init();
virtual void setup(int); virtual void setup(int);
void setup_pre_reverse(int eflag = 0, int vflag = 0);
void pair_init_check(const bool cdmessage=false); void pair_init_check(const bool cdmessage=false);
void bond_init_check(); void bond_init_check();
void kspace_init_check(); void kspace_init_check();

47
src/USER-OMP/pair_lj_cut_tip4p_cut_omp.cpp
View File

@ -127,9 +127,9 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce; double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj; double factor_coul,factor_lj;
double v[6],xH1[3],xH2[3]; double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz; double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
dbl3_t x1,x2; dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key; int i,j,ii,jj,jnum,itype,jtype,key;
@ -179,10 +179,14 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]); compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1; hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2; hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1; hneigh_thr[i].a = iH1;
} else { } else {
iH1 = hneigh_thr[i].a; iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b; iH2 = hneigh_thr[i].b;
@ -256,6 +260,9 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]); compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1; hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2; hneigh_thr[j].b = jH2;
@ -341,15 +348,14 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz; f[iH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[i].x,&x[iH1].x,xH1); xH1 = x[iH1];
domain->closest_image(&x[i].x,&x[iH2].x,xH2); xH2 = x[iH2];
v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = i; vlist[n++] = i;
@ -401,15 +407,14 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz; f[jH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[j].x,&x[jH1].x,xH1); xH1 = x[jH1];
domain->closest_image(&x[j].x,&x[jH2].x,xH2); xH2 = x[jH2];
v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = j; vlist[n++] = j;
@ -446,12 +451,10 @@ void PairLJCutTIP4PCutOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x; double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y; double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z; double delz1 = xH1.z - xO.z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x; double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y; double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z; double delz2 = xH2.z - xO.z;
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2); xM.x = xO.x + prefac * (delx1 + delx2);

46
src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
View File

@ -141,9 +141,9 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce; double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj; double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double v[6],xH1[3],xH2[3]; double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz; double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
dbl3_t x1,x2; dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,itable, key; int i,j,ii,jj,jnum,itype,jtype,itable, key;
@ -193,6 +193,9 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]); compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1; hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2; hneigh_thr[i].b = iH2;
@ -270,6 +273,9 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]); compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1; hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2; hneigh_thr[j].b = jH2;
@ -379,15 +385,14 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz; f[iH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[i].x,&x[iH1].x,xH1); xH1 = x[iH1];
domain->closest_image(&x[i].x,&x[iH2].x,xH2); xH2 = x[iH2];
v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = i; vlist[n++] = i;
@ -439,15 +444,14 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz; f[jH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[j].x,&x[jH1].x,xH1); xH1 = x[jH1];
domain->closest_image(&x[j].x,&x[jH2].x,xH2); xH2 = x[jH2];
v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = j; vlist[n++] = j;
@ -489,12 +493,10 @@ void PairLJCutTIP4PLongOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x; double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y; double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z; double delz1 = xH1.z - xO.z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x; double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y; double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z; double delz2 = xH2.z - xO.z;
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2); xM.x = xO.x + prefac * (delx1 + delx2);

46
src/USER-OMP/pair_lj_cut_tip4p_long_soft_omp.cpp
View File

@ -129,9 +129,9 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
double factor_coul,factor_lj; double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double denc, denlj, r4sig6; double denc, denlj, r4sig6;
double v[6],xH1[3],xH2[3]; double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz; double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
dbl3_t x1,x2; dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key; int i,j,ii,jj,jnum,itype,jtype,key;
@ -181,6 +181,9 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]); compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1; hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2; hneigh_thr[i].b = iH2;
@ -261,6 +264,9 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]); compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1; hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2; hneigh_thr[j].b = jH2;
@ -357,15 +363,14 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz; f[iH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[i].x,&x[iH1].x,xH1); xH1 = x[iH1];
domain->closest_image(&x[i].x,&x[iH2].x,xH2); xH2 = x[iH2];
v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = i; vlist[n++] = i;
@ -417,15 +422,14 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz; f[jH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[j].x,&x[jH1].x,xH1); xH1 = x[jH1];
domain->closest_image(&x[j].x,&x[jH2].x,xH2); xH2 = x[jH2];
v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = j; vlist[n++] = j;
@ -463,12 +467,10 @@ void PairLJCutTIP4PLongSoftOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x; double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y; double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z; double delz1 = xH1.z - xO.z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x; double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y; double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z; double delz2 = xH2.z - xO.z;
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2); xM.x = xO.x + prefac * (delx1 + delx2);

64
src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
View File

@ -737,8 +737,8 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,r2inv,forcecoul,forcelj,cforce; double r,r2inv,forcecoul,forcelj,cforce;
double factor_coul; double factor_coul;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; double fO[3],fH[3],fd[3],v[6];
dbl3_t x1,x2; dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
double rsq; double rsq;
@ -763,14 +763,20 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
itype = type[i]; itype = type[i];
if (itype == typeO) { if (itype == typeO) {
if (hneigh_thr[i].a < 0) { if (hneigh_thr[i].a < 0) {
hneigh_thr[i].a = iH1 = atom->map(atom->tag[i] + 1); iH1 = atom->map(atom->tag[i] + 1);
hneigh_thr[i].b = iH2 = atom->map(atom->tag[i] + 2); iH2 = atom->map(atom->tag[i] + 2);
hneigh_thr[i].t = 1;
if (iH1 == -1 || iH2 == -1) if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]); compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].a = iH1;
hneigh_thr[i].b = iH2;
hneigh_thr[i].t = 1;
} else { } else {
iH1 = hneigh_thr[i].a; iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b; iH2 = hneigh_thr[i].b;
@ -871,14 +877,20 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
if (itype == typeO || jtype == typeO) { if (itype == typeO || jtype == typeO) {
if (jtype == typeO) { if (jtype == typeO) {
if (hneigh_thr[j].a < 0) { if (hneigh_thr[j].a < 0) {
hneigh_thr[j].a = jH1 = atom->map(atom->tag[j] + 1); jH1 = atom->map(atom->tag[j] + 1);
hneigh_thr[j].b = jH2 = atom->map(atom->tag[j] + 2); jH2 = atom->map(atom->tag[j] + 2);
hneigh_thr[j].t = 1;
if (jH1 == -1 || jH2 == -1) if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]); compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].a = jH1;
hneigh_thr[j].b = jH2;
hneigh_thr[j].t = 1;
} else { } else {
jH1 = hneigh_thr[j].a; jH1 = hneigh_thr[j].a;
jH2 = hneigh_thr[j].b; jH2 = hneigh_thr[j].b;
@ -983,15 +995,14 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2][2] += fH[2]; f[iH2][2] += fH[2];
if (EVFLAG) { if (EVFLAG) {
domain->closest_image(&x[i].x,&x[iH1].x,xH1); xH1 = x[iH1];
domain->closest_image(&x[i].x,&x[iH2].x,xH2); xH2 = x[iH2];
v[0] = x[i].x*fO[0] + xH1.x*fH[0] + xH2.x*fH[0];
v[0] = x[i].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; v[1] = x[i].y*fO[1] + xH1.y*fH[1] + xH2.y*fH[1];
v[1] = x[i].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; v[2] = x[i].z*fO[2] + xH1.z*fH[2] + xH2.z*fH[2];
v[2] = x[i].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; v[3] = x[i].x*fO[1] + xH1.x*fH[1] + xH2.x*fH[1];
v[3] = x[i].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; v[4] = x[i].x*fO[2] + xH1.x*fH[2] + xH2.x*fH[2];
v[4] = x[i].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; v[5] = x[i].y*fO[2] + xH1.y*fH[2] + xH2.y*fH[2];
v[5] = x[i].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
} }
vlist[n++] = i; vlist[n++] = i;
vlist[n++] = iH1; vlist[n++] = iH1;
@ -1041,15 +1052,14 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2][2] += fH[2]; f[jH2][2] += fH[2];
if (EVFLAG) { if (EVFLAG) {
domain->closest_image(&x[j].x,&x[jH1].x,xH1); xH1 = x[jH1];
domain->closest_image(&x[j].x,&x[jH2].x,xH2); xH2 = x[jH2];
v[0] += x[j].x*fO[0] + xH1.x*fH[0] + xH2.x*fH[0];
v[0] += x[j].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; v[1] += x[j].y*fO[1] + xH1.y*fH[1] + xH2.y*fH[1];
v[1] += x[j].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; v[2] += x[j].z*fO[2] + xH1.z*fH[2] + xH2.z*fH[2];
v[2] += x[j].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; v[3] += x[j].x*fO[1] + xH1.x*fH[1] + xH2.x*fH[1];
v[3] += x[j].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; v[4] += x[j].x*fO[2] + xH1.x*fH[2] + xH2.x*fH[2];
v[4] += x[j].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; v[5] += x[j].y*fO[2] + xH1.y*fH[2] + xH2.y*fH[2];
v[5] += x[j].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
} }
vlist[n++] = j; vlist[n++] = j;
vlist[n++] = jH1; vlist[n++] = jH1;
@ -1981,12 +1991,10 @@ void PairLJLongTIP4PLongOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x; double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y; double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z; double delz1 = xH1.z - xO.z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x; double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y; double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z; double delz2 = xH2.z - xO.z;
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2); xM.x = xO.x + prefac * (delx1 + delx2);

47
src/USER-OMP/pair_tip4p_cut_omp.cpp
View File

@ -126,9 +126,9 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul;
double r,rsq,r2inv,forcecoul,cforce; double r,rsq,r2inv,forcecoul,cforce;
double factor_coul; double factor_coul;
double v[6],xH1[3],xH2[3]; double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz; double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
dbl3_t x1,x2; dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key; int i,j,ii,jj,jnum,itype,jtype,key;
@ -176,10 +176,14 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]); compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1; hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2; hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1; hneigh_thr[i].a = iH1;
} else { } else {
iH1 = hneigh_thr[i].a; iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b; iH2 = hneigh_thr[i].b;
@ -227,6 +231,9 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]); compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1; hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2; hneigh_thr[j].b = jH2;
@ -312,15 +319,14 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz; f[iH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[i].x,&x[iH1].x,xH1); xH1 = x[iH1];
domain->closest_image(&x[i].x,&x[iH2].x,xH2); xH2 = x[iH2];
v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = i; vlist[n++] = i;
@ -372,15 +378,14 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz; f[jH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[j].x,&x[jH1].x,xH1); xH1 = x[jH1];
domain->closest_image(&x[j].x,&x[jH2].x,xH2); xH2 = x[jH2];
v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = j; vlist[n++] = j;
@ -417,12 +422,10 @@ void PairTIP4PCutOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x; double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y; double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z; double delz1 = xH1.z - xO.z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x; double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y; double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z; double delz2 = xH2.z - xO.z;
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2); xM.x = xO.x + prefac * (delx1 + delx2);

48
src/USER-OMP/pair_tip4p_long_omp.cpp
View File

@ -141,9 +141,9 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,rsq,r2inv,forcecoul,cforce; double r,rsq,r2inv,forcecoul,cforce;
double factor_coul; double factor_coul;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double v[6],xH1[3],xH2[3]; double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz; double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
dbl3_t x1,x2; dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,itable, key; int i,j,ii,jj,jnum,itype,jtype,itable, key;
@ -191,10 +191,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]); compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1; hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2; hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1; hneigh_thr[i].a = iH1;
} else { } else {
iH1 = hneigh_thr[i].a; iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b; iH2 = hneigh_thr[i].b;
@ -242,10 +246,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]); compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1; hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2; hneigh_thr[j].b = jH2;
hneigh_thr[j].a = jH1; hneigh_thr[j].a = jH1;
} else { } else {
jH1 = hneigh_thr[j].a; jH1 = hneigh_thr[j].a;
jH2 = hneigh_thr[j].b; jH2 = hneigh_thr[j].b;
@ -351,15 +359,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz; f[iH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[i].x,&x[iH1].x,xH1); xH1 = x[iH1];
domain->closest_image(&x[i].x,&x[iH2].x,xH2); xH2 = x[iH2];
v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = i; vlist[n++] = i;
@ -411,15 +418,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz; f[jH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[j].x,&x[jH1].x,xH1); xH1 = x[jH1];
domain->closest_image(&x[j].x,&x[jH2].x,xH2); xH2 = x[jH2];
v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = j; vlist[n++] = j;
@ -461,12 +467,10 @@ void PairTIP4PLongOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x; double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y; double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z; double delz1 = xH1.z - xO.z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x; double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y; double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z; double delz2 = xH2.z - xO.z;
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2); xM.x = xO.x + prefac * (delx1 + delx2);

48
src/USER-OMP/pair_tip4p_long_soft_omp.cpp
View File

@ -128,9 +128,9 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,rsq,forcecoul,cforce; double r,rsq,forcecoul,cforce;
double factor_coul,denc; double factor_coul,denc;
double grij,expm2,prefactor,t,erfc; double grij,expm2,prefactor,t,erfc;
double v[6],xH1[3],xH2[3]; double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz; double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
dbl3_t x1,x2; dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key; int i,j,ii,jj,jnum,itype,jtype,key;
@ -178,10 +178,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]); compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1; hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2; hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1; hneigh_thr[i].a = iH1;
} else { } else {
iH1 = hneigh_thr[i].a; iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b; iH2 = hneigh_thr[i].b;
@ -229,10 +233,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing"); error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set jH1,jH2 to closest image to O
jH1 = domain->closest_image(j,jH1);
jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]); compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1; hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2; hneigh_thr[j].b = jH2;
hneigh_thr[j].a = jH1; hneigh_thr[j].a = jH1;
} else { } else {
jH1 = hneigh_thr[j].a; jH1 = hneigh_thr[j].a;
jH2 = hneigh_thr[j].b; jH2 = hneigh_thr[j].b;
@ -326,15 +334,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz; f[iH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[i].x,&x[iH1].x,xH1); xH1 = x[iH1];
domain->closest_image(&x[i].x,&x[iH2].x,xH2); xH2 = x[iH2];
v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = i; vlist[n++] = i;
@ -386,15 +393,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz; f[jH2].z += fHz;
if (VFLAG) { if (VFLAG) {
domain->closest_image(&x[j].x,&x[jH1].x,xH1); xH1 = x[jH1];
domain->closest_image(&x[j].x,&x[jH2].x,xH2); xH2 = x[jH2];
v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx; v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy; v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz; v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy; v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz; v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
} }
if (EVFLAG) { if (EVFLAG) {
vlist[n++] = j; vlist[n++] = j;
@ -432,12 +438,10 @@ void PairTIP4PLongSoftOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x; double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y; double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z; double delz1 = xH1.z - xO.z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x; double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y; double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z; double delz2 = xH2.z - xO.z;
domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5; const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2); xM.x = xO.x + prefac * (delx1 + delx2);

7
src/USER-OMP/pppm_disp_tip4p_omp.cpp
View File

@ -1852,17 +1852,20 @@ void PPPMDispTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
double delx1 = x[iH1].x - x[i].x; double delx1 = x[iH1].x - x[i].x;
double dely1 = x[iH1].y - x[i].y; double dely1 = x[iH1].y - x[i].y;
double delz1 = x[iH1].z - x[i].z; double delz1 = x[iH1].z - x[i].z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2].x - x[i].x; double delx2 = x[iH2].x - x[i].x;
double dely2 = x[iH2].y - x[i].y; double dely2 = x[iH2].y - x[i].y;
double delz2 = x[iH2].z - x[i].z; double delz2 = x[iH2].z - x[i].z;
domain->minimum_image(delx2,dely2,delz2);
xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2); xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2);
xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2); xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);

7
src/USER-OMP/pppm_tip4p_omp.cpp
View File

@ -742,17 +742,20 @@ void PPPMTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
// set iH1,iH2 to index of closest image to O
iH1 = domain->closest_image(i,iH1);
iH2 = domain->closest_image(i,iH2);
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
double delx1 = x[iH1].x - x[i].x; double delx1 = x[iH1].x - x[i].x;
double dely1 = x[iH1].y - x[i].y; double dely1 = x[iH1].y - x[i].y;
double delz1 = x[iH1].z - x[i].z; double delz1 = x[iH1].z - x[i].z;
domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2].x - x[i].x; double delx2 = x[iH2].x - x[i].x;
double dely2 = x[iH2].y - x[i].y; double dely2 = x[iH2].y - x[i].y;
double delz2 = x[iH2].z - x[i].z; double delz2 = x[iH2].z - x[i].z;
domain->minimum_image(delx2,dely2,delz2);
xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2); xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2);
xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2); xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);

7
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -440,13 +440,6 @@ void FixQEqReax::setup_pre_force_respa(int vflag, int ilevel)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixQEqReax::min_setup_pre_force(int vflag)
{
setup_pre_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixQEqReax::init_storage() void FixQEqReax::init_storage()
{ {
int NN; int NN;

1
src/USER-REAXC/fix_qeq_reax.h
View File

@ -48,7 +48,6 @@ class FixQEqReax : public Fix {
void setup_pre_force_respa(int, int); void setup_pre_force_respa(int, int);
void pre_force_respa(int, int, int); void pre_force_respa(int, int, int);
void min_setup_pre_force(int);
void min_pre_force(int); void min_pre_force(int);
int matvecs; int matvecs;

40
src/domain.cpp
View File

@ -938,28 +938,31 @@ void Domain::subbox_too_small_check(double thresh)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
minimum image convention minimum image convention in periodic dimensions
use 1/2 of box size as test use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed for triclinic, also add/subtract tilt factors in other dims as needed
changed "if" to "while" to enable distance to
far-away ghost atom returned by atom->map() to be wrapped back into box
could be problem for looking up atom IDs when cutoff > boxsize
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void Domain::minimum_image(double &dx, double &dy, double &dz) void Domain::minimum_image(double &dx, double &dy, double &dz)
{ {
if (triclinic == 0) { if (triclinic == 0) {
if (xperiodic) { if (xperiodic) {
if (fabs(dx) > xprd_half) { while (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd; if (dx < 0.0) dx += xprd;
else dx -= xprd; else dx -= xprd;
} }
} }
if (yperiodic) { if (yperiodic) {
if (fabs(dy) > yprd_half) { while (fabs(dy) > yprd_half) {
if (dy < 0.0) dy += yprd; if (dy < 0.0) dy += yprd;
else dy -= yprd; else dy -= yprd;
} }
} }
if (zperiodic) { if (zperiodic) {
if (fabs(dz) > zprd_half) { while (fabs(dz) > zprd_half) {
if (dz < 0.0) dz += zprd; if (dz < 0.0) dz += zprd;
else dz -= zprd; else dz -= zprd;
} }
@ -967,7 +970,7 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
} else { } else {
if (zperiodic) { if (zperiodic) {
if (fabs(dz) > zprd_half) { while (fabs(dz) > zprd_half) {
if (dz < 0.0) { if (dz < 0.0) {
dz += zprd; dz += zprd;
dy += yz; dy += yz;
@ -980,7 +983,7 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
} }
} }
if (yperiodic) { if (yperiodic) {
if (fabs(dy) > yprd_half) { while (fabs(dy) > yprd_half) {
if (dy < 0.0) { if (dy < 0.0) {
dy += yprd; dy += yprd;
dx += xy; dx += xy;
@ -991,7 +994,7 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
} }
} }
if (xperiodic) { if (xperiodic) {
if (fabs(dx) > xprd_half) { while (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd; if (dx < 0.0) dx += xprd;
else dx -= xprd; else dx -= xprd;
} }
@ -1000,28 +1003,31 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
minimum image convention minimum image convention in periodic dimensions
use 1/2 of box size as test use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed for triclinic, also add/subtract tilt factors in other dims as needed
changed "if" to "while" to enable distance to
far-away ghost atom returned by atom->map() to be wrapped back into box
could be problem for looking up atom IDs when cutoff > boxsize
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void Domain::minimum_image(double *delta) void Domain::minimum_image(double *delta)
{ {
if (triclinic == 0) { if (triclinic == 0) {
if (xperiodic) { if (xperiodic) {
if (fabs(delta[0]) > xprd_half) { while (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd; if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd; else delta[0] -= xprd;
} }
} }
if (yperiodic) { if (yperiodic) {
if (fabs(delta[1]) > yprd_half) { while (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) delta[1] += yprd; if (delta[1] < 0.0) delta[1] += yprd;
else delta[1] -= yprd; else delta[1] -= yprd;
} }
} }
if (zperiodic) { if (zperiodic) {
if (fabs(delta[2]) > zprd_half) { while (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) delta[2] += zprd; if (delta[2] < 0.0) delta[2] += zprd;
else delta[2] -= zprd; else delta[2] -= zprd;
} }
@ -1029,7 +1035,7 @@ void Domain::minimum_image(double *delta)
} else { } else {
if (zperiodic) { if (zperiodic) {
if (fabs(delta[2]) > zprd_half) { while (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) { if (delta[2] < 0.0) {
delta[2] += zprd; delta[2] += zprd;
delta[1] += yz; delta[1] += yz;
@ -1042,7 +1048,7 @@ void Domain::minimum_image(double *delta)
} }
} }
if (yperiodic) { if (yperiodic) {
if (fabs(delta[1]) > yprd_half) { while (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) { if (delta[1] < 0.0) {
delta[1] += yprd; delta[1] += yprd;
delta[0] += xy; delta[0] += xy;
@ -1053,7 +1059,7 @@ void Domain::minimum_image(double *delta)
} }
} }
if (xperiodic) { if (xperiodic) {
if (fabs(delta[0]) > xprd_half) { while (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd; if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd; else delta[0] -= xprd;
} }
@ -1092,12 +1098,16 @@ int Domain::closest_image(int i, int j)
closest = j; closest = j;
} }
} }
return closest; return closest;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
find and return Xj image = periodic image of Xj that is closest to Xi find and return Xj image = periodic image of Xj that is closest to Xi
for triclinic, add/subtract tilt factors in other dims as needed for triclinic, add/subtract tilt factors in other dims as needed
not currently used (Jan 2017):
used to be called by pair TIP4P styles but no longer,
due to use of other closest_image() method
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void Domain::closest_image(const double * const xi, const double * const xj, void Domain::closest_image(const double * const xi, const double * const xj,
@ -1984,7 +1994,7 @@ void Domain::subbox_corners()
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
compute 8 corner pts of any triclinic box with lo/hi in lamda coords compute 8 corner pts of any triclinic box with lo/hi in lamda coords
8 output conners are ordered with x changing fastest, then y, finally z 8 output corners are ordered with x changing fastest, then y, finally z
could be more efficient if just coded with xy,yz,xz explicitly could be more efficient if just coded with xy,yz,xz explicitly
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

73
src/finish.cpp
View File

@ -65,7 +65,8 @@ void Finish::end(int flag)
{ {
int i,m,nneigh,nneighfull; int i,m,nneigh,nneighfull;
int histo[10]; int histo[10];
int minflag,prdflag,tadflag,timeflag,fftflag,histoflag,neighflag; int minflag,prdflag,tadflag,hyperflag;
int timeflag,fftflag,histoflag,neighflag;
double time,tmp,ave,max,min; double time,tmp,ave,max,min;
double time_loop,time_other,cpu_loop; double time_loop,time_other,cpu_loop;
@ -86,9 +87,11 @@ void Finish::end(int flag)
// flag = 1 = dynamics or minimization // flag = 1 = dynamics or minimization
// flag = 2 = PRD // flag = 2 = PRD
// flag = 3 = TAD // flag = 3 = TAD
// flag = 4 = HYPER
// turn off neighflag for Kspace partition of verlet/split integrator // turn off neighflag for Kspace partition of verlet/split integrator
minflag = prdflag = tadflag = timeflag = fftflag = histoflag = neighflag = 0; minflag = prdflag = tadflag = hyperflag = 0;
timeflag = fftflag = histoflag = neighflag = 0;
time_loop = cpu_loop = time_other = 0.0; time_loop = cpu_loop = time_other = 0.0;
if (flag == 1) { if (flag == 1) {
@ -103,6 +106,7 @@ void Finish::end(int flag)
} }
if (flag == 2) prdflag = timeflag = histoflag = neighflag = 1; if (flag == 2) prdflag = timeflag = histoflag = neighflag = 1;
if (flag == 3) tadflag = histoflag = neighflag = 1; if (flag == 3) tadflag = histoflag = neighflag = 1;
if (flag == 4) hyperflag = timeflag = histoflag = neighflag = 1;
// loop stats // loop stats
@ -161,8 +165,10 @@ void Finish::end(int flag)
if (lmp->kokkos) { if (lmp->kokkos) {
const char fmt2[] = const char fmt2[] =
"%.1f%% CPU use with %d MPI tasks x %d OpenMP threads\n"; "%.1f%% CPU use with %d MPI tasks x %d OpenMP threads\n";
if (screen) fprintf(screen,fmt2,cpu_loop,nprocs,lmp->kokkos->num_threads); if (screen) fprintf(screen,fmt2,cpu_loop,nprocs,
if (logfile) fprintf(logfile,fmt2,cpu_loop,nprocs,lmp->kokkos->num_threads); lmp->kokkos->num_threads);
if (logfile) fprintf(logfile,fmt2,cpu_loop,nprocs,
lmp->kokkos->num_threads);
} else { } else {
#if defined(_OPENMP) #if defined(_OPENMP)
const char fmt2[] = const char fmt2[] =
@ -409,12 +415,61 @@ void Finish::end(int flag)
} }
} }
// HYPER stats using PAIR,BOND,KSPACE for dynamics,quench
if (hyperflag) {
if (me == 0) {
if (screen) fprintf(screen,"\nHyper stats:\n");
if (logfile) fprintf(logfile,"\nHyper stats:\n");
}
time = timer->get_wall(Timer::DYNAMICS);
MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
time = tmp/nprocs;
if (me == 0) {
if (screen)
fprintf(screen," Dynamics time (%%) = %g (%g)\n",
time,time/time_loop*100.0);
if (logfile)
fprintf(logfile," Dynamics time (%%) = %g (%g)\n",
time,time/time_loop*100.0);
}
time = timer->get_wall(Timer::QUENCH);
MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
time = tmp/nprocs;
if (me == 0) {
if (screen)
fprintf(screen," Quench time (%%) = %g (%g)\n",
time,time/time_loop*100.0);
if (logfile)
fprintf(logfile," Quench time (%%) = %g (%g)\n",
time,time/time_loop*100.0);
}
time = time_other;
MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
time = tmp/nprocs;
if (me == 0) {
if (screen)
fprintf(screen," Other time (%%) = %g (%g)\n",
time,time/time_loop*100.0);
if (logfile)
fprintf(logfile," Other time (%%) = %g (%g)\n",
time,time/time_loop*100.0);
}
}
// further timing breakdowns
if (timeflag && timer->has_normal()) { if (timeflag && timer->has_normal()) {
if (timer->has_full()) { if (timer->has_full()) {
const char hdr[] = "\nMPI task timing breakdown:\n" const char hdr[] = "\nMPI task timing breakdown:\n"
"Section | min time | avg time | max time |%varavg| %CPU | %total\n" "Section | min time | avg time | max time "
"-----------------------------------------------------------------------\n"; "|%varavg| %CPU | %total\n"
"-----------------------------------------------"
"------------------------\n";
if (me == 0) { if (me == 0) {
if (screen) fputs(hdr,screen); if (screen) fputs(hdr,screen);
if (logfile) fputs(hdr,logfile); if (logfile) fputs(hdr,logfile);
@ -516,8 +571,10 @@ void Finish::end(int flag)
if (lmp->kokkos && lmp->kokkos->ngpu > 0) if (lmp->kokkos && lmp->kokkos->ngpu > 0)
if (const char* env_clb = getenv("CUDA_LAUNCH_BLOCKING")) if (const char* env_clb = getenv("CUDA_LAUNCH_BLOCKING"))
if (!(strcmp(env_clb,"1") == 0)) { if (!(strcmp(env_clb,"1") == 0)) {
error->warning(FLERR,"Timing breakdown may not be accurate since GPU/CPU overlap is enabled. " error->warning(FLERR,"Timing breakdown may not be accurate "
"Using 'export CUDA_LAUNCH_BLOCKING=1' will give an accurate timing breakdown but will reduce performance"); "since GPU/CPU overlap is enabled\n"
"Using 'export CUDA_LAUNCH_BLOCKING=1' will give an "
"accurate timing breakdown but will reduce performance");
} }
// FFT timing statistics // FFT timing statistics

40
src/fix.h
View File

@ -113,6 +113,7 @@ class Fix : protected Pointers {
virtual void setup_pre_exchange() {} virtual void setup_pre_exchange() {}
virtual void setup_pre_neighbor() {} virtual void setup_pre_neighbor() {}
virtual void setup_pre_force(int) {} virtual void setup_pre_force(int) {}
virtual void setup_pre_reverse(int, int) {}
virtual void min_setup(int) {} virtual void min_setup(int) {}
virtual void initial_integrate(int) {} virtual void initial_integrate(int) {}
virtual void post_integrate() {} virtual void post_integrate() {}
@ -151,12 +152,10 @@ class Fix : protected Pointers {
virtual void post_force_respa(int, int, int) {} virtual void post_force_respa(int, int, int) {}
virtual void final_integrate_respa(int, int) {} virtual void final_integrate_respa(int, int) {}
virtual void min_setup_pre_exchange() {}
virtual void min_setup_pre_neighbor() {}
virtual void min_setup_pre_force(int) {}
virtual void min_pre_exchange() {} virtual void min_pre_exchange() {}
virtual void min_pre_neighbor() {} virtual void min_pre_neighbor() {}
virtual void min_pre_force(int) {} virtual void min_pre_force(int) {}
virtual void min_pre_reverse(int,int) {}
virtual void min_post_force(int) {} virtual void min_post_force(int) {}
virtual double min_energy(double *) {return 0.0;} virtual double min_energy(double *) {return 0.0;}
@ -241,23 +240,24 @@ namespace FixConst {
static const int PRE_EXCHANGE = 1<<2; static const int PRE_EXCHANGE = 1<<2;
static const int PRE_NEIGHBOR = 1<<3; static const int PRE_NEIGHBOR = 1<<3;
static const int PRE_FORCE = 1<<4; static const int PRE_FORCE = 1<<4;
static const int POST_FORCE = 1<<5; static const int PRE_REVERSE = 1<<5;
static const int FINAL_INTEGRATE = 1<<6; static const int POST_FORCE = 1<<6;
static const int END_OF_STEP = 1<<7; static const int FINAL_INTEGRATE = 1<<7;
static const int THERMO_ENERGY = 1<<8; static const int END_OF_STEP = 1<<8;
static const int INITIAL_INTEGRATE_RESPA = 1<<9; static const int POST_RUN = 1<<9;
static const int POST_INTEGRATE_RESPA = 1<<10; static const int THERMO_ENERGY = 1<<10;
static const int PRE_FORCE_RESPA = 1<<11; static const int INITIAL_INTEGRATE_RESPA = 1<<11;
static const int POST_FORCE_RESPA = 1<<12; static const int POST_INTEGRATE_RESPA = 1<<12;
static const int FINAL_INTEGRATE_RESPA = 1<<13; static const int PRE_FORCE_RESPA = 1<<13;
static const int MIN_PRE_EXCHANGE = 1<<14; static const int POST_FORCE_RESPA = 1<<14;
static const int MIN_PRE_NEIGHBOR = 1<<15; static const int FINAL_INTEGRATE_RESPA = 1<<15;
static const int MIN_PRE_FORCE = 1<<16; static const int MIN_PRE_EXCHANGE = 1<<16;
static const int MIN_POST_FORCE = 1<<17; static const int MIN_PRE_NEIGHBOR = 1<<17;
static const int MIN_ENERGY = 1<<18; static const int MIN_PRE_FORCE = 1<<18;
static const int POST_RUN = 1<<19; static const int MIN_PRE_REVERSE = 1<<19;
static const int PRE_REVERSE = 1<<20; static const int MIN_POST_FORCE = 1<<20;
static const int FIX_CONST_LAST = 1<<21; static const int MIN_ENERGY = 1<<21;
static const int FIX_CONST_LAST = 1<<22;
} }
} }

3
src/fix_restrain.cpp
View File

@ -45,7 +45,8 @@ enum{BOND,ANGLE,DIHEDRAL};
FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) : FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), Fix(lmp, narg, arg),
rstyle(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL), cos_target(NULL), sin_target(NULL) rstyle(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL),
cos_target(NULL), sin_target(NULL)
{ {
if (narg < 4) error->all(FLERR,"Illegal fix restrain command"); if (narg < 4) error->all(FLERR,"Illegal fix restrain command");

13
src/min.cpp
View File

@ -284,7 +284,7 @@ void Min::setup()
else force->kspace->compute_dummy(eflag,vflag); else force->kspace->compute_dummy(eflag,vflag);
} }
modify->pre_reverse(eflag,vflag); modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm(); if (force->newton) comm->reverse_comm();
// update per-atom minimization variables stored by pair styles // update per-atom minimization variables stored by pair styles
@ -365,7 +365,7 @@ void Min::setup_minimal(int flag)
else force->kspace->compute_dummy(eflag,vflag); else force->kspace->compute_dummy(eflag,vflag);
} }
modify->pre_reverse(eflag,vflag); modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm(); if (force->newton) comm->reverse_comm();
// update per-atom minimization variables stored by pair styles // update per-atom minimization variables stored by pair styles
@ -493,6 +493,11 @@ double Min::energy_force(int resetflag)
comm->borders(); comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM); timer->stamp(Timer::COMM);
if (modify->n_min_pre_neighbor) {
timer->stamp();
modify->min_pre_neighbor();
timer->stamp(Timer::MODIFY);
}
neighbor->build(); neighbor->build();
timer->stamp(Timer::NEIGH); timer->stamp(Timer::NEIGH);
} }
@ -525,8 +530,8 @@ double Min::energy_force(int resetflag)
timer->stamp(Timer::KSPACE); timer->stamp(Timer::KSPACE);
} }
if (modify->n_pre_reverse) { if (modify->n_min_pre_reverse) {
modify->pre_reverse(eflag,vflag); modify->min_pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY); timer->stamp(Timer::MODIFY);
} }

38
src/modify.cpp
View File

@ -47,7 +47,8 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
n_thermo_energy_atom = 0; n_thermo_energy_atom = 0;
n_initial_integrate_respa = n_post_integrate_respa = 0; n_initial_integrate_respa = n_post_integrate_respa = 0;
n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0; n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0;
n_min_pre_exchange = n_min_pre_force = n_min_post_force = n_min_energy = 0; n_min_pre_exchange = n_min_pre_force = n_min_pre_reverse = 0;
n_min_post_force = n_min_energy = 0;
fix = NULL; fix = NULL;
fmask = NULL; fmask = NULL;
@ -60,7 +61,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
list_pre_force_respa = list_post_force_respa = NULL; list_pre_force_respa = list_post_force_respa = NULL;
list_final_integrate_respa = NULL; list_final_integrate_respa = NULL;
list_min_pre_exchange = list_min_pre_neighbor = NULL; list_min_pre_exchange = list_min_pre_neighbor = NULL;
list_min_pre_force = list_min_post_force = NULL; list_min_pre_force = list_min_pre_reverse = list_min_post_force = NULL;
list_min_energy = NULL; list_min_energy = NULL;
end_of_step_every = NULL; end_of_step_every = NULL;
@ -136,6 +137,7 @@ Modify::~Modify()
delete [] list_min_pre_exchange; delete [] list_min_pre_exchange;
delete [] list_min_pre_neighbor; delete [] list_min_pre_neighbor;
delete [] list_min_pre_force; delete [] list_min_pre_force;
delete [] list_min_pre_reverse;
delete [] list_min_post_force; delete [] list_min_post_force;
delete [] list_min_energy; delete [] list_min_energy;
@ -188,6 +190,7 @@ void Modify::init()
list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange); list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange);
list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor); list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor);
list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force); list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force);
list_init(MIN_PRE_REVERSE,n_min_pre_reverse,list_min_pre_reverse);
list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force); list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
list_init(MIN_ENERGY,n_min_energy,list_min_energy); list_init(MIN_ENERGY,n_min_energy,list_min_energy);
@ -307,7 +310,7 @@ void Modify::setup_pre_exchange()
fix[list_pre_exchange[i]]->setup_pre_exchange(); fix[list_pre_exchange[i]]->setup_pre_exchange();
else if (update->whichflag == 2) else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_exchange; i++) for (int i = 0; i < n_min_pre_exchange; i++)
fix[list_min_pre_exchange[i]]->min_setup_pre_exchange(); fix[list_min_pre_exchange[i]]->setup_pre_exchange();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -322,7 +325,7 @@ void Modify::setup_pre_neighbor()
fix[list_pre_neighbor[i]]->setup_pre_neighbor(); fix[list_pre_neighbor[i]]->setup_pre_neighbor();
else if (update->whichflag == 2) else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_neighbor; i++) for (int i = 0; i < n_min_pre_neighbor; i++)
fix[list_min_pre_neighbor[i]]->min_setup_pre_neighbor(); fix[list_min_pre_neighbor[i]]->setup_pre_neighbor();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -337,7 +340,22 @@ void Modify::setup_pre_force(int vflag)
fix[list_pre_force[i]]->setup_pre_force(vflag); fix[list_pre_force[i]]->setup_pre_force(vflag);
else if (update->whichflag == 2) else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_force; i++) for (int i = 0; i < n_min_pre_force; i++)
fix[list_min_pre_force[i]]->min_setup_pre_force(vflag); fix[list_min_pre_force[i]]->setup_pre_force(vflag);
}
/* ----------------------------------------------------------------------
setup pre_reverse call, only for fixes that define pre_reverse
called from Verlet, RESPA, Min
------------------------------------------------------------------------- */
void Modify::setup_pre_reverse(int eflag, int vflag)
{
if (update->whichflag == 1)
for (int i = 0; i < n_pre_reverse; i++)
fix[list_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_reverse; i++)
fix[list_min_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -580,6 +598,16 @@ void Modify::min_pre_force(int vflag)
fix[list_min_pre_force[i]]->min_pre_force(vflag); fix[list_min_pre_force[i]]->min_pre_force(vflag);
} }
/* ----------------------------------------------------------------------
minimizer pre-reverse call, only for relevant fixes
------------------------------------------------------------------------- */
void Modify::min_pre_reverse(int eflag, int vflag)
{
for (int i = 0; i < n_min_pre_reverse; i++)
fix[list_min_pre_reverse[i]]->min_pre_reverse(eflag,vflag);
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
minimizer force adjustment call, only for relevant fixes minimizer force adjustment call, only for relevant fixes
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

6
src/modify.h
View File

@ -35,7 +35,7 @@ class Modify : protected Pointers {
int n_initial_integrate_respa,n_post_integrate_respa; int n_initial_integrate_respa,n_post_integrate_respa;
int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa; int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa;
int n_min_pre_exchange,n_min_pre_neighbor; int n_min_pre_exchange,n_min_pre_neighbor;
int n_min_pre_force,n_min_post_force,n_min_energy; int n_min_pre_force,n_min_pre_reverse,n_min_post_force,n_min_energy;
int restart_pbc_any; // 1 if any fix sets restart_pbc int restart_pbc_any; // 1 if any fix sets restart_pbc
int nfix_restart_global; // stored fix global info from restart file int nfix_restart_global; // stored fix global info from restart file
@ -54,6 +54,7 @@ class Modify : protected Pointers {
virtual void setup_pre_exchange(); virtual void setup_pre_exchange();
virtual void setup_pre_neighbor(); virtual void setup_pre_neighbor();
virtual void setup_pre_force(int); virtual void setup_pre_force(int);
virtual void setup_pre_reverse(int, int);
virtual void initial_integrate(int); virtual void initial_integrate(int);
virtual void post_integrate(); virtual void post_integrate();
virtual void pre_exchange(); virtual void pre_exchange();
@ -78,6 +79,7 @@ class Modify : protected Pointers {
virtual void min_pre_exchange(); virtual void min_pre_exchange();
virtual void min_pre_neighbor(); virtual void min_pre_neighbor();
virtual void min_pre_force(int); virtual void min_pre_force(int);
virtual void min_pre_reverse(int,int);
virtual void min_post_force(int); virtual void min_post_force(int);
virtual double min_energy(double *); virtual double min_energy(double *);
@ -124,7 +126,7 @@ class Modify : protected Pointers {
int *list_pre_force_respa,*list_post_force_respa; int *list_pre_force_respa,*list_post_force_respa;
int *list_final_integrate_respa; int *list_final_integrate_respa;
int *list_min_pre_exchange,*list_min_pre_neighbor; int *list_min_pre_exchange,*list_min_pre_neighbor;
int *list_min_pre_force,*list_min_post_force; int *list_min_pre_force,*list_min_pre_reverse,*list_min_post_force;
int *list_min_energy; int *list_min_energy;
int *end_of_step_every; int *end_of_step_every;

10
src/respa.cpp
View File

@ -44,8 +44,10 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg), Respa::Respa(LAMMPS *lmp, int narg, char **arg) :
step(NULL), loop(NULL), hybrid_level(NULL), hybrid_compute(NULL), newton(NULL), fix_respa(NULL) Integrate(lmp, narg, arg),
step(NULL), loop(NULL), hybrid_level(NULL), hybrid_compute(NULL),
newton(NULL), fix_respa(NULL)
{ {
nhybrid_styles = 0; nhybrid_styles = 0;
if (narg < 1) error->all(FLERR,"Illegal run_style respa command"); if (narg < 1) error->all(FLERR,"Illegal run_style respa command");
@ -473,7 +475,7 @@ void Respa::setup()
if (kspace_compute_flag) force->kspace->compute(eflag,vflag); if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
} }
modify->pre_reverse(eflag,vflag); modify->setup_pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm(); if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel); copy_f_flevel(ilevel);
} }
@ -549,7 +551,7 @@ void Respa::setup_minimal(int flag)
if (kspace_compute_flag) force->kspace->compute(eflag,vflag); if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
} }
modify->pre_reverse(eflag,vflag); modify->setup_pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm(); if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel); copy_f_flevel(ilevel);
} }

1
src/run.cpp
View File

@ -43,6 +43,7 @@ void Run::command(int narg, char **arg)
error->all(FLERR,"Run command before simulation box is defined"); error->all(FLERR,"Run command before simulation box is defined");
// ignore run command, if walltime limit was already reached // ignore run command, if walltime limit was already reached
if (timer->is_timeout()) return; if (timer->is_timeout()) return;
bigint nsteps_input = force->bnumeric(FLERR,arg[0]); bigint nsteps_input = force->bnumeric(FLERR,arg[0]);

4
src/verlet.cpp
View File

@ -144,7 +144,7 @@ void Verlet::setup()
else force->kspace->compute_dummy(eflag,vflag); else force->kspace->compute_dummy(eflag,vflag);
} }
modify->pre_reverse(eflag,vflag); modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm(); if (force->newton) comm->reverse_comm();
modify->setup(vflag); modify->setup(vflag);
@ -205,7 +205,7 @@ void Verlet::setup_minimal(int flag)
else force->kspace->compute_dummy(eflag,vflag); else force->kspace->compute_dummy(eflag,vflag);
} }
modify->pre_reverse(eflag,vflag); modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm(); if (force->newton) comm->reverse_comm();
modify->setup(vflag); modify->setup(vflag);