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

move UB bonds to angle amoeba

Browse Source
This commit is contained in:
Steve Plimpton
2022-03-08 13:28:11 -07:00
parent 7087dfc019
commit 64becd5642
6 changed files with 29612 additions and 3135 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3176
examples/amoeba/data.ubiquitin
View File

File diff suppressed because it is too large Load Diff

29069
examples/amoeba/data.ubiquitin.ub Normal file
View File

File diff suppressed because it is too large Load Diff

373
src/AMOEBA/angle_amoeba.cpp
View File

@ -35,6 +35,8 @@ using namespace MathConst;
AngleAmoeba::AngleAmoeba(LAMMPS *lmp) : Angle(lmp)
{
pflag = nullptr;
ubflag = nullptr;
theta0 = nullptr;
k2 = nullptr;
k3 = nullptr;
@ -46,6 +48,9 @@ AngleAmoeba::AngleAmoeba(LAMMPS *lmp) : Angle(lmp)
ba_k2 = nullptr;
ba_r1 = nullptr;
ba_r2 = nullptr;
ub_k = nullptr;
ub_r0 = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -58,8 +63,11 @@ AngleAmoeba::~AngleAmoeba()
memory->destroy(setflag);
memory->destroy(setflag_a);
memory->destroy(setflag_ba);
memory->destroy(setflag_ub);
memory->destroy(pflag);
memory->destroy(ubflag);
memory->destroy(theta0);
memory->destroy(k2);
memory->destroy(k3);
@ -71,6 +79,9 @@ AngleAmoeba::~AngleAmoeba()
memory->destroy(ba_k2);
memory->destroy(ba_r1);
memory->destroy(ba_r2);
memory->destroy(ub_k);
memory->destroy(ub_r0);
}
}
@ -78,7 +89,7 @@ AngleAmoeba::~AngleAmoeba()
void AngleAmoeba::compute(int eflag, int vflag)
{
int i1,i2,i3,n,type,tflag;
int i1,i2,i3,n,type,tflag,uflag;
double delx1,dely1,delz1,delx2,dely2,delz2;
double f1[3],f3[3];
double dtheta,dtheta2,dtheta3,dtheta4,dtheta5,dtheta6,de_angle;
@ -112,127 +123,17 @@ void AngleAmoeba::compute(int eflag, int vflag)
if (ba_k1[type] != 0.0)
tinker_bondangle(i1,i2,i3,type,eflag);
// Urey-Bradley H-H bond term within water molecules
uflag = ubflag[type];
if (uflag) tinker_urey_bradley(i1,i3,type,eflag);
}
}
/* ---------------------------------------------------------------------- */
void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag)
{
double delx1,dely1,delz1,delx2,dely2,delz2;
double rsq1,r1,rsq2,r2,c,s,dtheta;
double dr1,dr2,aa1,aa2,b1,b2;
double aa11,aa12,aa21,aa22;
double vx11,vx12,vy11,vy12,vz11,vz12,vx21,vx22,vy21,vy22,vz21,vz22;
double eangle,f1[3],f3[3];
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
// 1st bond
delx1 = x[i1][0] - x[i2][0];
dely1 = x[i1][1] - x[i2][1];
delz1 = x[i1][2] - x[i2][2];
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[i3][0] - x[i2][0];
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL;
s = 1.0/s;
dtheta = acos(c) - theta0[type];
// force & energy for bond-angle term
dr1 = r1 - ba_r1[type];
dr2 = r2 - ba_r2[type];
aa1 = s * dr1 * ba_k1[type];
aa2 = s * dr2 * ba_k2[type];
aa11 = aa1 * c / rsq1;
aa12 = -aa1 / (r1 * r2);
aa21 = aa2 * c / rsq1;
aa22 = -aa2 / (r1 * r2);
vx11 = (aa11 * delx1) + (aa12 * delx2);
vx12 = (aa21 * delx1) + (aa22 * delx2);
vy11 = (aa11 * dely1) + (aa12 * dely2);
vy12 = (aa21 * dely1) + (aa22 * dely2);
vz11 = (aa11 * delz1) + (aa12 * delz2);
vz12 = (aa21 * delz1) + (aa22 * delz2);
aa11 = aa1 * c / rsq2;
aa21 = aa2 * c / rsq2;
vx21 = (aa11 * delx2) + (aa12 * delx1);
vx22 = (aa21 * delx2) + (aa22 * delx1);
vy21 = (aa11 * dely2) + (aa12 * dely1);
vy22 = (aa21 * dely2) + (aa22 * dely1);
vz21 = (aa11 * delz2) + (aa12 * delz1);
vz22 = (aa21 * delz2) + (aa22 * delz1);
b1 = ba_k1[type] * dtheta / r1;
b2 = ba_k2[type] * dtheta / r2;
f1[0] -= vx11 + b1*delx1 + vx12;
f1[1] -= vy11 + b1*dely1 + vy12;
f1[2] -= vz11 + b1*delz1 + vz12;
f3[0] -= vx21 + b2*delx2 + vx22;
f3[1] -= vy21 + b2*dely2 + vy22;
f3[2] -= vz21 + b2*delz2 + vz22;
eangle = 0.0;
if (eflag) eangle = ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
// apply force to each of 3 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0];
f[i1][1] += f1[1];
f[i1][2] += f1[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= f1[0] + f3[0];
f[i2][1] -= f1[1] + f3[1];
f[i2][2] -= f1[2] + f3[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0];
f[i3][1] += f3[1];
f[i3][2] += f3[2];
}
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2);
}
/* ---------------------------------------------------------------------- */
void AngleAmoeba::tinker_angle(int i1, int i2, int i3, int type, int eflag)
{
double delx1,dely1,delz1,delx2,dely2,delz2;
@ -516,12 +417,178 @@ void AngleAmoeba::tinker_anglep(int i1, int i2, int i3, int type, int eflag)
/* ---------------------------------------------------------------------- */
void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag)
{
double delx1,dely1,delz1,delx2,dely2,delz2;
double rsq1,r1,rsq2,r2,c,s,dtheta;
double dr1,dr2,aa1,aa2,b1,b2;
double aa11,aa12,aa21,aa22;
double vx11,vx12,vy11,vy12,vz11,vz12,vx21,vx22,vy21,vy22,vz21,vz22;
double eangle,f1[3],f3[3];
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
// 1st bond
delx1 = x[i1][0] - x[i2][0];
dely1 = x[i1][1] - x[i2][1];
delz1 = x[i1][2] - x[i2][2];
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[i3][0] - x[i2][0];
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL;
s = 1.0/s;
dtheta = acos(c) - theta0[type];
// force & energy for bond-angle term
dr1 = r1 - ba_r1[type];
dr2 = r2 - ba_r2[type];
aa1 = s * dr1 * ba_k1[type];
aa2 = s * dr2 * ba_k2[type];
aa11 = aa1 * c / rsq1;
aa12 = -aa1 / (r1 * r2);
aa21 = aa2 * c / rsq1;
aa22 = -aa2 / (r1 * r2);
vx11 = (aa11 * delx1) + (aa12 * delx2);
vx12 = (aa21 * delx1) + (aa22 * delx2);
vy11 = (aa11 * dely1) + (aa12 * dely2);
vy12 = (aa21 * dely1) + (aa22 * dely2);
vz11 = (aa11 * delz1) + (aa12 * delz2);
vz12 = (aa21 * delz1) + (aa22 * delz2);
aa11 = aa1 * c / rsq2;
aa21 = aa2 * c / rsq2;
vx21 = (aa11 * delx2) + (aa12 * delx1);
vx22 = (aa21 * delx2) + (aa22 * delx1);
vy21 = (aa11 * dely2) + (aa12 * dely1);
vy22 = (aa21 * dely2) + (aa22 * dely1);
vz21 = (aa11 * delz2) + (aa12 * delz1);
vz22 = (aa21 * delz2) + (aa22 * delz1);
b1 = ba_k1[type] * dtheta / r1;
b2 = ba_k2[type] * dtheta / r2;
f1[0] -= vx11 + b1*delx1 + vx12;
f1[1] -= vy11 + b1*dely1 + vy12;
f1[2] -= vz11 + b1*delz1 + vz12;
f3[0] -= vx21 + b2*delx2 + vx22;
f3[1] -= vy21 + b2*dely2 + vy22;
f3[2] -= vz21 + b2*delz2 + vz22;
eangle = 0.0;
if (eflag) eangle = ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
// apply force to each of 3 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0];
f[i1][1] += f1[1];
f[i1][2] += f1[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= f1[0] + f3[0];
f[i2][1] -= f1[1] + f3[1];
f[i2][2] -= f1[2] + f3[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0];
f[i3][1] += f3[1];
f[i3][2] += f3[2];
}
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2);
}
/* ---------------------------------------------------------------------- */
void AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag)
{
double delx,dely,delz;
double rsq,r,dr,rk;
double fbond,ebond;
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
delx = x[i1][0] - x[i2][0];
dely = x[i1][1] - x[i2][1];
delz = x[i1][2] - x[i2][2];
rsq = delx*delx + dely*dely + delz*delz;
r = sqrt(rsq);
dr = r - ub_r0[type];
rk = ub_k[type] * dr;
// force & energy
if (r > 0.0) fbond = -2.0*rk/r;
else fbond = 0.0;
if (eflag) ebond = rk*dr;
// apply force to each of 2 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fbond;
f[i1][1] += dely*fbond;
f[i1][2] += delz*fbond;
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= delx*fbond;
f[i2][1] -= dely*fbond;
f[i2][2] -= delz*fbond;
}
// NOTE: there is no force on i2 for UB
//if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,ebond,f1,f3,
// delx1,dely1,delz1,delx2,dely2,delz2);
}
/* ---------------------------------------------------------------------- */
void AngleAmoeba::allocate()
{
allocated = 1;
int n = atom->nangletypes;
memory->create(pflag,n+1,"angle:pflag");
memory->create(ubflag,n+1,"angle:ubflag");
memory->create(theta0,n+1,"angle:theta0");
memory->create(k2,n+1,"angle:k2");
memory->create(k3,n+1,"angle:k3");
@ -534,11 +601,16 @@ void AngleAmoeba::allocate()
memory->create(ba_r1,n+1,"angle:ba_r1");
memory->create(ba_r2,n+1,"angle:ba_r2");
memory->create(setflag,n+1,"angle:setflag");
memory->create(setflag_a,n+1,"angle:setflag");
memory->create(setflag_ba,n+1,"angle:setflag");
memory->create(ub_k,n+1,"angle:ub_k");
memory->create(ub_r0,n+1,"angle:ub_r0");
for (int i = 1; i <= n; i++) setflag[i] = setflag_a[i] = setflag_ba[i] = 0;
memory->create(setflag,n+1,"angle:setflag");
memory->create(setflag_a,n+1,"angle:setflag_a");
memory->create(setflag_ba,n+1,"angle:setflag_ba");
memory->create(setflag_ub,n+1,"angle:setflag_ub");
for (int i = 1; i <= n; i++)
setflag[i] = setflag_a[i] = setflag_ba[i] = setflag_ub[i] = 0;
}
/* ----------------------------------------------------------------------
@ -572,6 +644,19 @@ void AngleAmoeba::coeff(int narg, char **arg)
count++;
}
} else if (strcmp(arg[1],"ub") == 0) {
if (narg != 4) error->all(FLERR,"Incorrect args for angle coefficients");
double ub_k_one = utils::numeric(FLERR,arg[2],false,lmp);
double ub_r0_one = utils::numeric(FLERR,arg[3],false,lmp);
for (int i = ilo; i <= ihi; i++) {
ub_k[i] = ub_k_one;
ub_r0[i] = ub_r0_one;
setflag_ub[i] = 1;
count++;
}
} else {
if (narg != 8) error->all(FLERR,"Incorrect args for angle coefficients");
@ -601,7 +686,8 @@ void AngleAmoeba::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
for (int i = ilo; i <= ihi; i++)
if (setflag_a[i] == 1 && setflag_ba[i] == 1) setflag[i] = 1;
if (setflag_a[i] == 1 && setflag_ba[i] == 1 && setflag_ub[i])
setflag[i] = 1;
}
/* ---------------------------------------------------------------------- */
@ -618,6 +704,8 @@ double AngleAmoeba::equilibrium_angle(int i)
void AngleAmoeba::write_restart(FILE *fp)
{
fwrite(&pflag[1],sizeof(int),atom->nangletypes,fp);
fwrite(&ubflag[1],sizeof(int),atom->nangletypes,fp);
fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
fwrite(&k2[1],sizeof(double),atom->nangletypes,fp);
fwrite(&k3[1],sizeof(double),atom->nangletypes,fp);
@ -629,6 +717,9 @@ void AngleAmoeba::write_restart(FILE *fp)
fwrite(&ba_k2[1],sizeof(double),atom->nangletypes,fp);
fwrite(&ba_r1[1],sizeof(double),atom->nangletypes,fp);
fwrite(&ba_r2[1],sizeof(double),atom->nangletypes,fp);
fwrite(&ub_k[1],sizeof(double),atom->nangletypes,fp);
fwrite(&ub_r0[1],sizeof(double),atom->nangletypes,fp);
}
/* ----------------------------------------------------------------------
@ -641,20 +732,34 @@ void AngleAmoeba::read_restart(FILE *fp)
if (comm->me == 0) {
utils::sfread(FLERR,&pflag[1],sizeof(int),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&ubflag[1],sizeof(int),atom->nangletypes,
fp,nullptr,error);
utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,
fp,nullptr,error);
utils::sfread(FLERR,&k2[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&k3[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&k4[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&k5[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&k6[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&ba_k1[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&ba_k2[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&ba_r1[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&ba_r2[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&ba_k1[1],sizeof(double),atom->nangletypes,
fp,nullptr,error);
utils::sfread(FLERR,&ba_k2[1],sizeof(double),atom->nangletypes,
fp,nullptr,error);
utils::sfread(FLERR,&ba_r1[1],sizeof(double),atom->nangletypes,
fp,nullptr,error);
utils::sfread(FLERR,&ba_r2[1],sizeof(double),atom->nangletypes,
fp,nullptr,error);
utils::sfread(FLERR,&ub_k[1],sizeof(double),atom->nangletypes,
fp,nullptr,error);
utils::sfread(FLERR,&ub_r0[1],sizeof(double),atom->nangletypes,
fp,nullptr,error);
}
MPI_Bcast(&pflag[1],atom->nangletypes,MPI_INT,0,world);
MPI_Bcast(&ubflag[1],atom->nangletypes,MPI_INT,0,world);
MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&k2[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&k3[1],atom->nangletypes,MPI_DOUBLE,0,world);
@ -667,6 +772,9 @@ void AngleAmoeba::read_restart(FILE *fp)
MPI_Bcast(&ba_r1[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&ba_r2[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&ub_k[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&ub_r0[1],atom->nangletypes,MPI_DOUBLE,0,world);
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
}
@ -677,12 +785,17 @@ void AngleAmoeba::read_restart(FILE *fp)
void AngleAmoeba::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp,"%d %d %g %g %g %g %g %g\n",
i,pflag[i],theta0[i]/MY_PI*180.0,k2[i],k3[i],k4[i],k5[i],k6[i]);
fprintf(fp,"%d %d %d %g %g %g %g %g %g\n",
i,pflag[i],ubflag[i],theta0[i]/MY_PI*180.0,
k2[i],k3[i],k4[i],k5[i],k6[i]);
fprintf(fp,"\nBondAngle Coeffs\n\n");
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp,"%d %g %g %g %g\n",i,ba_k1[i],ba_k2[i],ba_r1[i],ba_r2[i]);
fprintf(fp,"\nUreyBradley Coeffs\n\n");
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp,"%d %g %g\n",i,ub_k[i],ub_r0[i]);
}
/* ---------------------------------------------------------------------- */
@ -726,5 +839,7 @@ double AngleAmoeba::single(int type, int i1, int i2, int i3)
double dr2 = r2 - ba_r2[type];
energy += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
// NOTE: add UB term
return energy;
}

6
src/AMOEBA/angle_amoeba.h
View File

@ -37,14 +37,16 @@ class AngleAmoeba : public Angle {
double single(int, int, int, int);
protected:
int *pflag;
int *pflag, *ubflag;
double *theta0, *k2, *k3, *k4, *k5, *k6;
double *ba_k1, *ba_k2, *ba_r1, *ba_r2;
int *setflag_a, *setflag_ba;
double *ub_k, *ub_r0;
int *setflag_a, *setflag_ba, *setflag_ub;
void tinker_angle(int, int, int, int, int);
void tinker_anglep(int, int, int, int, int);
void tinker_bondangle(int, int, int, int, int);
void tinker_urey_bradley(int, int, int, int);
void allocate();
};

1
tools/tinker/data.py
View File

@ -372,6 +372,7 @@ skeywords = [["Masses","atom types"],
["PiTorsion Coeffs","pitorsion types"],
["BondBond Coeffs","angle types"],
["BondAngle Coeffs","angle types"],
["UreyBradley Coeffs","angle types"],
["MiddleBondTorsion Coeffs","dihedral types"],
["EndBondTorsion Coeffs","dihedral types"],
["AngleTorsion Coeffs","dihedral types"],

122
tools/tinker/tinker2lmp.py
View File

@ -206,6 +206,7 @@ class PRMfile:
def angles(self):
r2d = 180.0 / math.pi
ubflag = 0
params = []
iline = self.find_section("Angle Bending Parameters")
if iline < 0: return params
@ -236,17 +237,17 @@ class PRMfile:
lmp5 = self.angle_pentic * value1 * r2d*r2d*r2d
lmp6 = self.angle_sextic * value1 * r2d*r2d*r2d*r2d
option1 = (pflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
option1 = (pflag,ubflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
if len(words) >= 7:
value3 = float(words[6])
lmp1 = value3
option2 = (pflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
option2 = (pflag,ubflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
if len(words) == 8:
value4 = float(words[7])
lmp1 = value4
option3 = (pflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
option3 = (pflag,ubflag,lmp1,lmp2,lmp3,lmp4,lmp5,lmp6)
if not option2 and not option3:
params.append((class1,class2,class3,[option1]))
@ -272,7 +273,7 @@ class PRMfile:
bdict = {}
for m,bparams in enumerate(self.bondparams):
bdict[(bparams[0],bparams[1])] = bparams[2]
while iline < self.nlines:
words = self.lines[iline].split()
if len(words):
@ -372,7 +373,7 @@ class PRMfile:
iline += 1
return params
# convert PRMfile params to LAMMPS bond_style class2 bond params
# convert PRMfile params to LAMMPS angle_style amoeba UB params
# coeffs for K2,K3 = 0.0 since Urey-Bradley is simple harmonic
def ureybonds(self):
@ -394,7 +395,7 @@ class PRMfile:
lmp1 = value1
lmp2 = value2
params.append((class1,class2,class3,lmp1,lmp2,0.0,0.0))
params.append((class1,class2,class3,lmp1,lmp2))
iline += 1
return params
@ -715,6 +716,10 @@ for atom1 in id:
pitorsionlist.append((atom3,atom4,atom1,atom2,atom5,atom6))
# DEBUG - if uncommented, no PiTorsions appear in data file
pitorsionlist = []
# ----------------------------------------
# create lists of bond/angle/dihedral/improper types
# ----------------------------------------
@ -754,44 +759,6 @@ for atom1,atom2 in blist:
bparams.append((v1,v2,v3,v4))
flags[m] = len(bparams)
btype.append(flags[m])
# extend btype and bparams with Urey-Bradley H-H bond info
# loop over ublist of UB bonds
# convert prm.ureyparams to a dictionary for efficient searching
# key = (class1,class2,class3)
# value = (M,params) where M is index into prm.ureyparams
# also add the c1,c3 H-H bond to blist
id = xyz.id
type = xyz.type
classes = prm.classes
ubdict = {}
for m,params in enumerate(prm.ureyparams):
ubdict[(params[0],params[1],params[2])] = (m,params)
flags = len(prm.ureyparams)*[0]
for atom1,atom2,atom3 in ublist:
type1 = type[atom1-1]
type2 = type[atom2-1]
type3 = type[atom3-1]
c1 = classes[type1-1]
c2 = classes[type2-1]
c3 = classes[type3-1]
if (c1,c2,c3) in ubdict:
m,params = ubdict[(c1,c2,c3)]
blist.append((atom1,atom3))
elif (c3,c2,c1) in ubdict:
m,params = ubdict[(c3,c2,c1)]
blist.append((atom3,atom1))
if not flags[m]:
v1,v2,v3,v4 = params[3:]
bparams.append((v1,v2,v3,v4))
flags[m] = len(bparams)
btype.append(flags[m])
# generate atype = LAMMPS type of each angle
# generate aparams = LAMMPS params for each angle type
@ -905,19 +872,19 @@ for i,one in enumerate(alist):
error("Angle not found: %d %d %d: %d %d %d" % (atom1,atom2,atom3,c1,c2,c3))
if not flags[m]:
pflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
aparams.append((pflag,v1,v2,v3,v4,v5,v6))
pflag,ubflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
aparams.append((pflag,ubflag,v1,v2,v3,v4,v5,v6))
flags[m] = len(aparams)
atype.append(flags[m])
# augment the aparams with bond-angle cross terms from bondangleparams
# generate baparams = LAMMPS bond-angle params for each angle type
# sbdict = dictionary for angle tuples in bongangleparams
# badict = dictionary for angle tuples in bongangleparams
sbdict = {}
badict = {}
for v1,v2,v3,v4,v5,v6,v7 in prm.bondangleparams:
if (v1,v2,v3) in sbdict: continue
sbdict[(v1,v2,v3)] = (v4,v5,v6,v7)
if (v1,v2,v3) in badict: continue
badict[(v1,v2,v3)] = (v4,v5,v6,v7)
baparams = []
@ -931,16 +898,54 @@ for itype in range(len(aparams)):
c2 = classes[type2-1]
c3 = classes[type3-1]
if (c1,c2,c3) in sbdict:
n1,n2,r1,r2 = sbdict[(c1,c2,c3)]
elif (c3,c2,c1) in sbdict:
n1,n2,r1,r2 = sbdict[(c3,c2,c1)]
if (c1,c2,c3) in badict:
n1,n2,r1,r2 = badict[(c1,c2,c3)]
elif (c3,c2,c1) in badict:
n1,n2,r1,r2 = badict[(c3,c2,c1)]
else:
print "Bond-stretch angle triplet not found: %d %d %d" % (c1,c2,c3)
n1,n2,r1,r2 = 4*[0.0]
baparams.append((n1,n2,r1,r2))
# augment the aparams with Urey_Bradley terms from ureyparams
# generate ubparams = LAMMPS UB params for 1-3 bond in each angle type
# ubdict = dictionary for angle tuples in ureyparams
ubdict = {}
for v1,v2,v3,v4,v5 in prm.ureyparams:
if (v1,v2,v3) in ubdict: continue
ubdict[(v1,v2,v3)] = (v4,v5)
ubparams = []
for itype in range(len(aparams)):
iangle = atype.index(itype+1)
atom1,atom2,atom3 = alist[iangle]
type1 = type[atom1-1]
type2 = type[atom2-1]
type3 = type[atom3-1]
c1 = classes[type1-1]
c2 = classes[type2-1]
c3 = classes[type3-1]
# if UB settings exist for this angle type, set ubflag in aparams to 1
if (c1,c2,c3) in ubdict:
r1,r2 = ubdict[(c1,c2,c3)]
pflag,ubflag,v1,v2,v3,v4,v5,v6 = aparams[itype]
ubflag = 1
aparams[itype] = (pflag,ubflag,v1,v2,v3,v4,v5,v6)
elif (c3,c2,c1) in ubdict:
r1,r2 = ubdict[(c3,c2,c1)]
pflag,ubflag,v1,v2,v3,v4,v5,v6 = aparams[itype]
ubflag = 1
aparams[itype] = (pflag,ubflag,v1,v2,v3,v4,v5,v6)
else:
r1,r2 = 2*[0.0]
ubparams.append((r1,r2))
# generate dtype = LAMMPS type of each dihedral
# generate dparams = LAMMPS params for each dihedral type
# flags[i] = which LAMMPS dihedral type (1-N) the Tinker FF file dihedral I is
@ -1075,7 +1080,7 @@ for tmp1,tmp2,atom1,atom2,tmp3,tmp4 in pitorsionlist:
print "PRM pitor param",len(prm.pitorsionparams)
print "PiTor list",len(pitorsionlist)
print "Flags",flags
#print "Flags",flags
# ----------------------------------------
# assign each atom to a Tinker group
@ -1202,6 +1207,13 @@ if nangles:
lines.append(line+'\n')
d.sections["BondAngle Coeffs"] = lines
lines = []
for i,one in enumerate(ubparams):
strone = [str(single) for single in one]
line = "%d %s" % (i+1,' '.join(strone))
lines.append(line+'\n')
d.sections["UreyBradley Coeffs"] = lines
lines = []
for i,one in enumerate(alist):
line = "%d %d %d %d %d" % (i+1,atype[i],one[0],one[1],one[2])