git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10521 f3b2605a-c512-4ea7-a41b-209d697bcdaa

tools/msi2lmp/README

@@ -1,6 +1,11 @@
 Stephanie Teich-McGoldrick (Sandai) is the current maintainer
 of the msi2lmp tool.  She can be contacted at steichm at sandia.gov
 
+02 Aug 2013 Axel Kohlmeyer <akohlmey@gmail.com>
+
+Added rudimentary support for OPLS-AA based on
+input provided by jeff greathouse.
+
 18 Jul 2013 Axel Kohlmeyer <akohlmey@gmail.com>
 
 Added support for writing out image flags

tools/msi2lmp/src/GetParameters.c

@@ -77,7 +77,7 @@ void GetParameters()
       printf(" Unable to find vdw data for %s\n",atomtypes[i].potential);
       condexit(11);
     } else {
-      if (forcefield & FF_TYPE_CLASS1) {
+      if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
         if((ff_vdw.data[k].ff_param[0] != 0.0 ) &&
            (ff_vdw.data[k].ff_param[1] != 0.0)) {
           atomtypes[i].params[0] =
@@ -132,7 +132,7 @@ void GetParameters()
              potential_types[0],potential_types[1]);
       condexit(12);
     } else {
-      if (forcefield & FF_TYPE_CLASS1) {
+      if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
         bondtypes[i].params[0] = ff_bond.data[k].ff_param[1];
         bondtypes[i].params[1] = ff_bond.data[k].ff_param[0];
       } 
@@ -189,7 +189,7 @@ void GetParameters()
              potential_types[0],potential_types[1],potential_types[2]);
       condexit(13);
     } else {
-      if (forcefield & FF_TYPE_CLASS1) {
+      if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
         angletypes[i].params[0] = ff_ang.data[k].ff_param[1];
         angletypes[i].params[1] = ff_ang.data[k].ff_param[0];
       }
@@ -341,6 +341,10 @@ void GetParameters()
         }
         dihedraltypes[i].params[2] = ff_tor.data[k].ff_param[1];
       }
+      if (forcefield & FF_TYPE_OPLSAA) {
+        for (j=0; j < 4; j++)
+          dihedraltypes[i].params[j] = ff_tor.data[k].ff_param[j];
+      }
       if (forcefield & FF_TYPE_CLASS2) {
         for (j=0; j < 6; j++)
           dihedraltypes[i].params[j] = ff_tor.data[k].ff_param[j];

tools/msi2lmp/src/InitializeItems.c

@@ -31,7 +31,7 @@ void InitializeItems(void)
   /* BOND */
 
   ff_bond.number_of_members = 2;
-  if (forcefield & FF_TYPE_CLASS1) {
+  if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
     strcpy(ff_bond.keyword,"#quadratic_bond");
     ff_bond.number_of_parameters = 2;
   }
@@ -52,7 +52,7 @@ void InitializeItems(void)
   /* ANGLE */
 
   ff_ang.number_of_members = 3;
-  if (forcefield & FF_TYPE_CLASS1) {
+  if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
     strcpy(ff_ang.keyword,"#quadratic_angle");
     ff_ang.number_of_parameters = 2;
   }
@@ -65,7 +65,7 @@ void InitializeItems(void)
   /* TORSION */
 
   ff_tor.number_of_members = 4;
-  if (forcefield & FF_TYPE_CLASS1) {
+  if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
     strcpy(ff_tor.keyword,"#torsion_1");
     ff_tor.number_of_parameters = 3;
   } 
@@ -78,7 +78,7 @@ void InitializeItems(void)
   /* OOP */
 
   ff_oop.number_of_members = 4;
-  if (forcefield & FF_TYPE_CLASS1) {
+  if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
     strcpy(ff_oop.keyword,"#out_of_plane");
     ff_oop.number_of_parameters = 3;
   }

tools/msi2lmp/src/ReadCarFile.c

@@ -188,7 +188,7 @@ void ReadCarFile(void)
       atoms[k].molecule = m;
       atoms[k].no = k;
 
-      fscanf(CarF,"%s %lf %lf %lf %*s %d %s %s %f",
+      fscanf(CarF,"%s %lf %lf %lf %*s %d %s %s %lf",
              atoms[k].name,
              &(atoms[k].x[0]),
              &(atoms[k].x[1]),

tools/msi2lmp/src/ReadFrcFile.c

@@ -57,8 +57,9 @@ void ReadFrcFile(void)
             ff_vdw.keyword,ff_vdw.entries);
     fprintf(stderr," Item %s has %d entries\n",
             ff_bond.keyword,ff_bond.entries);
-    fprintf(stderr," Item %s has %d entries\n",
+    if (forcefield & FF_TYPE_CLASS1)
-            ff_morse.keyword,ff_morse.entries);
+        fprintf(stderr," Item %s has %d entries\n",
+                ff_morse.keyword,ff_morse.entries);
     fprintf(stderr," Item %s has %d entries\n",
             ff_ang.keyword,ff_ang.entries);
     if (forcefield & FF_TYPE_CLASS2) {

tools/msi2lmp/src/SearchAndFill.c

@@ -41,7 +41,7 @@ void SearchAndFill(struct FrcFieldItem *item)
   while (got_it == 0) {
     status = fgets( line, MAX_LINE_LENGTH, FrcF );
     if (status == NULL) {
-      fprintf(stderr," Unable to find forcefield keyword %s\n",item->keyword);
+      fprintf(stderr," Unable to find keyword '%s'\n",item->keyword);
       fprintf(stderr," Check consistency of forcefield name and class \n");
       fprintf(stderr," Exiting....\n");
       exit(1);

tools/msi2lmp/src/WriteDataFile.c

@@ -94,6 +94,7 @@ void WriteDataFile(char *nameroot)
   if (no_bond_types > 0) {
     m = 0;
     if (forcefield & FF_TYPE_CLASS1) m = 2;
+    if (forcefield & FF_TYPE_OPLSAA) m = 2;
     if (forcefield & FF_TYPE_CLASS2) m = 4;
 
     fprintf(DatF,"Bond Coeffs\n\n");
@@ -109,6 +110,7 @@ void WriteDataFile(char *nameroot)
   if (no_angle_types > 0) {
     m = 0;
     if (forcefield & FF_TYPE_CLASS1) m = 2;
+    if (forcefield & FF_TYPE_OPLSAA) m = 2;
     if (forcefield & FF_TYPE_CLASS2) m = 4;
 
     fprintf(DatF,"Angle Coeffs\n\n");
@@ -134,6 +136,17 @@ void WriteDataFile(char *nameroot)
                 (int) dihedraltypes[i].params[1],
                 (int) dihedraltypes[i].params[2]);
 
+    } else if (forcefield & FF_TYPE_OPLSAA) {
+
+      fprintf(DatF,"Dihedral Coeffs\n\n");
+
+      for (i=0; i < no_dihedral_types; i++) {
+        fprintf(DatF, "%3i",i+1);
+        for ( j = 0; j < 4; j++)
+          fprintf(DatF, " %10.4f",dihedraltypes[i].params[j]);
+
+        fputs("\n",DatF);
+      }
     } else if (forcefield & FF_TYPE_CLASS2) { 
 
       fprintf(DatF,"Dihedral Coeffs\n\n");
@@ -160,6 +173,17 @@ void WriteDataFile(char *nameroot)
       }
       fputs("\n",DatF);
     }
+  } else if (forcefield & FF_TYPE_OPLSAA) {
+    if (no_oop_types > 0) {
+      /* opls improper coeffs are like cvff: type K0 d(=-1) n(=2) */
+      fprintf(DatF,"Improper Coeffs\n\n");
+      for (i=0; i < no_oop_types; i++) {
+        fprintf(DatF,"%5i %10.4f %3i %3i\n",i+1,
+                ooptypes[i].params[0], (int) ooptypes[i].params[1],
+                (int) ooptypes[i].params[2]);
+      }
+      fputs("\n",DatF);
+    }
   } else if (forcefield & FF_TYPE_CLASS2) {
     if ((no_oop_types + no_angleangle_types) > 0) {
       fprintf(DatF,"Improper Coeffs\n\n");

tools/msi2lmp/src/msi2lmp.c

@@ -1,6 +1,8 @@
 /*
 *
-*  msi2lmp.exe  V3.8
+*  msi2lmp.exe  V3.9
+*
+*   v3.9 AK  - Rudimentary support for OPLS-AA
 *
 *   v3.6 KLA - Changes to output to either lammps 2001 (F90 version) or to
 *              lammps 2005 (C++ version)
@@ -317,15 +319,23 @@ int main (int argc, char *argv[])
     printf(" Forcefield file name: %s\n",FrcFileName);
   }
 
-  if (((forcefield & FF_TYPE_CLASS1) && (strstr(FrcFileName,"cff") != NULL)) ||
+  n = 0;
-      ((forcefield & FF_TYPE_CLASS2) &&
+  if (forcefield & FF_TYPE_CLASS1) {
-       ! ((strstr(FrcFileName,"cvff") == NULL)
+    if (strstr(FrcFileName,"cvff") != NULL) ++n;
-          || (strstr(FrcFileName,"clayff") == NULL)
+    if (strstr(FrcFileName,"clayff") != NULL) ++n;
-          || (strstr(FrcFileName,"compass") == NULL))) ||
+  } else if (forcefield & FF_TYPE_OPLSAA) {
-      ((forcefield & FF_TYPE_OPLSAA) &&
+    if (strstr(FrcFileName,"oplsaa") != NULL) ++n;
-       ! (strstr(FrcFileName,"opls") == NULL))) {
+  } else if (forcefield & FF_TYPE_CLASS2) {
-    fprintf(stderr," WARNING - forcefield name and class appear to\n");
+    if (strstr(FrcFileName,"pcff") != NULL) ++n;
-    fprintf(stderr,"           be inconsistent - Errors may result\n\n");
+    if (strstr(FrcFileName,"cff91") != NULL) ++n;
+    if (strstr(FrcFileName,"compass") != NULL) ++n;
+  }
+
+  if (n == 0) {
+    if (iflag > 0) fputs(" WARNING",stderr);
+    else           fputs(" Error  ",stderr);
+    
+    fputs("- forcefield name and class appear to be inconsistent\n\n",stderr);
     if (iflag == 0) return 7;
   }
 

tools/msi2lmp/test/data-compare.pl

@@ -826,6 +826,38 @@ sub syscompare {
     die "Number of ",$i," types does not match: $a vs $b" unless ($a == $b);
   }
 
+  topocompare($d1,$d2,'bond',2);
+  topocompare($d1,$d2,'angle',3);
+  topocompare($d1,$d2,'dihedral',4);
+  topocompare($d1,$d2,'improper',4);
+
+  coeffcompare($d1,$d2,'pair','atom');
+  coeffcompare($d1,$d2,'bond','bond');
+  coeffcompare($d1,$d2,'angle','angle');
+  coeffcompare($d1,$d2,'bondbond','angle');
+  coeffcompare($d1,$d2,'bondangle','angle');
+  coeffcompare($d1,$d2,'dihedral','dihedral');
+  coeffcompare($d1,$d2,'angleangletorsion','dihedral');
+  coeffcompare($d1,$d2,'bondbond13','dihedral');
+  coeffcompare($d1,$d2,'endbondtorsion','dihedral');
+  coeffcompare($d1,$d2,'middlebondtorsion','dihedral');
+  coeffcompare($d1,$d2,'improper','improper');
+  coeffcompare($d1,$d2,'angleangle','improper');
+
+  for ($i=0; $i < $d1->{natomtypes}; ++$i) {
+    $j = $i+1;
+    if (exists $d1->{mass}[$i]) {
+      $a = $d1->{mass}[$i];
+    } else {
+      die "No mass for atom type $j in data file 1";
+    }
+    if (exists $d2->{mass}[$i]) {
+      $a = $d2->{mass}[$i];
+    } else {
+      die "No mass for atom type $j in data file 2";
+    }
+  }
+
   # check box information
   die "Inconsistent box shape" if ($d1->{triclinic} != $d2->{triclinic});
 
@@ -868,38 +900,6 @@ sub syscompare {
     }
   }
 
-  for ($i=0; $i < $d1->{natomtypes}; ++$i) {
-    $j = $i+1;
-    if (exists $d1->{mass}[$i]) {
-      $a = $d1->{mass}[$i];
-    } else {
-      die "No mass for atom type $j in data file 1";
-    }
-    if (exists $d2->{mass}[$i]) {
-      $a = $d2->{mass}[$i];
-    } else {
-      die "No mass for atom type $j in data file 2";
-    }
-  }
-
-  topocompare($d1,$d2,'bond',2);
-  topocompare($d1,$d2,'angle',3);
-  topocompare($d1,$d2,'dihedral',4);
-  topocompare($d1,$d2,'improper',4);
-
-  coeffcompare($d1,$d2,'pair','atom');
-  coeffcompare($d1,$d2,'bond','bond');
-  coeffcompare($d1,$d2,'angle','angle');
-  coeffcompare($d1,$d2,'bondbond','angle');
-  coeffcompare($d1,$d2,'bondangle','angle');
-  coeffcompare($d1,$d2,'dihedral','dihedral');
-  coeffcompare($d1,$d2,'angleangletorsion','dihedral');
-  coeffcompare($d1,$d2,'bondbond13','dihedral');
-  coeffcompare($d1,$d2,'endbondtorsion','dihedral');
-  coeffcompare($d1,$d2,'middlebondtorsion','dihedral');
-  coeffcompare($d1,$d2,'improper','improper');
-  coeffcompare($d1,$d2,'angleangle','improper');
-
 }
 
 ########################################################################

tools/msi2lmp/test/runtests.sh

@@ -5,6 +5,24 @@ MSI2LMP=../src/msi2lmp.exe
 LAMMPS=../../../src/lmp_serial
 CHECKDATA=./data-compare.pl
 
+if [ ! -x $MSI2LMP ]
+then
+   echo "No executable $MSI2LMP"
+   exit 1
+fi
+
+if [ ! -d $BIOSYM_LIBRARY ]
+then
+   echo "No directory $BIOSYM_LIBRARY"
+   exit 1
+fi
+
+if [ ! -x $LAMMPS ]
+then
+   echo "No executable $LAMMPS"
+   exit 1
+fi
+
 verbose=1
 counter=0
 errors=0
@@ -41,6 +59,22 @@ do \
     counter=$(expr $counter + 4)
 done
 
+# OPLS-AA tests 
+for m in ethane
+do \
+    before=$errors
+    ${MSI2LMP} ${m}-oplsaa -c 0 -p ${verbose} -f oplsaa	-n	\
+        || errors=$(expr $errors + 1)
+    ${LAMMPS} -log none -screen none -in in.${m}-oplsaa		\
+        || errors=$(expr $errors + 1)
+    ${CHECKDATA} ${m}-oplsaa.data reference/${m}-oplsaa.data	\
+        || errors=$(expr $errors + 1)
+    ${CHECKDATA} ${m}-oplsaa.data2 reference/${m}-oplsaa.data2	\
+        || errors=$(expr $errors + 1)
+    [ $before -eq $errors ] && rm ${m}-oplsaa.data ${m}-oplsaa.data2 log.${m}-oplsaa
+    counter=$(expr $counter + 4)
+done
+
 # Class2 tests with compass
 for m in hydrogen ethane benzene naphthalene
 do \