diff --git a/tools/msi2lmp/README b/tools/msi2lmp/README index 924ca7024f..e3477a359f 100644 --- a/tools/msi2lmp/README +++ b/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 + +Added rudimentary support for OPLS-AA based on +input provided by jeff greathouse. + 18 Jul 2013 Axel Kohlmeyer Added support for writing out image flags diff --git a/tools/msi2lmp/src/GetParameters.c b/tools/msi2lmp/src/GetParameters.c index f4b83b926c..3fc6057a6f 100644 --- a/tools/msi2lmp/src/GetParameters.c +++ b/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]; diff --git a/tools/msi2lmp/src/InitializeItems.c b/tools/msi2lmp/src/InitializeItems.c index 954589ad04..4df9fd0f10 100644 --- a/tools/msi2lmp/src/InitializeItems.c +++ b/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; } diff --git a/tools/msi2lmp/src/ReadCarFile.c b/tools/msi2lmp/src/ReadCarFile.c index f07c9871b3..9a4f0e996b 100644 --- a/tools/msi2lmp/src/ReadCarFile.c +++ b/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]), diff --git a/tools/msi2lmp/src/ReadFrcFile.c b/tools/msi2lmp/src/ReadFrcFile.c index a64798add3..ce22288450 100644 --- a/tools/msi2lmp/src/ReadFrcFile.c +++ b/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", - ff_morse.keyword,ff_morse.entries); + if (forcefield & FF_TYPE_CLASS1) + 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) { diff --git a/tools/msi2lmp/src/SearchAndFill.c b/tools/msi2lmp/src/SearchAndFill.c index f3b3a37cbb..c9ea457d8e 100644 --- a/tools/msi2lmp/src/SearchAndFill.c +++ b/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); diff --git a/tools/msi2lmp/src/WriteDataFile.c b/tools/msi2lmp/src/WriteDataFile.c index ed64ec7d38..dbed331889 100644 --- a/tools/msi2lmp/src/WriteDataFile.c +++ b/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"); diff --git a/tools/msi2lmp/src/msi2lmp.c b/tools/msi2lmp/src/msi2lmp.c index f1f9af6a40..c53109faa3 100644 --- a/tools/msi2lmp/src/msi2lmp.c +++ b/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)) || - ((forcefield & FF_TYPE_CLASS2) && - ! ((strstr(FrcFileName,"cvff") == NULL) - || (strstr(FrcFileName,"clayff") == NULL) - || (strstr(FrcFileName,"compass") == NULL))) || - ((forcefield & FF_TYPE_OPLSAA) && - ! (strstr(FrcFileName,"opls") == NULL))) { - fprintf(stderr," WARNING - forcefield name and class appear to\n"); - fprintf(stderr," be inconsistent - Errors may result\n\n"); + n = 0; + if (forcefield & FF_TYPE_CLASS1) { + if (strstr(FrcFileName,"cvff") != NULL) ++n; + if (strstr(FrcFileName,"clayff") != NULL) ++n; + } else if (forcefield & FF_TYPE_OPLSAA) { + if (strstr(FrcFileName,"oplsaa") != NULL) ++n; + } else if (forcefield & FF_TYPE_CLASS2) { + if (strstr(FrcFileName,"pcff") != NULL) ++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; } diff --git a/tools/msi2lmp/test/data-compare.pl b/tools/msi2lmp/test/data-compare.pl index 63e40630fb..28d64b28f7 100755 --- a/tools/msi2lmp/test/data-compare.pl +++ b/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'); - } ######################################################################## diff --git a/tools/msi2lmp/test/runtests.sh b/tools/msi2lmp/test/runtests.sh index 6c0dd097e7..309de85064 100755 --- a/tools/msi2lmp/test/runtests.sh +++ b/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 \