make energy handling consistent for variable field

This commit is contained in:
Axel Kohlmeyer
2023-03-11 18:56:29 -05:00
parent 251fac2c60
commit a8c1359c54
4 changed files with 12 additions and 34 deletions

View File

@ -152,9 +152,9 @@ void FixEfieldTIP4P::post_force(int vflag)
}
if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR, "TIP4P hydrogen is missing");
if (iO == -1) error->one(FLERR, "TIP4P oxygen is missing");
if ((atom->type[iH1] != typeH) || (atom->type[iH2] != typeH))
if ((type[iH1] != typeH) || (type[iH2] != typeH))
error->one(FLERR, "TIP4P hydrogen has incorrect atom type");
if (atom->type[iO] != typeO) error->one(FLERR, "TIP4P oxygen has incorrect atom type");
if (type[iO] != typeO) error->one(FLERR, "TIP4P oxygen has incorrect atom type");
iH1 = domain->closest_image(i, iH1);
iH2 = domain->closest_image(i, iH2);
@ -399,9 +399,9 @@ void FixEfieldTIP4P::post_force(int vflag)
}
if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR, "TIP4P hydrogen is missing");
if (iO == -1) error->one(FLERR, "TIP4P oxygen is missing");
if ((atom->type[iH1] != typeH) || (atom->type[iH2] != typeH))
if ((type[iH1] != typeH) || (type[iH2] != typeH))
error->one(FLERR, "TIP4P hydrogen has incorrect atom type");
if (atom->type[iO] != typeO) error->one(FLERR, "TIP4P oxygen has incorrect atom type");
if (type[iO] != typeO) error->one(FLERR, "TIP4P oxygen has incorrect atom type");
iH1 = domain->closest_image(i, iH1);
iH2 = domain->closest_image(i, iH2);
@ -459,11 +459,7 @@ void FixEfieldTIP4P::post_force(int vflag)
}
if (i == iO) {
domain->unmap(xM, image[iO], unwrap);
if (estyle == ATOM)
fsum[0] += efield[0][3];
else
fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2];
if (estyle == ATOM) fsum[0] += efield[0][3];
fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;
@ -509,11 +505,7 @@ void FixEfieldTIP4P::post_force(int vflag)
// tally global force
domain->unmap(x[iH1], image[iH1], unwrap);
if (estyle == ATOM)
fsum[0] += efield[0][3];
else
fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2];
if (estyle == ATOM) fsum[0] += efield[0][3];
fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;
@ -559,11 +551,7 @@ void FixEfieldTIP4P::post_force(int vflag)
// tally global force
domain->unmap(x[iH2], image[iH2], unwrap);
if (estyle == ATOM)
fsum[0] += efield[0][3];
else
fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2];
if (estyle == ATOM) fsum[0] += efield[0][3];
fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;
@ -611,12 +599,7 @@ void FixEfieldTIP4P::post_force(int vflag)
f[i][2] += fz;
fsum[3] += fz;
if (estyle == ATOM) {
fsum[0] += efield[0][3];
} else {
domain->unmap(x[i], image[i], unwrap);
fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2];
}
if (estyle == ATOM) fsum[0] += efield[0][3];
}
}
}

View File

@ -34,7 +34,7 @@ class FixEfieldTIP4P : public FixEfield {
protected:
double alpha;
int typeO, typeH; // atom types for TIP4P molecule
int typeO, typeH; // atom types for TIP4P molecule
void find_M(double *, double *, double *, double *);
};
} // namespace LAMMPS_NS

View File

@ -387,12 +387,7 @@ void FixEfield::post_force(int vflag)
}
f[i][2] += fz;
fsum[3] += fz;
if (estyle == ATOM) {
fsum[0] += efield[0][3];
} else {
domain->unmap(x[i], image[i], unwrap);
fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2];
}
if (estyle == ATOM) fsum[0] += efield[0][3];
}
}

View File

@ -1,6 +1,6 @@
---
lammps_version: 8 Feb 2023
date_generated: Sat Mar 11 18:36:57 2023
date_generated: Sat Mar 11 18:52:03 2023
epsilon: 2e-13
skip_tests:
prerequisites: ! |
@ -24,7 +24,7 @@ post_commands: ! |
fix test all efield/tip4p v_xforce 0.0 v_zforce
input_file: in.fourmol
natoms: 29
global_scalar: -108.06823025980506
global_scalar: 0
global_vector: ! |-
3 9.596230365887404 0 -4.440892098500626e-16
run_pos: ! |2