Merge branch 'amoeba' into amoeba-gpu

examples/amoeba/in.ubiquitin

@@ -2,23 +2,34 @@
 
 units           real
 boundary        p p p
+atom_modify     sort 0 0.0
+
+atom_style      amoeba
 
-atom_style amoeba
-#atom_modify sort 1000 7.0
 bond_style      class2
 angle_style     amoeba
-dihedral_style  none
+dihedral_style  fourier
+improper_style  amoeba
 
 # per-atom properties required by AMOEBA or HIPPO
 
 fix             amtype all property/atom i_amtype ghost yes
+fix             pit all amoeba/pitorsion
+fix_modify      pit energy yes
+fix             bit all amoeba/bitorsion bitorsion.ubiquitin.data
+fix_modify      bit energy yes
+
 fix             extra all property/atom &
                 i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes
 fix             extra2 all property/atom i_polaxe
 
-read_data       data.ubiquitin fix amtype NULL "Tinker Types"
+#read_data       data.ubiquitin fix amtype NULL "Tinker Types"
+read_data       data.ubiquitin fix amtype NULL "Tinker Types" &
+                fix pit "pitorsion types" "PiTorsion Coeffs" &
+                fix pit pitorsions PiTorsions &
+                fix bit bitorsions BiTorsions
 
-pair_style      amoeba
+pair_style      amoeba exclude bitorsion
 pair_coeff      * * amoeba_ubiquitin.prm amoeba_ubiquitin.key
 
 special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
@@ -37,17 +48,8 @@ fix             fqxfer all store/state 0 fx fy fz
 compute virial all pressure NULL virial
 
 thermo_style    custom step temp epair ebond eangle edihed eimp &
-                emol etotal press c_virial[*]
+                emol etotal #press c_virial[*]
-
-# DEBUG: setup force components this way so can dump them
-
-fix             fhal all store/state 0 fx fy fz
-fix             frepulse all store/state 0 fx fy fz
-fix             fdisp all store/state 0 fx fy fz
-fix             fpolar all store/state 0 fx fy fz
-fix             fmpole all store/state 0 fx fy fz
-fix             fqxfer all store/state 0 fx fy fz
 
 # zero step run
 
-run             0
+run              0

examples/amoeba/in.water_box.amoeba

@@ -24,6 +24,8 @@ read_data       data.water_box.amoeba fix amtype NULL "Tinker Types"
 pair_style      amoeba
 pair_coeff      * * amoeba_water.prm amoeba_water.key
 
+special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
+
 # thermo output
 
 compute         virial all pressure NULL virial
@@ -31,15 +33,6 @@ compute         virial all pressure NULL virial
 thermo_style    custom step temp epair ebond eangle edihed eimp &
                 emol etotal press c_virial[*]
 
-# DEBUG: setup force components this way so can dump them
-
-fix             fhal all store/state 0 fx fy fz
-fix             frepulse all store/state 0 fx fy fz
-fix             fdisp all store/state 0 fx fy fz
-fix             fpolar all store/state 0 fx fy fz
-fix             fmpole all store/state 0 fx fy fz
-fix             fqxfer all store/state 0 fx fy fz
-
 # zero step run
 
 run             0

examples/amoeba/in.water_box.hippo

@@ -24,6 +24,8 @@ read_data       data.water_box.hippo fix amtype NULL "Tinker Types"
 pair_style      hippo
 pair_coeff      * * hippo_water.prm hippo_water.key
 
+special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
+
 # thermo output
 
 compute         virial all pressure NULL virial
@@ -31,15 +33,6 @@ compute         virial all pressure NULL virial
 thermo_style    custom step temp epair ebond eangle edihed eimp &
                 emol etotal press c_virial[*]
 
-# DEBUG: setup force components this way so can dump them
-
-fix             fhal all store/state 0 fx fy fz
-fix             frepulse all store/state 0 fx fy fz
-fix             fdisp all store/state 0 fx fy fz
-fix             fpolar all store/state 0 fx fy fz
-fix             fmpole all store/state 0 fx fy fz
-fix             fqxfer all store/state 0 fx fy fz
-
 # zero step run
 
 run             0

examples/amoeba/in.water_dimer.amoeba

@@ -24,6 +24,8 @@ read_data       data.water_dimer.amoeba fix amtype NULL "Tinker Types"
 pair_style      amoeba
 pair_coeff      * * amoeba_water.prm amoeba_water.key
 
+special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
+
 # thermo output
 
 compute         virial all pressure NULL virial
@@ -31,15 +33,6 @@ compute         virial all pressure NULL virial
 thermo_style    custom step temp epair ebond eangle edihed eimp &
                 emol etotal press c_virial[*]
 
-# DEBUG: setup force components this way so can dump them
-
-fix             fhal all store/state 0 fx fy fz
-fix             frepulse all store/state 0 fx fy fz
-fix             fdisp all store/state 0 fx fy fz
-fix             fpolar all store/state 0 fx fy fz
-fix             fmpole all store/state 0 fx fy fz
-fix             fqxfer all store/state 0 fx fy fz
-
 # zero step run
 
 run             0

examples/amoeba/in.water_dimer.hippo

@@ -24,6 +24,8 @@ read_data       data.water_dimer.hippo fix amtype NULL "Tinker Types"
 pair_style      hippo
 pair_coeff      * * hippo_water.prm hippo_water.key
 
+special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
+
 # thermo output
 
 compute         virial all pressure NULL virial
@@ -31,15 +33,6 @@ compute         virial all pressure NULL virial
 thermo_style    custom step temp epair ebond eangle edihed eimp &
                 emol etotal press c_virial[*]
 
-# DEBUG: setup force components this way so can dump them
-
-fix             fhal all store/state 0 fx fy fz
-fix             frepulse all store/state 0 fx fy fz
-fix             fdisp all store/state 0 fx fy fz
-fix             fpolar all store/state 0 fx fy fz
-fix             fmpole all store/state 0 fx fy fz
-fix             fqxfer all store/state 0 fx fy fz
-
 # zero step run
 
 run             0

examples/amoeba/in.water_hexamer.amoeba

@@ -24,6 +24,8 @@ read_data       data.water_hexamer.amoeba fix amtype NULL "Tinker Types"
 pair_style      amoeba
 pair_coeff      * * amoeba_water.prm amoeba_water.key
 
+special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
+
 # thermo output
 
 compute         virial all pressure NULL virial
@@ -31,15 +33,6 @@ compute         virial all pressure NULL virial
 thermo_style    custom step temp epair ebond eangle edihed eimp &
                 emol etotal press c_virial[*]
 
-# DEBUG: setup force components this way so can dump them
-
-fix             fhal all store/state 0 fx fy fz
-fix             frepulse all store/state 0 fx fy fz
-fix             fdisp all store/state 0 fx fy fz
-fix             fpolar all store/state 0 fx fy fz
-fix             fmpole all store/state 0 fx fy fz
-fix             fqxfer all store/state 0 fx fy fz
-
 # zero step run
 
 run             0

examples/amoeba/in.water_hexamer.hippo

@@ -24,6 +24,8 @@ read_data       data.water_hexamer.hippo fix amtype NULL "Tinker Types"
 pair_style      hippo
 pair_coeff      * * hippo_water.prm hippo_water.key
 
+special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
+
 # thermo output
 
 compute         virial all pressure NULL virial
@@ -31,15 +33,6 @@ compute         virial all pressure NULL virial
 thermo_style    custom step temp epair ebond eangle edihed eimp &
                 emol etotal press c_virial[*]
 
-# DEBUG: setup force components this way so can dump them
-
-fix             fhal all store/state 0 fx fy fz
-fix             frepulse all store/state 0 fx fy fz
-fix             fdisp all store/state 0 fx fy fz
-fix             fpolar all store/state 0 fx fy fz
-fix             fmpole all store/state 0 fx fy fz
-fix             fqxfer all store/state 0 fx fy fz
-
 # zero step run
 
 run             0

src/AMOEBA/amoeba_charge_transfer.cpp

@@ -57,6 +57,7 @@ void PairAmoeba::charge_transfer()
   // owned atoms
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   // neigh list
@@ -138,12 +139,12 @@ void PairAmoeba::charge_transfer()
 
       // increment the total charge transfer energy and derivatives
 
-      fqxfer[i][0] -= frcx;
+      f[i][0] -= frcx;
-      fqxfer[i][1] -= frcy;
+      f[i][1] -= frcy;
-      fqxfer[i][2] -= frcz;
+      f[i][2] -= frcz;
-      fqxfer[j][0] += frcx;
+      f[j][0] += frcx;
-      fqxfer[j][1] += frcy;
+      f[j][1] += frcy;
-      fqxfer[j][2] += frcz;
+      f[j][2] += frcz;
 
       // increment the internal virial tensor components
 

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
 
@@ -99,6 +99,7 @@ void PairAmoeba::dispersion_real()
   // owned atoms
   
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   // neigh list
@@ -202,12 +203,12 @@ void PairAmoeba::dispersion_real()
       dedx = de * xr;
       dedy = de * yr;
       dedz = de * zr;
-      fdisp[i][0] += dedx;
+      f[i][0] += dedx;
-      fdisp[i][1] += dedy;
+      f[i][1] += dedy;
-      fdisp[i][2] += dedz;
+      f[i][2] += dedz;
-      fdisp[j][0] -= dedx;
+      f[j][0] -= dedx;
-      fdisp[j][1] -= dedy;
+      f[j][1] -= dedy;
-      fdisp[j][2] -= dedz;
+      f[j][2] -= dedz;
 
       // increment the internal virial tensor components
 
@@ -270,6 +271,7 @@ void PairAmoeba::dispersion_kspace()
   // owned atoms
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   double volbox = domain->prd[0] * domain->prd[1] * domain->prd[2];
@@ -405,15 +407,14 @@ void PairAmoeba::dispersion_kspace()
     }
 
     fi = csix[iclass];
-    fdisp[m][0] += fi * (recip[0][0]*de1 + recip[0][1]*de2 + recip[0][2]*de3);
+    f[m][0] += fi * (recip[0][0]*de1 + recip[0][1]*de2 + recip[0][2]*de3);
-    fdisp[m][1] += fi * (recip[1][0]*de1 + recip[1][1]*de2 + recip[1][2]*de3);
+    f[m][1] += fi * (recip[1][0]*de1 + recip[1][1]*de2 + recip[1][2]*de3);
-    fdisp[m][2] += fi * (recip[2][0]*de1 + recip[2][1]*de2 + recip[2][2]*de3);
+    f[m][2] += fi * (recip[2][0]*de1 + recip[2][1]*de2 + recip[2][2]*de3);
   }
 
   // account for the energy and virial correction terms
 
   term = csixpr * aewald*aewald*aewald / denom0;
-
   
   if (me == 0) {
     edisp -= term;

src/AMOEBA/amoeba_hal.cpp

@@ -61,6 +61,7 @@ void PairAmoeba::hal()
   // owned atoms
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   // neigh list
@@ -158,31 +159,31 @@ void PairAmoeba::hal()
 		   "ghost comm is too short");
       
       if (i == iv) {
-        fhal[i][0] += dedx;
+        f[i][0] += dedx;
-        fhal[i][1] += dedy;
+        f[i][1] += dedy;
-        fhal[i][2] += dedz;
+        f[i][2] += dedz;
       } else {
-        fhal[i][0] += dedx*redi;
+        f[i][0] += dedx*redi;
-        fhal[i][1] += dedy*redi;
+        f[i][1] += dedy*redi;
-        fhal[i][2] += dedz*redi;
+        f[i][2] += dedz*redi;
-        fhal[iv][0] += dedx*rediv;
+        f[iv][0] += dedx*rediv;
-        fhal[iv][1] += dedy*rediv;
+        f[iv][1] += dedy*rediv;
-        fhal[iv][2] += dedz*rediv;
+        f[iv][2] += dedz*rediv;
       }
 
       if (j == jv) {
-        fhal[j][0] -= dedx;
+        f[j][0] -= dedx;
-        fhal[j][1] -= dedy;
+        f[j][1] -= dedy;
-        fhal[j][2] -= dedz;
+        f[j][2] -= dedz;
       } else {
         redj = kred[jclass];
         redjv = 1.0 - redj;
-        fhal[j][0] -= dedx*redj;
+        f[j][0] -= dedx*redj;
-        fhal[j][1] -= dedy*redj;
+        f[j][1] -= dedy*redj;
-        fhal[j][2] -= dedz*redj;
+        f[j][2] -= dedz*redj;
-        fhal[jv][0] -= dedx*redjv;
+        f[jv][0] -= dedx*redjv;
-        fhal[jv][1] -= dedy*redjv;
+        f[jv][1] -= dedy*redjv;
-        fhal[jv][2] -= dedz*redjv;
+        f[jv][2] -= dedz*redjv;
       }
 
       // increment the internal virial tensor components

src/AMOEBA/amoeba_induce.cpp

@@ -41,6 +41,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
@@ -586,7 +588,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
@@ -762,11 +764,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
 
@@ -1030,7 +1032,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++) {
@@ -1040,7 +1042,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
@@ -519,9 +519,9 @@ void PairAmoeba::multipole_real()
 
       // increment force-based gradient and torque on first site
 
-      fmpole[i][0] += frcx;
+      f[i][0] += frcx;
-      fmpole[i][1] += frcy;
+      f[i][1] += frcy;
-      fmpole[i][2] += frcz;
+      f[i][2] += frcz;
       tq[i][0] += ttmi[0];
       tq[i][1] += ttmi[1];
       tq[i][2] += ttmi[2];
@@ -529,9 +529,9 @@ void PairAmoeba::multipole_real()
       // increment force-based gradient and torque on second site
       // commenting out j parts for DEBUGGING
 
-      fmpole[j][0] -= frcx;
+      f[j][0] -= frcx;
-      fmpole[j][1] -= frcy;
+      f[j][1] -= frcy;
-      fmpole[j][2] -= frcz;
+      f[j][2] -= frcz;
       tq[j][0] += ttmk[0];
       tq[j][1] += ttmk[1];
       tq[j][2] += ttmk[2];
@@ -569,7 +569,7 @@ void PairAmoeba::multipole_real()
   // resolve site torques then increment forces and virial
 
   for (i = 0; i < nlocal; i++) {
-    torque2force(i,tq[i],fix,fiy,fiz,fmpole);
+    torque2force(i,tq[i],fix,fiy,fiz,f);  
 
     iz = zaxis2local[i];
     ix = xaxis2local[i];
@@ -648,6 +648,7 @@ void PairAmoeba::multipole_kspace()
   // owned atoms
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   double volbox = domain->prd[0] * domain->prd[1] * domain->prd[2];
@@ -817,12 +818,12 @@ void PairAmoeba::multipole_kspace()
     h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3;  // matvec?
     h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
     h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
-    fmpole[i][0] += h1;
+    f[i][0] += h1;
-    fmpole[i][1] += h2;
+    f[i][1] += h2;
-    fmpole[i][2] += h3;
+    f[i][2] += h3;
   }
   empole += 0.5*e;
-  //printf("mpole_force %g %g %g \n", fmpole[0][0], fmpole[0][1], fmpole[0][2]);
+  //printf("mpole_force %g %g %g \n", f[0][0], f[0][1], f[0][2]);
   
   // augment the permanent multipole virial contributions
 
@@ -860,7 +861,7 @@ void PairAmoeba::multipole_kspace()
       cmp[i][7]*cphi[i][4] + cmp[i][9]*cphi[i][8] - 
       cmp[i][7]*cphi[i][5] - cmp[i][8]*cphi[i][9];
 
-    torque2force(i,tem,fix,fiy,fiz,fmpole);
+    torque2force(i,tem,fix,fiy,fiz,f);
 
     iz = zaxis2local[i];
     ix = xaxis2local[i];

src/AMOEBA/amoeba_polar.cpp

@@ -68,6 +68,7 @@ void PairAmoeba::polar()
   // owned atoms
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   // set the energy unit conversion factor
@@ -80,11 +81,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
 
@@ -101,7 +102,7 @@ void PairAmoeba::polar()
     tep[1] = term * (diz*uix-dix*uiz);
     tep[2] = term * (dix*uiy-diy*uix);
 
-    torque2force(i,tep,fix,fiy,fiz,fpolar);
+    torque2force(i,tep,fix,fiy,fiz,f);
 
     iz = zaxis2local[i];
     ix = xaxis2local[i];
@@ -277,6 +278,7 @@ void PairAmoeba::polar_real()
   pval = atom->dvector[index_pval];
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
@@ -1138,12 +1140,12 @@ void PairAmoeba::polar_real()
       
       // increment force-based gradient on the interaction sites
       
-      fpolar[i][0] -= frcx;
+      f[i][0] -= frcx;
-      fpolar[i][1] -= frcy;
+      f[i][1] -= frcy;
-      fpolar[i][2] -= frcz;
+      f[i][2] -= frcz;
-      fpolar[j][0] += frcx;
+      f[j][0] += frcx;
-      fpolar[j][1] += frcy;
+      f[j][1] += frcy;
-      fpolar[j][2] += frcz;
+      f[j][2] += frcz;
 
       // increment the virial due to pairwise Cartesian forces
       
@@ -1197,10 +1199,8 @@ void PairAmoeba::polar_real()
       qiyz*dufld[i][3] - qixz*dufld[i][4] + 
       2.0*qixy*(dufld[i][0]-dufld[i][2]) + (qiyy-qixx)*dufld[i][1];
 
-    torque2force(i,tep,fix,fiy,fiz,fpolar);
+    torque2force(i,tep,fix,fiy,fiz,f);
-
+    
-    //if (i < 10) printf("i = %d: tep = %f %f %f\n", i, tep[0], tep[1], tep[2]);
-
     iz = zaxis2local[i];
     ix = xaxis2local[i];
     iy = yaxis2local[i];
@@ -1282,6 +1282,7 @@ void PairAmoeba::polar_kspace()
   // owned atoms
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   double volbox = domain->prd[0] * domain->prd[1] * domain->prd[2];
@@ -1517,9 +1518,9 @@ void PairAmoeba::polar_kspace()
     h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3;
     h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
     h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
-    fpolar[i][0] += h1;
+    f[i][0] += h1;
-    fpolar[i][1] += h2;
+    f[i][1] += h2;
-    fpolar[i][2] += h3;
+    f[i][2] += h3;
   }
 
   // set the potential to be the induced dipole average
@@ -1606,7 +1607,7 @@ void PairAmoeba::polar_kspace()
       2.0*(cmp[i][5]-cmp[i][4])*cphidp[i][7] + cmp[i][7]*cphidp[i][4] + 
       cmp[i][9]*cphidp[i][8] - cmp[i][7]*cphidp[i][5] - cmp[i][8]*cphidp[i][9];
     
-    torque2force(i,tep,fix,fiy,fiz,fpolar);
+    torque2force(i,tep,fix,fiy,fiz,f);
 
     iz = zaxis2local[i];
     ix = xaxis2local[i];
@@ -1671,9 +1672,9 @@ void PairAmoeba::polar_kspace()
           h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
           h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
 
-          fpolar[i][0] += copm[k+m+1]*h1;
+          f[i][0] += copm[k+m+1]*h1;
-          fpolar[i][1] += copm[k+m+1]*h2;
+          f[i][1] += copm[k+m+1]*h2;
-          fpolar[i][2] += copm[k+m+1]*h3;
+          f[i][2] += copm[k+m+1]*h3;
 
           for (j = 1; j < 4; j++) {
             cphid[j] = 0.0;
@@ -1761,9 +1762,9 @@ void PairAmoeba::polar_kspace()
         h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3;
         h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
         h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
-        fpolar[i][0] += h1;
+        f[i][0] += h1;
-        fpolar[i][1] += h2;
+        f[i][1] += h2;
-        fpolar[i][2] += h3;
+        f[i][2] += h3;
 
         for (j = 1; j < 4; j++) {
           cphid[j] = 0.0;
@@ -1837,9 +1838,9 @@ void PairAmoeba::polar_kspace()
         h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3;  // matvec
         h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
         h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
-        fpolar[i][0] += h1;
+        f[i][0] += h1;
-        fpolar[i][1] += h2;
+        f[i][1] += h2;
-        fpolar[i][2] += h3;
+        f[i][2] += h3;
 
         for (j = 1; j < 4; j++) {
           cphid[j] = 0.0;

src/AMOEBA/amoeba_repulsion.cpp

@@ -91,6 +91,7 @@ void PairAmoeba::repulsion()
   // owned atoms
 
   double **x = atom->x;
+  double **f = atom->f;
   int nlocal = atom->nlocal;
 
   // zero repulsion torque on owned + ghost atoms
@@ -337,18 +338,18 @@ void PairAmoeba::repulsion()
       
       // increment force-based gradient and torque on atom I
 
-      frepulse[i][0] += frcx;
+      f[i][0] += frcx;
-      frepulse[i][1] += frcy;
+      f[i][1] += frcy;
-      frepulse[i][2] += frcz;
+      f[i][2] += frcz;
       tq[i][0] += ttri[0];
       tq[i][1] += ttri[1];
       tq[i][2] += ttri[2];
 
       // increment force-based gradient and torque on atom J
 
-      frepulse[j][0] -= frcx;
+      f[j][0] -= frcx;
-      frepulse[j][1] -= frcy;
+      f[j][1] -= frcy;
-      frepulse[j][2] -= frcz;
+      f[j][2] -= frcz;
       tq[j][0] += ttrk[0];
       tq[j][1] += ttrk[1];
       tq[j][2] += ttrk[2];
@@ -389,7 +390,7 @@ void PairAmoeba::repulsion()
   // resolve site torques then increment forces and virial
 
   for (i = 0; i < nlocal; i++) {
-    torque2force(i,tq[i],fix,fiy,fiz,frepulse);
+    torque2force(i,tq[i],fix,fiy,fiz,f);
 
     iz = zaxis2local[i];
     ix = xaxis2local[i];

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"
@@ -35,12 +36,22 @@ using namespace MathConst;
 AngleAmoeba::AngleAmoeba(LAMMPS *lmp) : Angle(lmp) 
 {
   pflag = nullptr;
+  ubflag = nullptr;
+
   theta0 = nullptr;
   k2 = nullptr;
   k3 = nullptr;
   k4 = nullptr;
   k5 = nullptr;
   k6 = nullptr;
+
+  ba_k1 = nullptr;
+  ba_k2 = nullptr;
+  ba_r1 = nullptr;
+  ba_r2 = nullptr;
+
+  ub_k = nullptr;
+  ub_r0 = nullptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -51,14 +62,27 @@ AngleAmoeba::~AngleAmoeba()
 
   if (allocated) {
     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);
     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);
+
+    memory->destroy(ub_k);
+    memory->destroy(ub_r0);
   }
 }
 
@@ -66,22 +90,19 @@ 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 eangle,f1[3],f3[3];
+  double f1[3],f3[3];
   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);
-
-  double **x = atom->x;
-  double **f = atom->f;
   int **anglelist = neighbor->anglelist;
+  int **nspecial = atom->nspecial;
   int nanglelist = neighbor->nanglelist;
-  int nlocal = atom->nlocal;
+
-  int newton_bond = force->newton_bond;
+  ev_init(eflag,vflag);
 
   for (n = 0; n < nanglelist; n++) {
     i1 = anglelist[n][0];
@@ -90,101 +111,133 @@ void AngleAmoeba::compute(int eflag, int vflag)
     type = anglelist[n][3];
 
     // tflag = 0 for "angle", 1 for "anglep" in Tinker PRM file
-    // atom 2 must have 3 bond partners to invoke "anglep" option
+    // atom 2 must have exactly 3 bond partners to invoke anglep() variant
 
-    tflag = pflag[type];
+    if (enable_angle) {
+      tflag = pflag[type];
 
-    if (tflag && atom->num_bond[i2] == 3) {
+      if (tflag && nspecial[i2][0] == 3)
-      anglep(i1,i2,i3,type,eflag);
+        tinker_anglep(i1,i2,i3,type,eflag);
-      continue;
+      else
+        tinker_angle(i1,i2,i3,type,eflag);
+
+      // bondangle = bond-stretch cross term in Tinker
+
+      if (ba_k1[type] != 0.0)
+        tinker_bondangle(i1,i2,i3,type,eflag);
     }
 
-    // 1st bond
+    // Urey-Bradley H-H bond term within water molecules
 
-    delx1 = x[i1][0] - x[i2][0];
+    if (enable_urey) {
-    dely1 = x[i1][1] - x[i2][1];
+      uflag = ubflag[type];
-    delz1 = x[i1][2] - x[i2][2];
+      if (uflag) tinker_urey_bradley(i1,i3,type,eflag);
-
-    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;
-
-    // force & energy for angle term
-
-    dtheta = acos(c) - theta0[type];
-    dtheta2 = dtheta*dtheta;
-    dtheta3 = dtheta2*dtheta;
-    dtheta4 = dtheta3*dtheta;
-    dtheta5 = dtheta4*dtheta;
-    dtheta6 = dtheta5*dtheta;
-
-    de_angle = 2.0*k2[type]*dtheta + 3.0*k3[type]*dtheta2 +
-      4.0*k4[type]*dtheta3 + 5.0*k5[type]*dtheta4 + 6.0*k6[type]*dtheta5;
-
-    a = -de_angle*s;
-    a11 = a*c / rsq1;
-    a12 = -a / (r1*r2);
-    a22 = a*c / rsq2;
-
-    f1[0] = a11*delx1 + a12*delx2;
-    f1[1] = a11*dely1 + a12*dely2;
-    f1[2] = a11*delz1 + a12*delz2;
-
-    f3[0] = a22*delx2 + a12*delx1;
-    f3[1] = a22*dely2 + a12*dely1;
-    f3[2] = a22*delz2 + a12*delz1;
-
-    if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + 
-                 k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
-
-    // 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::anglep(int i1, int i2, int i3, int type, int eflag)
+void AngleAmoeba::tinker_angle(int i1, int i2, int i3, int type, int eflag)
+{
+  double delx1,dely1,delz1,delx2,dely2,delz2;
+  double eangle,f1[3],f3[3];
+  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;
+
+  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;
+  
+  // force & energy for angle term
+
+  dtheta = acos(c) - theta0[type];
+  dtheta2 = dtheta*dtheta;
+  dtheta3 = dtheta2*dtheta;
+  dtheta4 = dtheta3*dtheta;
+  dtheta5 = dtheta4*dtheta;
+  dtheta6 = dtheta5*dtheta;
+
+  de_angle = 2.0*k2[type]*dtheta + 3.0*k3[type]*dtheta2 +
+    4.0*k4[type]*dtheta3 + 5.0*k5[type]*dtheta4 + 6.0*k6[type]*dtheta5;
+
+  a = -de_angle*s;
+  a11 = a*c / rsq1;
+  a12 = -a / (r1*r2);
+  a22 = a*c / rsq2;
+  
+  f1[0] = a11*delx1 + a12*delx2;
+  f1[1] = a11*dely1 + a12*dely2;
+  f1[2] = a11*delz1 + a12*delz2;
+  
+  f3[0] = a22*delx2 + a12*delx1;
+  f3[1] = a22*dely2 + a12*dely1;
+  f3[2] = a22*delz2 + a12*delz1;
+
+  eangle = 0.0;
+  if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + 
+               k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
+
+  // 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_anglep(int i1, int i2, int i3, int type, int eflag)
 {
   int i4;
   tagint i1tag,i3tag,i4tag;
@@ -203,7 +256,7 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  tagint **bond_atom = atom->bond_atom;
+  tagint **special = atom->special;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
@@ -213,7 +266,7 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
   i3tag = atom->tag[i3];
 
   for (int ibond = 0; ibond < 3; ibond++) {
-    i4tag = bond_atom[i2][ibond];
+    i4tag = special[i2][ibond];
     if (i4tag != i1tag && i4tag != i3tag) break;
   }
 
@@ -287,18 +340,10 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
   deddt = 2.0*k2[type]*dtheta + 3.0*k3[type]*dtheta2 +
     4.0*k4[type]*dtheta3 + 5.0*k5[type]*dtheta4 + 6.0*k6[type]*dtheta5;
 
+  eangle = 0.0;
   if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + 
                k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
 
-  /*
-  printf("ANGLEP: %d %d %d %d: eng %g\n",
-         atom->tag[i1],
-         atom->tag[i2],
-         atom->tag[i3],
-         atom->tag[i4],
-         eangle);
-  */
-
   // chain rule terms for first derivative components
 
   terma = -deddt / (rap2*rm);
@@ -376,12 +421,175 @@ void AngleAmoeba::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;
+  }
+
+  if (evflag) ev_tally2(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
+}
+
+/* ---------------------------------------------------------------------- */
+
 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");
@@ -389,8 +597,21 @@ 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(ub_k,n+1,"angle:ub_k");
+  memory->create(ub_r0,n+1,"angle:ub_r0");
+
   memory->create(setflag,n+1,"angle:setflag");
-  for (int i = 1; i <= n; i++) setflag[i] = 0;
+  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;
 }
 
 /* ----------------------------------------------------------------------
@@ -399,6 +620,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;
@@ -406,32 +628,85 @@ 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 ba_k1_one = utils::numeric(FLERR,arg[2],false,lmp);
-  double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
+    double ba_k2_one = utils::numeric(FLERR,arg[3],false,lmp);
-  double k2_one = utils::numeric(FLERR,arg[3],false,lmp);
+    double ba_r1_one = utils::numeric(FLERR,arg[4],false,lmp);
-  double k3_one = utils::numeric(FLERR,arg[4],false,lmp);
+    double ba_r2_one = utils::numeric(FLERR,arg[5],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++) {
+      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 (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 != 9) error->all(FLERR,"Incorrect args for angle coefficients");
+
+    int pflag_one = utils::inumeric(FLERR,arg[1],false,lmp);
+    int ubflag_one = utils::inumeric(FLERR,arg[2],false,lmp);
+    double theta0_one = utils::numeric(FLERR,arg[3],false,lmp);
+    double k2_one = utils::numeric(FLERR,arg[4],false,lmp);
+    double k3_one = utils::numeric(FLERR,arg[5],false,lmp);
+    double k4_one = utils::numeric(FLERR,arg[6],false,lmp);
+    double k5_one = utils::numeric(FLERR,arg[7],false,lmp);
+    double k6_one = utils::numeric(FLERR,arg[8],false,lmp);
+
+    // convert theta0 from degrees to radians
   
-  for (int i = ilo; i <= ihi; i++) {
+    for (int i = ilo; i <= ihi; i++) {
-    pflag[i] = pflag_one;
+      pflag[i] = pflag_one;
-    theta0[i] = theta0_one/180.0 * MY_PI;
+      ubflag[i] = ubflag_one;
-    k2[i] = k2_one;
+      theta0[i] = theta0_one/180.0 * MY_PI;
-    k3[i] = k3_one;
+      k2[i] = k2_one;
-    k4[i] = k4_one;
+      k3[i] = k3_one;
-    k5[i] = k5_one;
+      k4[i] = k4_one;
-    k6[i] = k6_one;
+      k5[i] = k5_one;
-    count++;
+      k6[i] = k6_one;
+      setflag_a[i] = 1;
+      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_ub[i]) 
+      setflag[i] = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+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));
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -448,12 +723,22 @@ 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);
   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);
+
+  fwrite(&ub_k[1],sizeof(double),atom->nangletypes,fp);
+  fwrite(&ub_r0[1],sizeof(double),atom->nangletypes,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -466,15 +751,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,&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);
@@ -482,6 +786,14 @@ 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);
+
+  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;
 }
 
@@ -492,8 +804,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",
+    fprintf(fp,"%d %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]);
+            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]);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -533,5 +854,11 @@ 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;
+
+  // NOTE: add UB term
+
   return energy;
 }

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 *);
@@ -37,10 +38,18 @@ 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;
+  double *ub_k, *ub_r0;
+  int *setflag_a, *setflag_ba, *setflag_ub;
 
-  void anglep(int, int, int, int, int);
+  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);
+  void tinker_urey_bradley(int, int, int, int);
   void allocate();
 };
 

src/AMOEBA/fix_amoeba_bitorsion.h

@@ -0,0 +1,142 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+// clang-format off
+FixStyle(amoeba/bitorsion,FixAmoebaBiTorsion);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_AMOEBA_BITORSION_H
+#define LMP_FIX_AMOEBA_BITORSION_H
+
+#include "fix.h"
+namespace LAMMPS_NS {
+
+class FixAmoebaBiTorsion : public Fix {
+ public:
+  FixAmoebaBiTorsion(class LAMMPS *, int, char **);
+  ~FixAmoebaBiTorsion();
+  int setmask();
+  void init();
+  void setup(int);
+  void setup_pre_neighbor();
+  void setup_pre_reverse(int, int);
+  void min_setup(int);
+  void pre_neighbor();
+  void pre_reverse(int, int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  void min_post_force(int);
+  double compute_scalar();
+
+  void read_data_header(char *);
+  void read_data_section(char *, int, char *, tagint);
+  bigint read_data_skip_lines(char *);
+  void write_data_header(FILE *, int);
+  void write_data_section_size(int, int &, int &);
+  void write_data_section_pack(int, double **);
+  void write_data_section_keyword(int, FILE *);
+  void write_data_section(int, FILE *, int, double **, int);
+
+  void write_restart(FILE *);
+  void restart(char *);
+  int pack_restart(int, double *);
+  void unpack_restart(int, int);
+  int size_restart(int);
+  int maxsize_restart();
+
+  void grow_arrays(int);
+  void copy_arrays(int, int, int);
+  void set_arrays(int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
+
+  double memory_usage();
+
+ private:
+  int nprocs, me;
+  int eflag_caller;
+  int ilevel_respa;
+  int disable;
+  bigint nbitorsions;
+  double ebitorsion;
+  double onefifth;
+
+  // per-atom data for bitorsions stored with each owned atom
+
+  int *num_bitorsion;
+  int **bitorsion_type;
+  tagint **bitorsion_atom1, **bitorsion_atom2, **bitorsion_atom3;
+  tagint **bitorsion_atom4, **bitorsion_atom5;
+
+  // previous max atoms on this proc before grow() is called
+
+  int nmax_previous;
+
+  // list of all bitorsions to compute on this proc
+
+  int nbitorsion_list;
+  int max_bitorsion_list;
+  int **bitorsion_list;
+
+  // BiTorsion grid data
+
+  int ntypes;
+  int *nxgrid,*nygrid;
+  double ****btgrid;
+
+  // read BiTorsion grid data
+
+  void read_grid_data(char *);
+};
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+UNDOCUMENTED
+
+E: CMAP atoms %d %d %d %d %d missing on proc %d at step %ld
+
+UNDOCUMENTED
+
+E: Invalid CMAP crossterm_type
+
+UNDOCUMENTED
+
+E: Cannot open fix cmap file %s
+
+UNDOCUMENTED
+
+E: CMAP: atan2 function cannot take 2 zero arguments
+
+UNDOCUMENTED
+
+E: Invalid read data header line for fix cmap
+
+UNDOCUMENTED
+
+E: Incorrect %s format in data file
+
+UNDOCUMENTED
+
+E: Too many CMAP crossterms for one atom
+
+UNDOCUMENTED
+
+*/

src/AMOEBA/fix_amoeba_pitorsion.h

@@ -0,0 +1,135 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+// clang-format off
+FixStyle(amoeba/pitorsion,FixAmoebaPiTorsion);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_AMOEBA_PITORSION_H
+#define LMP_FIX_AMOEBA_PITORSION_H
+
+#include "fix.h"
+namespace LAMMPS_NS {
+
+class FixAmoebaPiTorsion : public Fix {
+ public:
+  FixAmoebaPiTorsion(class LAMMPS *, int, char **);
+  ~FixAmoebaPiTorsion();
+  int setmask();
+  void init();
+  void setup(int);
+  void setup_pre_neighbor();
+  void setup_pre_reverse(int, int);
+  void min_setup(int);
+  void pre_neighbor();
+  void pre_reverse(int, int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  void min_post_force(int);
+  double compute_scalar();
+
+  void read_data_header(char *);
+  void read_data_section(char *, int, char *, tagint);
+  bigint read_data_skip_lines(char *);
+  void write_data_header(FILE *, int);
+  void write_data_section_size(int, int &, int &);
+  void write_data_section_pack(int, double **);
+  void write_data_section_keyword(int, FILE *);
+  void write_data_section(int, FILE *, int, double **, int);
+
+  void write_restart(FILE *);
+  void restart(char *);
+  int pack_restart(int, double *);
+  void unpack_restart(int, int);
+  int size_restart(int);
+  int maxsize_restart();
+
+  void grow_arrays(int);
+  void copy_arrays(int, int, int);
+  void set_arrays(int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
+
+  double memory_usage();
+
+ private:
+  int nprocs, me;
+  int eflag_caller;
+  int ilevel_respa;
+  int disable;
+  bigint npitorsions;
+  int npitorsion_types;
+  double epitorsion;
+  double onesixth;
+
+  double *kpit;
+
+  // per-atom data for pitorsions stored with each owned atom
+
+  int *num_pitorsion;
+  int **pitorsion_type;
+  tagint **pitorsion_atom1, **pitorsion_atom2, **pitorsion_atom3;
+  tagint **pitorsion_atom4, **pitorsion_atom5, **pitorsion_atom6;
+
+  // previous max atoms on this proc before grow() is called
+
+  int nmax_previous;
+
+  // list of all pitorsions to compute on this proc
+
+  int npitorsion_list;
+  int max_pitorsion_list;
+  int **pitorsion_list;
+};
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+UNDOCUMENTED
+
+E: CMAP atoms %d %d %d %d %d missing on proc %d at step %ld
+
+UNDOCUMENTED
+
+E: Invalid CMAP crossterm_type
+
+UNDOCUMENTED
+
+E: Cannot open fix cmap file %s
+
+UNDOCUMENTED
+
+E: CMAP: atan2 function cannot take 2 zero arguments
+
+UNDOCUMENTED
+
+E: Invalid read data header line for fix cmap
+
+UNDOCUMENTED
+
+E: Incorrect %s format in data file
+
+UNDOCUMENTED
+
+E: Too many CMAP crossterms for one atom
+
+UNDOCUMENTED
+
+*/

src/AMOEBA/improper_amoeba.cpp

@@ -0,0 +1,328 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "improper_amoeba.h"
+
+#include <cmath>
+#include "atom.h"
+#include "comm.h"
+#include "update.h"
+#include "neighbor.h"
+#include "force.h"
+#include "pair.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define TOLERANCE 0.05
+#define SMALL     0.001
+
+/* ---------------------------------------------------------------------- */
+
+ImproperAmoeba::ImproperAmoeba(LAMMPS *lmp) : Improper(lmp)
+{
+  writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ImproperAmoeba::~ImproperAmoeba()
+{
+  if (allocated && !copymode) {
+    memory->destroy(setflag);
+    memory->destroy(k);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+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;
+  double rad2,rcd2,rdb2,dot,cc,ee;
+  double sine,angle;
+  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 fa[3],fb[3],fc[3],fd[3];
+
+  ev_init(eflag,vflag);
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int **improperlist = neighbor->improperlist;
+  int nimproperlist = neighbor->nimproperlist;
+  int nlocal = atom->nlocal;
+  int newton_bond = force->newton_bond;
+
+  // conversion factors for radians to degrees and vice versa
+
+  double rad2degree = 180.0/MY_PI;
+  double prefactor = 1.0 / (rad2degree*rad2degree);
+
+  for (n = 0; n < nimproperlist; n++) {
+
+    // in Tinker code, atom1 = D, atom2 = B, atom3 = A, atom4 = C
+    // for Alligner angle:
+    //   atoms A,C,D form a plane, B is out-of-plane
+    //   angle is between plane and the vector from D to B
+
+    id = improperlist[n][0];
+    ib = improperlist[n][1];
+    ia = improperlist[n][2];
+    ic = improperlist[n][3];
+    type = improperlist[n][4];
+
+    // coordinates of the atoms at trigonal center
+
+    xia = x[ia][0];
+    yia = x[ia][1];
+    zia = x[ia][2];
+    xib = x[ib][0];
+    yib = x[ib][1];
+    zib = x[ib][2];
+    xic = x[ic][0];
+    yic = x[ic][1];
+    zic = x[ic][2];
+    xid = x[id][0];
+    yid = x[id][1];
+    zid = x[id][2];
+
+    // compute the out-of-plane bending angle
+
+    xab = xia - xib;
+    yab = yia - yib;
+    zab = zia - zib;
+    xcb = xic - xib;
+    ycb = yic - yib;
+    zcb = zic - zib;
+    xdb = xid - xib;
+    ydb = yid - yib;
+    zdb = zid - zib;
+    xad = xia - xid;
+    yad = yia - yid;
+    zad = zia - zid;
+    xcd = xic - xid;
+    ycd = yic - yid;
+    zcd = zic - zid;
+    
+    // Allinger angle between A-C-D plane and D-B vector for D-B < AC
+      
+    rad2 = xad*xad + yad*yad + zad*zad;
+    rcd2 = xcd*xcd + ycd*ycd + zcd*zcd;
+    dot = xad*xcd + yad*ycd + zad*zcd;
+    cc = rad2*rcd2 - dot*dot;
+
+    // find the out-of-plane angle bending energy
+
+    ee = xdb*(yab*zcb-zab*ycb) + ydb*(zab*xcb-xab*zcb) + zdb*(xab*ycb-yab*xcb);
+    rdb2 = xdb*xdb + ydb*ydb + zdb*zdb;
+    if (rdb2 == 0.0 || cc == 0.0) continue;
+
+    sine = fabs(ee) / sqrt(cc*rdb2);
+    sine = MIN(1.0,sine);
+
+    // angle needs to be in degrees for Tinker formulas
+    // b/c opbend_3456 coeffs are in mixed units
+
+    angle = rad2degree * asin(sine);
+    dt = angle;
+    dt2 = dt * dt;
+    dt3 = dt2 * dt;
+    dt4 = dt2 * dt2;
+    e = prefactor * k[type] * dt2 * (1.0 + opbend_cubic*dt + opbend_quartic*dt2 + 
+                         opbend_pentic*dt3 + opbend_sextic*dt4);
+
+    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 || id < nlocal) {
+      f[id][0] += fd[0];
+      f[id][1] += fd[1];
+      f[id][2] += fd[2];
+    }
+
+    if (newton_bond || ib < nlocal) {
+      f[ib][0] += fb[0];
+      f[ib][1] += fb[1];
+      f[ib][2] += fb[2];
+    }
+
+    if (newton_bond || ia < nlocal) {
+      f[ia][0] += fa[0];
+      f[ia][1] += fa[1];
+      f[ia][2] += fa[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(id,ib,ia,ic,nlocal,newton_bond,e,fd,fa,fc,
+               xdb,ydb,zdb,xab,yab,zab,xic-xia,yic-yia,zic-zia);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ImproperAmoeba::allocate()
+{
+  allocated = 1;
+  int n = atom->nimpropertypes;
+
+  memory->create(k,n+1,"improper:k");
+
+  memory->create(setflag,n+1,"improper:setflag");
+  for (int i = 1; i <= n; i++) setflag[i] = 0;
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one type
+------------------------------------------------------------------------- */
+
+void ImproperAmoeba::coeff(int narg, char **arg)
+{
+  if (narg != 2) error->all(FLERR,"Incorrect args for improper coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi;
+  utils::bounds(FLERR,arg[0],1,atom->nimpropertypes,ilo,ihi,error);
+
+  double k_one = utils::numeric(FLERR,arg[1],false,lmp);
+
+  // convert chi from degrees to radians
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    k[i] = k_one;
+    setflag[i] = 1;
+    count++;
+  }
+
+  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);
+
+  // check if PairAmoeba disabled improper terms
+
+  int tmp;
+  int flag = *((int *) pair->extract("improper_flag",tmp));
+  disable = flag ? 0 : 1;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes out coeffs to restart file
+------------------------------------------------------------------------- */
+
+void ImproperAmoeba::write_restart(FILE *fp)
+{
+  fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads coeffs from restart file, bcasts them
+------------------------------------------------------------------------- */
+
+void ImproperAmoeba::read_restart(FILE *fp)
+{
+  allocate();
+
+  if (comm->me == 0)
+    utils::sfread(FLERR,&k[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
+  MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
+
+  for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void ImproperAmoeba::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->nimpropertypes; i++)
+    fprintf(fp,"%d %g\n",i,k[i]);
+}

src/AMOEBA/improper_amoeba.h

@@ -0,0 +1,62 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef IMPROPER_CLASS
+// clang-format off
+ImproperStyle(amoeba,ImproperAmoeba);
+// clang-format on
+#else
+
+#ifndef LMP_IMPROPER_AMOEBA_H
+#define LMP_IMPROPER_AMOEBA_H
+
+#include "improper.h"
+
+namespace LAMMPS_NS {
+
+class ImproperAmoeba : public Improper {
+ public:
+  ImproperAmoeba(class LAMMPS *);
+  ~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:
+  int disable;
+  double opbend_cubic,opbend_quartic,opbend_pentic,opbend_sextic;
+  double *k;
+
+  virtual void allocate();
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+W: Improper problem: %d %ld %d %d %d %d
+
+Conformation of the 4 listed improper atoms is extreme; you may want
+to check your simulation geometry.
+
+E: Incorrect args for improper coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/

src/AMOEBA/pair_amoeba.cpp

@@ -41,12 +41,14 @@ using namespace LAMMPS_NS;
 enum{INDUCE,RSD,SETUP_AMOEBA,SETUP_HIPPO,KMPOLE,AMGROUP,PVAL};  // forward comm
 enum{FIELD,ZRSD,TORQUE,UFLD};                                   // reverse comm
 enum{ARITHMETIC,GEOMETRIC,CUBIC_MEAN,R_MIN,SIGMA,DIAMETER,HARMONIC,HHG,W_H};
-enum{VDWL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
+enum{HAL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
 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);
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -247,9 +247,18 @@ void PairAmoeba::compute(int eflag, int vflag)
   evdwl = 0.0;
   ev_init(eflag,vflag);
 
-  // zero internal energy, force, virial terms
+  // zero energy/virial components
 
-  zero_energy_force_virial();
+  ehal = erepulse = edisp = epolar = empole = eqxfer = 0.0;
+
+  for (int i = 0; i < 6; i++) {
+    virhal[i] = 0.0;
+    virrepulse[i] = 0.0;
+    virdisp[i] = 0.0;
+    virpolar[i] = 0.0;
+    virmpole[i] = 0.0;
+    virqxfer[i] = 0.0;
+  }
 
   // grow local vectors and arrays if necessary
 
@@ -313,7 +322,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 +386,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 +406,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
@@ -414,15 +423,6 @@ void PairAmoeba::compute(int eflag, int vflag)
   nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
-  for (int i = 0; i < nall; i++) {
-    f[i][0] = fhal[i][0] + frepulse[i][0] + fdisp[i][0] + 
-      fmpole[i][0] + fpolar[i][0] + fqxfer[i][0];
-    f[i][1] = fhal[i][1] + frepulse[i][1] + fdisp[i][1] + 
-      fmpole[i][1] + fpolar[i][1] + fqxfer[i][1];
-    f[i][2] = fhal[i][2] + frepulse[i][2] + fdisp[i][2] + 
-      fmpole[i][2] + fpolar[i][2] + fqxfer[i][2];
-  }
-
   for (int i = 0; i < 6; i++)
     virial[i] = virhal[i] + virrepulse[i] + virdisp[i] + 
       virpolar[i] + virmpole[i] + virqxfer[i];
@@ -481,20 +481,24 @@ void PairAmoeba::compute(int eflag, int vflag)
     utils::logmesg(lmp,"\nAMEOBA/HIPPO timing info:\n");
     utils::logmesg(lmp,"  Init   time: {:.6g} {:.6g}\n",
                    time_init,time_init/time_amtot);
-    utils::logmesg(lmp,"  Hal    time: {:.6g} {:.6g}\n",
+    if (amoeba)
-                   time_hal,time_hal/time_amtot*100);
+      utils::logmesg(lmp,"  Hal    time: {:.6g} {:.6g}\n",
-    utils::logmesg(lmp,"  Repls  time: {:.6g} {:.6g}\n",
+                     time_hal,time_hal/time_amtot*100);
-                   time_repulse,time_repulse/time_amtot*100);
+    if (hippo)
-    utils::logmesg(lmp,"  Disp   time: {:.6g} {:.6g}\n",
+      utils::logmesg(lmp,"  Repls  time: {:.6g} {:.6g}\n",
-                   time_disp,time_disp/time_amtot*100);
+                     time_repulse,time_repulse/time_amtot*100);
+    if (hippo)
+      utils::logmesg(lmp,"  Disp   time: {:.6g} {:.6g}\n",
+                     time_disp,time_disp/time_amtot*100);
     utils::logmesg(lmp,"  Mpole  time: {:.6g} {:.6g}\n",
                    time_mpole,time_mpole/time_amtot*100);
     utils::logmesg(lmp,"  Induce time: {:.6g} {:.6g}\n",
                    time_induce,time_induce/time_amtot*100);
     utils::logmesg(lmp,"  Polar  time: {:.6g} {:.6g}\n",
                    time_polar,time_polar/time_amtot*100);
-    utils::logmesg(lmp,"  Qxfer  time: {:.6g} {:.6g}\n",
+    if (hippo)
-                   time_qxfer,time_qxfer/time_amtot*100);
+      utils::logmesg(lmp,"  Qxfer  time: {:.6g} {:.6g}\n",
+                     time_qxfer,time_qxfer/time_amtot*100);
     utils::logmesg(lmp,"  Total  time: {:.6g}\n\n",time_amtot);
   }
 }
@@ -522,49 +526,65 @@ 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 = bitorsion_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 = bitorsion_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;
+
+  if (narg < 2) error->all(FLERR,"Illegal pair_style command");
+
+  // 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],"bitorsion") == 0) bitorsion_flag = newvalue;
     else error->all(FLERR,"Illegal pair_style command");
-    iarg++;
   }
 }
 
@@ -677,14 +697,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");
     
@@ -719,19 +734,6 @@ void PairAmoeba::init_style()
   if (index_pval < 0 || !flag || cols) 
     error->all(FLERR,"Pair amoeba pval is not defined");
 
-  // check for fix store/state commands that stores force components
-
-  ifhal = modify->find_fix("fhal");
-  ifrepulse = modify->find_fix("frepulse");
-  ifdisp = modify->find_fix("fdisp");
-  ifpolar = modify->find_fix("fpolar");
-  ifmpole = modify->find_fix("fmpole");
-  ifqxfer = modify->find_fix("fqxfer");
-
-  if (ifhal < 0 || ifrepulse < 0 || ifdisp < 0 || 
-      ifpolar < 0 || ifmpole < 0 || ifqxfer < 0)
-    error->all(FLERR,"Pair amoeba fix store/state commands not defined");
-
   // -------------------------------------------------------------------
   // one-time initializations
   // can't do earlier b/c need all atoms to exist
@@ -982,7 +984,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");
   }
 }
 
@@ -993,32 +995,41 @@ void PairAmoeba::init_style()
 void PairAmoeba::print_settings(FILE *fp)
 {
   fprintf(fp,"AMOEBA/HIPPO force field settings\n");
-  choose(VDWL);
+
-  fprintf(fp,"  vdwl: cut %g taper %g vscale %g %g %g %g\n",
+  if (amoeba) {
-	  sqrt(off2),sqrt(cut2),
+    choose(HAL);
-	  special_hal[1],special_hal[2],special_hal[3],special_hal[4]);
+    fprintf(fp,"  hal: cut %g taper %g vscale %g %g %g %g\n",
+            sqrt(off2),sqrt(cut2),
+            special_hal[1],special_hal[2],special_hal[3],special_hal[4]);
+  }
+
+  if (hippo) {
+    choose(REPULSE);
+    fprintf(fp,"  repulsion: cut %g taper %g rscale %g %g %g %g\n",
+            sqrt(off2),sqrt(cut2),
+            special_repel[1],special_repel[2],special_repel[3],special_repel[4]);
+  }
   
-  choose(REPULSE);
+  if (hippo) {
-  fprintf(fp,"  repulsion: cut %g taper %g rscale %g %g %g %g\n",
+    choose(QFER);
-	  sqrt(off2),sqrt(cut2),
+    fprintf(fp,"  qxfer: cut %g taper %g mscale %g %g %g %g\n",
-	  special_repel[1],special_repel[2],special_repel[3],special_repel[4]);
+            sqrt(off2),sqrt(cut2),
+            special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
+  }
   
-  choose(QFER);
+  if (hippo) {
-  fprintf(fp,"  qxfer: cut %g taper %g mscale %g %g %g %g\n",
+    if (use_dewald) {
-	  sqrt(off2),sqrt(cut2),
+      choose(DISP_LONG);
-	  special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
+      fprintf(fp,"  dispersion: cut %g aewald %g bsorder %d "
-  
+              "FFT %d %d %d dspscale %g %g %g %g\n",
-  if (use_dewald) {
+              sqrt(off2),aewald,bsdorder,ndfft1,ndfft2,ndfft3,
-    choose(DISP_LONG);
+              special_disp[1],special_disp[2],special_disp[3],special_disp[4]);
-    fprintf(fp,"  dispersion: cut %g aewald %g bsorder %d "
+    } else {
-	    "FFT %d %d %d dspscale %g %g %g %g\n",
+      choose(DISP);
-	    sqrt(off2),aewald,bsdorder,ndfft1,ndfft2,ndfft3,
+      fprintf(fp,"  dispersion: cut %g aewald %g dspscale %g %g %g %g\n",
-	    special_disp[1],special_disp[2],special_disp[3],special_disp[4]);
+              sqrt(off2),aewald,
-  } else {
+              special_disp[1],special_disp[2],special_disp[3],special_disp[4]);
-    choose(DISP);
+    }
-    fprintf(fp,"  dispersion: cut %g aewald %g dspscale %g %g %g %g\n",
-	    sqrt(off2),aewald,
-	    special_disp[1],special_disp[2],special_disp[3],special_disp[4]);
   }
   
   if (use_ewald) {
@@ -1073,35 +1084,31 @@ double PairAmoeba::init_one(int i, int j)
 {
   double cutoff = 0.0;
 
-  if (hal_flag) {
+  if (amoeba) {
-    choose(VDWL);
+    choose(HAL);
     cutoff = MAX(cutoff,sqrt(off2));
   }
-  if (repulse_flag) {
+
+  if (hippo) {
     choose(REPULSE);
     cutoff = MAX(cutoff,sqrt(off2));
   }
-  if (disp_flag) {
+
+  if (hippo) {
     if (use_dewald) choose(DISP_LONG);
     else choose(DISP);
     cutoff = MAX(cutoff,sqrt(off2));
   }
-  if (mpole_flag) {
+
-    if (use_ewald) choose(MPOLE_LONG);
+  if (use_ewald) choose(MPOLE_LONG);
-    else choose(MPOLE);
+  else choose(MPOLE);
-    cutoff = MAX(cutoff,sqrt(off2));
+  cutoff = MAX(cutoff,sqrt(off2));
-  }
+
-  if (induce_flag) {
+  if (use_ewald) choose(POLAR_LONG);
-    if (use_ewald) choose(POLAR_LONG);
+  else choose(POLAR);
-    else choose(POLAR);
+  cutoff = MAX(cutoff,sqrt(off2));
-    cutoff = MAX(cutoff,sqrt(off2));
+
-  }
+  if (hippo) {
-  if (polar_flag) {
-    if (use_ewald) choose(POLAR_LONG);
-    else choose(POLAR);
-    cutoff = MAX(cutoff,sqrt(off2));
-  }
-  if (qxfer_flag) {
     choose(QFER);
     cutoff = MAX(cutoff,sqrt(off2));
   }
@@ -1927,7 +1934,7 @@ void PairAmoeba::choose(int which)
 
   // short-range only terms
   
-  if (which == VDWL) {
+  if (which == HAL) {
     off = vdwcut;
     cut = vdwtaper;
   } else if (which == REPULSE) {
@@ -2075,39 +2082,24 @@ void PairAmoeba::mix()
   }
 }
 
-/* ----------------------------------------------------------------------
+/* ---------------------------------------------------------------------- */
-   zero internal energy and virial terms
-------------------------------------------------------------------------- */
 
-void PairAmoeba::zero_energy_force_virial()
+void *PairAmoeba::extract(const char *str, int &dim)
 {
-  ehal = erepulse = edisp = epolar = empole = eqxfer = 0.0;
+  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,"bitorsion_flag") == 0) return (void *) &bitorsion_flag;
 
-  fhal = modify->fix[ifhal]->array_atom;
+  if (strcmp(str,"opbend_cubic") == 0) return (void *) &opbend_cubic;
-  frepulse = modify->fix[ifrepulse]->array_atom;
+  if (strcmp(str,"opbend_quartic") == 0) return (void *) &opbend_quartic;
-  fdisp = modify->fix[ifdisp]->array_atom;
+  if (strcmp(str,"opbend_pentic") == 0) return (void *) &opbend_pentic;
-  fpolar = modify->fix[ifpolar]->array_atom;
+  if (strcmp(str,"opbend_sextic") == 0) return (void *) &opbend_sextic;
-  fmpole = modify->fix[ifmpole]->array_atom;
+  return nullptr;
-  fqxfer = modify->fix[ifqxfer]->array_atom;
-
-  int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
-
-  memset(&fhal[0][0],0,3*nall*sizeof(double));
-  memset(&frepulse[0][0],0,3*nall*sizeof(double));
-  memset(&fdisp[0][0],0,3*nall*sizeof(double));
-  memset(&fpolar[0][0],0,3*nall*sizeof(double));
-  memset(&fmpole[0][0],0,3*nall*sizeof(double));
-  memset(&fqxfer[0][0],0,3*nall*sizeof(double));
-
-  for (int i = 0; i < 6; i++) {
-    virhal[i] = 0.0;
-    virrepulse[i] = 0.0;
-    virdisp[i] = 0.0;
-    virpolar[i] = 0.0;
-    virmpole[i] = 0.0;
-    virqxfer[i] = 0.0;
-  }
 }
 
 /* ----------------------------------------------------------------------

src/AMOEBA/pair_amoeba.h

@@ -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
@@ -73,23 +75,21 @@ 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,bitorsion_flag;
 
   // DEBUG timers
 
   double time_init,time_hal,time_repulse,time_disp;
   double time_mpole,time_induce,time_polar,time_qxfer;
   
-  // energy, force, and virial components
+  // energy/virial components
 
-  int ifhal,ifrepulse,ifdisp,ifpolar,ifmpole,ifqxfer;
   double ehal,erepulse,edisp,epolar,empole,eqxfer;
-  double **fhal,**frepulse,**fdisp,**fpolar,**fmpole,**fqxfer;
   double virhal[6],virrepulse[6],virdisp[6],virpolar[6],virmpole[6],virqxfer[6];
 
   // scalar values defined in force-field file

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

src/angle.cpp

@@ -179,6 +179,7 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
 /* ----------------------------------------------------------------------
    tally energy and virial into global and per-atom accumulators
    virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3
+   called by standard 3-body angles
 ------------------------------------------------------------------------- */
 
 void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
@@ -341,6 +342,7 @@ void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
 /* ----------------------------------------------------------------------
    tally energy and virial into global and per-atom accumulators
    virial = r1F1 + r2F2 + r3F3 + r4F4
+   called by AngleAmoeba for its 4-body angle term
 ------------------------------------------------------------------------- */
 
 void Angle::ev_tally4(int i, int j, int k, int m, int nlocal, int newton_bond,
@@ -437,6 +439,92 @@ void Angle::ev_tally4(int i, int j, int k, int m, int nlocal, int newton_bond,
   }
 }
 
+
+/* ----------------------------------------------------------------------
+   tally energy and virial into global and per-atom accumulators
+   called by AngleAmoeba for its 2-body Urey-Bradley H-H bond term
+   identical to Bond:ev_tally()
+------------------------------------------------------------------------- */
+
+void Angle::ev_tally2(int i, int j, int nlocal, int newton_bond,
+                      double ebond, double fbond,
+                      double delx, double dely, double delz)
+{
+  double ebondhalf,v[6];
+
+  if (eflag_either) {
+    if (eflag_global) {
+      if (newton_bond) energy += ebond;
+      else {
+        ebondhalf = 0.5*ebond;
+        if (i < nlocal) energy += ebondhalf;
+        if (j < nlocal) energy += ebondhalf;
+      }
+    }
+    if (eflag_atom) {
+      ebondhalf = 0.5*ebond;
+      if (newton_bond || i < nlocal) eatom[i] += ebondhalf;
+      if (newton_bond || j < nlocal) eatom[j] += ebondhalf;
+    }
+  }
+
+  if (vflag_either) {
+    v[0] = delx*delx*fbond;
+    v[1] = dely*dely*fbond;
+    v[2] = delz*delz*fbond;
+    v[3] = delx*dely*fbond;
+    v[4] = delx*delz*fbond;
+    v[5] = dely*delz*fbond;
+
+    if (vflag_global) {
+      if (newton_bond) {
+        virial[0] += v[0];
+        virial[1] += v[1];
+        virial[2] += v[2];
+        virial[3] += v[3];
+        virial[4] += v[4];
+        virial[5] += v[5];
+      } else {
+        if (i < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+        if (j < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+      }
+    }
+
+    if (vflag_atom) {
+      if (newton_bond || i < nlocal) {
+        vatom[i][0] += 0.5*v[0];
+        vatom[i][1] += 0.5*v[1];
+        vatom[i][2] += 0.5*v[2];
+        vatom[i][3] += 0.5*v[3];
+        vatom[i][4] += 0.5*v[4];
+        vatom[i][5] += 0.5*v[5];
+      }
+      if (newton_bond || j < nlocal) {
+        vatom[j][0] += 0.5*v[0];
+        vatom[j][1] += 0.5*v[1];
+        vatom[j][2] += 0.5*v[2];
+        vatom[j][3] += 0.5*v[3];
+        vatom[j][4] += 0.5*v[4];
+        vatom[j][5] += 0.5*v[5];
+      }
+    }
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 double Angle::memory_usage()

src/angle.h

@@ -78,6 +78,7 @@ class Angle : protected Pointers {
   void ev_tally(int, int, int, int, int, double, double *, double *, double, double, double, double,
                 double, double);
   void ev_tally4(int, int, int, int, int, int, double, double *, double *, double *, double *);
+  void ev_tally2(int, int, int, int, double, double, double, double, double);
 };
 
 }    // namespace LAMMPS_NS

src/modify.cpp

@@ -810,7 +810,8 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix)
 
   const char *exceptions[] =
     {"GPU", "OMP", "INTEL", "property/atom", "cmap", "cmap3", "rx",
-     "deprecated", "STORE/KIM", nullptr};
+     "deprecated", "STORE/KIM", "amoeba/pitorsion", "amoeba/bitorsion", 
+     nullptr};
 
   if (domain->box_exist == 0) {
     int m;

src/read_data.cpp

@@ -658,6 +658,13 @@ void ReadData::command(int narg, char **arg)
           error->all(FLERR,"Must define angle_style before BondAngle Coeffs");
         if (firstpass) anglecoeffs(2);
         else skip_lines(nangletypes);
+      } else if (strcmp(keyword,"UreyBradley Coeffs") == 0) {
+        if (atom->avec->angles_allow == 0)
+          error->all(FLERR,"Invalid data file section: UreyBradley Coeffs");
+        if (force->angle == nullptr)
+          error->all(FLERR,"Must define angle_style before UreyBradley Coeffs");
+        if (firstpass) anglecoeffs(3);
+        else skip_lines(nangletypes);
 
       } else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) {
         if (atom->avec->dihedrals_allow == 0)
@@ -722,7 +729,6 @@ void ReadData::command(int narg, char **arg)
             if (firstpass) fix(fix_index[i],keyword);
             else skip_lines(modify->fix[fix_index[i]]->
                             read_data_skip_lines(keyword));
-            parse_keyword(0);
             break;
           }
         if (i == nfix)
@@ -1852,6 +1858,7 @@ void ReadData::anglecoeffs(int which)
     if (which == 0) parse_coeffs(buf,nullptr,0,1,aoffset);
     else if (which == 1) parse_coeffs(buf,"bb",0,1,aoffset);
     else if (which == 2) parse_coeffs(buf,"ba",0,1,aoffset);
+    else if (which == 3) parse_coeffs(buf,"ub",0,1,aoffset);
     if (ncoeffarg == 0) error->all(FLERR,"Unexpected empty line in AngleCoeffs section");
     force->angle->coeff(ncoeffarg,coeffarg);
     buf = next + 1;

tools/tinker/data.py

@@ -0,0 +1,402 @@
+# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html
+# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+#
+# Copyright (2005) Sandia Corporation.  Under the terms of Contract
+# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+# certain rights in this software.  This software is distributed under 
+# the GNU General Public License.
+
+# data tool
+
+oneline = "Read, write, manipulate LAMMPS data files"
+
+docstr = """
+d = data("data.poly")            read a LAMMPS data file, can be gzipped
+d = data()			 create an empty data file
+
+d.map(1,"id",3,"x")              assign names to atom columns (1-N)
+
+coeffs = d.get("Pair Coeffs")    extract info from data file section
+q = d.get("Atoms",4)
+
+  1 arg = all columns returned as 2d array of floats
+  2 args = Nth column returned as vector of floats
+
+d.reorder("Atoms",1,3,2,4,5)     reorder columns (1-N) in a data file section
+
+  1,3,2,4,5 = new order of previous columns, can delete columns this way
+
+d.title = "My LAMMPS data file"	 set title of the data file
+d.headers["atoms"] = 1500        set a header value
+d.sections["Bonds"] = lines      set a section to list of lines (with newlines)
+d.delete("bonds")		 delete a keyword or section of data file
+d.delete("Bonds")
+d.replace("Atoms",5,vec)      	 replace Nth column of section with vector
+d.newxyz(dmp,1000)		 replace xyz in Atoms with xyz of snapshot N
+
+  newxyz assumes id,x,y,z are defined in both data and dump files
+    also replaces ix,iy,iz if they are defined
+  
+index,time,flag = d.iterator(0/1)          loop over single data file snapshot
+time,box,atoms,bonds,tris,lines = d.viz(index)   return list of viz objects
+
+  iterator() and viz() are compatible with equivalent dump calls
+  iterator() called with arg = 0 first time, with arg = 1 on subsequent calls
+    index = timestep index within dump object (only 0 for data file)
+    time = timestep value (only 0 for data file)
+    flag = -1 when iteration is done, 1 otherwise
+  viz() returns info for specified timestep index (must be 0)
+    time = 0
+    box = [xlo,ylo,zlo,xhi,yhi,zhi]
+    atoms = id,type,x,y,z for each atom as 2d array
+    bonds = id,type,x1,y1,z1,x2,y2,z2,t1,t2 for each bond as 2d array
+      NULL if bonds do not exist
+    tris = NULL
+    lines = NULL
+    
+d.write("data.new")             write a LAMMPS data file
+"""
+
+# History
+#   8/05, Steve Plimpton (SNL): original version
+#   11/07, added triclinic box support
+
+# ToDo list
+
+# Variables
+#   title = 1st line of data file
+#   names = dictionary with atom attributes as keys, col #s as values 
+#   headers = dictionary with header name as key, value or tuple as values
+#   sections = dictionary with section name as key, array of lines as values
+#   nselect = 1 = # of snapshots
+
+# Imports and external programs
+
+from os import popen
+
+try: tmp = PIZZA_GUNZIP
+except: PIZZA_GUNZIP = "gunzip"
+
+# Class definition
+
+class data:
+
+  # --------------------------------------------------------------------
+
+  def __init__(self,*list):
+    self.nselect = 1
+    
+    if len(list) == 0:
+      self.title = "LAMMPS data file"
+      self.names = {}
+      self.headers = {}
+      self.sections = {}
+      return
+
+    file = list[0]
+    if file[-3:] == ".gz": f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r')
+    else: f = open(file)
+
+    self.title = f.readline()
+    self.names = {}
+    
+    headers = {}
+    while 1:
+      line = f.readline()
+      line = line.strip()
+      if len(line) == 0:
+        continue
+      found = 0
+      for keyword in hkeywords:
+        if line.find(keyword) >= 0:
+	  found = 1
+	  words = line.split()
+	  if keyword == "xlo xhi" or keyword == "ylo yhi" or \
+	    keyword == "zlo zhi":
+	    headers[keyword] = (float(words[0]),float(words[1]))
+	  elif keyword == "xy xz yz":
+	    headers[keyword] = \
+              (float(words[0]),float(words[1]),float(words[2]))
+          else:
+	    headers[keyword] = int(words[0])
+      if not found:
+        break
+
+    sections = {}
+    while 1:
+      found = 0
+      for pair in skeywords:
+        keyword,length = pair[0],pair[1]
+        if keyword == line:
+	  found = 1
+          if not headers.has_key(length):
+            raise StandardError, \
+                  "data section %s has no matching header value" % line
+	  f.readline()
+          list = []
+          for i in xrange(headers[length]): list.append(f.readline())
+          sections[keyword] = list
+      if not found:
+        raise StandardError,"invalid section %s in data file" % line
+      f.readline()
+      line = f.readline()
+      if not line:
+        break
+      line = line.strip()
+      
+    f.close()
+    self.headers = headers
+    self.sections = sections
+
+  # --------------------------------------------------------------------
+  # assign names to atom columns
+
+  def map(self,*pairs):
+    if len(pairs) % 2 != 0:
+      raise StandardError, "data map() requires pairs of mappings"
+    for i in range(0,len(pairs),2):
+      j = i + 1
+      self.names[pairs[j]] = pairs[i]-1
+
+  # --------------------------------------------------------------------
+  # extract info from data file fields
+
+  def get(self,*list):
+    if len(list) == 1:
+      field = list[0]
+      array = []
+      lines = self.sections[field]
+      for line in lines:
+        words = line.split()
+        values = map(float,words)
+        array.append(values)
+      return array
+    elif len(list) == 2:
+      field = list[0]
+      n = list[1] - 1
+      vec = []
+      lines = self.sections[field]
+      for line in lines:
+        words = line.split()
+        vec.append(float(words[n]))
+      return vec
+    else:
+      raise StandardError, "invalid arguments for data.get()"
+
+  # --------------------------------------------------------------------
+  # reorder columns in a data file field
+
+  def reorder(self,name,*order):
+    n = len(order)
+    natoms = len(self.sections[name])
+    oldlines = self.sections[name]
+    newlines = natoms*[""]
+    for index in order:
+      for i in xrange(len(newlines)):
+        words = oldlines[i].split()
+        newlines[i] += words[index-1] + " "
+    for i in xrange(len(newlines)):
+      newlines[i] += "\n"
+    self.sections[name] = newlines
+
+  # --------------------------------------------------------------------
+  # replace a column of named section with vector of values
+
+  def replace(self,name,icol,vector):
+    lines = self.sections[name]
+    newlines = []
+    j = icol - 1
+    for i in xrange(len(lines)):
+      line = lines[i]
+      words = line.split()
+      words[j] = str(vector[i])
+      newline = ' '.join(words) + '\n'
+      newlines.append(newline)
+    self.sections[name] = newlines
+
+  # --------------------------------------------------------------------
+  # replace x,y,z in Atoms with x,y,z values from snapshot ntime of dump object
+  # assumes id,x,y,z are defined in both data and dump files
+  # also replaces ix,iy,iz if they are defined
+
+  def newxyz(self,dm,ntime):
+    nsnap = dm.findtime(ntime)
+
+    dm.sort(ntime)
+    x,y,z = dm.vecs(ntime,"x","y","z")
+
+    self.replace("Atoms",self.names['x']+1,x)
+    self.replace("Atoms",self.names['y']+1,y)
+    self.replace("Atoms",self.names['z']+1,z)
+    
+    if dm.names.has_key("ix") and self.names.has_key("ix"):
+      ix,iy,iz = dm.vecs(ntime,"ix","iy","iz")
+      self.replace("Atoms",self.names['ix']+1,ix)
+      self.replace("Atoms",self.names['iy']+1,iy)
+      self.replace("Atoms",self.names['iz']+1,iz)
+      
+  # --------------------------------------------------------------------
+  # delete header value or section from data file
+
+  def delete(self,keyword):
+
+    if self.headers.has_key(keyword): del self.headers[keyword]
+    elif self.sections.has_key(keyword): del self.sections[keyword]
+    else: raise StandardError, "keyword not found in data object"
+
+  # --------------------------------------------------------------------
+  # write out a LAMMPS data file
+
+  def write(self,file):
+    f = open(file,"w")
+    print >>f,self.title
+ 
+    # write any keywords in standard list hkeywords
+    #   in the order they are in hkeywords
+    # then write any extra keywords at end of header section
+   
+    for keyword in hkeywords:
+      if self.headers.has_key(keyword):
+        if keyword == "xlo xhi" or keyword == "ylo yhi" or \
+               keyword == "zlo zhi":
+	  pair = self.headers[keyword]
+	  print >>f,pair[0],pair[1],keyword
+        elif keyword == "xy xz yz":
+	  triple = self.headers[keyword]
+	  print >>f,triple[0],triple[1],triple[2],keyword
+        else:
+	  print >>f,self.headers[keyword],keyword
+
+    for keyword in self.headers.keys():
+      if keyword not in hkeywords:
+        print >>f,self.headers[keyword],keyword
+
+    # write any sections in standard list skeywords
+    #   in the order they are in skeywords
+    # then write any extra sections at end of file
+
+    for pair in skeywords:
+      keyword = pair[0]
+      if self.sections.has_key(keyword):
+        print >>f,"\n%s\n" % keyword
+        for line in self.sections[keyword]:
+	  print >>f,line,
+
+    skeyfirst = [pair[0] for pair in skeywords]
+    
+    for keyword in self.sections.keys():
+      if keyword not in skeyfirst:
+        print >>f,"\n%s\n" % keyword
+        for line in self.sections[keyword]:
+	  print >>f,line,
+
+    f.close()
+
+  # --------------------------------------------------------------------
+  # iterator called from other tools
+
+  def iterator(self,flag):
+    if flag == 0: return 0,0,1
+    return 0,0,-1
+
+  # --------------------------------------------------------------------
+  # time query from other tools
+
+  def findtime(self,n):
+    if n == 0: return 0
+    raise StandardError, "no step %d exists" % (n)
+   
+  # --------------------------------------------------------------------
+  # return list of atoms and bonds to viz for data object
+
+  def viz(self,isnap):
+    if isnap: raise StandardError, "cannot call data.viz() with isnap != 0"
+    
+    id = self.names["id"]
+    type = self.names["type"]
+    x = self.names["x"]
+    y = self.names["y"]
+    z = self.names["z"]
+    
+    xlohi = self.headers["xlo xhi"]
+    ylohi = self.headers["ylo yhi"]
+    zlohi = self.headers["zlo zhi"]
+    box = [xlohi[0],ylohi[0],zlohi[0],xlohi[1],ylohi[1],zlohi[1]]
+
+    # create atom list needed by viz from id,type,x,y,z
+
+    atoms = []
+    atomlines = self.sections["Atoms"]
+    for line in atomlines:
+      words = line.split()
+      atoms.append([int(words[id]),int(words[type]),
+                    float(words[x]),float(words[y]),float(words[z])])
+
+    # create list of current bond coords from list of bonds
+    # assumes atoms are sorted so can lookup up the 2 atoms in each bond
+
+    bonds = []
+    if self.sections.has_key("Bonds"):
+      bondlines = self.sections["Bonds"]
+      for line in bondlines:
+        words = line.split()
+        bid,btype   = int(words[0]),int(words[1])
+        atom1,atom2 = int(words[2]),int(words[3])
+        atom1words  = atomlines[atom1-1].split()
+        atom2words  = atomlines[atom2-1].split()
+        bonds.append([bid,btype,
+                      float(atom1words[x]),float(atom1words[y]),
+                      float(atom1words[z]),
+                      float(atom2words[x]),float(atom2words[y]),
+                      float(atom2words[z]),
+		      float(atom1words[type]),float(atom2words[type])])
+      
+    tris = []
+    lines = []
+    return 0,box,atoms,bonds,tris,lines
+
+  # --------------------------------------------------------------------
+  # return box size
+
+  def maxbox(self):
+    xlohi = self.headers["xlo xhi"]
+    ylohi = self.headers["ylo yhi"]
+    zlohi = self.headers["zlo zhi"]
+    return [xlohi[0],ylohi[0],zlohi[0],xlohi[1],ylohi[1],zlohi[1]]
+
+  # --------------------------------------------------------------------
+  # return number of atom types
+
+  def maxtype(self):
+    return self.headers["atom types"]
+
+# --------------------------------------------------------------------
+# standard data file keywords, both header and main sections
+
+hkeywords = ["atoms","ellipsoids","lines","triangles","bodies",
+             "bonds","angles","dihedrals","impropers",
+	     "atom types","bond types","angle types","dihedral types",
+	     "improper types",
+             "xlo xhi","ylo yhi","zlo zhi","xy xz yz"]
+
+skeywords = [["Masses","atom types"],
+             ["Atoms","atoms"],["Ellipsoids","ellipsoids"],
+             ["Lines","lines"],["Triangles","triangles"],["Bodies","bodies"],
+             ["Velocities","atoms"],
+             ["Bonds","bonds"],
+	     ["Angles","angles"],
+             ["Dihedrals","dihedrals"],
+	     ["Impropers","impropers"],
+             ["Pair Coeffs","atom types"],
+	     ["Bond Coeffs","bond types"],
+             ["Angle Coeffs","angle types"],
+	     ["Dihedral Coeffs","dihedral types"],
+	     ["Improper Coeffs","improper types"],
+             ["BondBond Coeffs","angle types"],
+             ["BondAngle Coeffs","angle types"],
+             ["MiddleBondTorsion Coeffs","dihedral types"],
+             ["EndBondTorsion Coeffs","dihedral types"],
+             ["AngleTorsion Coeffs","dihedral types"],
+             ["AngleAngleTorsion Coeffs","dihedral types"],
+             ["BondBond13 Coeffs","dihedral types"],
+             ["AngleAngle Coeffs","improper types"]]

tools/tinker/tinker2lmp.py

@@ -7,12 +7,13 @@
 #   -amoeba file = AMOEBA PRM force field file name (required, or hippo)
 #   -hippo file = HIPPO PRM force field file name (required, or amoeba)
 #   -data file = LAMMPS data file to output (required)
+#   -bitorsion file = LAMMPS fix bitorsion file to output (required if BiTorsions)
 #   -nopbc = non-periodic system (default)
 #   -pbc xhi yhi zhi = periodic system from 0 to hi in each dimension (optional)
 
+# Author: Steve Plimpton
+
 import sys,os,math
-path = os.environ["LAMMPS_PYTHON_TOOLS"]
-sys.path.append(path)
 from data import data
 
 BIG = 1.0e20
@@ -31,10 +32,11 @@ def error(txt=""):
     print "  -amoeba file"
     print "  -hippo file"
     print "  -data file"
+    print "  -bitorsion file"
     print "  -nopbc"
     print "  -pbc xhi yhi zhi"
   else: print "ERROR:",txt
-  sys.exit()
+  #sys.exit()
 
 # read and store values from a Tinker xyz file
 
@@ -105,8 +107,22 @@ class PRMfile:
     self.angleparams = self.angles()
     self.bondangleparams = self.bondangles()
     self.torsionparams = self.torsions()
+    self.opbendparams = self.opbend()
+    self.ureyparams = self.ureybonds()
+    self.pitorsionparams = self.pitorsions()
+    self.bitorsionparams = self.bitorsions()
     self.ntypes = len(self.masses)
 
+  # find a section in the PRM file
+  
+  def find_section(self,txt):
+    txt = "##  %s  ##" % txt
+    for iline,line in enumerate(self.lines):
+      if txt in line: return iline
+    return -1
+
+  # scalar params
+  
   def force_field_definition(self):
     iline = self.find_section("Force Field Definition")
     iline += 3
@@ -121,7 +137,9 @@ class PRMfile:
         elif words[0] == "angle-pentic": self.angle_pentic = float(words[1])
         elif words[0] == "angle-sextic": self.angle_sextic = float(words[1])
       iline += 1
-      
+
+  # atom classes and masses
+  
   def peratom(self):
     classes = []
     masses = []
@@ -139,6 +157,8 @@ class PRMfile:
       iline += 1
     return classes,masses
 
+  # polarization groups
+  
   def polarize(self):
     polgroup = []
     iline = self.find_section("Dipole Polarizability Parameters")
@@ -189,6 +209,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
@@ -219,17 +240,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]))
@@ -248,14 +269,14 @@ class PRMfile:
   
   def bondangles(self):
     params = []
-    iline = self.find_section("Stretch Bend Parameters")
+    iline = self.find_section("Stretch-Bend Parameters")
     if iline < 0: return params
     iline += 3
 
     bdict = {}
-    for m,params in enumerate(self.bondparams):
+    for m,bparams in enumerate(self.bondparams):
-      bdict[(params[0],params[1])] = params[2]
+      bdict[(bparams[0],bparams[1])] = bparams[2]
-    
+
     while iline < self.nlines:
       words = self.lines[iline].split()
       if len(words):
@@ -268,18 +289,33 @@ class PRMfile:
           value2 = float(words[5])
           lmp1 = value1
           lmp2 = value2
-          lmp3 = lmp4 = 0.0
+          
-          if (class2,class1) in bdict: lmp3 = bdict[(class2,class1)]
+          if (class1,class2) in bdict:
-          if (class1,class2) in bdict: lmp3 = bdict[(class1,class2)]
+            lmp3 = bdict[(class1,class2)]
-          if (class2,class3) in bdict: lmp4 = bdict[(class2,class3)]
+          elif (class2,class1) in bdict:
-          if (class3,class2) in bdict: lmp4 = bdict[(class3,class2)]
+            lmp3 = bdict[(class2,class1)]
-          if lmp3 == 0.0 or lmp4 == 0.0:
+          else:
-            print "Bond in BondAngle term not found",class1,class2,class3
+            error("1st bond in BondAngle term not found: %d %d %d" % \
-            sys.exit()
+                  (class1,class2,class3))
+            # NOTE: just for debugging
+            lmp3 = 0.0
+
+          if (class2,class3) in bdict:
+            lmp4 = bdict[(class2,class3)]
+          elif (class3,class2) in bdict:
+            lmp4 = bdict[(class3,class2)]
+          else:
+            error("2nd bond in BondAngle term not found: %d %d %d" % \
+                  (class1,class2,class3))
+            # NOTE: just for debugging
+            lmp4 = 0.0
+            
           params.append((class1,class2,class3,lmp1,lmp2,lmp3,lmp4))
       iline += 1
     return params
 
+  # dihedral interactions
+  
   def torsions(self):
     params = []
     iline = self.find_section("Torsional Parameters")
@@ -294,26 +330,136 @@ class PRMfile:
           class2 = int(words[2])
           class3 = int(words[3])
           class4 = int(words[4])
-          value1 = words[5]
+          
-          value2 = words[6]
+          if len(words) <= 5:
-          value3 = words[7]
+            error("torsion has no params: %d %d %d %d" % \
-          value4 = words[8]
+                  (class1,class2,class3,class4))
-          value5 = words[9]
+          if (len(words)-5) % 3: 
-          value6 = words[10]
+            error("torsion does not have triplets of params: %d %d %d %d" % \
-          value7 = words[11]
+                  (class1,class2,class3,class4))
-          value8 = words[12]
+
-          value9 = words[13]
+          mfourier = (len(words)-5) / 3
-          params.append((class1,class2,class3,class4,
+          oneparams = [class1,class2,class3,class4,mfourier]
-                         value1,value2,value3,value4,value5,
+
-                         value6,value7,value8,value9))
+          for iset in range(mfourier):
+            value1 = float(words[5 + iset*3 + 0])
+            value2 = float(words[5 + iset*3 + 1])
+            value3 = int(words[5 + iset*3 + 2])
+            lmp1 = value1/2.0
+            lmp2 = value3
+            lmp3 = value2
+            oneparams += [lmp1,lmp2,lmp3]
+            
+          params.append(oneparams)
+      iline += 1
+    return params
+
+  # improper or out-of-plane bend interactions
+  
+  def opbend(self):
+    params = []
+    iline = self.find_section("Out-of-Plane Bend Parameters")
+    if iline < 0: return params
+    iline += 3
+    while iline < self.nlines:
+      words = self.lines[iline].split()
+      if len(words):
+        if words[0].startswith("###########"): break
+        if words[0] == "opbend":
+          class1 = int(words[1])
+          class2 = int(words[2])
+          class3 = int(words[3])
+          class4 = int(words[4])
+          value1 = float(words[5])
+          lmp1 = value1
+          params.append((class1,class2,class3,class4,lmp1))
       iline += 1
     return params
     
-  def find_section(self,txt):
+  # convert PRMfile params to LAMMPS angle_style amoeba UB params
-    txt = "##  %s  ##" % txt
+  # coeffs for K2,K3 = 0.0 since Urey-Bradley is simple harmonic
-    for iline,line in enumerate(self.lines):
+  
-      if txt in line: return iline
+  def ureybonds(self):
-    return -1
+    params = []
+    iline = self.find_section("Urey-Bradley Parameters")
+    if iline < 0: return params
+    iline += 3
+
+    while iline < self.nlines:
+      words = self.lines[iline].split()
+      if len(words):
+        if words[0].startswith("###########"): break
+        if words[0] == "ureybrad":
+          class1 = int(words[1])
+          class2 = int(words[2])
+          class3 = int(words[3])
+          value1 = float(words[4])
+          value2 = float(words[5])
+          lmp1 = value1
+          lmp2 = value2
+            
+          params.append((class1,class2,class3,lmp1,lmp2))
+      iline += 1
+    return params
+
+  # PiTorsion params, will be read from data file by fix pitorsion
+
+  def pitorsions(self):
+    params = []
+    iline = self.find_section("Pi-Torsion Parameters")
+    if iline < 0: return params
+    iline += 3
+
+    while iline < self.nlines:
+      words = self.lines[iline].split()
+      if len(words):
+        if words[0].startswith("###########"): break
+        if words[0] == "pitors":
+          class1 = int(words[1])
+          class2 = int(words[2])
+          value1 = float(words[3])
+          lmp1 = value1
+            
+          params.append((class1,class2,lmp1))
+      iline += 1
+    return params
+
+  # BiTorsion params, will be read from data file by fix bitorsion
+
+  def bitorsions(self):
+    params = []
+    iline = self.find_section("Torsion-Torsion Parameters")
+    if iline < 0: return params
+    iline += 3
+
+    while iline < self.nlines:
+      words = self.lines[iline].split()
+      if len(words):
+        if words[0].startswith("###########"): break
+        if words[0] == "tortors":
+          class1 = int(words[1])
+          class2 = int(words[2])
+          class3 = int(words[3])
+          class4 = int(words[4])
+          class5 = int(words[5])
+          nx = int(words[6])
+          ny = int(words[7])
+          iline += 1
+          array = []
+          for iy in range(ny):
+            xrow = []
+            for ix in range(nx):
+              words = self.lines[iline].split()
+              xgrid = float(words[0])
+              ygrid = float(words[1])
+              value = float(words[2])
+              tuple3 = (xgrid,ygrid,value)
+              xrow.append(tuple3)
+              iline += 1
+            array.append(xrow)
+          params.append((class1,class2,class3,class4,class5,nx,ny,array))
+      iline += 1
+    return params
 
 # ----------------------------------------
 # main program
@@ -328,6 +474,7 @@ amoeba = hippo = 0
 xyzfile = ""
 prmfile = ""
 datafile = ""
+bitorsionfile = ""
 pbcflag = 0
 
 iarg = 0
@@ -350,6 +497,10 @@ while iarg < narg:
     if iarg + 2 > narg: error()
     datafile = args[iarg+1]
     iarg += 2
+  elif args[iarg] == "-bitorsion":
+    if iarg + 2 > narg: error()
+    bitorsionfile = args[iarg+1]
+    iarg += 2
   elif args[iarg] == "-nopbc":
     pbcflag = 0
     iarg += 1
@@ -455,7 +606,7 @@ for i in id:
         stack.append(k)
 
 # ----------------------------------------
-# create lists of bonds, angles, dihedrals
+# create lists of bonds, angles, dihedrals, impropers
 # ----------------------------------------
 
 # create blist = list of bonds
@@ -493,13 +644,18 @@ for atom2 in id:
 # generate topology via triple loop over neighbors of dihedral atom2
 #   double loop over bonds of atom2
 #   additional loop over bonds of atom3
-# avoid double counting by requiring atom1 < atom3
+# avoid double counting the reverse dihedral by use of ddict dictionary
+
+# NOTE: could just avoid double count by "if atom1 < atom4" as in bond, angle ?
+#       gives different list, but is it still identical ?
+#       would have to check 2 data files, write comparison Py script
 
 id = xyz.id
 type = xyz.type
 bonds = xyz.bonds
 
 dlist = []
+ddict = {}
 
 for atom2 in id:
   for atom1 in bonds[atom2-1]:
@@ -507,11 +663,125 @@ for atom2 in id:
       if atom3 == atom1: continue
       for atom4 in bonds[atom3-1]:
         if atom4 == atom2 or atom4 == atom1: continue
-        if atom1 < atom3:
+        if (atom4,atom3,atom2,atom1) in ddict: continue
-          dlist.append((atom1,atom2,atom3,atom4))
+        dlist.append((atom1,atom2,atom3,atom4))
+        ddict[(atom1,atom2,atom3,atom4)] = 1
+
+# create olist = list of out-of-plane impropers
+# generate topology by triple loop over bonds of center atom2
+# atom2 must have 3 or more bonds to be part of an improper
+# avoid double counting by requiring atom3 < atom4
+# this is since in Tinker the final 2 atoms in the improper are interchangeable
+
+id = xyz.id
+type = xyz.type
+bonds = xyz.bonds
+
+olist = []
+
+for atom2 in id:
+  if len(bonds[atom2-1]) < 3: continue
+  for atom1 in bonds[atom2-1]:
+    for atom3 in bonds[atom2-1]:
+      for atom4 in bonds[atom2-1]:
+        if atom1 == atom3: continue
+        if atom1 == atom4: continue
+        if atom3 >= atom4: continue
+        olist.append((atom1,atom2,atom3,atom4))
 
 # ----------------------------------------
-# create lists of bond/angle/dihedral types
+# create list of Urey-Bradley triplet matches
+# ----------------------------------------
+
+# scan list of angles to find triplets that match UB parameters
+# if match, add it to UB bond list
+
+type = xyz.type
+classes = prm.classes
+
+ublist = []
+
+ubdict = {}
+for m,params in enumerate(prm.ureyparams):
+  ubdict[(params[0],params[1],params[2])] = (m,params)
+
+for atom1,atom2,atom3 in alist:
+  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:
+    ublist.append((atom1,atom2,atom3))
+  elif (c3,c2,c1) in ubdict:
+    ublist.append((atom3,atom2,atom1))
+
+# create pitorslist = list of 6-body interactions
+# based on central bond, each bond atom is bonded to exactly 2 other atoms
+# avoid double counting by requiring atom1 < atom2
+# NOTE: need more info on how to order the 6 atoms for Tinker to compute on
+
+type = xyz.type
+classes = prm.classes
+bonds = xyz.bonds
+
+pitorsionlist = []
+
+for atom1 in id:
+  for atom2 in bonds[atom1-1]:
+    if atom1 < atom2:
+      if len(bonds[atom1-1]) != 3: continue
+      if len(bonds[atom2-1]) != 3: continue
+      
+      if bonds[atom1-1][0] == atom2:
+        atom3 = bonds[atom1-1][1]
+        atom4 = bonds[atom1-1][2]
+      elif bonds[atom1-1][1] == atom2:
+        atom3 = bonds[atom1-1][0]
+        atom4 = bonds[atom1-1][2]
+      elif bonds[atom1-1][2] == atom2:
+        atom3 = bonds[atom1-1][0]
+        atom4 = bonds[atom1-1][1]
+
+      if bonds[atom2-1][0] == atom1:
+        atom5 = bonds[atom2-1][1]
+        atom6 = bonds[atom2-1][2]
+      elif bonds[atom2-1][1] == atom1:
+        atom5 = bonds[atom2-1][0]
+        atom6 = bonds[atom2-1][2]
+      elif bonds[atom2-1][2] == atom1:
+        atom5 = bonds[atom2-1][0]
+        atom6 = bonds[atom2-1][1]
+
+      pitorsionlist.append((atom3,atom4,atom1,atom2,atom5,atom6))
+
+# create bitorslist = list of 5-body interactions
+# generate topology via double loop over neighbors of central atom3
+#   additional double loop over bonds of atom2 and bonds of atom4
+# avoid double counting the reverse bitorsion by use of btdict dictionary
+
+type = xyz.type
+classes = prm.classes
+bonds = xyz.bonds
+
+bitorsionlist = []
+btdict = {}
+
+for atom3 in id:
+  for atom2 in bonds[atom3-1]:
+    for atom4 in bonds[atom3-1]:
+      if atom2 == atom4: continue
+      for atom1 in bonds[atom2-1]:
+        for atom5 in bonds[atom4-1]:
+          if atom1 == atom3 or atom5 == atom3 or atom1 == atom5: continue
+          if (atom5,atom4,atom3,atom2,atom1) in btdict: continue
+          bitorsionlist.append((atom1,atom2,atom3,atom4,atom5))
+          btdict[(atom1,atom2,atom3,atom4,atom5)] = 1
+
+# ----------------------------------------
+# create lists of bond/angle/dihedral/improper types
 # ----------------------------------------
 
 # generate btype = LAMMPS type of each bond
@@ -542,19 +812,16 @@ for atom1,atom2 in blist:
   
   if (c1,c2) in bdict: m,params = bdict[(c1,c2)]
   elif (c2,c1) in bdict: m,params = bdict[(c2,c1)]
-  else:
+  else: error("bond not found: %d %d: %d %d" % (atom1,atom2,c1,c2))
-    print "Bond not found",atom1,atom2,c1,c2
-    sys.exit()
     
   if not flags[m]:
     v1,v2,v3,v4 = params[2:]
     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
-# generate baparams = LAMMPS bond-angle params for each angle type
 # flags[i] = which LAMMPS angle type (1-N) the Tinker FF file angle I is
 #        0 = none
 # Tinker FF file angle entries can have 1, 2, or 3 options
@@ -575,15 +842,11 @@ for m,params in enumerate(prm.angleparams):
   noptions += n
   
 flags = noptions*[0]
-#baflags = len(baprm)*[0]
 atype = []
 aparams = []
-baparams = []
 
-#DEBUG
+for i,one in enumerate(alist):
-opcount = 0
+  atom1,atom2,atom3 = one
-
-for atom1,atom2,atom3 in alist:
   type1 = type[atom1-1]
   type2 = type[atom2-1]
   type3 = type[atom3-1]
@@ -593,8 +856,19 @@ for atom1,atom2,atom3 in alist:
 
   if (c1,c2,c3) in adict or (c3,c2,c1) in adict:
     if (c1,c2,c3) in adict: m,params = adict[(c1,c2,c3)]
-    if (c3,c2,c1) in adict: m,params = adict[(c3,c2,c1)]
+
+    # IMPORTANT subtlety
+    # flip order of 3 atoms in alist if the angle
+    #   matches Angle Bending section of PRM file in reverse order
+    # necessary b/c BondAngle coeffs will be generated with r1,r2 params
+    #   from Bond Stretching section of PRM file
+    # since in general r1 != r2, the LAMMPS AngleAmoeba class requires
+    #   the 3 atoms in the angle be in the order that matches r1 and r2
     
+    if (c3,c2,c1) in adict:
+      m,params = adict[(c3,c2,c1)]
+      alist[i] = (atom3,atom2,atom1)
+      
     # params is a sequence of 1 or 2 or 3 options
     # which = which of 1,2,3 options this atom triplet matches
     # for which = 2 or 3, increment m to index correct position in flags
@@ -609,9 +883,6 @@ for atom1,atom2,atom3 in alist:
     #     option 2 if one of 2 additional bonds is to an H atom
     #     option 3 if both of 2 additional bonds is to an H atom
 
-    # DEBUG
-    debugflag = 0
-    
     if len(params[3]) == 1:
       which = 1
       
@@ -624,7 +895,7 @@ for atom1,atom2,atom3 in alist:
         print "  angle atom classes:",c1,c2,c3
         print "  Tinker FF file param options:",len(params[3])
         print "  Nbonds and hydrogen count:",nbonds,hcount
-        #sys.exit()      // NOTE: allow this for now
+        #sys.exit()      NOTE: allow this for now
 
       if hcount == 0: which = 1
       elif hcount == 1:
@@ -638,10 +909,6 @@ for atom1,atom2,atom3 in alist:
       print "  Nbonds and hydrogen count:",nbonds,hcount
       print "  which:",which,m
 
-      # DEBUG
-      debugflag = 1
-      opcount += 1
-      
     elif len(params[3]) == 3:
       nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass)
       
@@ -651,7 +918,7 @@ for atom1,atom2,atom3 in alist:
         print "  angle atom classes:",c1,c2,c3
         print "  Tinker FF file param options:",len(params[3])
         print "  Nbonds and hydrogen count:",nbonds,hcount
-        #sys.exit()     // NOTE: allow this for now
+        #sys.exit()     NOTE: allow this for now
         
       if hcount == 0: which = 1
       elif hcount == 1:
@@ -661,56 +928,90 @@ for atom1,atom2,atom3 in alist:
         which = 3
         m += 2
 
-      print "4-bond angle"
-      print "  angle atom IDs:",atom1,atom2,atom3
-      print "  angle atom classes:",c1,c2,c3
-      print "  Tinker FF file param options:",len(params[3])
-      print "  Nbonds and hydrogen count:",nbonds,hcount
-      print "  which:",which,m
-
-      # DEBUG
-      debugflag = 1
-      opcount += 1
-
   else:
-    print "Angle not found",atom1,atom2,atom3,c1,c2,c3
+    error("angle not found: %d %d %d: %d %d %d" % (atom1,atom2,atom3,c1,c2,c3))
-    sys.exit()
     
   if not flags[m]:
-    pflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
+    pflag,ubflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
-    # DEBUG single line
+    aparams.append((pflag,ubflag,v1,v2,v3,v4,v5,v6))
-    if debugflag: pflag = 2
-    aparams.append((pflag,v1,v2,v3,v4,v5,v6))
     flags[m] = len(aparams)
   atype.append(flags[m])
 
-print "OPTION angles",opcount
+# augment the aparams with bond-angle cross terms from bondangleparams
+# generate baparams = LAMMPS bond-angle params for each angle type
+# badict = dictionary for angle tuples in bongangleparams
 
-  # NOTE: baparams may need to be flipped if match is 3,2,1 instead of 1,2,3
+badict = {}
+for v1,v2,v3,v4,v5,v6,v7 in prm.bondangleparams:
+  if (v1,v2,v3) in badict: continue
+  badict[(v1,v2,v3)] = (v4,v5,v6,v7)
+
+baparams = []
+
+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 (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
   
-  # NOTE: mismatch between angle and bondangle params may not be handled right
+  if (c1,c2,c3) in ubdict:
-  # should be a new LAMMPS type if either angle or bondangle params do not match?
+    r1,r2 = ubdict[(c1,c2,c3)]
-  #for m,params in enumerate(baprm):
+    pflag,ubflag,v1,v2,v3,v4,v5,v6 = aparams[itype]
-  #  c1,c2,c3,v1,v2,v3,v4 = params
+    ubflag = 1
-  #  if (c1 == class1 and c2 == class2 and c3 == class3) or \
+    aparams[itype] = (pflag,ubflag,v1,v2,v3,v4,v5,v6)
-  #     (c1 == class3 and c2 == class2 and c3 == class1):
+  elif (c3,c2,c1) in ubdict:
-  #    found += 1
+    r1,r2 = ubdict[(c3,c2,c1)]
-  #    if baflags[m]:
+    pflag,ubflag,v1,v2,v3,v4,v5,v6 = aparams[itype]
-  #      continue
+    ubflag = 1
-  #      #atype.append(baflags[m])
+    aparams[itype] = (pflag,ubflag,v1,v2,v3,v4,v5,v6)
-  #    else:
+  else:
-  #      baparams.append((v1,v2,v3,v4))
+    r1,r2 = 2*[0.0]
-  #      baflags[m] = len(baparams)
+    
-  #      #atype.append(baflags[m])
+  ubparams.append((r1,r2))
-  #    break
-  #  if found != 1: print "Not found",atom1,atom2,atom3,class1,class2,class3
 
 # 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
 #        0 = none
 # convert prm.torsionparams to a dictionary for efficient searching
-# key = (class1,class2)
+# key = (class1,class2,class3,class4)
 # value = (M,params) where M is index into prm.torsionparams
 
 id = xyz.id
@@ -738,15 +1039,159 @@ for atom1,atom2,atom3,atom4 in dlist:
   if (c1,c2,c3,c4) in ddict: m,params = ddict[(c1,c2,c3,c4)]
   elif (c4,c3,c2,c1) in ddict: m,params = ddict[(c4,c3,c2,c1)]
   else:
-    print "Dihedral not found",atom1,atom2,atom3,atom4,c1,c2,c3,c4
+    error("dihedral not found: %d %d %d %d: %d %d %d %d" % \
-    sys.exit()
+          (atom1,atom2,atom3,atom4,c1,c2,c3,c4))
 
   if not flags[m]:
-    v1,v2,v3,v4,v5,v6,v7,v8,v9 = params[4:]
+    oneparams = params[4:]
-    dparams.append((v1,v2,v3,v4,v5,v6,v7,v8,v9))
+    dparams.append(oneparams)
     flags[m] = len(dparams)
   dtype.append(flags[m])
 
+# generate otype = LAMMPS type of each out-of-plane improper
+# generate oparams = LAMMPS params for each improper type
+# flags[i] = which LAMMPS improper type (1-N) the Tinker FF file improper I is
+#        0 = none
+# convert prm.opbendparams to a dictionary for efficient searching
+# key = (class1,class2)
+# value = (M,params) where M is index into prm.opbendparams
+
+id = xyz.id
+type = xyz.type
+classes = prm.classes
+
+odict = {}
+for m,params in enumerate(prm.opbendparams):
+  odict[(params[0],params[1])] = (m,params)
+
+flags = len(prm.opbendparams)*[0]
+otype = []
+oparams = []
+olist_reduced = []
+
+for atom1,atom2,atom3,atom4 in olist:
+  type1 = type[atom1-1]
+  type2 = type[atom2-1]
+  type3 = type[atom3-1]
+  type4 = type[atom4-1]
+  c1 = classes[type1-1]
+  c2 = classes[type2-1]
+  c3 = classes[type3-1]
+  c4 = classes[type4-1]
+
+  # 4-tuple is only an improper if matches an entry in PRM file
+  # olist_reduced = list of just these 4-tuples
+  
+  if (c1,c2) in odict:
+    m,params = odict[(c1,c2)]
+    olist_reduced.append((atom1,atom2,atom3,atom4))
+    
+    if not flags[m]:
+      oneparams = params[4:]
+      oparams.append(oneparams)
+      flags[m] = len(oparams)
+    otype.append(flags[m])
+
+# replace original olist with reduced version
+
+olist = olist_reduced
+
+# generate pitorsiontype = LAMMPS type of each pitorsion
+# generate pitorsionparams = LAMMPS params for each pitorsion type
+# flags[i] = which LAMMPS pitorsion type (1-N) the Ith Tinker PRM file ptors is
+#        0 = none
+# convert prm.pitorsionparams to a dictionary for efficient searching
+# key = (class1,class2)
+# value = (M,params) where M is index into prm.pitorsionparams
+
+id = xyz.id
+type = xyz.type
+classes = prm.classes
+
+pitdict = {}
+for m,params in enumerate(prm.pitorsionparams):
+  pitdict[(params[0],params[1])] = (m,params)
+         
+flags = len(prm.pitorsionparams)*[0]
+pitorsiontype = []
+pitorsionparams = []
+pitorsionlist_reduced = []
+
+for tmp1,tmp2,atom1,atom2,tmp3,tmp4 in pitorsionlist:
+  type1 = type[atom1-1]
+  type2 = type[atom2-1]
+  c1 = classes[type1-1]
+  c2 = classes[type2-1]
+
+  # 6-tuple is only a PiTorsion if central 2 atoms match an entry in PRM file
+  # pitorsionlist_reduced = list of just these 6-tuples
+  
+  if (c1,c2) in pitdict or (c2,c1) in pitdict:
+    if (c1,c2) in pitdict: m,params = pitdict[(c1,c2)]
+    else: m,params = pitdict[(c2,c1)]
+    pitorsionlist_reduced.append((tmp1,tmp2,atom1,atom2,tmp3,tmp4))
+
+    if not flags[m]:
+      v1 = params[2:]
+      pitorsionparams.append(v1)
+      flags[m] = len(pitorsionparams)
+    pitorsiontype.append(flags[m])
+    
+# replace original pitorsionlist with reduced version
+
+pitorsionlist = pitorsionlist_reduced
+
+# generate bitorsiontype = LAMMPS type of each bitorsion
+# generate bitorsionparams = LAMMPS params for each bitorsion type
+# flags[i] = which LAMMPS bitorsion type (1-N) the Ith Tinker PRM file btors is
+#        0 = none
+# convert prm.bitorsionparams to a dictionary for efficient searching
+# key = (class1,class2,class3,class4,class5)
+# value = (M,params) where M is index into prm.bitorsionparams
+
+id = xyz.id
+type = xyz.type
+classes = prm.classes
+
+bitdict = {}
+for m,params in enumerate(prm.bitorsionparams):
+  bitdict[(params[0],params[1],params[2],params[3],params[4])] = (m,params)
+
+flags = len(prm.bitorsionparams)*[0]
+bitorsiontype = []
+bitorsionparams = []
+bitorsionlist_reduced = []
+
+for atom1,atom2,atom3,atom4,atom5 in bitorsionlist:
+  type1 = type[atom1-1]
+  type2 = type[atom2-1]
+  type3 = type[atom3-1]
+  type4 = type[atom4-1]
+  type5 = type[atom5-1]
+  c1 = classes[type1-1]
+  c2 = classes[type2-1]
+  c3 = classes[type3-1]
+  c4 = classes[type4-1]
+  c5 = classes[type5-1]
+
+  # 5-tuple is only a BiTorsion if 5 atoms match an entry in PRM file
+  # bitorsionlist_reduced = list of just these 5-tuples
+  
+  if (c1,c2,c3,c4,c5) in bitdict or (c5,c4,c3,c2,c1) in bitdict:
+    if (c1,c2,c3,c4,c5) in bitdict: m,params = bitdict[(c1,c2,c3,c4,c5)]
+    else: m,params = bitdict[(c5,c4,c3,c2,c1)]
+    bitorsionlist_reduced.append((atom1,atom2,atom3,atom4,atom5))
+
+    if not flags[m]:
+      v1 = params[5:]
+      bitorsionparams.append(v1)
+      flags[m] = len(bitorsionparams)
+    bitorsiontype.append(flags[m])
+    
+# replace original bitorsionlist with reduced version
+
+bitorsionlist = bitorsionlist_reduced
+
 # ----------------------------------------
 # assign each atom to a Tinker group
 # NOTE: doing this inside LAMMPS now
@@ -800,6 +1245,9 @@ ttype = xyz.type
 nbonds = len(blist)
 nangles = len(alist)
 ndihedrals = len(dlist)
+nimpropers = len(olist)
+npitorsions = len(pitorsionlist)
+nbitorsions = len(bitorsionlist)
 
 # data file header values
 
@@ -863,17 +1311,19 @@ if nangles:
     lines.append(line+'\n')
   d.sections["Angle Coeffs"] = lines
 
-  #lines = []
+  lines = []
-  #for i,one in enumerate(aparams):
+  for i,one in enumerate(baparams):
-  #  line = "%d %g %g %g" % (i+1,0.0,0.0,0.0)
+    strone = [str(single) for single in one]
-  #  lines.append(line+'\n')
+    line = "%d %s" % (i+1,' '.join(strone))
-  #d.sections["BondBond Coeffs"] = lines
+    lines.append(line+'\n')
+  d.sections["BondAngle Coeffs"] = lines
 
-  #lines = []
+  lines = []
-  #for i,one in enumerate(aparams):
+  for i,one in enumerate(ubparams):
-  #  line = "%d %g %g %g %g" % (i+1,0.0,0.0,0.0,0.0)
+    strone = [str(single) for single in one]
-  #  lines.append(line+'\n')
+    line = "%d %s" % (i+1,' '.join(strone))
-  # d.sections["BondAngle Coeffs"] = lines
+    lines.append(line+'\n')
+  d.sections["UreyBradley Coeffs"] = lines
 
   lines = [] 
   for i,one in enumerate(alist):
@@ -887,7 +1337,8 @@ if ndihedrals:
   
   lines = []
   for i,one in enumerate(dparams):
-    line = "%d %s" % (i+1,' '.join(one))
+    strone = [str(single) for single in one]
+    line = "%d %s" % (i+1,' '.join(strone))
     lines.append(line+'\n')
   d.sections["Dihedral Coeffs"] = lines
 
@@ -897,6 +1348,71 @@ if ndihedrals:
     lines.append(line+'\n')
   d.sections["Dihedrals"] = lines
 
+if nimpropers:
+  d.headers["impropers"] = len(olist)
+  d.headers["improper types"] = len(oparams)
+  
+  lines = []
+  for i,one in enumerate(oparams):
+    strone = [str(single) for single in one]
+    line = "%d %s" % (i+1,' '.join(strone))
+    lines.append(line+'\n')
+  d.sections["Improper Coeffs"] = lines
+
+  lines = [] 
+  for i,one in enumerate(olist):
+    line = "%d %d %d %d %d %d" % (i+1,otype[i],one[0],one[1],one[2],one[3])
+    lines.append(line+'\n')
+  d.sections["Impropers"] = lines
+
+if npitorsions:
+  d.headers["pitorsions"] = len(pitorsionlist)
+  d.headers["pitorsion types"] = len(pitorsionparams)
+  
+  lines = []
+  for i,one in enumerate(pitorsionparams):
+    strone = [str(single) for single in one]
+    line = "%d %s" % (i+1,' '.join(strone))
+    lines.append(line+'\n')
+  d.sections["PiTorsion Coeffs"] = lines
+
+  lines = [] 
+  for i,one in enumerate(pitorsionlist):
+    line = "%d %d %d %d %d %d %d %d" % \
+           (i+1,pitorsiontype[i],one[0],one[1],one[2],one[3],one[4],one[5])
+    lines.append(line+'\n')
+  d.sections["PiTorsions"] = lines
+
+if nbitorsions:
+  d.headers["bitorsions"] = len(bitorsionlist)
+
+  # if there are bitorsions, then -bitorsion file must have been specified
+  
+  if not bitorsionfile:
+    error("no -bitorsion file was specified, but %d bitorsions exist" % \
+          nbitorsions)
+
+  fp = open(bitorsionfile,'w')
+  print >>fp,"Tinker BiTorsion parameter file for fix bitorsion\n"
+  print >>fp,"%d bitorsion types" % len(bitorsionparams)
+  itype = 0
+  for nx,ny,array in bitorsionparams:
+    itype += 1
+    print >>fp
+    print >>fp,itype,nx,ny
+    for ix in range(nx):
+      for iy in range(ny):
+        xgrid,ygrid,value = array[ix][iy]
+        print >>fp," ",xgrid,ygrid,value
+  fp.close()
+  
+  lines = [] 
+  for i,one in enumerate(bitorsionlist):
+    line = "%d %d %d %d %d %d %d" % \
+           (i+1,bitorsiontype[i],one[0],one[1],one[2],one[3],one[4])
+    lines.append(line+'\n')
+  d.sections["BiTorsions"] = lines
+
 d.write(datafile)
 
 # print stats to screen
@@ -907,9 +1423,15 @@ print "Tinker XYZ types =",len(tink2lmp)
 print "Tinker PRM types =",prm.ntypes
 #print "Tinker groups =",ngroups
 print "Nmol =",nmol
-print "Nbonds =",len(blist)
+print "Nbonds =",nbonds
-print "Nangles =",len(alist)
+print "Nangles =",nangles
-print "Ndihedrals =",len(dlist)
+print "Ndihedrals =",ndihedrals
+print "Nimpropers =",nimpropers
+print "Npitorsions =",npitorsions
+print "Nbitorsions =",nbitorsions
 print "Nbondtypes =",len(bparams)
 print "Nangletypes =",len(aparams)
 print "Ndihedraltypes =",len(dparams)
+print "Nimpropertypes =",len(oparams)
+print "Npitorsiontypes =",len(pitorsionparams)
+print "Nbitorsiontypes =",len(bitorsionparams)