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

changes to angle and improper amoeba terms

Browse Source
This commit is contained in:
Steve Plimpton
2021-12-15 11:44:38 -07:00
parent 161fdec540
commit 67f7e44688
6 changed files with 246 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files

146
src/AMOEBA/angle_amoeba.cpp
View File

@ -41,6 +41,11 @@ AngleAmoeba::AngleAmoeba(LAMMPS *lmp) : Angle(lmp)
k4 = nullptr;
k5 = nullptr;
k6 = nullptr;
ba_k1 = nullptr;
ba_k2 = nullptr;
ba_r1 = nullptr;
ba_r2 = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -51,6 +56,8 @@ AngleAmoeba::~AngleAmoeba()
if (allocated) {
memory->destroy(setflag);
memory->destroy(setflag_a);
memory->destroy(setflag_ba);
memory->destroy(pflag);
memory->destroy(theta0);
@ -59,6 +66,11 @@ AngleAmoeba::~AngleAmoeba()
memory->destroy(k4);
memory->destroy(k5);
memory->destroy(k6);
memory->destroy(ba_k1);
memory->destroy(ba_k2);
memory->destroy(ba_r1);
memory->destroy(ba_r2);
}
}
@ -72,6 +84,7 @@ void AngleAmoeba::compute(int eflag, int vflag)
double dtheta,dtheta2,dtheta3,dtheta4,dtheta5,dtheta6,de_angle;
double dr1,dr2,tk1,tk2,aa1,aa2,aa11,aa12,aa21,aa22;
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22,b1,b2;
double vx11,vx12,vy11,vy12,vz11,vz12,vx21,vx22,vy21,vy22,vz21,vz22;
eangle = 0.0;
ev_init(eflag,vflag);
@ -159,6 +172,50 @@ void AngleAmoeba::compute(int eflag, int vflag)
if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 +
k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
// force & energy for bond-angle term
// bond-stretch cross term in Tinker
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;
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) {
@ -382,8 +439,13 @@ void AngleAmoeba::allocate()
memory->create(k5,n+1,"angle:k5");
memory->create(k6,n+1,"angle:k6");
memory->create(ba_k1,n+1,"angle:ba_k1");
memory->create(ba_k2,n+1,"angle:ba_k2");
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");
for (int i = 1; i <= n; i++) setflag[i] = 0;
for (int i = 1; i <= n; i++) setflag[i] = setflag_a[i] = setflag_ba[i] = 0;
}
/* ----------------------------------------------------------------------
@ -392,6 +454,7 @@ void AngleAmoeba::allocate()
void AngleAmoeba::coeff(int narg, char **arg)
{
if (narg < 2) error->all(FLERR,"Incorrect args for angle coefficients");
if (!allocated) allocate();
int ilo,ihi;
@ -399,32 +462,52 @@ void AngleAmoeba::coeff(int narg, char **arg)
int count = 0;
if (narg != 8) error->all(FLERR,"Incorrect args for angle coefficients");
if (strcmp(arg[1],"ba") == 0) {
if (narg != 6) error->all(FLERR,"Incorrect args for angle coefficients");
int pflag_one = utils::inumeric(FLERR,arg[1],false,lmp);
double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
double k2_one = utils::numeric(FLERR,arg[3],false,lmp);
double k3_one = utils::numeric(FLERR,arg[4],false,lmp);
double k4_one = utils::numeric(FLERR,arg[5],false,lmp);
double k5_one = utils::numeric(FLERR,arg[6],false,lmp);
double k6_one = utils::numeric(FLERR,arg[7],false,lmp);
double ba_k1_one = utils::numeric(FLERR,arg[2],false,lmp);
double ba_k2_one = utils::numeric(FLERR,arg[3],false,lmp);
double ba_r1_one = utils::numeric(FLERR,arg[4],false,lmp);
double ba_r2_one = utils::numeric(FLERR,arg[5],false,lmp);
// convert theta0 from degrees to radians
for (int i = ilo; i <= ihi; i++) {
ba_k1[i] = ba_k1_one;
ba_k2[i] = ba_k2_one;
ba_r1[i] = ba_r1_one;
ba_r2[i] = ba_r2_one;
setflag_ba[i] = 1;
count++;
}
} else {
if (narg != 8) error->all(FLERR,"Incorrect args for angle coefficients");
int pflag_one = utils::inumeric(FLERR,arg[1],false,lmp);
double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
double k2_one = utils::numeric(FLERR,arg[3],false,lmp);
double k3_one = utils::numeric(FLERR,arg[4],false,lmp);
double k4_one = utils::numeric(FLERR,arg[5],false,lmp);
double k5_one = utils::numeric(FLERR,arg[6],false,lmp);
double k6_one = utils::numeric(FLERR,arg[7],false,lmp);
// convert theta0 from degrees to radians
for (int i = ilo; i <= ihi; i++) {
pflag[i] = pflag_one;
theta0[i] = theta0_one/180.0 * MY_PI;
k2[i] = k2_one;
k3[i] = k3_one;
k4[i] = k4_one;
k5[i] = k5_one;
k6[i] = k6_one;
count++;
for (int i = ilo; i <= ihi; i++) {
pflag[i] = pflag_one;
theta0[i] = theta0_one/180.0 * MY_PI;
k2[i] = k2_one;
k3[i] = k3_one;
k4[i] = k4_one;
k5[i] = k5_one;
k6[i] = k6_one;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
for (int i = ilo; i <= ihi; i++) setflag[i] = 1;
for (int i = ilo; i <= ihi; i++)
if (setflag_a[i] == 1 && setflag_ba[i] == 1) setflag[i] = 1;
}
/* ---------------------------------------------------------------------- */
@ -447,6 +530,11 @@ void AngleAmoeba::write_restart(FILE *fp)
fwrite(&k4[1],sizeof(double),atom->nangletypes,fp);
fwrite(&k5[1],sizeof(double),atom->nangletypes,fp);
fwrite(&k6[1],sizeof(double),atom->nangletypes,fp);
fwrite(&ba_k1[1],sizeof(double),atom->nangletypes,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);
}
/* ----------------------------------------------------------------------
@ -465,6 +553,11 @@ void AngleAmoeba::read_restart(FILE *fp)
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);
}
MPI_Bcast(&pflag[1],atom->nangletypes,MPI_INT,0,world);
@ -475,6 +568,11 @@ void AngleAmoeba::read_restart(FILE *fp)
MPI_Bcast(&k5[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&k6[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&ba_k1[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&ba_k2[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&ba_r1[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&ba_r2[1],atom->nangletypes,MPI_DOUBLE,0,world);
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
}
@ -487,6 +585,10 @@ 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,"\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]);
}
/* ---------------------------------------------------------------------- */
@ -526,5 +628,9 @@ double AngleAmoeba::single(int type, int i1, int i2, int i3)
double energy = k2[type]*dtheta2 + k3[type]*dtheta3 + k4[type]*dtheta4
+ k5[type]*dtheta5 + k6[type]*dtheta6;
double dr1 = r1 - ba_r1[type];
double dr2 = r2 - ba_r2[type];
energy += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
return energy;
}

2
src/AMOEBA/angle_amoeba.h
View File

@ -39,6 +39,8 @@ class AngleAmoeba : public Angle {
protected:
int *pflag;
double *theta0, *k2, *k3, *k4, *k5, *k6;
double *ba_k1, *ba_k2, *ba_r1, *ba_r2;
int *setflag_a, *setflag_ba;
void anglep(int, int, int, int, int);
void allocate();

133
src/AMOEBA/improper_amoeba.cpp
View File

@ -17,9 +17,10 @@
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "update.h"
#include "neighbor.h"
#include "force.h"
#include "update.h"
#include "pair.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
@ -56,7 +57,13 @@ void ImproperAmoeba::compute(int eflag, int vflag)
double xab,yab,zab,xcb,ycb,zcb,xdb,ydb,zdb,xad,yad,zad,xcd,ycd,zcd;
double rad2,rcd2,rdb2,dot,cc,ee;
double sine,angle;
double eimproper,f1[3],f2[3],f3[3],f4[3];
double dt,dt2,dt3,dt4,e;
double deddt,sign,dedcos,term;
double dccdxia,dccdyia,dccdzia,dccdxic,dccdyic,dccdzic;
double dccdxid,dccdyid,dccdzid;
double deedxia,deedyia,deedzia,deedxic,deedyic,deedzic;
double deedxid,deedyid,deedzid;
double eimproper,fa[3],fb[3],fc[3],fd[3];
eimproper = 0.0;
ev_init(eflag,vflag);
@ -129,48 +136,91 @@ void ImproperAmoeba::compute(int eflag, int vflag)
sine = fabs(ee) / sqrt(cc*rdb2);
sine = MIN(1.0,sine);
angle = asin(sine) * MY_PI/180.0;
angle = 180.0/MY_PI * asin(sine); // angle = degrees
//dt = angle;
//dt2 = dt * dt;
//dt3 = dt2 * dt;
//dt4 = dt2 * dt2;
//e = opbunit * force * dt2 * (1.0d0+copb*dt+qopb*dt2+popb*dt3+sopb*dt4);
//deddt = opbunit * force * dt * radian *
// (2.0d0 + 3.0d0*copb*dt + 4.0d0*qopb*dt2 + 5.0d0*popb*dt3 + 6.0d0*sopb*dt4);
//dedcos = -deddt * sign(1.0d0,ee) / sqrt(cc*rdb2-ee*ee);
dt = angle;
dt2 = dt * dt;
dt3 = dt2 * dt;
dt4 = dt2 * dt2;
double prefactor = 1.0 / (180.0/MY_PI) / (180.0/MY_PI);
e = prefactor * k[type] * dt2 * (1.0 + opbend_cubic*dt + opbend_quartic*dt2 +
opbend_pentic*dt3 + opbend_sextic*dt4);
eimproper += e;
deddt = k[type] * dt *
(2.0 + 3.0*opbend_cubic*dt + 4.0*opbend_quartic*dt2 +
5.0*opbend_pentic*dt3 + 6.0*opbend_sextic*dt4);
sign = (ee >= 0.0) ? 1.0 : -1.0;
dedcos = -deddt * sign / sqrt(cc*rdb2 - ee*ee);
// chain rule terms for first derivative components
term = ee / cc;
dccdxia = (xad*rcd2-xcd*dot) * term;
dccdyia = (yad*rcd2-ycd*dot) * term;
dccdzia = (zad*rcd2-zcd*dot) * term;
dccdxic = (xcd*rad2-xad*dot) * term;
dccdyic = (ycd*rad2-yad*dot) * term;
dccdzic = (zcd*rad2-zad*dot) * term;
dccdxid = -dccdxia - dccdxic;
dccdyid = -dccdyia - dccdyic;
dccdzid = -dccdzia - dccdzic;
term = ee / rdb2;
deedxia = ydb*zcb - zdb*ycb;
deedyia = zdb*xcb - xdb*zcb;
deedzia = xdb*ycb - ydb*xcb;
deedxic = yab*zdb - zab*ydb;
deedyic = zab*xdb - xab*zdb;
deedzic = xab*ydb - yab*xdb;
deedxid = ycb*zab - zcb*yab + xdb*term;
deedyid = zcb*xab - xcb*zab + ydb*term;
deedzid = xcb*yab - ycb*xab + zdb*term;
// compute first derivative components for this angle
fa[0] = dedcos * (dccdxia+deedxia);
fa[1] = dedcos * (dccdyia+deedyia);
fa[2] = dedcos * (dccdzia+deedzia);
fc[0] = dedcos * (dccdxic+deedxic);
fc[1] = dedcos * (dccdyic+deedyic);
fc[2] = dedcos * (dccdzic+deedzic);
fd[0] = dedcos * (dccdxid+deedxid);
fd[1] = dedcos * (dccdyid+deedyid);
fd[2] = dedcos * (dccdzid+deedzid);
fb[0] = -fa[0] - fc[0] - fd[0];
fb[1] = -fa[1] - fc[1] - fd[1];
fb[2] = -fa[1] - fc[2] - fd[2];
/*
// apply force to each of 4 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0];
f[i1][1] += f1[1];
f[i1][2] += f1[2];
if (newton_bond || id < nlocal) {
f[id][0] += fd[0];
f[id][1] += fd[1];
f[id][2] += fd[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] += f2[0];
f[i2][1] += f2[1];
f[i2][2] += f2[2];
if (newton_bond || ib < nlocal) {
f[ib][0] += fb[0];
f[ib][1] += fb[1];
f[ib][2] += fb[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0];
f[i3][1] += f3[1];
f[i3][2] += f3[2];
if (newton_bond || ia < nlocal) {
f[ia][0] += fa[0];
f[ia][1] += fa[1];
f[ia][2] += fa[2];
}
if (newton_bond || i4 < nlocal) {
f[i4][0] += f4[0];
f[i4][1] += f4[1];
f[i4][2] += f4[2];
if (newton_bond || ic < nlocal) {
f[ic][0] += fc[0];
f[ic][1] += fc[1];
f[ic][2] += fc[2];
}
if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
*/
ev_tally(id,ib,ia,ic,nlocal,newton_bond,e,fd,fa,fc,
xdb,ydb,zdb,xab,yab,zab,xic-xia,yic-yia,zic-zia);
}
}
@ -213,6 +263,27 @@ void ImproperAmoeba::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients");
}
/* ----------------------------------------------------------------------
set opbend higher-order term weights from PairAmoeba
------------------------------------------------------------------------- */
void ImproperAmoeba::init_style()
{
Pair *pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) error->all(FLERR,"Improper amoeba could not find pair amoega");
int dim;
opbend_cubic = *(double *) pair->extract("opbend_cubic",dim);
opbend_quartic = *(double *) pair->extract("opbend_quartic",dim);
opbend_pentic = *(double *) pair->extract("opbend_pentic",dim);
opbend_sextic = *(double *) pair->extract("opbend_sextic",dim);
printf("OPBEND %g %g %g %g\n",
opbend_cubic,opbend_quartic,opbend_pentic,opbend_sextic);
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */

2
src/AMOEBA/improper_amoeba.h
View File

@ -30,11 +30,13 @@ class ImproperAmoeba : public Improper {
~ImproperAmoeba();
void compute(int, int);
void coeff(int, char **);
void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
void write_data(FILE *);
protected:
double opbend_cubic,opbend_quartic,opbend_pentic,opbend_sextic;
double *k;
virtual void allocate();

12
src/AMOEBA/pair_amoeba.cpp
View File

@ -2107,6 +2107,18 @@ void PairAmoeba::zero_energy_force_virial()
}
}
/* ---------------------------------------------------------------------- */
void *PairAmoeba::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"opbend_cubic") == 0) return (void *) &opbend_cubic;
if (strcmp(str,"opbend_quartic") == 0) return (void *) &opbend_quartic;
if (strcmp(str,"opbend_pentic") == 0) return (void *) &opbend_pentic;
if (strcmp(str,"opbend_sextic") == 0) return (void *) &opbend_sextic;
return nullptr;
}
/* ----------------------------------------------------------------------
grow local vectors and arrays if necessary
keep them atom->nmax in length

2
src/AMOEBA/pair_amoeba.h
View File

@ -57,6 +57,8 @@ class PairAmoeba : public Pair {
void pack_reverse_grid(int, void *, int, int *);
void unpack_reverse_grid(int, void *, int, int *);
void *extract(const char *, int &);
protected:
int me,nprocs;
int nmax; // allocation for owned+ghost