use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class RotamerOptimizationTest method testPairEnergyElimination.
@Test
public void testPairEnergyElimination() {
// Load the test system.
load();
// Initialize Parallel Java
try {
String[] args = new String[0];
Comm.init(args);
} catch (Exception e) {
String message = String.format(" Exception starting up the Parallel Java communication layer.");
logger.log(Level.WARNING, message, e.toString());
message = String.format(" Skipping rotamer optimization test.");
logger.log(Level.WARNING, message, e.toString());
return;
}
// Run the optimization.
RotamerLibrary rLib = RotamerLibrary.getDefaultLibrary();
rLib.setLibrary(RotamerLibrary.ProteinLibrary.Richardson);
rLib.setUseOrigCoordsRotamer(useOriginalRotamers);
int counter = 1;
ArrayList<Residue> residueList = new ArrayList<Residue>();
Polymer[] polymers = molecularAssembly.getChains();
int nPolymers = polymers.length;
for (int p = 0; p < nPolymers; p++) {
Polymer polymer = polymers[p];
ArrayList<Residue> residues = polymer.getResidues();
for (int i = 0; i < endResID; i++) {
Residue residue = residues.get(i);
Rotamer[] rotamers = residue.getRotamers(rLib);
if (rotamers != null) {
int nrot = rotamers.length;
if (nrot == 1) {
RotamerLibrary.applyRotamer(residue, rotamers[0]);
}
if (counter >= startResID) {
residueList.add(residue);
}
}
counter++;
}
}
RotamerOptimization rotamerOptimization = new RotamerOptimization(molecularAssembly, forceFieldEnergy, null);
rotamerOptimization.setThreeBodyEnergy(useThreeBody);
rotamerOptimization.setUseGoldstein(useGoldstein);
rotamerOptimization.setPruning(pruningLevel);
rotamerOptimization.setEnergyRestartFile(restartFile);
rotamerOptimization.setResidues(residueList);
double energy;
int nRes = residueList.size();
if (doOverallOpt) {
rotamerOptimization.turnRotamerSingleEliminationOff();
energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
// System.out.println("The expected overall energy is: " + energy);
assertEquals(info + " Total Energy", expectedEnergy, energy, tolerance);
}
// ToDo: Test self-energy use for rotamer 2-body eliminations.
if (doSelfOpt) {
rotamerOptimization.turnRotamerSingleEliminationOff();
rotamerOptimization.setTestSelfEnergyEliminations(true);
energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
// System.out.println("The expected self energy is: " + energy);
assertEquals(info + " Self-Energy", expectedSelfEnergy, energy, tolerance);
// Check that optimized rotamers are equivalent to the lowest self-energy of each residue.
int[] optimum = rotamerOptimization.getOptimumRotamers();
// Loop over all residues
for (int i = 0; i < nRes; i++) {
Residue res = residueList.get(i);
Rotamer[] rotI = res.getRotamers(rLib);
int nRot = rotI.length;
int rotCounter = 0;
while (rotCounter < nRot && rotamerOptimization.checkPrunedSingles(i, rotCounter)) {
rotCounter++;
}
double lowEnergy = rotamerOptimization.getSelf(i, rotCounter);
int bestRot = rotCounter;
for (int ri = 1; ri < nRot; ri++) {
if (rotamerOptimization.checkPrunedSingles(i, ri)) {
continue;
} else {
double selfEnergy = rotamerOptimization.getSelf(i, ri);
if (selfEnergy < lowEnergy) {
lowEnergy = selfEnergy;
bestRot = ri;
}
}
}
assertEquals(String.format(" %s Self-Energy of residue %d", info, i), optimum[i], bestRot);
}
}
// ToDo: Test 2-body energy use for rotamer pair eliminations.
if (doPairOpt) {
rotamerOptimization.turnRotamerSingleEliminationOff();
rotamerOptimization.setTestPairEnergyEliminations(pairResidue);
energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
// System.out.println("The expected 2-body energy is: " + energy);
assertEquals(info + " Pair-Energy", expectedPairEnergy, energy, tolerance);
// Check that optimized rotamers are equivalent to the lowest 2-body energy sum for the "pairResidue".
int[] optimum = rotamerOptimization.getOptimumRotamers();
Residue resI = residueList.get(pairResidue);
Rotamer[] rotI = resI.getRotamers(rLib);
int ni = rotI.length;
double minEnergy = Double.POSITIVE_INFINITY;
int bestRotI = -1;
// Loop over the pairResidue rotamers to find its lowest energy rotamer.
for (int ri = 0; ri < ni; ri++) {
double energyForRi = 0.0;
if (rotamerOptimization.checkPrunedSingles(pairResidue, ri)) {
continue;
}
// Loop over residue J
for (int j = 0; j < nRes; j++) {
if (j == pairResidue) {
continue;
}
Residue resJ = residueList.get(j);
Rotamer[] rotJ = resJ.getRotamers(rLib);
int nRot = rotJ.length;
int rj = 0;
while (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, ri, j, rj)) {
if (++rj >= nRot) {
logger.warning("RJ is too large.");
}
}
double lowEnergy = rotamerOptimization.get2Body(pairResidue, ri, j, rj);
for (rj = 1; rj < nRot; rj++) {
if (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, ri, j, rj)) {
continue;
} else {
double pairEnergy = rotamerOptimization.get2Body(pairResidue, ri, j, rj);
if (pairEnergy < lowEnergy) {
lowEnergy = pairEnergy;
}
}
}
energyForRi += lowEnergy;
}
if (energyForRi < minEnergy) {
minEnergy = energyForRi;
bestRotI = ri;
}
}
assertEquals(String.format(" %s Best 2-body energy sum for residue %d is with rotamer %d at %10.4f.", info, pairResidue, bestRotI, minEnergy), optimum[pairResidue], bestRotI);
// Given the minimum energy rotamer for "pairResidue" is "bestRotI", we can check selected rotamers for all other residues.
for (int j = 0; j < nRes; j++) {
if (j == pairResidue) {
continue;
}
Residue resJ = residueList.get(j);
Rotamer[] rotJ = resJ.getRotamers(rLib);
int nRotJ = rotJ.length;
int rotCounter = 0;
while (rotamerOptimization.checkPrunedPairs(pairResidue, bestRotI, j, rotCounter) && rotCounter < nRotJ) {
rotCounter++;
}
double lowEnergy = rotamerOptimization.get2Body(pairResidue, bestRotI, j, rotCounter);
int bestRotJ = rotCounter;
for (int rj = 1; rj < nRotJ; rj++) {
if (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, bestRotI, j, rj)) {
continue;
} else {
double pairEnergy = rotamerOptimization.get2Body(pairResidue, bestRotI, j, rj);
if (pairEnergy < lowEnergy) {
lowEnergy = pairEnergy;
bestRotJ = rj;
}
}
}
if (bestRotJ != optimum[j]) {
// Check if 2-body energies are equal.
if (lowEnergy == rotamerOptimization.get2Body(pairResidue, bestRotI, j, optimum[j])) {
logger.warning(String.format(" Identical 2-body energies for %s: resi %d-%d, resj %d, best rotamer J %d, optimum J %d, 2-body energy (both) %10.6f", info, pairResidue, bestRotI, j, bestRotJ, optimum[j], lowEnergy));
} else {
assertEquals(String.format(" %s Pair-Energy of residue (%d,%d) with residue %d", info, pairResidue, bestRotI, j), optimum[j], bestRotJ);
}
}
}
}
// ToDo: Test 3-Body use for rotamer pair eliminations.
if (doTripleOpt) {
rotamerOptimization.turnRotamerSingleEliminationOff();
rotamerOptimization.setTestTripleEnergyEliminations(tripleResidue1, tripleResidue2);
try {
energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
assertEquals(info + " Triple-Energy", expectedTripleEnergy, energy, tolerance);
} catch (Exception e) {
e.fillInStackTrace();
e.printStackTrace();
logger.log(java.util.logging.Level.INFO, "Error in doTripleOpt", e);
}
// Check that optimized rotamers are equivalent to the lowest 3-body energy of each residue with the tripleResidue1 and 2.
int[] optimum = rotamerOptimization.getOptimumRotamers();
// fix residue 1 and gets its rotamers
Residue resI = residueList.get(tripleResidue1);
Rotamer[] rotI = resI.getRotamers(rLib);
int ni = rotI.length;
// fix residue 2 and get its rotamers
Residue resJ = residueList.get(tripleResidue2);
Rotamer[] rotJ = resJ.getRotamers(rLib);
int nj = rotJ.length;
double minEnergyIJ = Double.POSITIVE_INFINITY;
int bestRotI = -1;
int bestRotJ = -1;
for (int ri = 0; ri < ni; ri++) {
// loop through rot I
if (rotamerOptimization.check(tripleResidue1, ri)) {
continue;
}
for (int rj = 0; rj < nj; rj++) {
// loop through rot J
if (rotamerOptimization.checkPrunedSingles(tripleResidue2, rj) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, tripleResidue2, rj)) {
continue;
}
double currentEnergy = 0.0;
for (int k = 0; k < nRes; k++) {
// loop through all other residues
if (k == tripleResidue1 || k == tripleResidue2) {
continue;
}
Residue resK = residueList.get(k);
Rotamer[] rotK = resK.getRotamers(rLib);
int nk = rotK.length;
int rkStart = 0;
while (rotamerOptimization.checkPrunedSingles(k, rkStart) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, k, rkStart) || rotamerOptimization.checkPrunedPairs(tripleResidue2, rj, k, rkStart)) {
if (++rkStart >= nk) {
logger.warning("RJ is too large.");
}
}
double lowEnergy = rotamerOptimization.get3Body(tripleResidue1, ri, tripleResidue2, rj, k, rkStart);
for (int rk = rkStart; rk < nk; rk++) {
if (rotamerOptimization.checkPrunedSingles(k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue2, rj, k, rk)) {
continue;
} else {
double tripleEnergy = rotamerOptimization.get3Body(tripleResidue1, ri, tripleResidue2, rj, k, rk);
if (tripleEnergy < lowEnergy) {
lowEnergy = tripleEnergy;
}
}
}
// adds lowest energy conformation of residue k to that of the rotamer I
currentEnergy += lowEnergy;
}
if (currentEnergy < minEnergyIJ) {
minEnergyIJ = currentEnergy;
bestRotI = ri;
bestRotJ = rj;
}
}
}
assertEquals(String.format(" %s Best three-body energy sum for residue %d is with rotamer %d at %10.4f.", info, tripleResidue1, bestRotI, minEnergyIJ), optimum[tripleResidue1], bestRotI);
assertEquals(String.format(" %s Best three-body energy sum for residue %d is with rotamer %d at %10.4f.", info, tripleResidue2, bestRotJ, minEnergyIJ), optimum[tripleResidue2], bestRotJ);
// loop over the residues to find the best rotamer per residue given bestRotI and bestRotJ
for (int k = 0; k < nRes; k++) {
if (k == tripleResidue1 || k == tripleResidue2) {
continue;
}
Residue resK = residueList.get(k);
Rotamer[] rotK = resK.getRotamers(rLib);
int nk = rotK.length;
int rotCounter = 0;
while (rotamerOptimization.checkPrunedPairs(tripleResidue1, bestRotI, k, rotCounter) && rotamerOptimization.checkPrunedPairs(tripleResidue2, bestRotJ, k, rotCounter) && rotCounter < nk) {
rotCounter++;
}
double lowEnergy = rotamerOptimization.get3Body(tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k, rotCounter);
int bestRotK = rotCounter;
for (int rk = 1; rk < nk; rk++) {
if (rotamerOptimization.checkPrunedSingles(k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue1, bestRotI, k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue2, bestRotJ, k, rk)) {
continue;
} else {
double tripleEnergy = rotamerOptimization.get3Body(tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k, rk);
if (tripleEnergy < lowEnergy) {
lowEnergy = tripleEnergy;
bestRotK = rk;
}
}
}
assertEquals(String.format(" %s Triple-Energy of residue (%d,%d) and residue (%d,%d) with residue %d", info, tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k), optimum[k], bestRotK);
}
}
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class ExtendedSystem method populate.
public void populate(List<String> residueIDs) {
// Locate the Residue identified by the given resid.
Polymer[] polymers = mola.getChains();
for (String token : residueIDs) {
char chainID = token.charAt(0);
int resNum = Integer.parseInt(token.substring(1));
Residue target = null;
for (Polymer p : polymers) {
char pid = p.getChainID().charValue();
if (pid == chainID) {
for (Residue res : p.getResidues()) {
if (res.getResidueNumber() == resNum) {
target = res;
break;
}
}
if (target != null) {
break;
}
}
}
if (target == null) {
logger.severe("Couldn't find target residue " + token);
}
MultiResidue titrating = TitrationUtils.titratingMultiresidueFactory(mola, target);
TitrationESV esv = new TitrationESV(this, titrating);
this.addVariable(esv);
}
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class TitrationUtils method chooseTitratables.
/**
* Identify titratable residues and choose them all.
*/
public static List<Residue> chooseTitratables(MolecularAssembly searchMe) {
List<Residue> chosen = new ArrayList<>();
Polymer[] polymers = searchMe.getChains();
for (int i = 0; i < polymers.length; i++) {
ArrayList<Residue> residues = polymers[i].getResidues();
for (int j = 0; j < residues.size(); j++) {
Residue res = residues.get(j);
Titration[] avail = Titration.multiLookup(res);
if (avail != null) {
chosen.add(residues.get(j));
}
}
}
return chosen;
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class TitrationUtils method propagateInactiveResidues.
/**
* Copies atomic coordinates from each active residue to its inactive
* counterparts. Inactive hydrogen coordinates are updated by geometry with
* the propagated heavies.
*/
public static void propagateInactiveResidues(MultiResidue multiRes, boolean propagateDynamics) {
// Propagate all atom coordinates from active residues to their inactive counterparts.
Residue active = multiRes.getActive();
String activeResName = active.getName();
List<Residue> inactives = multiRes.getInactive();
for (Atom activeAtom : active.getAtomList()) {
String activeName = activeAtom.getName();
for (Residue inactive : inactives) {
Atom inactiveAtom = (Atom) inactive.getAtomNode(activeName);
if (inactiveAtom != null) {
// Propagate position and gradient.
double[] activeXYZ = activeAtom.getXYZ(null);
inactiveAtom.setXYZ(activeXYZ);
double[] grad = new double[3];
activeAtom.getXYZGradient(grad);
inactiveAtom.setXYZGradient(grad[0], grad[1], grad[2]);
if (propagateDynamics) {
// Propagate velocity, acceleration, and previous acceleration.
double[] activeVelocity = new double[3];
activeAtom.getVelocity(activeVelocity);
inactiveAtom.setVelocity(activeVelocity);
double[] activeAccel = new double[3];
activeAtom.getAcceleration(activeAccel);
inactiveAtom.setAcceleration(activeAccel);
double[] activePrevAcc = new double[3];
activeAtom.getPreviousAcceleration(activePrevAcc);
inactiveAtom.setPreviousAcceleration(activePrevAcc);
}
} else {
if (activeName.equals("C") || activeName.equals("O") || activeName.equals("N") || activeName.equals("CA") || activeName.equals("H") || activeName.equals("HA")) {
// Backbone atoms aren't supposed to exist in inactive multiResidue components; so no problem.
} else if (isTitratableHydrogen(activeAtom)) {
/**
* i.e. ((activeResName.equals("LYS") &&
* activeName.equals("HZ3")) ||
* (activeResName.equals("TYR") &&
* activeName.equals("HH")) ||
* (activeResName.equals("CYS") &&
* activeName.equals("HG")) ||
* (activeResName.equals("HIS") &&
* (activeName.equals("HD1") ||
* activeName.equals("HE2"))) ||
* (activeResName.equals("HID") &&
* activeName.equals("HD1")) ||
* (activeResName.equals("HIE") &&
* activeName.equals("HE2")) ||
* (activeResName.equals("ASH") &&
* activeName.equals("HD2")) ||
* (activeResName.equals("GLH") &&
* activeName.equals("HE2")))
*/
// These titratable protons are handled below; so no problem.
} else {
// Now we have a problem.
logger.warning(format("Couldn't propagate inactive MultiResidue atom: %s: %s, %s", multiRes, activeName, activeAtom));
}
}
}
}
rebuildStrandedProtons(multiRes);
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class TitrationUtils method chooseTitratables.
public static List<Residue> chooseTitratables(List<String> crIDs, MolecularAssembly searchMe) {
List<Residue> chosen = new ArrayList<>();
for (String crID : crIDs) {
char chain = crID.charAt(0);
int num = Integer.parseInt(crID.substring(1));
boolean found = false;
List<Residue> allRes = searchMe.getResidueList();
for (Residue res : allRes) {
if (res.getChainID() == chain && res.getResidueNumber() == num) {
chosen.add(res);
found = true;
break;
}
}
if (!found) {
logger.severe(format("Couldn't find residue for crID %c,%d.", chain, num));
}
}
return chosen;
}
Aggregations