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

angles issue with angle vs anglep

Browse Source
This commit is contained in:
Steve Plimpton
2021-10-19 16:10:51 -06:00
parent 7eba439388
commit a2f62ae2db
2 changed files with 6 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
src/AMOEBA/angle_amoeba.cpp
View File

@ -79,6 +79,8 @@ void AngleAmoeba::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
int **anglelist = neighbor->anglelist;
int **nspecial = atom->nspecial;
int nanglelist = neighbor->nanglelist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
@ -90,11 +92,11 @@ void AngleAmoeba::compute(int eflag, int vflag)
type = anglelist[n][3];
// tflag = 0 for "angle", 1 for "anglep" in Tinker PRM file
// atom 2 must have 3 bond partners to invoke "anglep" option
// atom 2 must have 3 bond partners to invoke anglep() variant
tflag = pflag[type];
if (tflag && atom->num_bond[i2] == 3) {
if (tflag && nspecial[i2][0] == 3) {
anglep(i1,i2,i3,type,eflag);
continue;
}
@ -203,7 +205,7 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
double **x = atom->x;
double **f = atom->f;
tagint **bond_atom = atom->bond_atom;
tagint **special = atom->special;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
@ -213,7 +215,7 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
i3tag = atom->tag[i3];
for (int ibond = 0; ibond < 3; ibond++) {
i4tag = bond_atom[i2][ibond];
i4tag = special[i2][ibond];
if (i4tag != i1tag && i4tag != i3tag) break;
}
@ -290,15 +292,6 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 +
k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
/*
printf("ANGLEP: %d %d %d %d: eng %g\n",
atom->tag[i1],
atom->tag[i2],
atom->tag[i3],
atom->tag[i4],
eangle);
*/
// chain rule terms for first derivative components
terma = -deddt / (rap2*rm);

18
tools/tinker2lmp.py
View File

@ -580,9 +580,6 @@ atype = []
aparams = []
baparams = []
#DEBUG
opcount = 0
for atom1,atom2,atom3 in alist:
type1 = type[atom1-1]
type2 = type[atom2-1]
@ -609,9 +606,6 @@ for atom1,atom2,atom3 in alist:
# option 2 if one of 2 additional bonds is to an H atom
# option 3 if both of 2 additional bonds is to an H atom
# DEBUG
debugflag = 0
if len(params[3]) == 1:
which = 1
@ -638,10 +632,6 @@ for atom1,atom2,atom3 in alist:
print " Nbonds and hydrogen count:",nbonds,hcount
print " which:",which,m
# DEBUG
debugflag = 1
opcount += 1
elif len(params[3]) == 3:
nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass)
@ -668,24 +658,16 @@ for atom1,atom2,atom3 in alist:
print " Nbonds and hydrogen count:",nbonds,hcount
print " which:",which,m
# DEBUG
debugflag = 1
opcount += 1
else:
print "Angle not found",atom1,atom2,atom3,c1,c2,c3
sys.exit()
if not flags[m]:
pflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
# DEBUG single line
if debugflag: pflag = 2
aparams.append((pflag,v1,v2,v3,v4,v5,v6))
flags[m] = len(aparams)
atype.append(flags[m])
print "OPTION angles",opcount
# NOTE: baparams may need to be flipped if match is 3,2,1 instead of 1,2,3
# NOTE: mismatch between angle and bondangle params may not be handled right