first working version of fix pitorsion

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

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
 

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
 

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

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
 

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)
+      if (tflag && nspecial[i2][0] == 3)
-      tinker_anglep(i1,i2,i3,type,eflag);
+        tinker_anglep(i1,i2,i3,type,eflag);
-    else
+      else
-      tinker_angle(i1,i2,i3,type,eflag);
+        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)
+      if (ba_k1[type] != 0.0)
-      tinker_bondangle(i1,i2,i3,type,eflag);
+        tinker_bondangle(i1,i2,i3,type,eflag);
+    }
 
     // Urey-Bradley H-H bond term within water molecules
 
-    uflag = ubflag[type];
+    if (enable_urey) {
-    if (uflag) tinker_urey_bradley(i1,i3,type,eflag);
+      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];

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);

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;
+  int nlist,list[6];
-  double **f = atom->f;
+  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];
+    ib = pitorsion_list[n][1];
-    ic = pitorsion_list[n][0];
+    ic = pitorsion_list[n][2];
-    id = pitorsion_list[n][0];
+    id = pitorsion_list[n][3];
-    ie = pitorsion_list[n][0];
+    ie = pitorsion_list[n][4];
-    ig = pitorsion_list[n][0];
+    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;
+    if (ib < nlocal) {
-    f[ia][1] += dedyia;
+      epitorsion += engfraction;
-    f[ia][2] += dedzia;
+      f[ib][0] += dedxib;
-    f[ib][0] += dedxib;
+      f[ib][1] += dedyib;
-    f[ib][1] += dedyib;
+      f[ib][2] += dedzib;
-    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;
 
-    // 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;
+    if (id < nlocal) {
-    vyterm = dedyid + dedyia + dedyib;
+      epitorsion += engfraction;
-    vzterm = dedzid + dedzia + dedzib;
+      f[id][0] += dedxid;
-    vxx = xdc*vxterm + xcp*dedxip - xqd*dedxiq;
+      f[id][1] += dedyid;
-    vyx = ydc*vxterm + ycp*dedxip - yqd*dedxiq;
+      f[id][2] += dedzid;
-    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;
 
-    virial[0] += vxx;
+    if (ie < nlocal) {
-    virial[1] += vyy;
+      epitorsion += engfraction;
-    virial[2] += vzz;
+      f[ie][0] += dedxie;
-    virial[3] += vyx;
+      f[ie][1] += dedyie;
-    virial[4] += vzx;
+      f[ie][2] += dedzie;
-    virial[5] += vzy;
+    }
+
+    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(atom4)) >= 0) {
-        if ((m = atom->map(atom1)) >= 0) {
+        if (num_pitorsion[m] == PITORSIONMAX)
-          if (num_pitorsion[m] == PITORSIONMAX)
+          error->one(FLERR,"Too many PiTorsions for one atom");
-            error->one(FLERR,"Too many PiTorsions for one atom");
+        pitorsion_type[m][num_pitorsion[m]] = itype;
-          pitorsion_type[m][num_pitorsion[m]] = itype;
+        pitorsion_atom1[m][num_pitorsion[m]] = atom1;
-          pitorsion_atom1[m][num_pitorsion[m]] = atom1;
+        pitorsion_atom2[m][num_pitorsion[m]] = atom2;
-          pitorsion_atom2[m][num_pitorsion[m]] = atom2;
+        pitorsion_atom3[m][num_pitorsion[m]] = atom3;
-          pitorsion_atom3[m][num_pitorsion[m]] = atom3;
+        pitorsion_atom4[m][num_pitorsion[m]] = atom4;
-          pitorsion_atom4[m][num_pitorsion[m]] = atom4;
+        pitorsion_atom5[m][num_pitorsion[m]] = atom5;
-          pitorsion_atom5[m][num_pitorsion[m]] = atom5;
+        pitorsion_atom6[m][num_pitorsion[m]] = atom6;
-          pitorsion_atom6[m][num_pitorsion[m]] = atom6;
+        num_pitorsion[m]++;
-          num_pitorsion[m]++;
-        }
       }
-
+      
-      if (newton_bond == 0) {
+      if ((m = atom->map(atom5)) >= 0) {
-        if ((m = atom->map(atom2)) >= 0) {
+        if (num_pitorsion[m] == PITORSIONMAX)
-          if (num_pitorsion[m] == PITORSIONMAX)
+          error->one(FLERR,"Too many PiTorsions for one atom");
-            error->one(FLERR,"Too many PiTorsions for one atom");
+        pitorsion_type[m][num_pitorsion[m]] = itype;
-          pitorsion_type[m][num_pitorsion[m]] = itype;
+        pitorsion_atom1[m][num_pitorsion[m]] = atom1;
-          pitorsion_atom1[m][num_pitorsion[m]] = atom1;
+        pitorsion_atom2[m][num_pitorsion[m]] = atom2;
-          pitorsion_atom2[m][num_pitorsion[m]] = atom2;
+        pitorsion_atom3[m][num_pitorsion[m]] = atom3;
-          pitorsion_atom3[m][num_pitorsion[m]] = atom3;
+        pitorsion_atom4[m][num_pitorsion[m]] = atom4;
-          pitorsion_atom4[m][num_pitorsion[m]] = atom4;
+        pitorsion_atom5[m][num_pitorsion[m]] = atom5;
-          pitorsion_atom5[m][num_pitorsion[m]] = atom5;
+        pitorsion_atom6[m][num_pitorsion[m]] = atom6;
-          pitorsion_atom6[m][num_pitorsion[m]] = atom6;
+        num_pitorsion[m]++;
-          num_pitorsion[m]++;
-        }
       }
-
+      
-      if (newton_bond == 0) {
+      if ((m = atom->map(atom6)) >= 0) {
-        if ((m = atom->map(atom4)) >= 0) {
+        if (num_pitorsion[m] == PITORSIONMAX)
-          if (num_pitorsion[m] == PITORSIONMAX)
+          error->one(FLERR,"Too many PiTorsions for one atom");
-            error->one(FLERR,"Too many PiTorsions for one atom");
+        pitorsion_type[m][num_pitorsion[m]] = itype;
-          pitorsion_type[m][num_pitorsion[m]] = itype;
+        pitorsion_atom1[m][num_pitorsion[m]] = atom1;
-          pitorsion_atom1[m][num_pitorsion[m]] = atom1;
+        pitorsion_atom2[m][num_pitorsion[m]] = atom2;
-          pitorsion_atom2[m][num_pitorsion[m]] = atom2;
+        pitorsion_atom3[m][num_pitorsion[m]] = atom3;
-          pitorsion_atom3[m][num_pitorsion[m]] = atom3;
+        pitorsion_atom4[m][num_pitorsion[m]] = atom4;
-          pitorsion_atom4[m][num_pitorsion[m]] = atom4;
+        pitorsion_atom5[m][num_pitorsion[m]] = atom5;
-          pitorsion_atom5[m][num_pitorsion[m]] = atom5;
+        pitorsion_atom6[m][num_pitorsion[m]] = atom6;
-          pitorsion_atom6[m][num_pitorsion[m]] = atom6;
+        num_pitorsion[m]++;
-          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) {

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;
 

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;
 }
 
 /* ----------------------------------------------------------------------

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;
 

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 (me == 0 && UIND_DEBUG) fclose(fp_uind);
-    if (uind_flag) 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 = 
+  hal_flag = repulse_flag = qxfer_flag = 1;
-    polar_flag = mpole_flag = qxfer_flag = 1;
+  disp_rspace_flag = disp_kspace_flag = 1;
-  rspace_flag = 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
+  if (narg && (strcmp(arg[0],"include") == 0)) {
-  // none turns all FF components off
+    newvalue = 1;
-  
+    hal_flag = repulse_flag = qxfer_flag = 0;
-  int iarg = 0;
+    disp_rspace_flag = disp_kspace_flag = 0;
-  while (iarg < narg) {
+    polar_rspace_flag = polar_kspace_flag = 0;
-    if (strcmp(arg[iarg],"none") == 0) {
+    mpole_rspace_flag = mpole_kspace_flag = 0;
-      hal_flag = repulse_flag = disp_flag = induce_flag = 
+    bond_flag = angle_flag = dihedral_flag = improper_flag = 0;
-        polar_flag = mpole_flag = qxfer_flag = 0;
+    urey_flag = pitorsion_flag = xtorsion_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;
 
-    // DEBUG flags
+  // exclude only specified FF components
-    
+
-    else if (strcmp(arg[iarg],"uind") == 0) uind_flag = 1;
+  } 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 (polar_rspace_flag || polar_kspace_flag) {
-    if (use_ewald) choose(POLAR_LONG);
-    else choose(POLAR);
-    cutoff = MAX(cutoff,sqrt(off2));
-  }
-  if (polar_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;

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 hal_flag,repulse_flag,qxfer_flag;
-  int rspace_flag,kspace_flag;
+  int disp_rspace_flag,disp_kspace_flag;
-
+  int polar_rspace_flag,polar_kspace_flag;
-  // DEBUG flags
+  int mpole_rspace_flag,mpole_kspace_flag;
-
+  int bond_flag,angle_flag,dihedral_flag,improper_flag;
-  int uind_flag;
+  int urey_flag,pitorsion_flag,xtorsion_flag;
 
   // DEBUG timers
 

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
 ------------------------------------------------------------------------- */

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();
 };

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);

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();
 };