use of ffx.potential.parameters.BondType in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method addRestraintBonds.
/**
* Adds restraint bonds, if any.
*/
private void addRestraintBonds() {
List<RestraintBond> restraintBonds = super.getRestraintBonds();
if (restraintBonds != null && !restraintBonds.isEmpty()) {
// OpenMM's HarmonicBondForce class uses k, not 1/2*k as does FFX.
double kParameterConversion = BondType.units * 2.0 * OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
for (RestraintBond rbond : super.getRestraintBonds()) {
// Is not a valid substitute for a targeting computer.
PointerByReference theForce;
BondType bType = rbond.bondType;
BondType.BondFunction funct = bType.bondFunction;
if (restraintForces.containsKey(funct)) {
theForce = restraintForces.get(funct);
} else {
theForce = OpenMM_CustomBondForce_create(funct.toMathematicalForm());
OpenMM_CustomBondForce_addPerBondParameter(theForce, "k");
OpenMM_CustomBondForce_addPerBondParameter(theForce, "r0");
if (funct.hasFlatBottom()) {
OpenMM_CustomBondForce_addPerBondParameter(theForce, "fb");
}
// Wholly untested code.
switch(funct) {
case QUARTIC:
case FLAT_BOTTOM_QUARTIC:
OpenMM_CustomBondForce_addGlobalParameter(theForce, "cubic", BondType.cubic / OpenMM_NmPerAngstrom);
OpenMM_CustomBondForce_addGlobalParameter(theForce, "quartic", BondType.quartic / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
break;
default:
break;
}
OpenMM_System_addForce(system, theForce);
}
double forceConst = bType.forceConstant * kParameterConversion;
double equilDist = bType.distance * OpenMM_NmPerAngstrom;
Atom[] ats = rbond.getAtomArray();
int at1 = ats[0].getXyzIndex() - 1;
int at2 = ats[1].getXyzIndex() - 1;
PointerByReference bondParams = OpenMM_DoubleArray_create(0);
OpenMM_DoubleArray_append(bondParams, forceConst);
OpenMM_DoubleArray_append(bondParams, equilDist);
if (funct.hasFlatBottom()) {
OpenMM_DoubleArray_append(bondParams, bType.flatBottomRadius * OpenMM_NmPerAngstrom);
}
OpenMM_CustomBondForce_addBond(theForce, at1, at2, bondParams);
OpenMM_DoubleArray_destroy(bondParams);
}
}
}
use of ffx.potential.parameters.BondType in project ffx by mjschnie.
the class PDBFilter method buildDisulfideBonds.
/**
* Assign parameters to disulfide bonds.
*
* @param ssBondList
*/
private void buildDisulfideBonds(List<Bond> ssBondList) {
StringBuilder sb = new StringBuilder(" Disulfide Bonds:");
for (Bond bond : ssBondList) {
Atom a1 = bond.getAtom(0);
Atom a2 = bond.getAtom(1);
int[] c = new int[2];
c[0] = a1.getAtomType().atomClass;
c[1] = a2.getAtomType().atomClass;
String key = BondType.sortKey(c);
BondType bondType = forceField.getBondType(key);
if (bondType == null) {
logger.severe(format("No BondType for key: %s\n %s\n %s", key, a1, a2));
} else {
bond.setBondType(bondType);
}
double d = VectorMath.dist(a1.getXYZ(null), a2.getXYZ(null));
Polymer c1 = activeMolecularAssembly.getChain(a1.getSegID());
Polymer c2 = activeMolecularAssembly.getChain(a2.getSegID());
Residue r1 = c1.getResidue(a1.getResidueNumber());
Residue r2 = c2.getResidue(a2.getResidueNumber());
sb.append(format("\n S-S distance of %6.2f for %s and %s.", d, r1.toString(), r2.toString()));
bondList.add(bond);
}
if (ssBondList.size() > 0) {
logger.info(sb.toString());
}
}
use of ffx.potential.parameters.BondType in project ffx by mjschnie.
the class XYZFilter method readFile.
/**
* {@inheritDoc}
*
* Parse the XYZ File
*/
@Override
public boolean readFile() {
File xyzFile = activeMolecularAssembly.getFile();
if (forceField == null) {
logger.warning(format(" No force field is associated with %s.", xyzFile.toString()));
return false;
}
try {
FileReader fr = new FileReader(xyzFile);
BufferedReader br = new BufferedReader(fr);
String data = br.readLine();
// Read blank lines at the top of the file
while (data != null && data.trim().equals("")) {
data = br.readLine();
}
if (data == null) {
return false;
}
String[] tokens = data.trim().split(" +", 2);
int numberOfAtoms = Integer.parseInt(tokens[0]);
if (numberOfAtoms < 1) {
return false;
}
if (tokens.length == 2) {
getActiveMolecularSystem().setName(tokens[1]);
}
logger.info(format("\n Opening %s with %d atoms\n", xyzFile.getName(), numberOfAtoms));
// The header line is reasonable. Check for periodic box dimensions.
br.mark(10000);
data = br.readLine();
if (!readPBC(data, activeMolecularAssembly)) {
br.reset();
}
// Prepare to parse atom lines.
HashMap<Integer, Integer> labelHash = new HashMap<>();
int[] label = new int[numberOfAtoms];
int[][] bonds = new int[numberOfAtoms][8];
double[][] d = new double[numberOfAtoms][3];
boolean renumber = false;
atomList = new ArrayList<>();
// Loop over the expected number of atoms.
for (int i = 0; i < numberOfAtoms; i++) {
if (!br.ready()) {
return false;
}
data = br.readLine();
if (data == null) {
logger.warning(format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
return false;
}
tokens = data.trim().split(" +");
if (tokens == null || tokens.length < 6) {
logger.warning(format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
return false;
}
// Valid number of tokens, so try to parse this line.
label[i] = Integer.parseInt(tokens[0]);
// Check for valid atom numbering, or flag for re-numbering.
if (label[i] != i + 1) {
renumber = true;
}
String atomName = tokens[1];
d[i][0] = Double.parseDouble(tokens[2]);
d[i][1] = Double.parseDouble(tokens[3]);
d[i][2] = Double.parseDouble(tokens[4]);
int type = Integer.parseInt(tokens[5]);
AtomType atomType = forceField.getAtomType(Integer.toString(type));
if (atomType == null) {
StringBuilder message = new StringBuilder("Check atom type ");
message.append(type).append(" for Atom ").append(i + 1);
message.append(" in ").append(activeMolecularAssembly.getFile().getName());
logger.warning(message.toString());
return false;
}
Atom a = new Atom(i + 1, atomName, atomType, d[i]);
atomList.add(a);
// Bond Data
int numberOfBonds = tokens.length - 6;
for (int b = 0; b < 8; b++) {
if (b < numberOfBonds) {
int bond = Integer.parseInt(tokens[6 + b]);
bonds[i][b] = bond;
} else {
bonds[i][b] = 0;
}
}
}
// Check if this is an archive.
if (br.ready()) {
// Read past blank lines between archive files
data = br.readLine().trim();
while (data != null && data.equals("")) {
data = br.readLine().trim();
}
if (data != null) {
tokens = data.split(" +", 2);
if (tokens != null && tokens.length > 0) {
try {
int archiveNumberOfAtoms = Integer.parseInt(tokens[0]);
if (archiveNumberOfAtoms == numberOfAtoms) {
setType(FileType.ARC);
}
} catch (NumberFormatException e) {
tokens = null;
}
}
}
}
br.close();
fr.close();
// Try to renumber
if (renumber) {
for (int i = 0; i < numberOfAtoms; i++) {
if (labelHash.containsKey(label[i])) {
logger.warning(format(" Two atoms have the same index: %d.", label[i]));
return false;
}
labelHash.put(label[i], i + 1);
}
for (int i = 0; i < numberOfAtoms; i++) {
int j = -1;
while (j < 3 && bonds[i][++j] > 0) {
bonds[i][j] = labelHash.get(bonds[i][j]);
}
}
}
bondList = new ArrayList<>();
int[] c = new int[2];
for (int i = 1; i <= numberOfAtoms; i++) {
int a1 = i;
int j = -1;
while (j < 7 && bonds[i - 1][++j] > 0) {
int a2 = bonds[i - 1][j];
if (a1 < a2) {
if (a1 > numberOfAtoms || a1 < 1 || a2 > numberOfAtoms || a2 < 1) {
logger.warning(format(" Check the bond between %d and %d in %s.", a1, a2, activeMolecularAssembly.getFile().getName()));
return false;
}
// Check for bidirectional connection
boolean bidirectional = false;
int k = -1;
while (k < 7 && bonds[a2 - 1][++k] > 0) {
int a3 = bonds[a2 - 1][k];
if (a3 == a1) {
bidirectional = true;
break;
}
}
if (!bidirectional) {
logger.warning(format(" Check the bond between %d and %d in %s.", a1, a2, activeMolecularAssembly.getFile().getName()));
return false;
}
Atom atom1 = atomList.get(a1 - 1);
Atom atom2 = atomList.get(a2 - 1);
if (atom1 == null || atom2 == null) {
logger.warning(format(" Check the bond between %d and %d in %s.", a1, a2, activeMolecularAssembly.getFile().getName()));
return false;
}
Bond bond = new Bond(atom1, atom2);
c[0] = atom1.getAtomType().atomClass;
c[1] = atom2.getAtomType().atomClass;
String key = BondType.sortKey(c);
BondType bondType = forceField.getBondType(key);
if (bondType == null) {
logger.severe(format(" No BondType for key %s", key));
} else {
bond.setBondType(bondType);
}
bondList.add(bond);
}
}
}
return true;
} catch (IOException e) {
logger.severe(e.toString());
}
return false;
}
use of ffx.potential.parameters.BondType in project ffx by mjschnie.
the class ForceFieldFilter method parseBond.
private void parseBond(String input, String[] tokens) {
if (tokens.length < 5) {
logger.log(Level.WARNING, "Invalid BOND type:\n{0}", input);
return;
}
try {
int[] atomClasses = new int[2];
atomClasses[0] = Integer.parseInt(tokens[1]);
atomClasses[1] = Integer.parseInt(tokens[2]);
double forceConstant = Double.parseDouble(tokens[3]);
double distance = Double.parseDouble(tokens[4]);
String forceFieldName = forceField.toString().toUpperCase();
BondType bondType;
if (forceFieldName.contains("OPLS") || forceFieldName.contains("AMBER")) {
bondType = new BondType(atomClasses, forceConstant, distance, BondType.BondFunction.HARMONIC);
} else {
bondType = new BondType(atomClasses, forceConstant, distance, BondType.BondFunction.QUARTIC);
}
forceField.addForceFieldType(bondType);
} catch (NumberFormatException e) {
String message = "Exception parsing BOND type:\n" + input + "\n";
logger.log(Level.SEVERE, message, e);
}
}
Aggregations