use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class PDBFilter method writeSIFTFile.
public boolean writeSIFTFile(File saveFile, boolean append, String[] resAndScore) {
if (saveFile == null) {
return false;
}
if (vdwH) {
logger.info(" Printing hydrogens to van der Waals centers instead of nuclear locations.");
}
if (nSymOp != 0) {
logger.info(String.format(" Printing atoms with symmetry operator %s", activeMolecularAssembly.getCrystal().spaceGroup.getSymOp(nSymOp).toString()));
}
/**
* Create StringBuilders for ATOM, ANISOU and TER records that can be
* reused.
*/
StringBuilder sb = new StringBuilder("ATOM ");
StringBuilder anisouSB = new StringBuilder("ANISOU");
StringBuilder terSB = new StringBuilder("TER ");
for (int i = 6; i < 80; i++) {
sb.append(' ');
anisouSB.append(' ');
terSB.append(' ');
}
FileWriter fw;
BufferedWriter bw;
try {
File newFile = saveFile;
if (!append && !noVersioning) {
newFile = version(saveFile);
}
activeMolecularAssembly.setFile(newFile);
activeMolecularAssembly.setName(newFile.getName());
if (logWrites) {
logger.log(Level.INFO, " Saving {0}", newFile.getName());
}
fw = new FileWriter(newFile, append);
bw = new BufferedWriter(fw);
// =============================================================================
// The CRYST1 record presents the unit cell parameters, space group, and Z
// value. If the structure was not determined by crystallographic means, CRYST1
// simply provides the unitary values, with an appropriate REMARK.
//
// 7 - 15 Real(9.3) a a (Angstroms).
// 16 - 24 Real(9.3) b b (Angstroms).
// 25 - 33 Real(9.3) c c (Angstroms).
// 34 - 40 Real(7.2) alpha alpha (degrees).
// 41 - 47 Real(7.2) beta beta (degrees).
// 48 - 54 Real(7.2) gamma gamma (degrees).
// 56 - 66 LString sGroup Space group.
// 67 - 70 Integer z Z value.
// =============================================================================
Crystal crystal = activeMolecularAssembly.getCrystal();
if (crystal != null && !crystal.aperiodic()) {
Crystal c = crystal.getUnitCell();
if (!listMode) {
bw.write(format("CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %10s\n", c.a, c.b, c.c, c.alpha, c.beta, c.gamma, padRight(c.spaceGroup.pdbName, 10)));
} else {
listOutput.add(format("CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %10s", c.a, c.b, c.c, c.alpha, c.beta, c.gamma, padRight(c.spaceGroup.pdbName, 10)));
}
}
// =============================================================================
// The SSBOND record identifies each disulfide bond in protein and polypeptide
// structures by identifying the two residues involved in the bond.
// The disulfide bond distance is included after the symmetry operations at
// the end of the SSBOND record.
//
// 8 - 10 Integer serNum Serial number.
// 12 - 14 LString(3) "CYS" Residue name.
// 16 Character chainID1 Chain identifier.
// 18 - 21 Integer seqNum1 Residue sequence number.
// 22 AChar icode1 Insertion code.
// 26 - 28 LString(3) "CYS" Residue name.
// 30 Character chainID2 Chain identifier.
// 32 - 35 Integer seqNum2 Residue sequence number.
// 36 AChar icode2 Insertion code.
// 60 - 65 SymOP sym1 Symmetry oper for 1st resid
// 67 - 72 SymOP sym2 Symmetry oper for 2nd resid
// 74 – 78 Real(5.2) Length Disulfide bond distance
//
// If SG of cysteine is disordered then there are possible alternate linkages.
// wwPDB practice is to put together all possible SSBOND records. This is
// problematic because the alternate location identifier is not specified in
// the SSBOND record.
// =============================================================================
int serNum = 1;
Polymer[] polymers = activeMolecularAssembly.getChains();
if (polymers != null) {
for (Polymer polymer : polymers) {
ArrayList<Residue> residues = polymer.getResidues();
for (Residue residue : residues) {
if (residue.getName().equalsIgnoreCase("CYS")) {
List<Atom> cysAtoms = residue.getAtomList();
Atom SG1 = null;
for (Atom atom : cysAtoms) {
if (atom.getName().equalsIgnoreCase("SG")) {
SG1 = atom;
break;
}
}
List<Bond> bonds = SG1.getBonds();
for (Bond bond : bonds) {
Atom SG2 = bond.get1_2(SG1);
if (SG2.getName().equalsIgnoreCase("SG")) {
if (SG1.getIndex() < SG2.getIndex()) {
bond.energy(false);
if (!listMode) {
bw.write(format("SSBOND %3d CYS %1s %4s CYS %1s %4s %36s %5.2f\n", serNum++, SG1.getChainID().toString(), Hybrid36.encode(4, SG1.getResidueNumber()), SG2.getChainID().toString(), Hybrid36.encode(4, SG2.getResidueNumber()), "", bond.getValue()));
} else {
listOutput.add(format("SSBOND %3d CYS %1s %4s CYS %1s %4s %36s %5.2f\n", serNum++, SG1.getChainID().toString(), Hybrid36.encode(4, SG1.getResidueNumber()), SG2.getChainID().toString(), Hybrid36.encode(4, SG2.getResidueNumber()), "", bond.getValue()));
}
}
}
}
}
}
}
}
// =============================================================================
//
// 7 - 11 Integer serial Atom serial number.
// 13 - 16 Atom name Atom name.
// 17 Character altLoc Alternate location indicator.
// 18 - 20 Residue name resName Residue name.
// 22 Character chainID Chain identifier.
// 23 - 26 Integer resSeq Residue sequence number.
// 27 AChar iCode Code for insertion of residues.
// 31 - 38 Real(8.3) x Orthogonal coordinates for X in Angstroms.
// 39 - 46 Real(8.3) y Orthogonal coordinates for Y in Angstroms.
// 47 - 54 Real(8.3) z Orthogonal coordinates for Z in Angstroms.
// 55 - 60 Real(6.2) occupancy Occupancy.
// 61 - 66 Real(6.2) tempFactor Temperature factor.
// 77 - 78 LString(2) element Element symbol, right-justified.
// 79 - 80 LString(2) charge Charge on the atom.
// =============================================================================
// 1 2 3 4 5 6 7
// 123456789012345678901234567890123456789012345678901234567890123456789012345678
// ATOM 1 N ILE A 16 60.614 71.140 -10.592 1.00 7.38 N
// ATOM 2 CA ILE A 16 60.793 72.149 -9.511 1.00 6.91 C
MolecularAssembly[] molecularAssemblies = this.getMolecularAssemblys();
int serial = 1;
// Loop over biomolecular chains
if (polymers != null) {
for (Polymer polymer : polymers) {
currentSegID = polymer.getName();
currentChainID = polymer.getChainID();
sb.setCharAt(21, currentChainID);
// Loop over residues
ArrayList<Residue> residues = polymer.getResidues();
for (Residue residue : residues) {
String resName = residue.getName();
if (resName.length() > 3) {
resName = resName.substring(0, 3);
}
int resID = residue.getResidueNumber();
int i = 0;
String[] entries = null;
for (; i < resAndScore.length; i++) {
entries = resAndScore[i].split("\\t");
if (!entries[0].equals(entries[0].replaceAll("\\D+", ""))) {
String[] subEntries = entries[0].split("[^0-9]");
entries[0] = subEntries[0];
}
if (entries[0].equals(String.valueOf(resID)) && !".".equals(entries[1])) {
break;
}
}
sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
// Loop over atoms
ArrayList<Atom> residueAtoms = residue.getAtomList();
boolean altLocFound = false;
for (Atom atom : residueAtoms) {
if (i != resAndScore.length) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, entries[1]);
} else {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
}
Character altLoc = atom.getAltLoc();
if (altLoc != null && !altLoc.equals(' ')) {
altLocFound = true;
}
}
// Write out alternate conformers
if (altLocFound) {
for (int ma = 1; ma < molecularAssemblies.length; ma++) {
MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
Polymer altPolymer = altMolecularAssembly.getPolymer(currentChainID, currentSegID, false);
Residue altResidue = altPolymer.getResidue(resName, resID, false);
residueAtoms = altResidue.getAtomList();
for (Atom atom : residueAtoms) {
if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
if (i != resAndScore.length) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, entries[1]);
} else {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
}
}
}
}
}
}
terSB.replace(6, 11, String.format("%5s", Hybrid36.encode(5, serial++)));
terSB.replace(12, 16, " ");
terSB.replace(16, 26, sb.substring(16, 26));
if (!listMode) {
bw.write(terSB.toString());
bw.newLine();
} else {
listOutput.add(terSB.toString());
}
}
}
sb.replace(0, 6, "HETATM");
sb.setCharAt(21, 'A');
int resID = 1;
Polymer polymer = activeMolecularAssembly.getPolymer('A', "A", false);
if (polymer != null) {
ArrayList<Residue> residues = polymer.getResidues();
for (Residue residue : residues) {
int resID2 = residue.getResidueNumber();
if (resID2 >= resID) {
resID = resID2 + 1;
}
}
}
/**
* Loop over molecules, ions and then water.
*/
ArrayList<Molecule> molecules = activeMolecularAssembly.getMolecules();
for (int i = 0; i < molecules.size(); i++) {
Molecule molecule = (Molecule) molecules.get(i);
Character chainID = molecule.getChainID();
sb.setCharAt(21, chainID);
String resName = molecule.getResidueName();
if (resName.length() > 3) {
resName = resName.substring(0, 3);
}
sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
ArrayList<Atom> moleculeAtoms = molecule.getAtomList();
boolean altLocFound = false;
for (Atom atom : moleculeAtoms) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
Character altLoc = atom.getAltLoc();
if (altLoc != null && !altLoc.equals(' ')) {
altLocFound = true;
}
}
// Write out alternate conformers
if (altLocFound) {
for (int ma = 1; ma < molecularAssemblies.length; ma++) {
MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
MSNode altmolecule = altMolecularAssembly.getMolecules().get(i);
moleculeAtoms = altmolecule.getAtomList();
for (Atom atom : moleculeAtoms) {
if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
}
}
}
}
resID++;
}
ArrayList<MSNode> ions = activeMolecularAssembly.getIons();
for (int i = 0; i < ions.size(); i++) {
Molecule ion = (Molecule) ions.get(i);
Character chainID = ion.getChainID();
sb.setCharAt(21, chainID);
String resName = ion.getResidueName();
if (resName.length() > 3) {
resName = resName.substring(0, 3);
}
sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
ArrayList<Atom> ionAtoms = ion.getAtomList();
boolean altLocFound = false;
for (Atom atom : ionAtoms) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
Character altLoc = atom.getAltLoc();
if (altLoc != null && !altLoc.equals(' ')) {
altLocFound = true;
}
}
// Write out alternate conformers
if (altLocFound) {
for (int ma = 1; ma < molecularAssemblies.length; ma++) {
MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
MSNode altion = altMolecularAssembly.getIons().get(i);
ionAtoms = altion.getAtomList();
for (Atom atom : ionAtoms) {
if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
}
}
}
}
resID++;
}
ArrayList<MSNode> waters = activeMolecularAssembly.getWaters();
for (int i = 0; i < waters.size(); i++) {
Molecule water = (Molecule) waters.get(i);
Character chainID = water.getChainID();
sb.setCharAt(21, chainID);
String resName = water.getResidueName();
if (resName.length() > 3) {
resName = resName.substring(0, 3);
}
sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
ArrayList<Atom> waterAtoms = water.getAtomList();
boolean altLocFound = false;
for (Atom atom : waterAtoms) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
Character altLoc = atom.getAltLoc();
if (altLoc != null && !altLoc.equals(' ')) {
altLocFound = true;
}
}
// Write out alternate conformers
if (altLocFound) {
for (int ma = 1; ma < molecularAssemblies.length; ma++) {
MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
MSNode altwater = altMolecularAssembly.getWaters().get(i);
waterAtoms = altwater.getAtomList();
for (Atom atom : waterAtoms) {
if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
}
}
}
}
resID++;
}
if (!listMode) {
bw.write("END");
bw.newLine();
} else {
listOutput.add("END");
}
bw.close();
} catch (Exception e) {
String message = "Exception writing to file: " + saveFile.toString();
logger.log(Level.WARNING, message, e);
return false;
}
return true;
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class PDBFilter method assignAtomTypes.
/**
* Assign force field atoms types to common chemistries using "biotype"
* records.
*/
private void assignAtomTypes() {
/**
* Create a list to store bonds defined by PDB atom names.
*/
bondList = new ArrayList<>();
/**
* To Do: Look for cyclic peptides and disulfides.
*/
Polymer[] polymers = activeMolecularAssembly.getChains();
/**
* Loop over chains.
*/
if (polymers != null) {
logger.info(format("\n Assigning atom types for %d chains.", polymers.length));
for (Polymer polymer : polymers) {
List<Residue> residues = polymer.getResidues();
int numberOfResidues = residues.size();
/**
* Check if all residues are known amino acids.
*/
boolean isProtein = true;
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
Residue residue = residues.get(residueNumber);
String name = residue.getName().toUpperCase();
boolean aa = false;
for (AminoAcid3 amino : aminoAcidList) {
if (amino.toString().equalsIgnoreCase(name)) {
aa = true;
checkHydrogenAtomNames(residue);
break;
}
}
// Check for a patch.
if (!aa) {
logger.info(" Checking for non-standard amino acid patch " + name);
HashMap<String, AtomType> types = forceField.getAtomTypes(name);
if (types.isEmpty()) {
isProtein = false;
break;
} else {
logger.info(" Patch found for non-standard amino acid " + name);
}
}
}
/**
* If all the residues in this chain have known amino acids
* names, then attempt to assign atom types.
*/
if (isProtein) {
try {
logger.info(format(" Amino acid chain %s", polymer.getName()));
double dist = properties.getDouble("chainbreak", 3.0);
// Detect main chain breaks!
List<List<Residue>> subChains = findChainBreaks(residues, dist);
for (List<Residue> subChain : subChains) {
assignAminoAcidAtomTypes(subChain);
}
} catch (MissingHeavyAtomException missingHeavyAtomException) {
logger.severe(missingHeavyAtomException.toString());
} catch (MissingAtomTypeException missingAtomTypeException) {
logger.severe(missingAtomTypeException.toString());
}
continue;
}
/**
* Check if all residues have known nucleic acids names.
*/
boolean isNucleicAcid = true;
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
Residue residue = residues.get(residueNumber);
String name = residue.getName().toUpperCase();
/**
* Convert 1 and 2-character nucleic acid names to
* 3-character names.
*/
if (name.length() == 1) {
if (name.equals("A")) {
name = NucleicAcid3.ADE.toString();
} else if (name.equals("C")) {
name = NucleicAcid3.CYT.toString();
} else if (name.equals("G")) {
name = NucleicAcid3.GUA.toString();
} else if (name.equals("T")) {
name = NucleicAcid3.THY.toString();
} else if (name.equals("U")) {
name = NucleicAcid3.URI.toString();
}
} else if (name.length() == 2) {
if (name.equals("YG")) {
name = NucleicAcid3.YYG.toString();
}
}
residue.setName(name);
NucleicAcid3 nucleicAcid = null;
for (NucleicAcid3 nucleic : nucleicAcidList) {
String nuc3 = nucleic.toString();
nuc3 = nuc3.substring(nuc3.length() - 3);
if (nuc3.equalsIgnoreCase(name)) {
nucleicAcid = nucleic;
break;
}
}
if (nucleicAcid == null) {
logger.info(format("Nucleic acid was not recognized %s.", name));
isNucleicAcid = false;
break;
}
}
/**
* If all the residues in this chain have known nucleic acids
* names, then attempt to assign atom types.
*/
if (isNucleicAcid) {
try {
logger.info(format(" Nucleic acid chain %s", polymer.getName()));
assignNucleicAcidAtomTypes(residues, forceField, bondList);
} catch (MissingHeavyAtomException | MissingAtomTypeException e) {
logger.severe(e.toString());
}
}
}
}
// Assign ion atom types.
ArrayList<MSNode> ions = activeMolecularAssembly.getIons();
if (ions != null && ions.size() > 0) {
logger.info(format(" Assigning atom types for %d ions.", ions.size()));
for (MSNode m : ions) {
Molecule ion = (Molecule) m;
String name = ion.getResidueName().toUpperCase();
HetAtoms hetatm = HetAtoms.valueOf(name);
Atom atom = ion.getAtomList().get(0);
if (ion.getAtomList().size() != 1) {
logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
}
try {
switch(hetatm) {
case NA:
atom.setAtomType(findAtomType(2003));
break;
case K:
atom.setAtomType(findAtomType(2004));
break;
case MG:
case MG2:
atom.setAtomType(findAtomType(2005));
break;
case CA:
case CA2:
atom.setAtomType(findAtomType(2006));
break;
case CL:
atom.setAtomType(findAtomType(2007));
break;
case ZN:
case ZN2:
atom.setAtomType(findAtomType(2008));
break;
case BR:
atom.setAtomType(findAtomType(2009));
break;
default:
logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
}
} catch (Exception e) {
String message = "Error assigning atom types.";
logger.log(Level.SEVERE, message, e);
}
}
}
// Assign water atom types.
ArrayList<MSNode> water = activeMolecularAssembly.getWaters();
if (water != null && water.size() > 0) {
logger.info(format(" Assigning atom types for %d waters.", water.size()));
for (MSNode m : water) {
Molecule wat = (Molecule) m;
try {
Atom O = buildHeavy(wat, "O", null, 2001);
Atom H1 = buildHydrogen(wat, "H1", O, 0.96e0, null, 109.5e0, null, 120.0e0, 0, 2002);
H1.setHetero(true);
Atom H2 = buildHydrogen(wat, "H2", O, 0.96e0, H1, 109.5e0, null, 120.0e0, 0, 2002);
H2.setHetero(true);
} catch (Exception e) {
String message = "Error assigning atom types to a water.";
logger.log(Level.SEVERE, message, e);
}
}
}
// Assign small molecule atom types.
ArrayList<Molecule> molecules = activeMolecularAssembly.getMolecules();
for (MSNode m : molecules) {
Molecule molecule = (Molecule) m;
String moleculeName = molecule.getResidueName();
logger.info(" Attempting to patch " + moleculeName);
ArrayList<Atom> moleculeAtoms = molecule.getAtomList();
boolean patched = true;
HashMap<String, AtomType> types = forceField.getAtomTypes(moleculeName);
/**
* Assign atom types for all known atoms.
*/
for (Atom atom : moleculeAtoms) {
String atomName = atom.getName().toUpperCase();
AtomType atomType = types.get(atomName);
if (atomType == null) {
logger.info(" No atom type was found for " + atomName + " of " + moleculeName + ".");
patched = false;
break;
} else {
logger.fine(" " + atom.toString() + " -> " + atomType.toString());
atom.setAtomType(atomType);
types.remove(atomName);
}
}
/**
* Create missing hydrogen atoms. Check for missing heavy atoms.
*/
if (patched && !types.isEmpty()) {
for (AtomType type : types.values()) {
if (type.atomicNumber != 1) {
logger.info(" Missing heavy atom " + type.name);
patched = false;
break;
}
}
}
// Create bonds between known atoms.
if (patched) {
for (Atom atom : moleculeAtoms) {
String atomName = atom.getName();
String[] bonds = forceField.getBonds(moleculeName, atomName);
if (bonds != null) {
for (String name : bonds) {
Atom atom2 = molecule.getAtom(name);
if (atom2 != null && !atom.isBonded(atom2)) {
buildBond(atom, atom2);
}
}
}
}
}
// Create missing hydrogen atoms.
if (patched && !types.isEmpty()) {
// Create a hashmap of the molecule's atoms
HashMap<String, Atom> atomMap = new HashMap<String, Atom>();
for (Atom atom : moleculeAtoms) {
atomMap.put(atom.getName().toUpperCase(), atom);
}
for (String atomName : types.keySet()) {
AtomType type = types.get(atomName);
String[] bonds = forceField.getBonds(moleculeName, atomName.toUpperCase());
if (bonds == null || bonds.length != 1) {
patched = false;
logger.info(" Check biotype for hydrogen " + type.name + ".");
break;
}
// Get the heavy atom the hydrogen is bonded to.
Atom ia = atomMap.get(bonds[0].toUpperCase());
Atom hydrogen = new Atom(0, atomName, ia.getAltLoc(), new double[3], ia.getResidueName(), ia.getResidueNumber(), ia.getChainID(), ia.getOccupancy(), ia.getTempFactor(), ia.getSegID());
logger.fine(" Created hydrogen " + atomName + ".");
hydrogen.setAtomType(type);
hydrogen.setHetero(true);
molecule.addMSNode(hydrogen);
int valence = ia.getAtomType().valence;
List<Bond> aBonds = ia.getBonds();
int numBonds = aBonds.size();
/**
* Try to find the following configuration: ib-ia-ic
*/
Atom ib = null;
Atom ic = null;
Atom id = null;
if (numBonds > 0) {
Bond bond = aBonds.get(0);
ib = bond.get1_2(ia);
}
if (numBonds > 1) {
Bond bond = aBonds.get(1);
ic = bond.get1_2(ia);
}
if (numBonds > 2) {
Bond bond = aBonds.get(2);
id = bond.get1_2(ia);
}
/**
* Building the hydrogens depends on hybridization and the
* locations of other bonded atoms.
*/
logger.fine(" Bonding " + atomName + " to " + ia.getName() + " (" + numBonds + " of " + valence + ").");
switch(valence) {
case 4:
switch(numBonds) {
case 3:
// Find the average coordinates of atoms ib, ic and id.
double[] b = ib.getXYZ(null);
double[] c = ib.getXYZ(null);
double[] d = ib.getXYZ(null);
double[] a = new double[3];
a[0] = (b[0] + c[0] + d[0]) / 3.0;
a[1] = (b[1] + c[1] + d[1]) / 3.0;
a[2] = (b[2] + c[2] + d[2]) / 3.0;
// Place the hydrogen at chiral position #1.
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 1);
double[] e1 = new double[3];
hydrogen.getXYZ(e1);
double[] ret = new double[3];
diff(a, e1, ret);
double l1 = r(ret);
// Place the hydrogen at chiral position #2.
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, -1);
double[] e2 = new double[3];
hydrogen.getXYZ(e2);
diff(a, e2, ret);
double l2 = r(ret);
// Revert to #1 if it is farther from the average.
if (l1 > l2) {
hydrogen.setXYZ(e1);
}
break;
case 2:
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 0);
break;
case 1:
intxyz(hydrogen, ia, 1.0, ib, 109.5, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 3:
switch(numBonds) {
case 2:
intxyz(hydrogen, ia, 1.0, ib, 120.0, ic, 0.0, 0);
break;
case 1:
intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 2:
switch(numBonds) {
case 1:
intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 1:
switch(numBonds) {
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
if (!patched) {
break;
} else {
buildBond(ia, hydrogen);
}
}
}
if (!patched) {
logger.log(Level.WARNING, format(" Deleting unrecognized molecule %s.", m.toString()));
activeMolecularAssembly.deleteMolecule((Molecule) m);
} else {
logger.info(" Patch for " + moleculeName + " succeeded.");
}
}
resolvePolymerLinks(molecules);
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class XYZFilter method writeFile.
/**
* {@inheritDoc}
*/
@Override
public boolean writeFile(File saveFile, boolean append) {
if (saveFile == null) {
return false;
}
try {
File newFile = saveFile;
if (!append) {
newFile = version(saveFile);
}
activeMolecularAssembly.setFile(newFile);
activeMolecularAssembly.setName(newFile.getName());
FileWriter fw;
if (append && !newFile.exists()) {
fw = new FileWriter(newFile);
} else {
fw = new FileWriter(newFile, append);
}
BufferedWriter bw = new BufferedWriter(fw);
// XYZ File First Line
int numberOfAtoms = activeMolecularAssembly.getAtomList().size();
String output = format("%7d %s\n", numberOfAtoms, activeMolecularAssembly.toString());
bw.write(output);
Crystal crystal = activeMolecularAssembly.getCrystal();
if (!crystal.aperiodic()) {
Crystal uc = crystal.getUnitCell();
String params = String.format("%14.8f%14.8f%14.8f%14.8f%14.8f%14.8f\n", uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
bw.write(params);
}
Atom a2;
StringBuilder line;
StringBuilder[] lines = new StringBuilder[numberOfAtoms];
// XYZ File Atom Lines
ArrayList<Atom> atoms = activeMolecularAssembly.getAtomList();
Vector3d offset = activeMolecularAssembly.getOffset();
for (Atom a : atoms) {
if (vdwH) {
line = new StringBuilder(String.format("%7d %3s%14.8f%14.8f%14.8f%6d", a.getIndex(), a.getAtomType().name, a.getRedX() - offset.x, a.getRedY() - offset.y, a.getRedZ() - offset.z, a.getType()));
} else {
line = new StringBuilder(String.format("%7d %3s%14.8f%14.8f%14.8f%6d", a.getIndex(), a.getAtomType().name, a.getX() - offset.x, a.getY() - offset.y, a.getZ() - offset.z, a.getType()));
}
for (Bond b : a.getBonds()) {
a2 = b.get1_2(a);
line.append(format("%8d", a2.getIndex()));
}
lines[a.getIndex() - 1] = line.append("\n");
}
try {
for (int i = 0; i < numberOfAtoms; i++) {
bw.write(lines[i].toString());
}
} catch (IOException e) {
String message = format(" Their was an unexpected error writing to %s.", getActiveMolecularSystem().toString());
logger.log(Level.WARNING, message, e);
return false;
}
bw.close();
fw.close();
} catch (IOException e) {
String message = format(" Their was an unexpected error writing to %s.", getActiveMolecularSystem().toString());
logger.log(Level.WARNING, message, e);
return false;
}
return true;
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class BiojavaFilter method assignAtomTypes.
/**
* Assign force field atoms types to common chemistries using "biotype"
* records.
*/
public void assignAtomTypes() {
/**
* Create a new List to store bonds determined based on PDB atom names.
*/
bondList = new ArrayList<>();
/**
* To Do: Look for cyclic peptides and disulfides.
*/
Polymer[] polymers = activeMolecularAssembly.getChains();
/**
* Loop over chains.
*/
if (polymers != null) {
logger.info(format("\n Assigning atom types for %d chains.", polymers.length));
for (Polymer polymer : polymers) {
List<Residue> residues = polymer.getResidues();
int numberOfResidues = residues.size();
/**
* Check if all residues are known amino acids.
*/
boolean isProtein = true;
if (!residues.isEmpty()) {
// renameNTerminusHydrogens(residues.get(0)); Not safe to use until it distinguishes between true N-termini and N-terminal residues in general.
}
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
Residue residue = residues.get(residueNumber);
String name = residue.getName().toUpperCase();
boolean aa = false;
for (AminoAcid3 amino : aminoAcidList) {
if (amino.toString().equalsIgnoreCase(name)) {
aa = true;
renameNonstandardHydrogens(residue);
break;
}
}
// Check for a patch.
if (!aa) {
HashMap<String, AtomType> types = forceField.getAtomTypes(name);
if (types.isEmpty()) {
isProtein = false;
break;
} else {
logger.info(" Patch found for non-standard amino acid " + name);
}
}
}
/**
* If all the residues in this chain have known amino acids
* names, then attempt to assign atom types.
*/
if (isProtein) {
try {
logger.info(format(" Amino acid chain %s", polymer.getName()));
double dist = properties.getDouble("chainbreak", 3.0);
// Detect main chain breaks!
List<List<Residue>> subChains = findChainBreaks(residues, dist);
for (List<Residue> subChain : subChains) {
assignAminoAcidAtomTypes(subChain);
}
} catch (MissingHeavyAtomException missingHeavyAtomException) {
logger.severe(missingHeavyAtomException.toString());
} catch (MissingAtomTypeException missingAtomTypeException) {
logger.severe(missingAtomTypeException.toString());
}
continue;
}
/**
* Check if all residues have known nucleic acids names.
*/
boolean isNucleicAcid = true;
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
Residue residue = residues.get(residueNumber);
String name = residue.getName().toUpperCase();
/**
* Convert 1 and 2-character nucleic acid names to
* 3-character names.
*/
if (name.length() == 1) {
if (name.equals("A")) {
name = NucleicAcid3.ADE.toString();
} else if (name.equals("C")) {
name = NucleicAcid3.CYT.toString();
} else if (name.equals("G")) {
name = NucleicAcid3.GUA.toString();
} else if (name.equals("T")) {
name = NucleicAcid3.THY.toString();
} else if (name.equals("U")) {
name = NucleicAcid3.URI.toString();
}
} else if (name.length() == 2) {
if (name.equals("YG")) {
name = NucleicAcid3.YYG.toString();
}
}
residue.setName(name);
NucleicAcid3 nucleicAcid = null;
for (NucleicAcid3 nucleic : nucleicAcidList) {
String nuc3 = nucleic.toString();
nuc3 = nuc3.substring(nuc3.length() - 3);
if (nuc3.equalsIgnoreCase(name)) {
nucleicAcid = nucleic;
break;
}
}
if (nucleicAcid == null) {
logger.info(format("Nucleic acid was not recognized %s.", name));
isNucleicAcid = false;
break;
}
}
/**
* If all the residues in this chain have known nucleic acids
* names, then attempt to assign atom types.
*/
if (isNucleicAcid) {
try {
logger.info(format(" Nucleic acid chain %s", polymer.getName()));
assignNucleicAcidAtomTypes(residues);
} catch (MissingHeavyAtomException missingHeavyAtomException) {
logger.severe(missingHeavyAtomException.toString());
} catch (MissingAtomTypeException missingAtomTypeException) {
logger.severe(missingAtomTypeException.toString());
}
}
}
}
// Assign ion atom types.
ArrayList<MSNode> ions = activeMolecularAssembly.getIons();
if (ions != null && ions.size() > 0) {
logger.info(format(" Assigning atom types for %d ions.", ions.size()));
for (MSNode m : ions) {
Molecule ion = (Molecule) m;
String name = ion.getResidueName().toUpperCase();
HetAtoms hetatm = HetAtoms.valueOf(name);
Atom atom = ion.getAtomList().get(0);
if (ion.getAtomList().size() != 1) {
logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
}
try {
switch(hetatm) {
case NA:
atom.setAtomType(findAtomType(2003));
break;
case K:
atom.setAtomType(findAtomType(2004));
break;
case MG:
case MG2:
atom.setAtomType(findAtomType(2005));
break;
case CA:
case CA2:
atom.setAtomType(findAtomType(2006));
break;
case CL:
atom.setAtomType(findAtomType(2007));
break;
case ZN:
case ZN2:
atom.setAtomType(findAtomType(2008));
break;
case BR:
atom.setAtomType(findAtomType(2009));
break;
default:
logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
}
} catch (Exception e) {
String message = "Error assigning atom types.";
logger.log(Level.SEVERE, message, e);
}
}
}
// Assign water atom types.
ArrayList<MSNode> water = activeMolecularAssembly.getWaters();
if (water != null && water.size() > 0) {
logger.info(format(" Assigning atom types for %d waters.", water.size()));
for (MSNode m : water) {
Molecule wat = (Molecule) m;
try {
Atom O = buildHeavy(wat, "O", null, 2001);
Atom H1 = buildHydrogen(wat, "H1", O, 0.96e0, null, 109.5e0, null, 120.0e0, 0, 2002);
H1.setHetero(true);
Atom H2 = buildHydrogen(wat, "H2", O, 0.96e0, H1, 109.5e0, null, 120.0e0, 0, 2002);
H2.setHetero(true);
} catch (Exception e) {
String message = "Error assigning atom types to a water.";
logger.log(Level.SEVERE, message, e);
}
}
}
// Assign small molecule atom types.
ArrayList<Molecule> molecules = activeMolecularAssembly.getMolecules();
for (MSNode m : molecules) {
Molecule molecule = (Molecule) m;
String moleculeName = molecule.getResidueName();
logger.info(" Attempting to patch " + moleculeName);
ArrayList<Atom> moleculeAtoms = molecule.getAtomList();
boolean patched = true;
HashMap<String, AtomType> types = forceField.getAtomTypes(moleculeName);
/**
* Assign atom types for all known atoms.
*/
for (Atom atom : moleculeAtoms) {
String atomName = atom.getName().toUpperCase();
AtomType atomType = types.get(atomName);
if (atomType == null) {
logger.info(" No atom type was found for " + atomName + " of " + moleculeName + ".");
patched = false;
break;
} else {
atom.setAtomType(atomType);
types.remove(atomName);
}
}
/**
* Create missing hydrogen atoms. Check for missing heavy atoms.
*/
if (patched && !types.isEmpty()) {
for (AtomType type : types.values()) {
if (type.atomicNumber != 1) {
logger.info(" Missing heavy atom " + type.name);
patched = false;
break;
}
}
}
// Create bonds between known atoms.
if (patched) {
for (Atom atom : moleculeAtoms) {
String atomName = atom.getName();
String[] bonds = forceField.getBonds(moleculeName, atomName);
if (bonds != null) {
for (String name : bonds) {
Atom atom2 = molecule.getAtom(name);
if (atom2 != null && !atom.isBonded(atom2)) {
buildBond(atom, atom2);
}
}
}
}
}
// Create missing hydrogen atoms.
if (patched && !types.isEmpty()) {
// Create a hashmap of the molecule's atoms
HashMap<String, Atom> atomMap = new HashMap<String, Atom>();
for (Atom atom : moleculeAtoms) {
atomMap.put(atom.getName().toUpperCase(), atom);
}
for (String atomName : types.keySet()) {
AtomType type = types.get(atomName);
String[] bonds = forceField.getBonds(moleculeName, atomName.toUpperCase());
if (bonds == null || bonds.length != 1) {
patched = false;
logger.info(" Check biotype for hydrogen " + type.name + ".");
break;
}
// Get the heavy atom the hydrogen is bonded to.
Atom ia = atomMap.get(bonds[0].toUpperCase());
Atom hydrogen = new Atom(0, atomName, ia.getAltLoc(), new double[3], ia.getResidueName(), ia.getResidueNumber(), ia.getChainID(), ia.getOccupancy(), ia.getTempFactor(), ia.getSegID());
logger.fine(" Created hydrogen " + atomName + ".");
hydrogen.setAtomType(type);
hydrogen.setHetero(true);
molecule.addMSNode(hydrogen);
int valence = ia.getAtomType().valence;
List<Bond> aBonds = ia.getBonds();
int numBonds = aBonds.size();
/**
* Try to find the following configuration: ib-ia-ic
*/
Atom ib = null;
Atom ic = null;
Atom id = null;
if (numBonds > 0) {
Bond bond = aBonds.get(0);
ib = bond.get1_2(ia);
}
if (numBonds > 1) {
Bond bond = aBonds.get(1);
ic = bond.get1_2(ia);
}
if (numBonds > 2) {
Bond bond = aBonds.get(2);
id = bond.get1_2(ia);
}
/**
* Building the hydrogens depends on hybridization and the
* locations of other bonded atoms.
*/
logger.fine(" Bonding " + atomName + " to " + ia.getName() + " (" + numBonds + " of " + valence + ").");
switch(valence) {
case 4:
switch(numBonds) {
case 3:
// Find the average coordinates of atoms ib, ic and id.
double[] b = ib.getXYZ(null);
double[] c = ib.getXYZ(null);
double[] d = ib.getXYZ(null);
double[] a = new double[3];
a[0] = (b[0] + c[0] + d[0]) / 3.0;
a[1] = (b[1] + c[1] + d[1]) / 3.0;
a[2] = (b[2] + c[2] + d[2]) / 3.0;
// Place the hydrogen at chiral position #1.
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 1);
double[] e1 = new double[3];
hydrogen.getXYZ(e1);
double[] ret = new double[3];
diff(a, e1, ret);
double l1 = r(ret);
// Place the hydrogen at chiral position #2.
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, -1);
double[] e2 = new double[3];
hydrogen.getXYZ(e2);
diff(a, e2, ret);
double l2 = r(ret);
// Revert to #1 if it is farther from the average.
if (l1 > l2) {
hydrogen.setXYZ(e1);
}
break;
case 2:
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 0);
break;
case 1:
intxyz(hydrogen, ia, 1.0, ib, 109.5, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 3:
switch(numBonds) {
case 2:
intxyz(hydrogen, ia, 1.0, ib, 120.0, ic, 0.0, 0);
break;
case 1:
intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 2:
switch(numBonds) {
case 1:
intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 1:
switch(numBonds) {
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
if (!patched) {
break;
} else {
buildBond(ia, hydrogen);
}
}
}
if (!patched) {
logger.log(Level.WARNING, format(" Deleting unrecognized molecule %s.", m.toString()));
activeMolecularAssembly.deleteMolecule((Molecule) m);
} else {
logger.info(" Patch for " + moleculeName + " succeeded.");
}
}
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class MultipoleType method multipoleTypeFactory.
public static MultipoleType multipoleTypeFactory(Atom atom, ForceField forceField) {
AtomType atomType = atom.getAtomType();
if (atomType == null) {
String message = " Multipoles can only be assigned to atoms that have been typed.";
logger.severe(message);
return null;
}
PolarizeType polarizeType = forceField.getPolarizeType(atomType.getKey());
if (polarizeType != null) {
atom.setPolarizeType(polarizeType);
} else {
String message = " No polarization type was found for " + atom.toString();
logger.fine(message);
double polarizability = 0.0;
double thole = 0.0;
int[] polarizationGroup = null;
polarizeType = new PolarizeType(atomType.type, polarizability, thole, polarizationGroup);
forceField.addForceFieldType(polarizeType);
atom.setPolarizeType(polarizeType);
}
String key;
// No reference atoms.
key = atomType.getKey() + " 0 0";
MultipoleType multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
atom.setAxisAtoms(null);
return multipoleType;
}
// No bonds.
List<Bond> bonds = atom.getBonds();
if (bonds == null || bonds.size() < 1) {
String message = "Multipoles can only be assigned after bonded relationships are defined.\n";
logger.severe(message);
}
// 1 reference atom.
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
key = atomType.getKey() + " " + atom2.getAtomType().getKey() + " 0";
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
atom.setAxisAtoms(atom2);
return multipoleType;
}
}
// 2 reference atoms.
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
String key2 = atom2.getAtomType().getKey();
for (Bond b2 : bonds) {
if (b == b2) {
continue;
}
Atom atom3 = b2.get1_2(atom);
String key3 = atom3.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
atom.setAxisAtoms(atom2, atom3);
return multipoleType;
}
}
}
/**
* 3 reference atoms.
*/
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
String key2 = atom2.getAtomType().getKey();
for (Bond b2 : bonds) {
if (b == b2) {
continue;
}
Atom atom3 = b2.get1_2(atom);
String key3 = atom3.getAtomType().getKey();
for (Bond b3 : bonds) {
if (b == b3 || b2 == b3) {
continue;
}
Atom atom4 = b3.get1_2(atom);
String key4 = atom4.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
atom.setAxisAtoms(atom2, atom3, atom4);
return multipoleType;
}
}
List<Angle> angles = atom.getAngles();
for (Angle angle : angles) {
Atom atom4 = angle.get1_3(atom);
if (atom4 != null) {
String key4 = atom4.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
atom.setAxisAtoms(atom2, atom3, atom4);
return multipoleType;
}
}
}
}
}
/**
* Revert to a 2 reference atom definition that may include a 1-3 site.
* For example a hydrogen on water.
*/
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
String key2 = atom2.getAtomType().getKey();
List<Angle> angles = atom.getAngles();
for (Angle angle : angles) {
Atom atom3 = angle.get1_3(atom);
if (atom3 != null) {
String key3 = atom3.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
atom.setAxisAtoms(atom2, atom3);
return multipoleType;
}
for (Angle angle2 : angles) {
Atom atom4 = angle2.get1_3(atom);
if (atom4 != null && atom4 != atom3) {
String key4 = atom4.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
atom.setAxisAtoms(atom2, atom3, atom4);
return multipoleType;
}
}
}
}
}
}
return null;
}
Aggregations