package name.mjw.jamber.IO.OpenMM;
import com.google.common.collect.HashMultimap;
import com.google.common.collect.Multimap;
import name.mjw.jamber.IO.AMBER.AngleType;
import name.mjw.jamber.IO.AMBER.Atom;
import name.mjw.jamber.IO.AMBER.AtomType;
import name.mjw.jamber.IO.AMBER.BondType;
import name.mjw.jamber.IO.AMBER.Connection;
import name.mjw.jamber.IO.AMBER.DihedralType;
import name.mjw.jamber.IO.AMBER.Lib;
import name.mjw.jamber.IO.AMBER.ParameterStore;
import name.mjw.jamber.IO.AMBER.Residue;
import name.mjw.jamber.util.ChemicalSymbol;
import nu.xom.*;
import org.apache.log4j.Logger;
import java.io.IOException;
import java.io.OutputStream;
import java.text.DecimalFormat;
import java.util.*;
import java.util.Map.Entry;
/**
* Representation of OpenMM's XML format as found in their Python interface
* simtk/openmm/app/data/*.xml .
* <p>
*
* Contains the ability to convert AMBER {@link ParameterStore} (parm) and
* {@link Lib} files into a set of XML formats suitable for use by OpenMM's
* Python wrappers.
*
* @author mjw
*
*/
class OpenMMXML {
private final Logger LOG = Logger.getLogger(OpenMMXML.class);
private final Double ANGSTROM_TO_NANOMETER = 0.1;
private final Double KCAL_TO_KJOULES = 4.184;
/**
* XML needed by simtk.openmm.app.forcefield(). Contains Force field terms
* and parameters.
*/
private Document ffDoc;
/**
* XML needed by simtk.openmm.app.topology.loadBondDefinitions() to parse a
* PDB file; it contains a list of bonds. This is then used by
* simtk.openmm.app.topology.createStandardBonds().
*/
private Document residuesDoc;
/**
* XML needed by simtk.openmm.app.modeller.loadHydrogenDefinitions().
* Contains a list of hydrogen atoms and which parent atoms they are bonded
* to.
*/
private Document hydrogensDoc;
/**
* XML needed by simtk.openmm.app.modeller.GBSAOBCForce.
*/
private Document GBSA_OBCDoc;
/**
* Parameters passed in
*/
private final ParameterStore parameterStore;
/**
* Topologies passed in
*/
private final Lib lib;
private final Calendar instance = Calendar.getInstance();
private final long timeOfDay = instance.getTimeInMillis();
/**
* The name attribute that links an Atom element to Atom elements within a
* Residue element and Atom elements within the NonbondedForce element,
* needs to be unique across *ALL* XML files passed to the ForceField
* constructor. Hence, use the time of day as a name prefix. Not the best
* solution, but it is functional.
*/
private final String uid = Long.toString(timeOfDay) + "_";
/**
* Store only the atomTypes found in residues, so that only parameters
* containing those are harvested from the Parm objects
*/
private HashSet<String> requiredAtomTypes = new HashSet<>();
/**
* Harvest all externalAtomTypes found in residues. Needed because extra
* parameters will be sought for these.
*/
private HashSet<String> requiredExternalAtomTypes = new HashSet<>();
/**
* Generate an OpenMMXML object from an AMBER parm file and library file.
*
* @param parameterStore
* Parm object
* @param lib
* Lib object
*
*/
OpenMMXML(ParameterStore parameterStore, Lib lib) {
this.lib = lib;
this.parameterStore = parameterStore;
buildFFDoc();
buildHydrogensDoc();
buildResiduesDoc();
buildGBSA_OBC_Doc();
// removeDuplicateTypes();
}
/**
* Within the AMBER force field, an atomtype, the two letter label given to
* an atom in a topology (e.g. CT), does NOT completely define its
* interaction; more information is required. Specifically, the mapping
* between an atomtype and its charge is NOT an injective (one-to-one)
* function. This is due to the nature of the RESP protocol; the atom's
* charge is actually a function of its molecular environment.
* <p>
*
* Hence for an atomtype to be completely unique, it must have the SAME
* atomtype AND charge. Within OpenMM's parlance, this uniqueness is termed
* a "type", and the two letter label is termed "class".
* <p>
*
* During the construction the {@link OpenMMXML#ffDoc} XML, it is assumed
* every single atom within the passed {@link OpenMMXML#lib} object is
* unique. This may not be always true if there exits multiple duplicate
* pairs of identical atomtype/charges, i.e. OpenMM types. This method
* removes and resolves these types.
*
*
*/
@SuppressWarnings("unused")
private void removeDuplicateTypes() {
Multimap<String, String> classChargeToTypeMultiMap = HashMultimap
.create();
/*
* 1) Get all types present i.e.
*
* <ForceField> <AtomTypes> <Type name="0" class="nc" element="N"
* mass="14.01"/> <Type name="1" class="cc" element="C" mass="12.01"/>
* </> </>
*
* will return: 0, 1
*/
Nodes atomTypes;
atomTypes = ffDoc.query("//AtomTypes/Type/@name");
/*
* 2) For each type, look up its class attribute in //AtomTypes/Type and
* charge attribute from //NonbondedForce/Atom and add it to the
* MultiMap, classChargeToTypeMultiMap, object.
*
* MultiMap.keySet() can be used to determine unique type and
* MultiMap.key() can be used to determine which are duplicate types
*
* Key Value
*
* CT -0.0048 --> [1366600275087_682] HC 0.0627 --> [1366600275087_776,
* 1366600275087_777, 1366600275087_778] CT -0.4121 -->
* [1366600275087_239, 1366600275087_243] HC 0.0621 -->
* [1366600275087_282, 1366600275087_281]
*/
Nodes atomCharge;
Nodes atomClass;
for (int i = 0; i < atomTypes.size(); i++) {
String currentAtomType = atomTypes.get(i).getValue();
atomClass = ffDoc.query("//AtomTypes/Type[ @name='"
+ currentAtomType + "' ]/@class");
atomCharge = ffDoc.query("//NonbondedForce/Atom[ @type='"
+ currentAtomType + "']/@charge");
String uniqueClassChargeKey = atomClass.get(0).getValue() + " "
+ atomCharge.get(0).getValue();
LOG.debug("uniqueClassChargeKey:currentAtomType \t"
+ uniqueClassChargeKey + ":" + currentAtomType);
classChargeToTypeMultiMap
.put(uniqueClassChargeKey, currentAtomType);
}
LOG.debug("Total number of AtomTypes: "
+ classChargeToTypeMultiMap.keys().size());
LOG.debug("Total number of unique AtomTypes: "
+ classChargeToTypeMultiMap.keySet().size());
/*
* 3) Build a map, classChargeToTypeMultiMap, linking from all types to
* the first instance of the duplicated type.
*/
Map<String, String> typeToUniqueType = new HashMap<>();
for (String uniqueKey : classChargeToTypeMultiMap.keySet()) {
Collection<String> values = classChargeToTypeMultiMap
.get(uniqueKey);
// Get first instance of the duplicated type and use it as value
String zeroValue = (String) values.toArray()[0];
for (String value : values) {
typeToUniqueType.put(value, zeroValue);
}
}
LOG.debug(typeToUniqueType.size());
for (Entry<String, String> entry : typeToUniqueType.entrySet()) {
LOG.debug(entry.getKey() + " --> " + entry.getValue());
}
/*
* 4) Map the following:
*
* //AtomTypes/Type[@name="old"]/@name=typeToUniqueTyp.getKey(old)
* //Residues/Residue[@type="old"]/@type=typeToUniqueType.getKey(old)
* //NonbondedForce/Atom[@type="old"]/@type=typeToUniqueType.getKey(old)
*/
/*
* 5) Remove duplicate entries from //AtomTypes/
*/
}
private void buildHydrogensDoc() {
// Build the basic XML structure
generateHydrogensDocSkeleton();
// Navigate to the parent Residues element
Element element = hydrogensDoc.getRootElement();
for (Element item : generateHydrogenAtomInformation()) {
element.appendChild(item);
}
}
private void buildResiduesDoc() {
// Build the basic XML structure
generateResiduesDocSkeleton();
// Navigate to the parent Residues element
Element element = residuesDoc.getRootElement();
for (Element item : generateResidueElementsForResiduesDoc()) {
element.appendChild(item);
}
}
private void buildFFDoc() {
// Build the basic XML structure
generateFFDocSkeleton();
// Navigate to the AtomType element
Element element = (Element) ffDoc.query("//AtomTypes").get(0);
// Populate AtomType element with children
// <Type name="1360511932278_0" class="na" element="N" mass="14.01"/>
for (Element item : generateAtomElements()) {
element.appendChild(item);
}
// Navigate to the Residues element
element = (Element) ffDoc.query("//Residues").get(0);
// Populate Residues element with children
// <Residue name="KTN">
// <Atom name="N" type="1360512194219_0"/>
// <Atom name="C" type="1360512194219_1"/>
// ...
// <Bond from="0" to="1"/>
// <Bond from="0" to="4"/>
for (Element item : generateResidueElements()) {
element.appendChild(item);
}
// Navigate to the HarmonicBondForce element
element = (Element) ffDoc.query("//HarmonicBondForce").get(0);
// Populate HarmonicBondForce element with children
// <Bond class1="c3" class2="c3" length="0.1535" k="253634.1"/>
for (Element item : generateHarmonicBondForceElements()) {
element.appendChild(item);
}
// Navigate to the HarmonicAngleForce element
element = (Element) ffDoc.query("//HarmonicAngleForce").get(0);
// Populate HarmonicAngleForce element with children
// <Angle class1="c3" class2="c3" class3="c3" angle="1.93085775148"
// k="528.941"/>
for (Element item : generateHarmonicAngleForceElements()) {
element.appendChild(item);
}
// Navigate to the PeriodicTorsionForce element
element = (Element) ffDoc.query("//PeriodicTorsionForce").get(0);
// Populate PeriodicTorsionForce element with children
// <Proper class1="n" class2="c3" class3="c" class4="n" periodicity1="1"
// periodicity2="2" k1="7.1128" k2="8.368" phase1="3.14159265359"
// phase2="3.14159265359"/>
// ...
// <Improper class1="ca" class2="c" class3="ca" class4="c3"
// periodicity1="2" k1="4.6024" phase1="3.14159265359"/>
for (Element item : generatePeriodicTorsionForceElements()) {
element.appendChild(item);
}
// Navigate to the NonbondedForce element
element = (Element) ffDoc.query("//NonbondedForce").get(0);
// Populate NonbondedForce element with children
// <Atom type="1360512194219_0" charge="-0.2039" sigma="0.324999852378"
// epsilon="0.71128"/>
for (Element item : generateNonbondedForceElements()) {
element.appendChild(item);
}
}
private void buildGBSA_OBC_Doc() {
// Build the basic XML structure
generateGBSA_OBC_DocSkeleton();
// Navigate to the GBSAOBCForce element
Element element = (Element) GBSA_OBCDoc.query("//GBSAOBCForce").get(0);
for (Element item : generateGBSAOBCForceElements()) {
element.appendChild(item);
}
}
/**
* Write object to an output stream.
*/
private void toXMLStream(Document doc, OutputStream out) {
try {
Serializer serializer = new Serializer(out, "ISO-8859-1");
serializer.setIndent(2);
serializer.setMaxLength(120);
serializer.write(doc);
} catch (IOException ex) {
System.err.println(ex);
}
}
/**
* Sysout dump of the object
*/
private void toXML(Document doc) {
try {
Serializer serializer = new Serializer(System.out, "ISO-8859-1");
serializer.setIndent(2);
serializer.setMaxLength(120);
serializer.write(doc);
} catch (IOException ex) {
System.err.println(ex);
}
}
/**
* Create the bare bones OpenMM FF XML Document that will be populated
* later.
*/
private void generateFFDocSkeleton() {
// Root node
Element ForceFieldElement = new Element("ForceField");
// AtomTypes
Element AtomTypesElement = new Element("AtomTypes");
// Residues
Element ResiduesElement = new Element("Residues");
// Harmonic Bond
Element HarmonicBondForceElement = new Element("HarmonicBondForce");
// Harmonic angle
Element HarmonicAngleForceElement = new Element("HarmonicAngleForce");
// Torsions
Element PeriodicTorsionForceElement = new Element(
"PeriodicTorsionForce");
// Non bonded
Attribute coulomb14scaleAttribute = new Attribute("coulomb14scale",
"0.833333");
Attribute lj14scaleAttribute = new Attribute("lj14scale", "0.5");
Element NonbondedForceElement = new Element("NonbondedForce");
NonbondedForceElement.addAttribute(coulomb14scaleAttribute);
NonbondedForceElement.addAttribute(lj14scaleAttribute);
// Append all elements to the root node
ForceFieldElement.appendChild(AtomTypesElement);
ForceFieldElement.appendChild(ResiduesElement);
ForceFieldElement.appendChild(HarmonicBondForceElement);
ForceFieldElement.appendChild(HarmonicAngleForceElement);
ForceFieldElement.appendChild(PeriodicTorsionForceElement);
ForceFieldElement.appendChild(NonbondedForceElement);
ffDoc = new Document(ForceFieldElement);
}
/**
* Create the bare bones OpenMM residues XML Document that will be populated
* later.
*/
private void generateResiduesDocSkeleton() {
// Root node
Element ResiduesElement = new Element("Residues");
residuesDoc = new Document(ResiduesElement);
}
/**
* Create the bare bones OpenMM hydrogens XML Document that will be
* populated later.
*/
private void generateHydrogensDocSkeleton() {
// Root node
Element ResiduesElement = new Element("Residues");
hydrogensDoc = new Document(ResiduesElement);
}
/**
* Create the bare bones OpenMM GBSA_OBC XML Document that will be populated
* later.
*/
private void generateGBSA_OBC_DocSkeleton() {
// Root node
Element ForceFieldElement = new Element("ForceField");
// GBSAOBC
Element GBSAOBCForceElement = new Element("GBSAOBCForce");
ForceFieldElement.appendChild(GBSAOBCForceElement);
GBSA_OBCDoc = new Document(ForceFieldElement);
}
/**
* Creates the Atom elements.
*
* @return List of Atom elements
*/
private ArrayList<Element> generateAtomElements() {
ArrayList<Element> elements = new ArrayList<>();
// Enumerate ALL atoms found within the lib file
for (int i = 0; i < lib.getAtoms().size(); i++) {
String CurrentAtomType = lib.getAtomByIndex(i).getType();
// Store this atomType for later
requiredAtomTypes.add(CurrentAtomType);
Element TypeElement = new Element("Type");
Attribute AtomTypeName = new Attribute("name", uid
+ Integer.toString(i));
Attribute AtomTypeClass = new Attribute("class", CurrentAtomType);
// Map the Z value to the Chemical Symbol
String chemicalElement = ChemicalSymbol.zToChemicalSymbol(lib
.getAtomByIndex(i).getZ());
Attribute AtomTypeElement = new Attribute("element",
chemicalElement);
AtomType atomtype = parameterStore.getAtomByType(CurrentAtomType);
Attribute AtomTypeMass = new Attribute("mass", atomtype.getMass()
.toString());
TypeElement.addAttribute(AtomTypeName);
TypeElement.addAttribute(AtomTypeClass);
TypeElement.addAttribute(AtomTypeElement);
TypeElement.addAttribute(AtomTypeMass);
elements.add(TypeElement);
LOG.debug(TypeElement.toXML());
}
LOG.debug("requiredAtomTypes" + requiredAtomTypes);
return elements;
}
/**
* Creates the Residue elements.
*
* @return List of Residue elements
*/
private ArrayList<Element> generateResidueElements() {
ArrayList<Element> elements = new ArrayList<>();
for (Residue residue : lib.getResidues()) {
for (Atom atom : residue.getExternalBonds()) {
requiredExternalAtomTypes.add(atom.getType());
}
}
LOG.debug("requiredExternalAtomTypes " + requiredExternalAtomTypes);
for (Residue residue : lib.getResidues()) {
Element ResidueElement = new Element("Residue");
Attribute ResidueName = new Attribute("name", residue.getName());
ResidueElement.addAttribute(ResidueName);
for (Atom atom : residue.getAtoms()) {
Element ResidueAtom = new Element("Atom");
Attribute ResidueAtomName = new Attribute("name",
atom.getName());
ResidueAtom.addAttribute(ResidueAtomName);
// Find the corresponding type in the AtomTypes element
// i.e. position in Lib.atoms<>
Attribute ResidueAtomType = new Attribute("type", uid
+ lib.getAtomByIndex(atom).toString());
ResidueAtom.addAttribute(ResidueAtomType);
// <Atom name="O" type="tip3p-O"/>
ResidueElement.appendChild(ResidueAtom);
}
for (Connection connection : residue.getConnections()) {
Element ResidueBond = new Element("Bond");
Attribute ResidueBondFrom = new Attribute("from", connection
.getI().getIndexWithinCurrentResidue().toString());
Attribute ResidueBondTo = new Attribute("to", connection.getJ()
.getIndexWithinCurrentResidue().toString());
ResidueBond.addAttribute(ResidueBondFrom);
ResidueBond.addAttribute(ResidueBondTo);
// <Bond from="0" to="1"/>
ResidueElement.appendChild(ResidueBond);
}
for (Atom atom : residue.getExternalBonds()) {
Element ResidueExternalBond = new Element("ExternalBond");
Attribute ResidueExternalBondFrom = new Attribute("from", atom
.getIndexWithinCurrentResidue().toString());
ResidueExternalBond.addAttribute(ResidueExternalBondFrom);
ResidueElement.appendChild(ResidueExternalBond);
}
LOG.debug(ResidueElement.toXML());
elements.add(ResidueElement);
}
return elements;
}
/**
* Creates the Harmonic bond elements.
*
* @return List of Harmonic bond elements
*/
private ArrayList<Element> generateHarmonicBondForceElements() {
ArrayList<Element> elements = new ArrayList<>();
// Filter for bonds that contain both atom types
ArrayList<BondType> requiredBondTypes = new ArrayList<>();
for (BondType bondType : parameterStore.getBondTypes()) {
if (requiredAtomTypes.contains(bondType.getAtomI().getName())) {
if (requiredAtomTypes.contains(bondType.getAtomB().getName())) {
requiredBondTypes.add(bondType);
}
}
}
// Filter for bonds that contain any requiredExternalAtomTypes
for (BondType bondType : parameterStore.getBondTypes()) {
if (requiredExternalAtomTypes.contains(bondType.getAtomI()
.getName())) {
requiredBondTypes.add(bondType);
}
// TODO may need to check both atomTypes
}
// Create the XML
for (BondType bondType : requiredBondTypes) {
// <Bond class1="OW" class2="HW" length="0.09572" k="462750.4"/>
Element BondElement = new Element("Bond");
Attribute bondClass1Element = new Attribute("class1", bondType
.getAtomI().toString());
Attribute bondClass2Element = new Attribute("class2", bondType
.getAtomB().toString());
Double length = bondType.getEquilibriumLength();
length = length * ANGSTROM_TO_NANOMETER;
DecimalFormat df = new DecimalFormat("#.####");
Attribute bondLengthElement = new Attribute("length",
df.format(length));
// Initially in kcal/(mol Angstrom**2), but needs to be in kJ/(mole
// nanometer**2)
Double k = bondType.getForceConstant();
// Two factor is because OpenMM
k = k * KCAL_TO_KJOULES * 2 * 100;
df = new DecimalFormat("######.#");
Attribute kElement = new Attribute("k", df.format(k));
BondElement.addAttribute(bondClass1Element);
BondElement.addAttribute(bondClass2Element);
BondElement.addAttribute(bondLengthElement);
BondElement.addAttribute(kElement);
LOG.debug(BondElement.toXML());
elements.add(BondElement);
}
return elements;
}
/**
* Creates the Harmonic angle elements.
*
* @return List of Harmonic angle elements
*/
private ArrayList<Element> generateHarmonicAngleForceElements() {
ArrayList<Element> elements = new ArrayList<>();
// Filter for bonds that contain all atom types
ArrayList<AngleType> requiredAngleTypes = new ArrayList<>();
for (AngleType angleType : parameterStore.getAngleTypes()) {
if (requiredAtomTypes.contains(angleType.getAtomI().getName())) {
if (requiredAtomTypes.contains(angleType.getAtomJ().getName())) {
if (requiredAtomTypes.contains(angleType.getAtomK()
.getName())) {
requiredAngleTypes.add(angleType);
}
}
}
}
// Filter for angles that contain any requiredExternalAtomTypes
for (AngleType angleType : parameterStore.getAngleTypes()) {
if (requiredExternalAtomTypes.contains(angleType.getAtomJ()
.getName())) {
requiredAngleTypes.add(angleType);
}
}
// Create the XML
for (AngleType angleType : requiredAngleTypes) {
Element AngleElement = new Element("Angle");
Attribute AngleClass1Element = new Attribute("class1", angleType
.getAtomI().toString());
Attribute AngleClass2Element = new Attribute("class2", angleType
.getAtomJ().toString());
Attribute AngleClass3Element = new Attribute("class3", angleType
.getAtomK().toString());
// Convert to radians from degrees to radians
Double angle = angleType.getEquilibriumAngle();
angle = (angle / 360) * 2 * Math.PI;
DecimalFormat df = new DecimalFormat("#.###########");
Attribute AngleAngleElement = new Attribute("angle",
df.format(angle));
// Convert from kcal/(mol radian**2) to kJoules/(mol radian**2);
Double k = angleType.getForceConstant();
// 2 is from AMBER convention
k = k * KCAL_TO_KJOULES * 2;
df = new DecimalFormat("###.###");
Attribute AnglekElement = new Attribute("k", df.format(k));
AngleElement.addAttribute(AngleClass1Element);
AngleElement.addAttribute(AngleClass2Element);
AngleElement.addAttribute(AngleClass3Element);
AngleElement.addAttribute(AngleAngleElement);
AngleElement.addAttribute(AnglekElement);
LOG.debug(AngleElement.toXML());
elements.add(AngleElement);
}
return elements;
}
/**
* Creates the proper and improper elements.
*
* @return List of proper and improper elements
*/
private ArrayList<Element> generatePeriodicTorsionForceElements() {
ArrayList<Element> elements = new ArrayList<>();
// Filter for dihedrals that contain all atom types
ArrayList<DihedralType> requiredDihedralType = new ArrayList<>();
for (DihedralType dihedralType : parameterStore.getDihedralTypes()) {
if (requiredAtomTypes.contains(dihedralType.getAtomI().getName())) {
if (requiredAtomTypes.contains(dihedralType.getAtomJ()
.getName())) {
if (requiredAtomTypes.contains(dihedralType.getAtomK()
.getName())) {
if (requiredAtomTypes.contains(dihedralType.getAtomL()
.getName())) {
requiredDihedralType.add(dihedralType);
}
}
}
}
}
// Filter for any wildcarded proper dihedrals
for (DihedralType dihedralType : parameterStore
.getProperDihedralTypes()) {
if (requiredAtomTypes.contains(dihedralType.getAtomJ().getName())) {
if (requiredAtomTypes.contains(dihedralType.getAtomK()
.getName())) {
if (dihedralType.getAtomI().getName().contains("X")
&& dihedralType.getAtomL().getName().contains("X")) {
requiredDihedralType.add(dihedralType);
}
}
}
}
// Filter for any wildcarded improper dihedrals
for (DihedralType dihedralType : parameterStore
.getImproperDihedralTypes()) {
if (requiredAtomTypes.contains(dihedralType.getAtomK().getName())) {
if (requiredAtomTypes.contains(dihedralType.getAtomL()
.getName())) {
dihedralType.reorderForOpenMM();
requiredDihedralType.add(dihedralType);
}
}
}
// Filter for dihedrals that contain any requiredExternalAtomTypes
for (DihedralType dihedralType : parameterStore.getDihedralTypes()) {
if (requiredExternalAtomTypes.contains(dihedralType.getAtomJ()
.getName())) {
if (requiredExternalAtomTypes.contains(dihedralType.getAtomK()
.getName())) {
dihedralType.reorderForOpenMM();
requiredDihedralType.add(dihedralType);
}
}
}
// Write the XML
// Proper
for (DihedralType dihedral : requiredDihedralType) {
if (dihedral instanceof name.mjw.jamber.IO.AMBER.ProperDihedralType) {
Element ProperElement = new Element("Proper");
Attribute ProperClass1Element = new Attribute("class1",
dihedral.getAtomI().toString());
Attribute ProperClass2Element = new Attribute("class2",
dihedral.getAtomJ().toString());
Attribute ProperClass3Element = new Attribute("class3",
dihedral.getAtomK().toString());
Attribute ProperClass4Element = new Attribute("class4",
dihedral.getAtomL().toString());
ProperElement.addAttribute(ProperClass1Element);
ProperElement.addAttribute(ProperClass2Element);
ProperElement.addAttribute(ProperClass3Element);
ProperElement.addAttribute(ProperClass4Element);
Integer counter = 1;
for (Integer period : dihedral.getPeriodicity()) {
String periodicityName = "periodicity" + counter.toString();
Attribute ProperPeriodicity = new Attribute(
periodicityName, period.toString());
ProperElement.addAttribute(ProperPeriodicity);
counter = counter + 1;
}
counter = 1;
for (Double height : dihedral.getBarrierHeight()) {
Double convertedHeight = (height * KCAL_TO_KJOULES)
/ dihedral.getMultiplicity();
String barrierHeightName = "k" + counter.toString();
DecimalFormat df = new DecimalFormat("#.#####");
Attribute barrierHeightAttribute = new Attribute(
barrierHeightName, df.format(convertedHeight));
ProperElement.addAttribute(barrierHeightAttribute);
counter = counter + 1;
}
counter = 1;
for (Double phase : dihedral.getPhase()) {
phase = (phase / 360) * 2 * Math.PI;
String phaseName = "phase" + counter.toString();
DecimalFormat df = new DecimalFormat("#.###########");
Attribute phaseAttribute = new Attribute(phaseName,
df.format(phase));
ProperElement.addAttribute(phaseAttribute);
counter = counter + 1;
}
LOG.debug(ProperElement.toXML());
elements.add(ProperElement);
}
}
// ImproperProper
for (DihedralType dihedral : requiredDihedralType) {
if (dihedral instanceof name.mjw.jamber.IO.AMBER.ImproperDihedralType) {
Element ImproperElement = new Element("Improper");
Attribute ImproperClass1Element = new Attribute("class1",
dihedral.getAtomI().toString());
Attribute ImproperClass2Element = new Attribute("class2",
dihedral.getAtomJ().toString());
Attribute ImproperClass3Element = new Attribute("class3",
dihedral.getAtomK().toString());
Attribute ImproperClass4Element = new Attribute("class4",
dihedral.getAtomL().toString());
ImproperElement.addAttribute(ImproperClass1Element);
ImproperElement.addAttribute(ImproperClass2Element);
ImproperElement.addAttribute(ImproperClass3Element);
ImproperElement.addAttribute(ImproperClass4Element);
Integer counter = 1;
for (Integer period : dihedral.getPeriodicity()) {
String periodicityName = "periodicity" + counter.toString();
Attribute ProperPeriodicity = new Attribute(
periodicityName, period.toString());
ImproperElement.addAttribute(ProperPeriodicity);
counter = counter + 1;
}
counter = 1;
for (Double height : dihedral.getBarrierHeight()) {
Double convertedHeight = (height * KCAL_TO_KJOULES);
String barrierHeightName = "k" + counter.toString();
DecimalFormat df = new DecimalFormat("#.#####");
Attribute barrierHeightAttribute = new Attribute(
barrierHeightName, df.format(convertedHeight));
ImproperElement.addAttribute(barrierHeightAttribute);
counter = counter + 1;
}
counter = 1;
for (Double phase : dihedral.getPhase()) {
phase = (phase / 360) * 2 * Math.PI;
String phaseName = "phase" + counter.toString();
DecimalFormat df = new DecimalFormat("#.###########");
Attribute phaseAttribute = new Attribute(phaseName,
df.format(phase));
ImproperElement.addAttribute(phaseAttribute);
counter = counter + 1;
}
LOG.debug(ImproperElement.toXML());
elements.add(ImproperElement);
}
}
return elements;
}
/**
* Creates the non-bonded elements.
*
* @return List of non-bonded elements
*/
private ArrayList<Element> generateNonbondedForceElements() {
ArrayList<Element> elements = new ArrayList<>();
// <Atom type="1960" charge="1.0" sigma="0.526699322165"
// epsilon="0.00071128"/>
for (int i = 0; i < lib.getAtoms().size(); i++) {
Atom atom = lib.getAtomByIndex(i);
Element AtomElement = new Element("Atom");
Attribute NonbondedForceElementTypeAttribute = new Attribute(
"type", uid + Integer.toString(i));
AtomElement.addAttribute(NonbondedForceElementTypeAttribute);
DecimalFormat df = new DecimalFormat("#.####");
Attribute NonbondedForceElementChargeAttribute = new Attribute(
"charge", df.format(atom.getCharge()));
AtomElement.addAttribute(NonbondedForceElementChargeAttribute);
// Now look up the atomtype from the parm file
AtomType atomType = parameterStore.getAtomByType(atom.getType());
/*
* Convert values; getVdwRadius is in the form of r_min and in
* Angstroms, this needs to be converted to sigma and into nm.
*
* See section 14.6.1 of the theory OpenMM user's manual
*/
Double sigma = (atomType.getVdwRadius() * 2 * ANGSTROM_TO_NANOMETER)
/ (Math.pow(2.0, (1.0 / 6)));
df = new DecimalFormat("#.############");
Attribute NonbondedForceElementSigmaAttribute = new Attribute(
"sigma", df.format(sigma));
// Convert values; epsilon is in kcal/mol, but should be in kJ
Double epsilon = atomType.getVdwWellDepth() * KCAL_TO_KJOULES;
df = new DecimalFormat("#.#######");
Attribute NonbondedForceElementEpsilonAttribute = new Attribute(
"epsilon", df.format(epsilon));
AtomElement.addAttribute(NonbondedForceElementSigmaAttribute);
AtomElement.addAttribute(NonbondedForceElementEpsilonAttribute);
LOG.debug(AtomElement.toXML());
elements.add(AtomElement);
}
return elements;
}
/**
* Used by Modeller.addHydrogens in the simtk/openmm/app/data/hydrogens.xml
* to add missing hydrogen atoms on a residue when building from a PDB that
* does not have them.
*
* @return
*/
private ArrayList<Element> generateHydrogenAtomInformation() {
ArrayList<Element> elements = new ArrayList<>();
for (Residue residue : lib.getResidues()) {
Element ResidueElement = new Element("Residue");
Attribute ResidueName = new Attribute("name", residue.getName());
ResidueElement.addAttribute(ResidueName);
// H atom information for Modeller.addHydrogens()
for (Atom atom : residue.getAtoms()) {
// Select only hydrogens
if (atom.getZ().equals(1)) {
Element ResidueH = new Element("H");
Attribute ResidueHName = new Attribute("name",
atom.getName());
ResidueH.addAttribute(ResidueHName);
// Look over connections in the residue
for (Connection connection : residue.getConnections()) {
// TODO super fragile; we are assuming the 2nd atom is
// always H
if (connection.getJ().equals(atom)) {
Attribute ResidueHParent = new Attribute("parent",
connection.getI().getName());
ResidueH.addAttribute(ResidueHParent);
}
}
ResidueElement.appendChild(ResidueH);
}
}
LOG.debug(ResidueElement.toXML());
elements.add(ResidueElement);
}
return elements;
}
private ArrayList<Element> generateResidueElementsForResiduesDoc() {
ArrayList<Element> elements = new ArrayList<>();
for (Residue residue : lib.getResidues()) {
for (Atom atom : residue.getExternalBonds()) {
requiredExternalAtomTypes.add(atom.getType());
}
}
LOG.debug("requiredExternalAtomTypes " + requiredExternalAtomTypes);
for (Residue residue : lib.getResidues()) {
Element ResidueElement = new Element("Residue");
Attribute ResidueName = new Attribute("name", residue.getName());
ResidueElement.addAttribute(ResidueName);
/*
* Use the AMBER atom type to decide how it will connect to the next
* residue.
*/
for (Atom atom : residue.getExternalBonds()) {
if (atom.getType().equals("N")) {
Element ResidueBond = new Element("Bond");
Attribute ResidueBondFrom = new Attribute("from", "-C");
Attribute ResidueBondTo = new Attribute("to",
atom.getName());
ResidueBond.addAttribute(ResidueBondFrom);
ResidueBond.addAttribute(ResidueBondTo);
ResidueElement.appendChild(ResidueBond);
}
if (atom.getType().equals("C")) {
Element ResidueBond = new Element("Bond");
Attribute ResidueBondTo = new Attribute("to", "+N");
Attribute ResidueBondFrom = new Attribute("from",
atom.getName());
ResidueBond.addAttribute(ResidueBondFrom);
ResidueBond.addAttribute(ResidueBondTo);
ResidueElement.appendChild(ResidueBond);
}
}
for (Connection connection : residue.getConnections()) {
Element ResidueBond = new Element("Bond");
Attribute ResidueBondFrom = new Attribute("from", connection
.getI().getName());
Attribute ResidueBondTo = new Attribute("to", connection.getJ()
.getName());
ResidueBond.addAttribute(ResidueBondFrom);
ResidueBond.addAttribute(ResidueBondTo);
// <Bond from="0" to="1"/>
ResidueElement.appendChild(ResidueBond);
}
LOG.debug(ResidueElement.toXML());
elements.add(ResidueElement);
}
return elements;
}
/**
* Add in GBSA-OBC solvation parameters:
* http://dx.doi.org/10.1002/prot.20033 This uses model II and mimics LeaP's
* parameter assignment algorithm. See
* amber12/AmberTools/src/leap/src/unitio.c:5071
*
* "H(N)-modified Bondi radii (mbondi2)"
*
*/
private ArrayList<Element> generateGBSAOBCForceElements() {
ArrayList<Element> GBSAOBCForceElements = new ArrayList<>();
DecimalFormat df = new DecimalFormat("#.####");
for (int i = 0; i < lib.getAtoms().size(); i++) {
Element AtomElement = new Element("Atom");
Atom atom = lib.getAtomByIndex(i);
Attribute GBSAOBCTypeAttribute = new Attribute("type", uid
+ Integer.toString(i));
Attribute GBSAOBCChargeAttribute = new Attribute("charge",
df.format(atom.getCharge()));
Attribute GBSAOBCRadiusAttribute = new Attribute("radius",
df.format(getGBSAOBCIIRadius(atom)));
Attribute GBSAOBCScaleAttribute = new Attribute("scale",
df.format(getGBSAOBCIIScale(atom)));
AtomElement.addAttribute(GBSAOBCTypeAttribute);
AtomElement.addAttribute(GBSAOBCChargeAttribute);
AtomElement.addAttribute(GBSAOBCRadiusAttribute);
AtomElement.addAttribute(GBSAOBCScaleAttribute);
GBSAOBCForceElements.add(AtomElement);
}
return GBSAOBCForceElements;
}
private Double getGBSAOBCIIRadius(Atom atom) {
Double radius;
switch (atom.getZ()) {
case 1:
radius = 1.2;
// TODO set to 1.3 if its neighbour is Nitrogen
break;
case 6:
radius = 2.2;
// TODO set to 1.7 if type = C{1,2,3}
break;
case 7:
radius = 1.55;
break;
case 8:
radius = 1.5;
break;
case 9:
radius = 1.5;
break;
case 14:
radius = 2.1;
break;
case 15:
radius = 1.85;
break;
case 16:
radius = 1.8;
break;
case 17:
radius = 1.7;
break;
default:
radius = 1.5;
break;
}
radius = (radius * ANGSTROM_TO_NANOMETER);
return radius;
}
private Double getGBSAOBCIIScale(Atom atom) {
Double scale;
switch (atom.getZ()) {
case 1:
scale = 0.85;
break;
case 6:
scale = 0.72;
break;
case 7:
scale = 0.79;
break;
case 8:
scale = 0.85;
break;
case 9:
scale = 0.88;
break;
case 15:
scale = 0.86;
break;
case 16:
scale = 0.96;
break;
default:
scale = 0.8;
break;
}
return scale;
}
/**
* Sysout dump of the of the XML needed by XML needed by
* simtk.openmm.app.forcefield(). Contains Force field terms and parameters.
*/
public void toFFXML() {
toXML(ffDoc);
}
/**
* Sysout dump of the of the XML needed by
* simtk.openmm.app.Topology.loadBondDefinitions() to parse a PDB file.
*/
public void toResiduesXML() {
toXML(residuesDoc);
}
/**
* Sysout dump of the of the XML needed by
* simtk.openmm.app.Modeller.loadHydrogenDefinitions(). Contains a list of
* hydrogen atoms and which atom they are bonded to.
*/
public void toHydrogensXML() {
toXML(hydrogensDoc);
}
/**
* Write to an output stream, the XML needed by simtk.openmm.app.forcefield.
* Contains Force field terms and parameters.
*/
public void toFFXMLOutputStream(OutputStream os) {
toXMLStream(ffDoc, os);
}
/**
* Write to an output stream, the XML needed by simtk.openmm.app.modeller to
* parse a PDB file.
*/
public void toResiduesXMLOutputStream(OutputStream os) {
toXMLStream(residuesDoc, os);
}
/**
* Write to an output stream, the XML needed by
* simtk.openmm.app.modeller.addHydrogens(). Contains a list of hydrogen
* atoms and which atom they are bonded to.
*/
public void toHydrogensXMLOutputStream(OutputStream os) {
toXMLStream(hydrogensDoc, os);
}
/**
* Write to an output stream, the XML needed by
* simtk.openmm.app.modeller.GBSAOBCForce
*/
public void toGBSA_OBC_XMLOutputStream(OutputStream os) {
toXMLStream(GBSA_OBCDoc, os);
}
}