initial version of fix bitorsion

src/AMOEBA/fix_bitorsion.h

@@ -0,0 +1,136 @@
+/* -*- 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(bitorsion,FixBiTorsion);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_BITORSION_H
+#define LMP_FIX_BITORSION_H
+
+#include "fix.h"
+namespace LAMMPS_NS {
+
+class FixBiTorsion : public Fix {
+ public:
+  FixBiTorsion(class LAMMPS *, int, char **);
+  ~FixBiTorsion();
+  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;
+
+  // read BiTorsion grid data
+
+  void read_grid_map(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_pitorsion.cpp

@@ -32,10 +32,11 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-#define PITORSIONMAX 6
+#define PITORSIONMAX 6   // max # of PiTorsion terms stored by one atom
 #define LISTDELTA 8196
 #define LB_FACTOR 1.5
 
+
 /* ---------------------------------------------------------------------- */
 
 FixPiTorsion::FixPiTorsion(LAMMPS *lmp, int narg, char **arg) :
@@ -116,7 +117,7 @@ FixPiTorsion::~FixPiTorsion()
   memory->destroy(pitorsion_atom5);
   memory->destroy(pitorsion_atom6);
 
-  // compute list
+  // local list of bitorsions to compute
 
   memory->destroy(pitorsion_list);
 
@@ -249,13 +250,6 @@ void FixPiTorsion::pre_neighbor()
       atom5 = domain->closest_image(i,atom5);
       atom6 = domain->closest_image(i,atom6);
 
-      /*
-      if (atom1 >= nlocal || atom2 >= nlocal || atom3 >= nlocal ||
-          atom4 >= nlocal || atom5 >= nlocal || atom6 >= nlocal)
-        printf("ATOM 123456: %d %d %d %d %d %d\n",
-               atom1,atom2,atom3,atom4,atom5,atom6);
-      */
-
       if (i <= atom1 && i <= atom2 && i <= atom3 &&
           i <= atom4 && i <= atom5 && i <= atom6) {
         if (npitorsion_list == max_pitorsion_list) {
@@ -437,7 +431,7 @@ void FixPiTorsion::post_force(int vflag)
     dphi2 = 2.0 * (cosine2*s2 - sine2*c2);
 
     // calculate pi-system torsion energy and master chain rule term
-    // NOTE: remove ptornunit if 1.0 ?
+    // NOTE: remove ptorunit if 1.0 ?
 
     double ptorunit = 1.0;
     e = ptorunit * v2 * phi2;
@@ -898,10 +892,12 @@ void FixPiTorsion::write_data_section(int mth, FILE *fp,
   if (mth == 0) {
     for (int i = 0; i < n; i++)
       fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
-              " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT "\n",
+              " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT 
+              " " TAGINT_FORMAT "\n",
               index+i,(int) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
               (tagint) ubuf(buf[i][2]).i,(tagint) ubuf(buf[i][3]).i,
-              (tagint) ubuf(buf[i][4]).i,(tagint) ubuf(buf[i][5]).i);
+              (tagint) ubuf(buf[i][4]).i,(tagint) ubuf(buf[i][5]).i,
+              (tagint) ubuf(buf[i][6]).i);
 
   // PiTorsion Coeffs section
 

src/AMOEBA/fix_pitorsion.h

@@ -66,7 +66,7 @@ class FixPiTorsion : public Fix {
 
  private:
   int nprocs, me;
-  int newton_bond, eflag_caller;
+  int eflag_caller;
   int ilevel_respa;
   int disable;
   bigint npitorsions;

src/AMOEBA/pair_amoeba.cpp

@@ -533,7 +533,7 @@ void PairAmoeba::settings(int narg, char **arg)
   polar_rspace_flag = polar_kspace_flag = 1;
   mpole_rspace_flag = mpole_kspace_flag = 1;
   bond_flag = angle_flag = dihedral_flag = improper_flag = 1;
-  urey_flag = pitorsion_flag = xtorsion_flag = 1;
+  urey_flag = pitorsion_flag = bitorsion_flag = 1;
 
   int newvalue = -1;
 
@@ -546,7 +546,7 @@ void PairAmoeba::settings(int narg, char **arg)
     polar_rspace_flag = polar_kspace_flag = 0;
     mpole_rspace_flag = mpole_kspace_flag = 0;
     bond_flag = angle_flag = dihedral_flag = improper_flag = 0;
-    urey_flag = pitorsion_flag = xtorsion_flag = 0;
+    urey_flag = pitorsion_flag = bitorsion_flag = 0;
 
   // exclude only specified FF components
 
@@ -583,7 +583,7 @@ void PairAmoeba::settings(int narg, char **arg)
     else if (strcmp(arg[iarg],"improper") == 0) improper_flag = newvalue;
     else if (strcmp(arg[iarg],"urey") == 0) urey_flag = newvalue;
     else if (strcmp(arg[iarg],"pitorsion") == 0) pitorsion_flag = newvalue;
-    else if (strcmp(arg[iarg],"xtorsion") == 0) xtorsion_flag = newvalue;
+    else if (strcmp(arg[iarg],"bitorsion") == 0) bitorsion_flag = newvalue;
     else error->all(FLERR,"Illegal pair_style command");
   }
 }
@@ -2090,6 +2090,7 @@ void *PairAmoeba::extract(const char *str, int &dim)
   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;
 
   if (strcmp(str,"opbend_cubic") == 0) return (void *) &opbend_cubic;
   if (strcmp(str,"opbend_quartic") == 0) return (void *) &opbend_quartic;

src/AMOEBA/pair_amoeba.h

@@ -80,7 +80,7 @@ class PairAmoeba : public Pair {
   int polar_rspace_flag,polar_kspace_flag;
   int mpole_rspace_flag,mpole_kspace_flag;
   int bond_flag,angle_flag,dihedral_flag,improper_flag;
-  int urey_flag,pitorsion_flag,xtorsion_flag;
+  int urey_flag,pitorsion_flag,bitorsion_flag;
 
   // DEBUG timers