From 0658844e04ad83f814a43e08308d6ad827732d2c Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 10 Mar 2022 16:00:05 -0700 Subject: [PATCH] first working version of fix pitorsion --- examples/amoeba/in.ubiquitin | 1 + src/AMOEBA/amoeba_dispersion.cpp | 4 +- src/AMOEBA/amoeba_induce.cpp | 12 +- src/AMOEBA/amoeba_multipole.cpp | 4 +- src/AMOEBA/amoeba_polar.cpp | 4 +- src/AMOEBA/angle_amoeba.cpp | 41 +++- src/AMOEBA/angle_amoeba.h | 3 + src/AMOEBA/fix_pitorsion.cpp | 305 ++++++++++++++---------- src/AMOEBA/fix_pitorsion.h | 2 + src/AMOEBA/improper_amoeba.cpp | 12 +- src/AMOEBA/improper_amoeba.h | 1 + src/AMOEBA/pair_amoeba.cpp | 131 +++++----- src/AMOEBA/pair_amoeba.h | 12 +- src/CLASS2/bond_class2.cpp | 20 +- src/CLASS2/bond_class2.h | 2 + src/EXTRA-MOLECULE/dihedral_fourier.cpp | 21 +- src/EXTRA-MOLECULE/dihedral_fourier.h | 2 + 17 files changed, 354 insertions(+), 223 deletions(-) diff --git a/examples/amoeba/in.ubiquitin b/examples/amoeba/in.ubiquitin index 367030ba74..6f7d605af0 100644 --- a/examples/amoeba/in.ubiquitin +++ b/examples/amoeba/in.ubiquitin @@ -15,6 +15,7 @@ improper_style amoeba fix amtype all property/atom i_amtype ghost yes fix pitorsion all pitorsion +fix_modify pitorsion energy yes fix extra all property/atom & i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes diff --git a/src/AMOEBA/amoeba_dispersion.cpp b/src/AMOEBA/amoeba_dispersion.cpp index bf9b3be11a..d19bceb14a 100644 --- a/src/AMOEBA/amoeba_dispersion.cpp +++ b/src/AMOEBA/amoeba_dispersion.cpp @@ -44,11 +44,11 @@ void PairAmoeba::dispersion() // compute the real space portion of the Ewald summation - if (rspace_flag) dispersion_real(); + if (disp_rspace_flag) dispersion_real(); // compute the reciprocal space part of the Ewald summation - if (kspace_flag) dispersion_kspace(); + if (disp_kspace_flag) dispersion_kspace(); // compute the self-energy portion of the Ewald summation diff --git a/src/AMOEBA/amoeba_induce.cpp b/src/AMOEBA/amoeba_induce.cpp index c8f361053c..e644f7759a 100644 --- a/src/AMOEBA/amoeba_induce.cpp +++ b/src/AMOEBA/amoeba_induce.cpp @@ -40,6 +40,8 @@ enum{GORDON1,GORDON2}; #define DEBYE 4.80321 // conversion factor from q-Angs (real units) to Debye +#define UIND_DEBUG 0 + /* ---------------------------------------------------------------------- induce = induced dipole moments via pre-conditioned CG solver adapted from Tinker induce0a() routine @@ -539,7 +541,7 @@ void PairAmoeba::induce() // DEBUG output to dump file - if (uind_flag) + if (UIND_DEBUG) dump6(fp_uind,"id uindx uindy uindz uinpx uinpy uinpz",DEBYE,uind,uinp); // deallocation of arrays @@ -715,11 +717,11 @@ void PairAmoeba::ufield0c(double **field, double **fieldp) // get the reciprocal space part of the mutual field - if (kspace_flag) umutual1(field,fieldp); + if (polar_kspace_flag) umutual1(field,fieldp); // get the real space portion of the mutual field - if (rspace_flag) umutual2b(field,fieldp); + if (polar_rspace_flag) umutual2b(field,fieldp); // add the self-energy portion of the mutual field @@ -983,7 +985,7 @@ void PairAmoeba::dfield0c(double **field, double **fieldp) // get the reciprocal space part of the permanent field - if (kspace_flag) udirect1(field); + if (polar_kspace_flag) udirect1(field); for (i = 0; i < nlocal; i++) { for (j = 0; j < 3; j++) { @@ -993,7 +995,7 @@ void PairAmoeba::dfield0c(double **field, double **fieldp) // get the real space portion of the permanent field - if (rspace_flag) udirect2b(field,fieldp); + if (polar_rspace_flag) udirect2b(field,fieldp); // get the self-energy portion of the permanent field diff --git a/src/AMOEBA/amoeba_multipole.cpp b/src/AMOEBA/amoeba_multipole.cpp index c06f07d70c..894a1df835 100644 --- a/src/AMOEBA/amoeba_multipole.cpp +++ b/src/AMOEBA/amoeba_multipole.cpp @@ -85,12 +85,12 @@ void PairAmoeba::multipole() // compute the real space part of the Ewald summation - if (rspace_flag) multipole_real(); + if (mpole_rspace_flag) multipole_real(); // compute the reciprocal space part of the Ewald summation //printf("zero check %e \n",empole); - if (kspace_flag) multipole_kspace(); + if (mpole_kspace_flag) multipole_kspace(); //printf("kspace energy %e \n",empole); // compute the Ewald self-energy term over all the atoms diff --git a/src/AMOEBA/amoeba_polar.cpp b/src/AMOEBA/amoeba_polar.cpp index 1503243220..83a07caa2e 100644 --- a/src/AMOEBA/amoeba_polar.cpp +++ b/src/AMOEBA/amoeba_polar.cpp @@ -80,11 +80,11 @@ void PairAmoeba::polar() // compute the real space part of the dipole interactions - if (rspace_flag) polar_real(); + if (polar_rspace_flag) polar_real(); // compute the reciprocal space part of dipole interactions - if (kspace_flag) polar_kspace(); + if (polar_kspace_flag) polar_kspace(); // compute the Ewald self-energy torque and virial terms diff --git a/src/AMOEBA/angle_amoeba.cpp b/src/AMOEBA/angle_amoeba.cpp index e81bbff00f..216a79ff83 100644 --- a/src/AMOEBA/angle_amoeba.cpp +++ b/src/AMOEBA/angle_amoeba.cpp @@ -21,6 +21,7 @@ #include "domain.h" #include "comm.h" #include "force.h" +#include "pair.h" #include "math_const.h" #include "memory.h" #include "error.h" @@ -112,22 +113,26 @@ void AngleAmoeba::compute(int eflag, int vflag) // tflag = 0 for "angle", 1 for "anglep" in Tinker PRM file // atom 2 must have exactly 3 bond partners to invoke anglep() variant - tflag = pflag[type]; + if (enable_angle) { + tflag = pflag[type]; - if (tflag && nspecial[i2][0] == 3) - tinker_anglep(i1,i2,i3,type,eflag); - else - tinker_angle(i1,i2,i3,type,eflag); + if (tflag && nspecial[i2][0] == 3) + tinker_anglep(i1,i2,i3,type,eflag); + else + tinker_angle(i1,i2,i3,type,eflag); - // bondangle = bond-stretch cross term in Tinker + // bondangle = bond-stretch cross term in Tinker - if (ba_k1[type] != 0.0) - tinker_bondangle(i1,i2,i3,type,eflag); + 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); + if (enable_urey) { + uflag = ubflag[type]; + if (uflag) tinker_urey_bradley(i1,i3,type,eflag); + } } } @@ -690,6 +695,22 @@ void AngleAmoeba::coeff(int narg, char **arg) /* ---------------------------------------------------------------------- */ +void AngleAmoeba::init_style() +{ + // check if PairAmoeba disabled angle or Urey-Bradley terms + + Pair *pair = force->pair_match("amoeba",1,0); + + if (!pair) enable_angle = enable_urey = 1; + else { + int tmp; + enable_angle = *((int *) pair->extract("angle_flag",tmp)); + enable_urey = *((int *) pair->extract("urey_flag",tmp)); + } +} + +/* ---------------------------------------------------------------------- */ + double AngleAmoeba::equilibrium_angle(int i) { return theta0[i]; diff --git a/src/AMOEBA/angle_amoeba.h b/src/AMOEBA/angle_amoeba.h index 69a782b48b..05d8041286 100644 --- a/src/AMOEBA/angle_amoeba.h +++ b/src/AMOEBA/angle_amoeba.h @@ -30,6 +30,7 @@ class AngleAmoeba : public Angle { virtual ~AngleAmoeba(); void compute(int, int); void coeff(int, char **); + void init_style(); double equilibrium_angle(int); void write_restart(FILE *); void read_restart(FILE *); @@ -43,6 +44,8 @@ class AngleAmoeba : public Angle { double *ub_k, *ub_r0; int *setflag_a, *setflag_ba, *setflag_ub; + int enable_angle,enable_urey; + void tinker_angle(int, int, int, int, int); void tinker_anglep(int, int, int, int, int); void tinker_bondangle(int, int, int, int, int); diff --git a/src/AMOEBA/fix_pitorsion.cpp b/src/AMOEBA/fix_pitorsion.cpp index ebfb6b37ff..b474c34202 100644 --- a/src/AMOEBA/fix_pitorsion.cpp +++ b/src/AMOEBA/fix_pitorsion.cpp @@ -22,6 +22,7 @@ #include "respa.h" #include "domain.h" #include "force.h" +#include "pair.h" #include "comm.h" #include "math_const.h" #include "memory.h" @@ -145,6 +146,21 @@ void FixPiTorsion::init() ilevel_respa = ((Respa *) update->integrate)->nlevels-1; if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); } + + // check if PairAmoeba disabled pitorsion terms + + Pair *pair = force->pair_match("amoeba",1,0); + + if (!pair) disable = 0; + else { + int tmp; + int flag = *((int *) pair->extract("pitorsion_flag",tmp)); + disable = flag ? 0 : 1; + } + + // constant + + onesixth = 1.0/6.0; } /* --------------------------------------------------------------------- */ @@ -185,7 +201,9 @@ void FixPiTorsion::min_setup(int vflag) } /* ---------------------------------------------------------------------- - store local neighbor list for newton_bond ON + create local list of pitorsions + if one or more atoms in pitorsion are on this proc, + this proc lists the pitorsion exactly once ------------------------------------------------------------------------- */ void FixPiTorsion::pre_neighbor() @@ -231,6 +249,13 @@ void FixPiTorsion::pre_neighbor() atom5 = domain->closest_image(i,atom5); atom6 = domain->closest_image(i,atom6); + /* + if (atom1 >= nlocal || atom2 >= nlocal || atom3 >= nlocal || + atom4 >= nlocal || atom5 >= nlocal || atom6 >= nlocal) + printf("ATOM 123456: %d %d %d %d %d %d\n", + atom1,atom2,atom3,atom4,atom5,atom6); + */ + if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4 && i <= atom5 && i <= atom6) { if (npitorsion_list == max_pitorsion_list) { @@ -261,13 +286,15 @@ void FixPiTorsion::pre_reverse(int eflag, int /*vflag*/) } /* ---------------------------------------------------------------------- - compute PiTorsion terms as if newton_bond = OFF, even if actually ON + compute PiTorsion terms ------------------------------------------------------------------------- */ void FixPiTorsion::post_force(int vflag) { + if (disable) return; + int ia,ib,ic,id,ie,ig,ptype; - double e,dedphi; + double e,dedphi,engfraction; double xt,yt,zt,rt2; double xu,yu,zu,ru2; double xtu,ytu,ztu; @@ -307,18 +334,24 @@ void FixPiTorsion::post_force(int vflag) double vxx,vyy,vzz; double vyx,vzx,vzy; - double **x = atom->x; - double **f = atom->f; + int nlist,list[6]; + double v[6]; epitorsion = 0.0; + int eflag = eflag_caller; + ev_init(eflag,vflag); + + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; for (int n = 0; n < npitorsion_list; n++) { ia = pitorsion_list[n][0]; - ib = pitorsion_list[n][0]; - ic = pitorsion_list[n][0]; - id = pitorsion_list[n][0]; - ie = pitorsion_list[n][0]; - ig = pitorsion_list[n][0]; + ib = pitorsion_list[n][1]; + ic = pitorsion_list[n][2]; + id = pitorsion_list[n][3]; + ie = pitorsion_list[n][4]; + ig = pitorsion_list[n][5]; ptype = pitorsion_list[n][6]; // compute the value of the pi-system torsion angle @@ -404,12 +437,16 @@ void FixPiTorsion::post_force(int vflag) dphi2 = 2.0 * (cosine2*s2 - sine2*c2); // calculate pi-system torsion energy and master chain rule term + // NOTE: remove ptornunit if 1.0 ? - // NOTE: is ptornunit 1.0 ? double ptorunit = 1.0; e = ptorunit * v2 * phi2; dedphi = ptorunit * v2 * dphi2; + // fraction of energy for each atom + + engfraction = e * onesixth; + // chain rule terms for first derivative components xdp = xid - xip; @@ -461,47 +498,73 @@ void FixPiTorsion::post_force(int vflag) dedyid = dedyid + dedyiq - dedyia - dedyib; dedzid = dedzid + dedziq - dedzia - dedzib; - // increment the total pi-system torsion energy and gradient + // forces and energy - epitorsion += e; + if (ia < nlocal) { + epitorsion += engfraction; + f[ia][0] += dedxia; + f[ia][1] += dedyia; + f[ia][2] += dedzia; + } - f[ia][0] += dedxia; - f[ia][1] += dedyia; - f[ia][2] += dedzia; - f[ib][0] += dedxib; - f[ib][1] += dedyib; - f[ib][2] += dedzib; - f[ic][0] += dedxic; - f[ic][1] += dedyic; - f[ic][2] += dedzic; - f[id][0] += dedxid; - f[id][1] += dedyid; - f[id][2] += dedzid; - f[ie][0] += dedxie; - f[ie][1] += dedyie; - f[ie][2] += dedzie; - f[ig][0] += dedxig; - f[ig][1] += dedyig; - f[ig][2] += dedzig; + if (ib < nlocal) { + epitorsion += engfraction; + f[ib][0] += dedxib; + f[ib][1] += dedyib; + f[ib][2] += dedzib; + } - // increment the internal virial tensor components + if (ic < nlocal) { + epitorsion += engfraction; + f[ic][0] += dedxic; + f[ic][1] += dedyic; + f[ic][2] += dedzic; + } - vxterm = dedxid + dedxia + dedxib; - vyterm = dedyid + dedyia + dedyib; - vzterm = dedzid + dedzia + dedzib; - vxx = xdc*vxterm + xcp*dedxip - xqd*dedxiq; - vyx = ydc*vxterm + ycp*dedxip - yqd*dedxiq; - vzx = zdc*vxterm + zcp*dedxip - zqd*dedxiq; - vyy = ydc*vyterm + ycp*dedyip - yqd*dedyiq; - vzy = zdc*vyterm + zcp*dedyip - zqd*dedyiq; - vzz = zdc*vzterm + zcp*dedzip - zqd*dedziq; + if (id < nlocal) { + epitorsion += engfraction; + f[id][0] += dedxid; + f[id][1] += dedyid; + f[id][2] += dedzid; + } - virial[0] += vxx; - virial[1] += vyy; - virial[2] += vzz; - virial[3] += vyx; - virial[4] += vzx; - virial[5] += vzy; + if (ie < nlocal) { + epitorsion += engfraction; + f[ie][0] += dedxie; + f[ie][1] += dedyie; + f[ie][2] += dedzie; + } + + if (ig < nlocal) { + epitorsion += engfraction; + f[ig][0] += dedxig; + f[ig][1] += dedyig; + f[ig][2] += dedzig; + } + + // virial tensor components + + if (evflag) { + nlist = 0; + if (ia < nlocal) list[nlist++] = ia; + if (ib < nlocal) list[nlist++] = ib; + if (ic < nlocal) list[nlist++] = ic; + if (id < nlocal) list[nlist++] = id; + if (ie < nlocal) list[nlist++] = ie; + if (ig < nlocal) list[nlist++] = ig; + + vxterm = dedxid + dedxia + dedxib; + vyterm = dedyid + dedyia + dedyib; + vzterm = dedzid + dedzia + dedzib; + v[0] = xdc*vxterm + xcp*dedxip - xqd*dedxiq; + v[1] = ydc*vyterm + ycp*dedyip - yqd*dedyiq; + v[2] = zdc*vzterm + zcp*dedzip - zqd*dedziq; + v[3] = ydc*vxterm + ycp*dedxip - yqd*dedxiq; + v[4] = zdc*vxterm + zcp*dedxip - zqd*dedxiq; + v[5] = zdc*vyterm + zcp*dedyip - zqd*dedyiq; + + ev_tally(nlist,list,6.0,e,v); + } } } @@ -543,17 +606,11 @@ void FixPiTorsion::read_data_header(char *line) } else if (strstr(line,"pitorsion types")) { sscanf(line,"%d",&npitorsion_types); } else error->all(FLERR,"Invalid read data header line for fix pitorsion"); - - // didn't set in constructor because this fix could be defined - // before newton command - - newton_bond = force->newton_bond; } /* ---------------------------------------------------------------------- unpack N lines in buf from section of data file labeled by keyword id_offset is applied to atomID fields if multiple data files are read - store PiTorsion interactions as if newton_bond = OFF, even if actually ON ------------------------------------------------------------------------- */ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf, @@ -583,7 +640,7 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf, for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; - sscanf(buf,"%d %g",&itype,&value); + sscanf(buf,"%d %lg",&itype,&value); if (itype <= 0 || itype > npitorsion_types) error->all(FLERR,"Incorrect args for pitorsion coefficients"); kpit[itype] = value; @@ -593,7 +650,7 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf, // loop over lines of PiTorsions // tokenize the line into values - // add PiTorsion to one or more of my atoms, depending on newton_bond + // each atom in the PiTorsion stores it if (which == 0) { @@ -623,8 +680,32 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf, atom5 += id_offset; atom6 += id_offset; - // atom3 owns the pitorsion when newton_bond on - + if ((m = atom->map(atom1)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + + if ((m = atom->map(atom2)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + if ((m = atom->map(atom3)) >= 0) { if (num_pitorsion[m] == PITORSIONMAX) error->one(FLERR,"Too many PiTorsions for one atom"); @@ -638,81 +719,45 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf, num_pitorsion[m]++; } - if (newton_bond == 0) { - if ((m = atom->map(atom1)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } + if ((m = atom->map(atom4)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; } - - if (newton_bond == 0) { - if ((m = atom->map(atom2)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } + + if ((m = atom->map(atom5)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; } - - if (newton_bond == 0) { - if ((m = atom->map(atom4)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } + + if ((m = atom->map(atom6)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; } - - if (newton_bond == 0) { - if ((m = atom->map(atom5)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } - } - - if (newton_bond == 0) { - if ((m = atom->map(atom6)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } - } - + buf = next + 1; } } @@ -750,7 +795,7 @@ void FixPiTorsion::write_data_section_size(int mth, int &nx, int &ny) // PiTorsions section // nx = # of PiTorsions owned by my local atoms - // if newton_bond off, atom only owns PiTorsion if it is atom3 + // only atom3 owns the PiTorsion in this context // ny = 7 columns = type + 6 atom IDs if (mth == 0) { diff --git a/src/AMOEBA/fix_pitorsion.h b/src/AMOEBA/fix_pitorsion.h index 7e3258c66a..718f6407d7 100644 --- a/src/AMOEBA/fix_pitorsion.h +++ b/src/AMOEBA/fix_pitorsion.h @@ -68,9 +68,11 @@ class FixPiTorsion : public Fix { int nprocs, me; int newton_bond, eflag_caller; int ilevel_respa; + int disable; bigint npitorsions; int npitorsion_types; double epitorsion; + double onesixth; double *kpit; diff --git a/src/AMOEBA/improper_amoeba.cpp b/src/AMOEBA/improper_amoeba.cpp index 0f6bb6705c..93a9b412cd 100644 --- a/src/AMOEBA/improper_amoeba.cpp +++ b/src/AMOEBA/improper_amoeba.cpp @@ -52,6 +52,8 @@ ImproperAmoeba::~ImproperAmoeba() void ImproperAmoeba::compute(int eflag, int vflag) { + if (disable) return; + int ia,ib,ic,id,n,type; double xia,yia,zia,xib,yib,zib,xic,yic,zic,xid,yid,zid; double xab,yab,zab,xcb,ycb,zcb,xdb,ydb,zdb,xad,yad,zad,xcd,ycd,zcd; @@ -63,9 +65,8 @@ void ImproperAmoeba::compute(int eflag, int vflag) 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]; + double fa[3],fb[3],fc[3],fd[3]; - eimproper = 0.0; ev_init(eflag,vflag); double **x = atom->x; @@ -152,7 +153,6 @@ void ImproperAmoeba::compute(int eflag, int vflag) dt4 = dt2 * dt2; 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 + @@ -285,6 +285,12 @@ void ImproperAmoeba::init_style() opbend_quartic = *(double *) pair->extract("opbend_quartic",dim); opbend_pentic = *(double *) pair->extract("opbend_pentic",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; } /* ---------------------------------------------------------------------- diff --git a/src/AMOEBA/improper_amoeba.h b/src/AMOEBA/improper_amoeba.h index b25a60ac47..2c0cbaab6d 100644 --- a/src/AMOEBA/improper_amoeba.h +++ b/src/AMOEBA/improper_amoeba.h @@ -36,6 +36,7 @@ class ImproperAmoeba : public Improper { void write_data(FILE *); protected: + int disable; double opbend_cubic,opbend_quartic,opbend_pentic,opbend_sextic; double *k; diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp index 474a764f22..446ca490a1 100644 --- a/src/AMOEBA/pair_amoeba.cpp +++ b/src/AMOEBA/pair_amoeba.cpp @@ -46,7 +46,9 @@ enum{MPOLE_GRID,POLAR_GRID,POLAR_GRIDC,DISP_GRID,INDUCE_GRID,INDUCE_GRIDC}; enum{MUTUAL,OPT,TCG,DIRECT}; enum{GEAR,ASPC,LSQR}; -#define DELTASTACK 1 // change this when debugged +#define DELTASTACK 16 + +#define UIND_DEBUG 0 // also in amoeba_induce.cpp /* ---------------------------------------------------------------------- */ @@ -226,9 +228,7 @@ PairAmoeba::~PairAmoeba() // DEBUG - if (me == 0) { - if (uind_flag) fclose(fp_uind); - } + if (me == 0 && UIND_DEBUG) fclose(fp_uind); } /* ---------------------------------------------------------------------- */ @@ -313,7 +313,7 @@ void PairAmoeba::compute(int eflag, int vflag) add_onefive_neighbors(); if (amoeba) find_hydrogen_neighbors(); find_multipole_neighbors(); - if (induce_flag && poltyp == MUTUAL && pcgprec) precond_neigh(); + if (poltyp == MUTUAL && pcgprec) precond_neigh(); } // reset KSpace recip matrix if box size/shape change dynamically @@ -377,18 +377,18 @@ void PairAmoeba::compute(int eflag, int vflag) // Ewald dispersion, pairwise and long range - if (hippo && disp_flag) dispersion(); + if (hippo && (disp_rspace_flag || disp_kspace_flag)) dispersion(); time4 = MPI_Wtime(); // multipole, pairwise and long range - if (mpole_flag) multipole(); + if (mpole_rspace_flag || mpole_kspace_flag) multipole(); time5 = MPI_Wtime(); // induced dipoles, interative CG relaxation // communicate induce() output values needed by ghost atoms - if (induce_flag) { + if (polar_rspace_flag || polar_kspace_flag) { induce(); cfstyle = INDUCE; comm->forward_comm_pair(this); @@ -397,7 +397,7 @@ void PairAmoeba::compute(int eflag, int vflag) // dipoles, pairwise and long range - if (polar_flag) polar(); + if (polar_rspace_flag || polar_kspace_flag) polar(); time7 = MPI_Wtime(); // charge transfer, pairwise @@ -522,49 +522,63 @@ void PairAmoeba::allocate() void PairAmoeba::settings(int narg, char **arg) { - // for now, allow turn on/off of individual FF terms - - // if (narg) error->all(FLERR,"Illegal pair_style command"); - // turn on all FF components by default - hal_flag = repulse_flag = disp_flag = induce_flag = - polar_flag = mpole_flag = qxfer_flag = 1; - rspace_flag = kspace_flag = 1; + hal_flag = repulse_flag = qxfer_flag = 1; + disp_rspace_flag = disp_kspace_flag = 1; + polar_rspace_flag = polar_kspace_flag = 1; + mpole_rspace_flag = mpole_kspace_flag = 1; + bond_flag = angle_flag = dihedral_flag = improper_flag = 1; + urey_flag = pitorsion_flag = xtorsion_flag = 1; - // turn off debug output by default + int newvalue = -1; - uind_flag = 0; + // include only specified FF components - // process optional args - // none turns all FF components off - - int iarg = 0; - while (iarg < narg) { - if (strcmp(arg[iarg],"none") == 0) { - hal_flag = repulse_flag = disp_flag = induce_flag = - polar_flag = mpole_flag = qxfer_flag = 0; - rspace_flag = kspace_flag = 0; - } - else if (strcmp(arg[iarg],"hal") == 0) hal_flag = 1; - else if (strcmp(arg[iarg],"repulse") == 0) repulse_flag = 1; - else if (strcmp(arg[iarg],"disp") == 0) disp_flag = 1; - else if (strcmp(arg[iarg],"induce") == 0) induce_flag = 1; - else if (strcmp(arg[iarg],"polar") == 0) polar_flag = 1; - else if (strcmp(arg[iarg],"mpole") == 0) mpole_flag = 1; - else if (strcmp(arg[iarg],"qxfer") == 0) qxfer_flag = 1; - - else if (strcmp(arg[iarg],"rspace") == 0) rspace_flag = 1; - else if (strcmp(arg[iarg],"kspace") == 0) kspace_flag = 1; - else if (strcmp(arg[iarg],"no-rspace") == 0) rspace_flag = 0; - else if (strcmp(arg[iarg],"no-kspace") == 0) kspace_flag = 0; + if (narg && (strcmp(arg[0],"include") == 0)) { + newvalue = 1; + hal_flag = repulse_flag = qxfer_flag = 0; + disp_rspace_flag = disp_kspace_flag = 0; + polar_rspace_flag = polar_kspace_flag = 0; + mpole_rspace_flag = mpole_kspace_flag = 0; + bond_flag = angle_flag = dihedral_flag = improper_flag = 0; + urey_flag = pitorsion_flag = xtorsion_flag = 0; - // DEBUG flags - - else if (strcmp(arg[iarg],"uind") == 0) uind_flag = 1; - + // exclude only specified FF components + + } else if (narg && (strcmp(arg[0],"exclude") == 0)) { + newvalue = 0; + + } else if (narg) error->all(FLERR,"Illegal pair_style command"); + + if (narg == 0) return; + + // toggle components to include or exclude + + for (int iarg = 1; iarg < narg; iarg++) { + if (strcmp(arg[iarg],"hal") == 0) hal_flag = newvalue; + else if (strcmp(arg[iarg],"repulse") == 0) repulse_flag = newvalue; + else if (strcmp(arg[iarg],"qxfer") == 0) qxfer_flag = newvalue; + else if (strcmp(arg[iarg],"disp") == 0) + disp_rspace_flag = disp_kspace_flag = newvalue; + else if (strcmp(arg[iarg],"disp/rspace") == 0) disp_rspace_flag = newvalue; + else if (strcmp(arg[iarg],"disp/kspace") == 0) disp_rspace_flag = newvalue; + else if (strcmp(arg[iarg],"polar") == 0) + polar_rspace_flag = polar_kspace_flag = newvalue; + else if (strcmp(arg[iarg],"polar/rspace") == 0) polar_rspace_flag = newvalue; + else if (strcmp(arg[iarg],"polar/kspace") == 0) polar_rspace_flag = newvalue; + else if (strcmp(arg[iarg],"mpole") == 0) + mpole_rspace_flag = mpole_kspace_flag = newvalue; + else if (strcmp(arg[iarg],"mpole/rspace") == 0) mpole_rspace_flag = newvalue; + else if (strcmp(arg[iarg],"mpole/kspace") == 0) mpole_rspace_flag = newvalue; + else if (strcmp(arg[iarg],"bond") == 0) bond_flag = newvalue; + else if (strcmp(arg[iarg],"angle") == 0) angle_flag = newvalue; + else if (strcmp(arg[iarg],"dihedral") == 0) dihedral_flag = newvalue; + else if (strcmp(arg[iarg],"improper") == 0) improper_flag = newvalue; + else if (strcmp(arg[iarg],"urey") == 0) urey_flag = newvalue; + else if (strcmp(arg[iarg],"pitorsion") == 0) pitorsion_flag = newvalue; + else if (strcmp(arg[iarg],"xtorsion") == 0) xtorsion_flag = newvalue; else error->all(FLERR,"Illegal pair_style command"); - iarg++; } } @@ -677,14 +691,9 @@ void PairAmoeba::init_style() //if (!force->special_onefive) // error->all(FLERR,"Pair style amoeba/hippo requires special_bonds one/five be set"); - // b/c induce computes fopt,foptp used by polar - - if (kspace_flag && optorder && polar_flag && !induce_flag) - error->all(FLERR,"Pair amoeba with optorder requires induce and polar together"); - // b/c polar uses mutipole virial terms - if (kspace_flag && apewald == aeewald && polar_flag && !mpole_flag) + if (apewald == aeewald && polar_kspace_flag && !mpole_kspace_flag) error->all(FLERR, "Pair amoeba with apewald = aeewald requires mpole and polar together"); @@ -979,7 +988,7 @@ void PairAmoeba::init_style() if (me == 0) { char fname[32]; sprintf(fname,"tmp.uind.kspace.%d",nprocs); - if (uind_flag) fp_uind = fopen(fname,"w"); + if (UIND_DEBUG) fp_uind = fopen(fname,"w"); } } @@ -1078,22 +1087,17 @@ double PairAmoeba::init_one(int i, int j) choose(REPULSE); cutoff = MAX(cutoff,sqrt(off2)); } - if (disp_flag) { + if (disp_rspace_flag || disp_kspace_flag) { if (use_dewald) choose(DISP_LONG); else choose(DISP); cutoff = MAX(cutoff,sqrt(off2)); } - if (mpole_flag) { + if (mpole_rspace_flag || mpole_kspace_flag) { if (use_ewald) choose(MPOLE_LONG); else choose(MPOLE); cutoff = MAX(cutoff,sqrt(off2)); } - if (induce_flag) { - if (use_ewald) choose(POLAR_LONG); - else choose(POLAR); - cutoff = MAX(cutoff,sqrt(off2)); - } - if (polar_flag) { + if (polar_rspace_flag || polar_kspace_flag) { if (use_ewald) choose(POLAR_LONG); else choose(POLAR); cutoff = MAX(cutoff,sqrt(off2)); @@ -2112,6 +2116,13 @@ void PairAmoeba::zero_energy_force_virial() void *PairAmoeba::extract(const char *str, int &dim) { dim = 0; + if (strcmp(str,"bond_flag") == 0) return (void *) &bond_flag; + if (strcmp(str,"angle_flag") == 0) return (void *) &angle_flag; + if (strcmp(str,"dihedral_flag") == 0) return (void *) &dihedral_flag; + if (strcmp(str,"improper_flag") == 0) return (void *) &improper_flag; + if (strcmp(str,"urey_flag") == 0) return (void *) &urey_flag; + if (strcmp(str,"pitorsion_flag") == 0) return (void *) &pitorsion_flag; + 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; diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h index 607c261695..342e0c5205 100644 --- a/src/AMOEBA/pair_amoeba.h +++ b/src/AMOEBA/pair_amoeba.h @@ -75,12 +75,12 @@ class PairAmoeba : public Pair { // turn on/off components of force field - int hal_flag,repulse_flag,disp_flag,induce_flag,polar_flag,mpole_flag,qxfer_flag; - int rspace_flag,kspace_flag; - - // DEBUG flags - - int uind_flag; + int hal_flag,repulse_flag,qxfer_flag; + int disp_rspace_flag,disp_kspace_flag; + int polar_rspace_flag,polar_kspace_flag; + int mpole_rspace_flag,mpole_kspace_flag; + int bond_flag,angle_flag,dihedral_flag,improper_flag; + int urey_flag,pitorsion_flag,xtorsion_flag; // DEBUG timers diff --git a/src/CLASS2/bond_class2.cpp b/src/CLASS2/bond_class2.cpp index 1a4c525ac6..d056b2aea9 100644 --- a/src/CLASS2/bond_class2.cpp +++ b/src/CLASS2/bond_class2.cpp @@ -24,10 +24,10 @@ #include "neighbor.h" #include "comm.h" #include "force.h" +#include "pair.h" #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -53,6 +53,8 @@ BondClass2::~BondClass2() void BondClass2::compute(int eflag, int vflag) { + if (disable) return; + int i1,i2,n,type; double delx,dely,delz,ebond,fbond; double rsq,r,dr,dr2,dr3,dr4,de_bond; @@ -155,6 +157,22 @@ void BondClass2::coeff(int narg, char **arg) if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients"); } +/* ---------------------------------------------------------------------- */ + +void BondClass2::init_style() +{ + // check if PairAmoeba disabled bond terms + + Pair *pair = force->pair_match("amoeba",1,0); + + if (!pair) disable = 0; + else { + int tmp; + int flag = *((int *) pair->extract("bond_flag",tmp)); + disable = flag ? 0 : 1; + } +} + /* ---------------------------------------------------------------------- return an equilbrium bond length ------------------------------------------------------------------------- */ diff --git a/src/CLASS2/bond_class2.h b/src/CLASS2/bond_class2.h index 75c1d52c1f..d9f7cafedd 100644 --- a/src/CLASS2/bond_class2.h +++ b/src/CLASS2/bond_class2.h @@ -30,6 +30,7 @@ class BondClass2 : public Bond { virtual ~BondClass2(); virtual void compute(int, int); virtual void coeff(int, char **); + void init_style(); double equilibrium_distance(int); void write_restart(FILE *); virtual void read_restart(FILE *); @@ -39,6 +40,7 @@ class BondClass2 : public Bond { protected: double *r0, *k2, *k3, *k4; + int disable; void allocate(); }; diff --git a/src/EXTRA-MOLECULE/dihedral_fourier.cpp b/src/EXTRA-MOLECULE/dihedral_fourier.cpp index b88c866e7c..1dcbdf4563 100644 --- a/src/EXTRA-MOLECULE/dihedral_fourier.cpp +++ b/src/EXTRA-MOLECULE/dihedral_fourier.cpp @@ -24,12 +24,12 @@ #include "comm.h" #include "neighbor.h" #include "force.h" +#include "pair.h" #include "update.h" #include "math_const.h" #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; using namespace MathConst; @@ -70,6 +70,8 @@ DihedralFourier::~DihedralFourier() void DihedralFourier::compute(int eflag, int vflag) { + if (disable) return; + int i1,i2,i3,i4,i,j,m,n,type; double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double edihedral,f1[3],f2[3],f3[3],f4[3]; @@ -327,13 +329,28 @@ void DihedralFourier::coeff(int narg, char **arg) if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); } +/* ---------------------------------------------------------------------- */ + +void DihedralFourier::init_style() +{ + // check if PairAmoeba disabled dihedral terms + + Pair *pair = force->pair_match("amoeba",1,0); + + if (!pair) disable = 0; + else { + int tmp; + int flag = *((int *) pair->extract("dihedral_flag",tmp)); + disable = flag ? 0 : 1; + } +} + /* ---------------------------------------------------------------------- proc 0 writes out coeffs to restart file ------------------------------------------------------------------------- */ void DihedralFourier::write_restart(FILE *fp) { - fwrite(&nterms[1],sizeof(int),atom->ndihedraltypes,fp); for (int i = 1; i <= atom->ndihedraltypes; i++) { fwrite(k[i],sizeof(double),nterms[i],fp); diff --git a/src/EXTRA-MOLECULE/dihedral_fourier.h b/src/EXTRA-MOLECULE/dihedral_fourier.h index 557d4b5277..208cf2099f 100644 --- a/src/EXTRA-MOLECULE/dihedral_fourier.h +++ b/src/EXTRA-MOLECULE/dihedral_fourier.h @@ -30,6 +30,7 @@ class DihedralFourier : public Dihedral { virtual ~DihedralFourier(); virtual void compute(int, int); void coeff(int, char **); + void init_style(); void write_restart(FILE *); void read_restart(FILE *); void write_data(FILE *); @@ -39,6 +40,7 @@ class DihedralFourier : public Dihedral { int **multiplicity; int *nterms; int implicit, weightflag; + int disable; void allocate(); };