turn off terms for both amoeba and hippo

This commit is contained in:
Steve Plimpton
2022-04-08 16:51:08 -06:00
parent 2111797ed8
commit d0af0fa456
10 changed files with 45 additions and 30 deletions

View File

@ -21,7 +21,7 @@ read_data data.water_dimer.amoeba fix amtype NULL "Tinker Types"
# force field # force field
pair_style amoeba pair_style amoeba include angle
pair_coeff * * amoeba_water.prm amoeba_water.key pair_coeff * * amoeba_water.prm amoeba_water.key
special_bonds lj/coul 0.5 0.5 0.5 one/five yes special_bonds lj/coul 0.5 0.5 0.5 one/five yes
@ -35,10 +35,10 @@ thermo_style custom step temp epair ebond eangle edihed eimp &
# zero step run # zero step run
#run 0 run 0
# dynamics # dynamics
thermo 1 #thermo 1
fix 1 all nve #fix 1 all nve
run 100 #run 100

View File

@ -21,7 +21,7 @@ read_data data.water_dimer.hippo fix amtype NULL "Tinker Types"
# force field # force field
pair_style hippo pair_style hippo include angle
pair_coeff * * hippo_water.prm hippo_water.key pair_coeff * * hippo_water.prm hippo_water.key
special_bonds lj/coul 0.5 0.5 0.5 one/five yes special_bonds lj/coul 0.5 0.5 0.5 one/five yes
@ -35,10 +35,10 @@ thermo_style custom step temp epair ebond eangle edihed eimp &
# zero step run # zero step run
#run 0 run 0
# dynamics # dynamics
thermo 1 #thermo 1
fix 1 all nve #fix 1 all nve
run 100 #run 100

View File

@ -812,6 +812,7 @@ void PairAmoeba::multipole_kspace()
f[i][2] -= h3; f[i][2] -= h3;
} }
empole += 0.5*e; empole += 0.5*e;
//printf("mpole_force %g %g %g \n", f[0][0], f[0][1], f[0][2]);
// augment the permanent multipole virial contributions // augment the permanent multipole virial contributions

View File

@ -697,9 +697,11 @@ void AngleAmoeba::coeff(int narg, char **arg)
void AngleAmoeba::init_style() void AngleAmoeba::init_style()
{ {
// check if PairAmoeba disabled angle or Urey-Bradley terms // check if PairAmoeba or PairHippo disabled angle or Urey-Bradley terms
Pair *pair = force->pair_match("amoeba",1,0); Pair *pair = NULL;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) enable_angle = enable_urey = 1; if (!pair) enable_angle = enable_urey = 1;
else { else {

View File

@ -185,9 +185,11 @@ void FixAmoebaBiTorsion::init()
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
} }
// check if PairAmoeba disabled bitorsion terms // check if PairAmoeba or PairHippo disabled bitorsion terms
Pair *pair = force->pair_match("amoeba",1,0); Pair *pair = NULL;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) disable = 0; if (!pair) disable = 0;
else { else {

View File

@ -151,9 +151,11 @@ void FixAmoebaPiTorsion::init()
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
} }
// check if PairAmoeba disabled pitorsion terms // check if PairAmoeba or PairHippo disabled pitorsion terms
Pair *pair = force->pair_match("amoeba",1,0); Pair *pair = NULL;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) disable = 0; if (!pair) disable = 0;
else { else {

View File

@ -276,21 +276,24 @@ void ImproperAmoeba::coeff(int narg, char **arg)
void ImproperAmoeba::init_style() void ImproperAmoeba::init_style()
{ {
Pair *pair = force->pair_match("amoeba",1,0); // check if PairAmoeba disabled improper terms
Pair *pair = NULL;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("hippo",1,0); if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) error->all(FLERR,"Improper amoeba could not find pair amoega"); if (!pair) error->all(FLERR,"Improper amoeba could not find pair amoeba/hippo");
int tmp;
int flag = *((int *) pair->extract("improper_flag",tmp));
disable = flag ? 0 : 1;
// also extract opbend params
int dim; int dim;
opbend_cubic = *(double *) pair->extract("opbend_cubic",dim); opbend_cubic = *(double *) pair->extract("opbend_cubic",dim);
opbend_quartic = *(double *) pair->extract("opbend_quartic",dim); opbend_quartic = *(double *) pair->extract("opbend_quartic",dim);
opbend_pentic = *(double *) pair->extract("opbend_pentic",dim); opbend_pentic = *(double *) pair->extract("opbend_pentic",dim);
opbend_sextic = *(double *) pair->extract("opbend_sextic",dim); opbend_sextic = *(double *) pair->extract("opbend_sextic",dim);
// check if PairAmoeba disabled improper terms
int tmp;
int flag = *((int *) pair->extract("improper_flag",tmp));
disable = flag ? 0 : 1;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

View File

@ -161,9 +161,11 @@ void BondClass2::coeff(int narg, char **arg)
void BondClass2::init_style() void BondClass2::init_style()
{ {
// check if PairAmoeba disabled bond terms // check if PairAmoeba or PairHippo disabled bond terms
Pair *pair = force->pair_match("amoeba",1,0); Pair *pair = NULL;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) disable = 0; if (!pair) disable = 0;
else { else {

View File

@ -333,9 +333,11 @@ void DihedralFourier::coeff(int narg, char **arg)
void DihedralFourier::init_style() void DihedralFourier::init_style()
{ {
// check if PairAmoeba disabled dihedral terms // check if PairAmoeba or PairHippo disabled dihedral terms
Pair *pair = force->pair_match("amoeba",1,0); Pair *pair = NULL;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) disable = 0; if (!pair) disable = 0;
else { else {

View File

@ -860,12 +860,13 @@ for i,one in enumerate(alist):
# IMPORTANT subtlety # IMPORTANT subtlety
# flip order of 3 atoms in alist if the angle # flip order of 3 atoms in alist if the angle
# matches Angle Bending section of PRM file in reverse order # matches Angle Bending section of PRM file in reverse order
# no need to flip if c1 = c3
# necessary b/c BondAngle coeffs will be generated with r1,r2 params # necessary b/c BondAngle coeffs will be generated with r1,r2 params
# from Bond Stretching section of PRM file # from Bond Stretching section of PRM file
# since in general r1 != r2, the LAMMPS AngleAmoeba class requires # since in general r1 != r2, the LAMMPS AngleAmoeba class requires
# the 3 atoms in the angle be in the order that matches r1 and r2 # the 3 atoms in the angle be in the order that matches r1 and r2
if (c3,c2,c1) in adict: if c1 != c3 and (c3,c2,c1) in adict:
m,params = adict[(c3,c2,c1)] m,params = adict[(c3,c2,c1)]
alist[i] = (atom3,atom2,atom1) alist[i] = (atom3,atom2,atom1)