reading of bitorsion data file

examples/amoeba/data.ubiquitin

@@ -14,6 +14,7 @@ LAMMPS data file created from Tinker ubiquitin.xyz and amoeba_ubiquitin.prm file
 0 41.91 ylo yhi
 0 41.91 zlo zhi
 6 pitorsion types
+106 bitorsions
 106 pitorsions
 
 Masses
@@ -26348,6 +26349,80 @@ UreyBradley Coeffs
 110 0.0 0.0
 111 -7.6 1.5537
 
+BiTorsions
+
+1 1 3 20 21 22 37
+2 1 22 37 38 39 56
+3 1 39 56 57 58 76
+4 1 58 76 77 78 92
+5 1 78 92 93 94 114
+6 1 94 114 115 116 128
+7 1 116 128 129 130 147
+8 1 130 147 148 149 161
+9 2 149 161 162 163 168
+10 1 163 168 169 170 190
+11 1 170 190 191 192 204
+12 1 192 204 205 206 223
+13 1 206 223 224 225 237
+14 1 225 237 238 239 256
+15 1 239 256 257 258 271
+16 1 258 271 272 273 287
+17 3 273 287 288 289 302
+18 1 304 316 317 318 327
+19 1 318 327 328 329 339
+20 1 329 339 340 341 353
+21 1 341 353 354 355 372
+22 1 355 372 373 374 387
+23 1 374 387 388 389 401
+24 1 389 401 402 403 417
+25 1 403 417 418 419 439
+26 1 419 439 440 441 449
+27 1 441 449 450 451 471
+28 1 451 471 472 473 490
+29 1 473 490 491 492 507
+30 1 492 507 508 509 519
+31 1 509 519 520 521 541
+32 1 521 541 542 543 556
+33 2 543 556 557 558 563
+34 3 558 563 564 565 582
+35 1 598 610 611 612 622
+36 1 612 622 623 624 639
+37 1 624 639 640 641 656
+38 1 641 656 657 658 680
+39 1 658 680 681 682 699
+40 1 682 699 700 701 718
+41 1 701 718 719 720 738
+42 1 720 738 739 740 748
+43 2 740 748 749 750 755
+44 1 750 755 756 757 777
+45 1 757 777 778 779 794
+46 1 779 794 795 796 813
+47 1 796 813 814 815 828
+48 1 815 828 829 830 840
+49 2 830 840 841 842 847
+50 1 842 847 848 849 871
+51 1 849 871 872 873 885
+52 1 873 885 886 887 904
+53 1 887 904 905 906 915
+54 1 906 915 916 917 927
+55 1 917 927 928 929 948
+56 1 929 948 949 950 962
+57 1 950 962 963 964 981
+58 1 964 981 982 983 998
+59 1 983 998 999 1000 1020
+60 1 1000 1020 1021 1022 1035
+61 1 1022 1035 1036 1037 1046
+62 1 1037 1046 1047 1048 1060
+63 1 1048 1060 1061 1062 1079
+64 1 1062 1079 1080 1081 1097
+65 1 1081 1097 1098 1099 1116
+66 1 1099 1116 1117 1118 1132
+67 1 1118 1132 1133 1134 1151
+68 1 1134 1151 1152 1153 1175
+69 1 1153 1175 1176 1177 1194
+70 1 1177 1194 1195 1196 1218
+71 2 1196 1218 1219 1220 1225
+
 PiTorsions
 
 1 1 2 4 3 20 21 24

examples/amoeba/in.ubiquitin

@@ -16,6 +16,8 @@ improper_style  amoeba
 fix             amtype all property/atom i_amtype ghost yes
 fix             pitorsion all pitorsion
 fix_modify      pitorsion energy yes
+fix             bitorsion all bitorsion bitorsion.ubiquitin.data
+fix_modify      bitorsion energy yes
 
 fix             extra all property/atom &
                 i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes
@@ -26,7 +28,7 @@ read_data       data.ubiquitin fix amtype NULL "Tinker Types" &
                 fix pitorsion "pitorsion types" "PiTorsion Coeffs" &
                 fix pitorsion pitorsions PiTorsions
 
-pair_style      amoeba
+pair_style      amoeba include bitorsion
 pair_coeff      * * amoeba_ubiquitin.prm amoeba_ubiquitin.key
 
 special_bonds   lj/coul 0.5 0.5 0.5 one/five yes

src/AMOEBA/fix_bitorsion.cpp

@@ -35,6 +35,7 @@ using namespace MathConst;
 #define BITORSIONMAX 6   // max # of BiTorsion terms stored by one atom
 #define LISTDELTA 10000
 #define LB_FACTOR 1.5
+#define MAXLINE 1024
 
 // NOTE: extra until figure things out
 
@@ -106,7 +107,7 @@ FixBiTorsion::FixBiTorsion(LAMMPS *lmp, int narg, char **arg) :
 
   // read and setup BiTorsion grid data
 
-  read_grid_map(arg[3]);
+  read_grid_data(arg[3]);
 
   // perform initial allocation of atom-based arrays
 
@@ -130,6 +131,12 @@ FixBiTorsion::FixBiTorsion(LAMMPS *lmp, int narg, char **arg) :
   nbitorsion_list = 0;
   max_bitorsion_list = 0;
   bitorsion_list = nullptr;
+
+  // BiTorsion grid data
+
+  ntypes = 0;
+  nxgrid,nygrid = nullptr;
+  btgrid = nullptr;
 }
 
 /* --------------------------------------------------------------------- */
@@ -154,6 +161,14 @@ FixBiTorsion::~FixBiTorsion()
   // local list of bitorsions to compute
 
   memory->destroy(bitorsion_list);
+
+  // BiTorsion grid data
+
+  delete [] nxgrid;
+  delete [] nygrid;
+  for (int itype = 0; itype < ntypes; itype++)
+    memory->destroy(btgrid[itype]);
+  delete [] btgrid;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -679,47 +694,92 @@ double FixBiTorsion::compute_scalar()
 
 // ----------------------------------------------------------------------
 // ----------------------------------------------------------------------
-// methods to read BiTorsion potential file, perform interpolation
+// methods to read BiTorsion grid file, perform interpolation
 // ----------------------------------------------------------------------
 // ----------------------------------------------------------------------
 
-void FixBiTorsion::read_grid_map(char *bitorsion_file)
+void FixBiTorsion::read_grid_data(char *bitorsion_file)
 {
-  int MAXLINE = 1024;
+  char line[MAXLINE];
-  char linebuf[MAXLINE];
+  char *eof;
-  char *chunk,*line;
-  int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter;
 
   FILE *fp = nullptr;
   if (me == 0) {
     fp = utils::open_potential(bitorsion_file,lmp,nullptr);
     if (fp == nullptr)
-      error->one(FLERR,"Cannot open fix cmap file {}: {}",
+      error->one(FLERR,"Cannot open fix bitorsion file {}: {}",
                  bitorsion_file, utils::getsyserror());
 
+    eof = fgets(line,MAXLINE,fp);
+    eof = fgets(line,MAXLINE,fp);
+    eof = fgets(line,MAXLINE,fp);
+    if (eof == nullptr) error->one(FLERR,"Unexpected end of fix bitorsion file");
+
+    sscanf(line,"%d",&ntypes);
   }
 
-  int done;
+  MPI_Bcast(&ntypes,1,MPI_INT,0,world);
+  if (ntypes == 0) error->all(FLERR,"Fix bitorsion file has no types");
 
-  while (!done) {
+  btgrid = new double***[ntypes];
-    // only read on rank 0 and broadcast to all other ranks
+  nxgrid = new int[ntypes];
-    if (me == 0)
+  nygrid = new int[ntypes];
-      done = (fgets(linebuf,MAXLINE,fp) == nullptr);
 
-    MPI_Bcast(&done,1,MPI_INT,0,world);
+  // read one array for each BiTorsion type from file
-    if (done) continue;
 
-    MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world);
+  int tmp,nx,ny;
+  double xgrid,ygrid,value;
 
-    // remove leading whitespace
+  for (int itype = 0; itype < ntypes; itype++) {
-    line = linebuf;
+    if (me == 0) {
-    while (line && (*line == ' ' || *line == '\t' || *line == '\r')) ++line;
+      eof = fgets(line,MAXLINE,fp);
+      eof = fgets(line,MAXLINE,fp);
+      if (eof == nullptr) 
+        error->one(FLERR,"Unexpected end of fix bitorsion file");
+      sscanf(line,"%d %d %d",&tmp,&nx,&ny);
+    }
 
-    // skip if empty line or comment
+    MPI_Bcast(&nx,1,MPI_INT,0,world);
-    if (!line || *line =='\n' || *line == '\0' || *line == '#') continue;
+    MPI_Bcast(&ny,1,MPI_INT,0,world);
+    nxgrid[itype] = nx;
+    nygrid[itype] = ny;
+
+    memory->create(btgrid[itype],nx,ny,3,"bitorsion:btgrid");
+
+    // NOTE: should read this chunk of lines with utils in single read
+
+    if (me == 0) {
+      for (int iy = 0; iy < ny; iy++) {
+        for (int ix = 0; ix < nx; ix++) {
+          eof = fgets(line,MAXLINE,fp);
+          if (eof == nullptr) 
+            error->one(FLERR,"Unexpected end of fix bitorsion file");
+          sscanf(line,"%lg %lg %lg",&xgrid,&ygrid,&value);
+          btgrid[itype][ix][iy][0] = xgrid;
+          btgrid[itype][ix][iy][1] = ygrid;
+          btgrid[itype][ix][iy][2] = value;
+        }
+      }
+    }
+
+    MPI_Bcast(&btgrid[itype][0][0][0],nx*ny*3,MPI_DOUBLE,0,world);
   }
 
   if (me == 0) fclose(fp);
+
+  // DEBUG
+
+  for (int i = 0; i < ntypes; i++) {
+    printf("ITYPE %d NXY %d %d\n",i+1,nxgrid[i],nygrid[i]);
+    for (int iy = 0; iy < ny; iy++) {
+      for (int ix = 0; ix < nx; ix++) {
+        printf("  IXY %d %d, values %g %g %g\n",ix+1,iy+1,
+               btgrid[i][ix][iy][0],
+               btgrid[i][ix][iy][1],
+               btgrid[i][ix][iy][2]);
+      }
+    }
+  }
 }
 
 // ----------------------------------------------------------------------

src/AMOEBA/fix_bitorsion.h

@@ -90,9 +90,15 @@ class FixBiTorsion : public Fix {
   int max_bitorsion_list;
   int **bitorsion_list;
 
+  // BiTorsion grid data
+
+  int ntypes;
+  int *nxgrid,*nygrid;
+  double ****btgrid;
+
   // read BiTorsion grid data
 
-  void read_grid_map(char *);
+  void read_grid_data(char *);
 };
 }    // namespace LAMMPS_NS
 

src/modify.cpp

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