From 239916df00c051dbae36ee619b47dacfbe28b2a8 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Sun, 4 Dec 2016 15:12:57 +0000 Subject: [PATCH] Fixed - FIELD file import was not adding terms correctly at all. Imported FIELD files now have energy units converted correctly. --- src/plugins/io_dlpoly/field_funcs.cpp | 184 +++++++++++++------------- 1 file changed, 93 insertions(+), 91 deletions(-) diff --git a/src/plugins/io_dlpoly/field_funcs.cpp b/src/plugins/io_dlpoly/field_funcs.cpp index fcb2726b4..4390d8540 100644 --- a/src/plugins/io_dlpoly/field_funcs.cpp +++ b/src/plugins/io_dlpoly/field_funcs.cpp @@ -116,43 +116,48 @@ bool DLPExpressionPlugin::importData() // First line is header information - use as FF name QString line; - Forcefield* newFF; + Forcefield* newFF; if (!fileParser_.readLine(line)) { - newFF = createForcefield(QString("unnamed field")); - }else{ - newFF = createForcefield(line); - } - // Create a new forcefield + newFF = createForcefield(QString("unnamed field")); + }else{ + newFF = createForcefield(line); + } // Next meaningful line is energy units if (!fileParser_.parseLine(Parser::StripComments+Parser::SkipBlanks) ) return false; - QString word = fileParser_.argc(0).toLower(); - if (word == "units"){ - QString unitString = fileParser_.argc(1).toLower(); - Prefs::EnergyUnit unit = Prefs::energyUnit(unitString); - if (unit == Prefs::nEnergyUnits) - { - Messenger::error("FIELD file energy unit (%s) is not compatible with Aten.", qPrintable(unitString)); - return false; - } - } - - int nMols; + QString word = fileParser_.argc(0).toLower(); + if (word == "units") + { + QString unitString = fileParser_.argc(1).toLower(); + Prefs::EnergyUnit unit = Prefs::energyUnit(unitString); + if (unit == Prefs::nEnergyUnits) + { + Messenger::error("FIELD file energy unit (%s) is not compatible with Aten.", qPrintable(unitString)); + return false; + } + newFF->setEnergyUnit(unit); + } + // Next is number of molecule types described or multipolar data section in the file if (!fileParser_.parseLine(Parser::StripComments+Parser::SkipBlanks)) return false; - word = fileParser_.argc(0).toLower(); - - if ((word == "mult")||(word=="multipolar")) { - if (!fileParser_.parseLine(Parser::StripComments+Parser::SkipBlanks)) return false; - //more stuff may need to be read from here but last line needs to be a reading - } - if (fileParser_.argc(0).toLower() == "molecules"){ - nMols = fileParser_.argi(1); - Messenger::print("Number of molecule types specified in FIELD file : %i", nMols); - } else { - Messenger::error("Didn't find 'molecules' directive where expected."); - return false; + word = fileParser_.argc(0).toLower(); + + if ((word == "mult")||(word=="multipolar")) + { + if (!fileParser_.parseLine(Parser::StripComments+Parser::SkipBlanks)) return false; + //more stuff may need to be read from here but last line needs to be a reading + } + + int nMols; + if (fileParser_.argc(0).toLower() == "molecules") + { + nMols = fileParser_.argi(1); + Messenger::print("Number of molecule types specified in FIELD file : %i", nMols); + } else { + Messenger::error("Didn't find 'molecules' directive where expected."); + return false; } + // Loop over molecule types for (int mol=0; molfindBond(atomNames.at(ii), atomNames.at(jj))) - { - Messenger::print("Bond %s-%s has already been created in the forcefield. Skipped...", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj))); - continue; - } + if (newFF->findBond(atomNames.at(ii), atomNames.at(jj))) continue; // Create new definition - ffb = newFF->addBond(bf, atomNames.at(ii), atomNames.at(jj)); - if (bf == BondFunctions::Harmonic) + if (formString == "harm") { + ffb = newFF->addBond(BondFunctions::Harmonic, atomNames.at(ii), atomNames.at(jj)); ffb->setParameter(BondFunctions::HarmonicK, fileParser_.argd(3)); ffb->setParameter(BondFunctions::HarmonicEq, fileParser_.argd(4)); } - else if (bf == BondFunctions::Morse) + else if (formString == "morse") { + ffb = newFF->addBond(BondFunctions::Morse, atomNames.at(ii), atomNames.at(jj)); ffb->setParameter(BondFunctions::MorseD, fileParser_.argd(3)); ffb->setParameter(BondFunctions::MorseK, fileParser_.argd(5)); ffb->setParameter(BondFunctions::MorseEq, fileParser_.argd(4)); } + else + { + Messenger::print("Functional form of bond term (%s) is not present in Aten.", qPrintable(formString)); + continue; + } + + Messenger::print("Created bond term %s-%s", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj))); } } else if (keyword == "constraints") { Messenger::print("Found 'constraints' block..."); - for (n=0; nfindBond(atomNames.at(ii), atomNames.at(jj))) - { - Messenger::print("Constraint bond %s-%s has already been created in the forcefield. Skipped...", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj))); - continue; - } + if (newFF->findBond(atomNames.at(ii), atomNames.at(jj))) continue; // Create new definition newFF->addBond(BondFunctions::Constraint, atomNames.at(ii), atomNames.at(jj)); ffb->setParameter(BondFunctions::HarmonicK, 1000.0); ffb->setParameter(BondFunctions::HarmonicEq, fileParser_.argd(2)); + + Messenger::print("Created bond constraint term %s-%s", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj))); } } else if (keyword == "angles") { Messenger::print("Found 'angles' block..."); - for (n=0; nfindAngle(atomNames.at(ii), atomNames.at(jj), atomNames.at(kk))) - { - Messenger::print("Angle %s-%s-%s has already been created in the forcefield. Skipped...", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj)), qPrintable(atomNames.at(kk))); - continue; - } + if (newFF->findAngle(atomNames.at(ii), atomNames.at(jj), atomNames.at(kk))) continue; // Create new definition - ffb = newFF->addAngle(af, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk)); - if (af == AngleFunctions::Harmonic) + if (formString == "harm") { + ffb = newFF->addAngle(AngleFunctions::Harmonic, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk)); ffb->setParameter(AngleFunctions::HarmonicK, fileParser_.argd(4)); ffb->setParameter(AngleFunctions::HarmonicEq, fileParser_.argd(5)); } - else if (af == AngleFunctions::Cosine) + else if (formString == "cos") { + ffb = newFF->addAngle(AngleFunctions::Cosine, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk)); ffb->setParameter(AngleFunctions::CosineK, fileParser_.argd(4)); ffb->setParameter(AngleFunctions::CosineN, fileParser_.argd(6)); ffb->setParameter(AngleFunctions::CosineEq, fileParser_.argd(5)); } - else if (af == AngleFunctions::HarmonicCosine) + else if (formString == "hcos") { + ffb = newFF->addAngle(AngleFunctions::HarmonicCosine, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk)); ffb->setParameter(AngleFunctions::HarmonicCosineK, fileParser_.argd(4)); ffb->setParameter(AngleFunctions::HarmonicCosineEq, fileParser_.argd(5)); } + else + { + Messenger::print("Functional form of angle term (%s) is not present in Aten.", qPrintable(formString)); + continue; + } + + Messenger::print("Created angle term %s-%s-%s", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj)), qPrintable(atomNames.at(kk))); } } else if (keyword == "dihedrals") { Messenger::print("Found 'dihedrals' block..."); - for (n=0; nfindTorsion(atomNames.at(ii), atomNames.at(jj), atomNames.at(kk), atomNames.at(ll))) - { - Messenger::print("Torsion %s-%s-%s-%s has already been created in the forcefield. Skipped...", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj)), qPrintable(atomNames.at(kk)), qPrintable(atomNames.at(ll))); - continue; - } + if (newFF->findTorsion(atomNames.at(ii), atomNames.at(jj), atomNames.at(kk), atomNames.at(ll))) continue; // Create new definition - ffb = newFF->addTorsion(tf, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk), atomNames.at(ll)); - if (tf == TorsionFunctions::Cosine) + if (formString == "cos") { + ffb = newFF->addTorsion(TorsionFunctions::Cosine, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk), atomNames.at(ll)); ffb->setParameter(TorsionFunctions::CosineK, fileParser_.argd(5)); ffb->setParameter(TorsionFunctions::CosineN, fileParser_.argd(7)); ffb->setParameter(TorsionFunctions::CosineEq, fileParser_.argd(6)); + ffb->setElecScale(fileParser_.argd(8)); + ffb->setVdwScale(fileParser_.argd(9)); } - else if (tf == TorsionFunctions::Cos3) + else if (formString == "cos3") { + ffb = newFF->addTorsion(TorsionFunctions::Cos3, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk), atomNames.at(ll)); ffb->setParameter(TorsionFunctions::Cos3K1, fileParser_.argd(5)); ffb->setParameter(TorsionFunctions::Cos3K2, fileParser_.argd(6)); ffb->setParameter(TorsionFunctions::Cos3K3, fileParser_.argd(7)); + ffb->setElecScale(fileParser_.argd(8)); + ffb->setVdwScale(fileParser_.argd(9)); } - else if (tf == TorsionFunctions::Cos3C) + else if (formString == "opls") { + ffb = newFF->addTorsion(TorsionFunctions::Cos3C, atomNames.at(ii), atomNames.at(jj), atomNames.at(kk), atomNames.at(ll)); ffb->setParameter(TorsionFunctions::Cos3CK0, fileParser_.argd(5)); ffb->setParameter(TorsionFunctions::Cos3CK1, fileParser_.argd(6)); ffb->setParameter(TorsionFunctions::Cos3CK2, fileParser_.argd(7)); - ffb->setParameter(TorsionFunctions::Cos3CK3, fileParser_.argd(7)); + ffb->setParameter(TorsionFunctions::Cos3CK3, fileParser_.argd(8)); } - if (tf != TorsionFunctions::nTorsionFunctions) + else { - ffb->setElecScale(fileParser_.argd(8)); - ffb->setVdwScale(fileParser_.argd(9)); + Messenger::print("Functional form of torsion term (%s) is not present in Aten.", qPrintable(formString)); + continue; } + + Messenger::print("Created torsion term %s-%s-%s-%s", qPrintable(atomNames.at(ii)), qPrintable(atomNames.at(jj)), qPrintable(atomNames.at(kk)), qPrintable(atomNames.at(ll))); } } else if (keyword == "finish") break; @@ -394,7 +394,6 @@ bool DLPExpressionPlugin::importData() // Read in each line of data, searching for those where the first atomtype is equal to the second if (!fileParser_.parseLine()) return false; QString nameI = fileParser_.argc(0); -// readLine(names[1], names[2], form, data[1], data[2], data[3], data[4], data[5]); if (nameI != fileParser_.argc(1)) continue; // We may have added multiple types of the same name earlier, so search the whole list explicitly @@ -415,6 +414,9 @@ bool DLPExpressionPlugin::importData() } } + // Finally, convert energy based on unit read in earlier + newFF->convertParameters(); + return true; }