Example 26 with Rotamer
use of ffx.potential.bonded.Rotamer in project ffx by mjschnie.
the class RosenbluthChiAllMove method engage_controlAll.
/**
* For validation. Performs Monte Carlo chi moves WITHOUT biasing. Give ALL
* CHIs a random theta simultaneously. Accept on the vanilla Metropolis
* criterion.
*/
private boolean engage_controlAll() {
report.append(String.format(" Rosenbluth Control Move: %4d %s\n", moveNumber, target));
double origEnergy = totalEnergy();
double[] origChi = RotamerLibrary.measureRotamer(target, false);
double[] newChi = new double[origChi.length];
System.arraycopy(origChi, 0, newChi, 0, origChi.length);
for (int i = 0; i < origChi.length; i++) {
if (doChi[i]) {
double theta = rand.nextDouble(360.0) - 180;
newChi[i] = theta;
}
}
proposedChis = newChi;
Rotamer newState = createRotamer(target, newChi);
RotamerLibrary.applyRotamer(target, newState);
finalEnergy = totalEnergy();
if (this.finalEnergy < CATASTROPHE_THRESHOLD) {
report.append("\nWARNING: Multipole catastrophe in CBMC.\n");
report.append(" Discarding move.\n");
target.revertState(origState);
updateAll();
Wn = -1.0;
Wo = 1000;
logger.info(report.toString());
return false;
}
double dU = finalEnergy - origEnergy;
double criterion = FastMath.exp(-beta * dU);
double rng = rand.nextDouble();
report.append(String.format(" move (thetas): "));
for (int i = 0; i < newChi.length; i++) {
report.append(String.format("%7.2f ", newChi[i]));
}
report.append(String.format("\n"));
report.append(String.format(" orig, final, dU: %.2g %.2g %.2g\n", origEnergy, finalEnergy, dU));
report.append(String.format(" crit, rng: %.2g %.2g\n", criterion, rng));
if (rng < criterion) {
accepted = true;
numAccepted++;
report.append(String.format(" Accepted! %5d NewEnergy: %.4f Chi:", numAccepted, finalEnergy));
for (int k = 0; k < proposedChis.length; k++) {
report.append(String.format(" %7.2f", proposedChis[k]));
}
report.append(String.format("\n"));
updateAll();
if (!noSnaps) {
PDBFilter writer = new PDBFilter(mola.getFile(), mola, null, null);
String filename = FilenameUtils.removeExtension(mola.getFile().toString());
filename = mola.getFile().getAbsolutePath();
if (!filename.contains("_mc")) {
filename = FilenameUtils.removeExtension(filename) + "_mc.pdb";
}
File file = new File(filename);
writer.writeFile(file, false);
}
} else {
accepted = false;
report.append(String.format(" Denied. %5d NewEnergy: %.4f Chi:", numAccepted, origEnergy));
for (int k = 0; k < origChi.length; k++) {
report.append(String.format(" %7.2f", origChi[k]));
}
report.append(String.format("\n"));
target.revertState(origState);
}
updateAll();
endTime = System.nanoTime();
double took = (endTime - startTime) * NS_TO_MS;
if (logTimings) {
report.append(String.format(" Timing (ms): %.2f", took));
}
logger.info(report.toString());
return accepted;
}Example 27 with Rotamer
use of ffx.potential.bonded.Rotamer in project ffx by mjschnie.
the class RotamerOptimizationTest method testSelfEnergyElimination.
@Test
public void testSelfEnergyElimination() {
// 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.turnRotamerPairEliminationOff();
rotamerOptimization.setTestOverallOpt(true);
energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
// System.out.println("The expected overall energy is: " + energy);
assertEquals(info + " Total Energy", expectedEnergy, energy, tolerance);
}
if (doSelfOpt) {
rotamerOptimization.turnRotamerPairEliminationOff();
rotamerOptimization.setTestSelfEnergyEliminations(true);
energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
// System.out.println("The expected self 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);
}
}
if (doPairOpt) {
rotamerOptimization.turnRotamerPairEliminationOff();
rotamerOptimization.setTestPairEnergyEliminations(pairResidue);
energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
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;
}
}
}
assertEquals(String.format(" %s Pair-Energy of residue (%d,%d) with residue %d", info, pairResidue, bestRotI, j), optimum[j], bestRotJ);
}
}
// Test 3-Body Energy Eliminations.
if (doTripleOpt) {
rotamerOptimization.turnRotamerPairEliminationOff();
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);
}
}
}Example 28 with Rotamer
use of ffx.potential.bonded.Rotamer 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);
}
}
}Example 29 with Rotamer
use of ffx.potential.bonded.Rotamer in project ffx by mjschnie.
the class RotamerOptimization method slidingWindowOptimization.
private double slidingWindowOptimization(ArrayList<Residue> residueList, int windowSize, int increment, boolean revert, double distance, Direction direction) {
long beginTime = -System.nanoTime();
boolean incrementTruncated = false;
boolean firstWindowSaved = false;
int counter = 1;
int windowEnd;
int nOptimize = residueList.size();
if (nOptimize < windowSize) {
windowSize = nOptimize;
logger.info(" Window size too small for given residue range; truncating window size.");
}
switch(direction) {
case BACKWARD:
ArrayList<Residue> temp = new ArrayList<>();
for (int i = nOptimize - 1; i >= 0; i--) {
temp.add(residueList.get(i));
}
residueList = temp;
// Fall through into the FORWARD case.
case FORWARD:
for (int windowStart = 0; windowStart + (windowSize - 1) < nOptimize; windowStart += increment) {
long windowTime = -System.nanoTime();
windowEnd = windowStart + (windowSize - 1);
logIfMaster(format("\n Iteration %d of the sliding window.\n", counter++));
Residue firstResidue = residueList.get(windowStart);
Residue lastResidue = residueList.get(windowEnd);
if (firstResidue != lastResidue) {
logIfMaster(format(" Residues %s ... %s", firstResidue.toString(), lastResidue.toString()));
} else {
logIfMaster(format(" Residue %s", firstResidue.toString()));
}
ArrayList<Residue> currentWindow = new ArrayList<>();
// Not filled if useForcedResidues == false.
ArrayList<Residue> onlyRotameric = new ArrayList<>();
for (int i = windowStart; i <= windowEnd; i++) {
Residue residue = residueList.get(i);
if (useForcedResidues && residue.getRotamers(library) != null) {
onlyRotameric.add(residue);
}
currentWindow.add(residueList.get(i));
}
if (distance > 0) {
for (int i = windowStart; i <= windowEnd; i++) {
Residue residuei = residueList.get(i);
int indexI = allResiduesList.indexOf(residuei);
int lengthRi;
if (checkIfForced(residuei)) {
lengthRi = 1;
} else {
lengthRi = residuei.getRotamers(library).length;
}
for (int ri = 0; ri < lengthRi; ri++) {
for (int j = 0; j < numResidues; j++) {
Residue residuej = allResiduesArray[j];
Rotamer[] rotamersj = residuej.getRotamers(library);
if (currentWindow.contains(residuej) || rotamersj == null) {
continue;
}
int lengthRj = rotamersj.length;
for (int rj = 0; rj < lengthRj; rj++) {
double rotamerSeparation = get2BodyDistance(indexI, ri, j, rj);
// if (distanceMatrix[indexI][ri][j][rj] <= distance) {
if (rotamerSeparation <= distance) {
if (!currentWindow.contains(residuej)) {
logIfMaster(format(" Adding residue %s at distance %16.8f Ang from %s %d.", residuej.toFormattedString(false, true), rotamerSeparation, residuei.toFormattedString(false, true), ri));
currentWindow.add(residuej);
if (useForcedResidues) {
onlyRotameric.add(residuej);
}
}
break;
}
}
}
}
}
}
/**
* If the window starts with a nucleic acid, and there is a
* 5' NA residue, ensure that 5' NA residue has been
* included in the window. Otherwise, that previous residue
* may not have had a chance to be flexible about its sugar
* pucker.
*
* If window size is greater than increment, however, this
* has already been handled. Additionally, do not perform
* this for the first window (counter is already incremented
* by the time this check is performed, so first window's
* counter will be 2). Furthermore, do not include Residues
* with null Rotamer lists (this breaks things).
*
* The issue: if window size = increment, the last NA
* residue in each window will not have flexibility about
* its sugar pucker, because its self-energy includes the
* O3' (i) to P (i+1) bond, so it must remain in the
* original sugar pucker to meet the i+1 residue. However,
* this problem can be solved by ensuring that final residue
* is included in the next window, where it will have
* flexibility about its sugar pucker.
*
* If you are running successive sliding window jobs on the
* same file, I would suggest starting the next job on the
* last residue of the previous job, unless you know your
* settings will include it.
*/
if (counter > 2 && windowSize <= increment && firstResidue.getResidueType() == NA) {
Residue prevResidue = firstResidue.getPreviousResidue();
if (prevResidue != null && prevResidue.getResidueType() == NA && !currentWindow.contains(prevResidue) && prevResidue.getRotamers(library) != null) {
logIfMaster(format(" Adding nucleic acid residue 5' of window start %s to give it flexibility about its sugar pucker.", prevResidue.toString()));
currentWindow.add(prevResidue);
if (useForcedResidues) {
onlyRotameric.add(prevResidue);
}
}
}
if (useForcedResidues) {
sortResidues(onlyRotameric);
} else {
sortResidues(currentWindow);
}
if (revert) {
ResidueState[] coordinates = ResidueState.storeAllCoordinates(currentWindow);
// If x has not yet been constructed, construct it.
if (x == null) {
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
x = new double[nAtoms * 3];
}
double startingEnergy = Double.NaN;
try {
startingEnergy = currentEnergy(currentWindow);
} catch (ArithmeticException ex) {
logger.severe(String.format(" Exception %s in calculating starting energy of a window; FFX shutting down", ex.toString()));
}
if (useForcedResidues) {
if (onlyRotameric.size() < 1) {
logger.info(" Window has no rotameric residues.");
ResidueState.revertAllCoordinates(currentWindow, coordinates);
} else {
globalOptimization(onlyRotameric);
double finalEnergy = Double.NaN;
try {
finalEnergy = currentEnergy(currentWindow);
} catch (ArithmeticException ex) {
logger.severe(String.format(" Exception %s in calculating final energy of a window; FFX shutting down", ex.toString()));
}
if (startingEnergy <= finalEnergy) {
logger.warning("Optimization did not yield a better energy. Reverting to orginal coordinates.");
ResidueState.revertAllCoordinates(currentWindow, coordinates);
}
}
} else {
globalOptimization(currentWindow);
double finalEnergy = Double.NaN;
try {
finalEnergy = currentEnergy(currentWindow);
} catch (ArithmeticException ex) {
logger.severe(String.format(" Exception %s in calculating final energy of a window; FFX shutting down", ex.toString()));
}
if (startingEnergy <= finalEnergy) {
logger.warning("Optimization did not yield a better energy. Reverting to orginal coordinates.");
ResidueState.revertAllCoordinates(currentWindow, coordinates);
}
}
} else if (useForcedResidues) {
if (onlyRotameric.size() < 1) {
logger.info(" Window has no rotameric residues.");
} else {
globalOptimization(onlyRotameric);
}
} else {
globalOptimization(currentWindow);
}
if (!incrementTruncated) {
if (windowStart + (windowSize - 1) + increment > nOptimize - 1) {
increment = nOptimize - windowStart - windowSize;
if (increment == 0) {
break;
}
logger.warning(" Increment truncated in order to optimize entire residue range.");
incrementTruncated = true;
}
}
if (master && printFiles) {
File file = molecularAssembly.getFile();
if (firstWindowSaved) {
file.delete();
}
// Don't write a file if its the final iteration.
if (windowStart + windowSize == nOptimize) {
continue;
}
// String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
// File tempFile = new File(filename + ".win");
// PDBFilter windowFilter = new PDBFilter(new File(filename + ".win"), molecularAssembly, null, null);
PDBFilter windowFilter = new PDBFilter(file, molecularAssembly, null, null);
// StringBuilder header = new StringBuilder(format("Iteration %d of the sliding window\n", counter - 1));
try {
windowFilter.writeFile(file, false);
if (firstResidue != lastResidue) {
logger.info(format(" File with residues %s ... %s in window written to.", firstResidue.toString(), lastResidue.toString()));
} else {
logger.info(format(" File with residue %s in window written to.", firstResidue.toString()));
}
} catch (Exception e) {
logger.warning(format("Exception writing to file: %s", file.getName()));
}
firstWindowSaved = true;
}
long currentTime = System.nanoTime();
windowTime += currentTime;
logIfMaster(format(" Time elapsed for this iteration: %11.3f sec", windowTime * 1.0E-9));
logIfMaster(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
/*for (Residue residue : residueList) {
if (residue instanceof MultiResidue) {
((MultiResidue) residue).setDefaultResidue();
residue.reInitOriginalAtomList();
}
}*/
}
break;
default:
// No default case.
break;
}
return 0.0;
}
Also used :
ResidueState(ffx.potential.bonded.ResidueState)
ArrayList(java.util.ArrayList)
RotamerLibrary.applyRotamer(ffx.potential.bonded.RotamerLibrary.applyRotamer)
Rotamer(ffx.potential.bonded.Rotamer)
Atom(ffx.potential.bonded.Atom)
IOException(java.io.IOException)
NACorrectionException(ffx.potential.bonded.NACorrectionException)
Residue(ffx.potential.bonded.Residue)
MultiResidue(ffx.potential.bonded.MultiResidue)
File(java.io.File)
PDBFilter(ffx.potential.parsers.PDBFilter)
Example 30 with Rotamer
use of ffx.potential.bonded.Rotamer in project ffx by mjschnie.
the class RotamerOptimization method loadEnergyRestart.
private int loadEnergyRestart(File restartFile, Residue[] residues, int boxIteration, int[] cellIndices) {
try {
int nResidues = residues.length;
Path path = Paths.get(restartFile.getCanonicalPath());
List<String> lines = Files.readAllLines(path, StandardCharsets.UTF_8);
List<String> linesThisBox = new ArrayList<>();
try {
backboneEnergy = computeBackboneEnergy(residues);
} catch (ArithmeticException ex) {
logger.severe(String.format(" Exception %s in calculating backbone energy; FFX shutting down.", ex.toString()));
}
if (usingBoxOptimization && boxIteration >= 0) {
boolean foundBox = false;
for (int i = 0; i < lines.size(); i++) {
String line = lines.get(i);
if (line.startsWith("Box")) {
String[] tok = line.replaceAll("Box", "").replaceAll(":", ",").replaceAll(" ", "").split(",");
int readIteration = Integer.parseInt(tok[0]);
int readCellIndexX = Integer.parseInt(tok[1]);
int readCellIndexY = Integer.parseInt(tok[2]);
int readCellIndexZ = Integer.parseInt(tok[3]);
if (readIteration == boxIteration && readCellIndexX == cellIndices[0] && readCellIndexY == cellIndices[1] && readCellIndexZ == cellIndices[2]) {
foundBox = true;
for (int j = i + 1; j < lines.size(); j++) {
String l = lines.get(j);
if (l.startsWith("Box")) {
break;
}
linesThisBox.add(l);
}
break;
}
}
}
if (!foundBox) {
logIfMaster(format(" Didn't find restart energies for Box %d: %d,%d,%d", boxIteration, cellIndices[0], cellIndices[1], cellIndices[2]));
return 0;
} else if (linesThisBox.size() == 0) {
return 0;
} else {
lines = linesThisBox;
}
}
List<String> singleLines = new ArrayList<>();
List<String> pairLines = new ArrayList<>();
List<String> tripleLines = new ArrayList<>();
for (String line : lines) {
String[] tok = line.split("\\s");
if (tok[0].startsWith("Self")) {
singleLines.add(line);
} else if (tok[0].startsWith("Pair")) {
pairLines.add(line);
} else if (tok[0].startsWith("Triple")) {
tripleLines.add(line);
}
}
int loaded = 0;
if (tripleLines.size() > 0) {
loaded = 3;
} else if (pairLines.size() > 0) {
loaded = 2;
} else if (singleLines.size() > 0) {
loaded = 1;
} else {
logger.warning(format(" Empty or unreadable energy restart file: %s.", restartFile.getCanonicalPath()));
}
if (loaded >= 1) {
selfEnergyMap.clear();
// allocate selfEnergy array and create self jobs
HashMap<String, Integer> reverseJobMapSingles = new HashMap<>();
int singleJobIndex = 0;
selfEnergy = new double[nResidues][];
for (int i = 0; i < nResidues; i++) {
Residue resi = residues[i];
Rotamer[] roti = resi.getRotamers(library);
selfEnergy[i] = new double[roti.length];
for (int ri = 0; ri < roti.length; ri++) {
Integer[] selfJob = { i, ri };
if (decomposeOriginal && ri != 0) {
continue;
}
selfEnergyMap.put(singleJobIndex, selfJob);
String revKey = format("%d %d", i, ri);
reverseJobMapSingles.put(revKey, singleJobIndex);
singleJobIndex++;
}
}
// fill in self-energies from file while removing the corresponding jobs from selfEnergyMap
for (String line : singleLines) {
try {
String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
int i;
if (tok[1].contains("-")) {
i = nameToNumber(tok[1], residues);
} else {
i = Integer.parseInt(tok[1]);
}
int ri = Integer.parseInt(tok[2]);
double energy = Double.parseDouble(tok[3]);
try {
selfEnergy[i][ri] = energy;
if (verbose) {
logIfMaster(format(" From restart file: Self energy %3d (%8s,%2d): %s", i, residues[i].toFormattedString(false, true), ri, formatEnergy(energy)));
}
} catch (Exception e) {
if (verbose) {
logIfMaster(format(" Restart file out-of-bounds index: %s", line));
}
}
// remove that job from the pool
String revKey = format("%d %d", i, ri);
Integer[] ret = selfEnergyMap.remove(reverseJobMapSingles.get(revKey));
if (ret == null) {
// logIfMaster(format("(sdl %d) Restart file contained unnecessary value for %s", BOXNUM, revKey));
}
} catch (NumberFormatException ex) {
logger.log(Level.WARNING, format(" Unparsable line in energy restart file: \n%s", line), ex);
}
}
logIfMaster(" Loaded self energies from restart file.");
// Pre-Prune if self-energy is Double.NaN.
for (int i = 0; i < residues.length; i++) {
Residue residue = residues[i];
Rotamer[] rotamers = residue.getRotamers(library);
int nrot = rotamers.length;
for (int ri = 0; ri < nrot; ri++) {
if (!check(i, ri) && Double.isNaN(getSelf(i, ri))) {
logIfMaster(format(" Rotamer (%7s,%2d) self-energy %12.4f pre-pruned since energy is NaN.", residue, ri, getSelf(i, ri)));
eliminateRotamer(residues, i, ri, false);
}
}
}
// prune singles
if (pruneClashes) {
pruneSingleClashes(residues);
}
}
/**
* Remap to sequential integer keys.
*/
condenseEnergyMap(selfEnergyMap);
if (loaded >= 2) {
if (selfEnergyMap.size() > 0) {
logIfMaster(" Double-check that parameters match original run due to missing self-energies.");
}
twoBodyEnergyMap.clear();
// allocated twoBodyEnergy array and create pair jobs
HashMap<String, Integer> reverseJobMapPairs = new HashMap<>();
int pairJobIndex = 0;
twoBodyEnergy = new double[nResidues][][][];
for (int i = 0; i < nResidues; i++) {
Residue resi = residues[i];
Rotamer[] roti = resi.getRotamers(library);
twoBodyEnergy[i] = new double[roti.length][][];
for (int ri = 0; ri < roti.length; ri++) {
if (check(i, ri)) {
continue;
}
twoBodyEnergy[i][ri] = new double[nResidues][];
for (int j = i + 1; j < nResidues; j++) {
Residue resj = residues[j];
Rotamer[] rotj = resj.getRotamers(library);
twoBodyEnergy[i][ri][j] = new double[rotj.length];
for (int rj = 0; rj < rotj.length; rj++) {
if (check(j, rj) || check(i, ri, j, rj)) {
continue;
}
Integer[] pairJob = { i, ri, j, rj };
if (decomposeOriginal && (ri != 0 || rj != 0)) {
continue;
}
twoBodyEnergyMap.put(pairJobIndex, pairJob);
String revKey = format("%d %d %d %d", i, ri, j, rj);
reverseJobMapPairs.put(revKey, pairJobIndex);
pairJobIndex++;
}
}
}
}
// fill in pair-energies from file while removing the corresponding jobs from twoBodyEnergyMap
for (String line : pairLines) {
try {
String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
int i;
if (tok[1].contains("-")) {
i = nameToNumber(tok[1], residues);
} else {
i = Integer.parseInt(tok[1]);
}
int ri = Integer.parseInt(tok[2]);
int j;
if (tok[3].contains("-")) {
j = nameToNumber(tok[3], residues);
} else {
j = Integer.parseInt(tok[3]);
}
int rj = Integer.parseInt(tok[4]);
double energy = Double.parseDouble(tok[5]);
try {
twoBodyEnergy[i][ri][j][rj] = energy;
if (verbose) {
logIfMaster(format(" From restart file: Pair energy [(%8s,%2d),(%8s,%2d)]: %12.4f", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, energy));
}
} catch (Exception e) {
if (verbose) {
logIfMaster(format(" Restart file out-of-bounds index: %s", line));
}
}
// remove that job from the pool
String revKey = format("%d %d %d %d", i, ri, j, rj);
Integer[] ret = twoBodyEnergyMap.remove(reverseJobMapPairs.get(revKey));
} catch (NumberFormatException ex) {
logger.log(Level.WARNING, format("Unparsable line in energy restart file: \n%s", line), ex);
}
}
logIfMaster(" Loaded 2-body energies from restart file.");
// Loop over first residue.
for (int i = 0; i < nResidues - 1; i++) {
Residue resi = residues[i];
Rotamer[] roti = resi.getRotamers(library);
int ni = roti.length;
// Loop over second residue.
for (int j = i + 1; j < nResidues; j++) {
Residue resj = residues[j];
Rotamer[] rotj = resj.getRotamers(library);
int nj = rotj.length;
// Loop over the rotamers for residue i.
for (int ri = 0; ri < ni; ri++) {
if (!validRotamer(residues, i, ri)) {
continue;
}
// Loop over rotamers for residue j.
for (int rj = 0; rj < nj; rj++) {
if (!validRotamer(residues, j, rj) || check(i, ri, j, rj)) {
continue;
}
if (!check(i, ri, j, rj) && Double.isNaN(get2Body(i, ri, j, rj))) {
logIfMaster(format(" Rotamer Pair (%7s,%2d) (%7s,%2d) 2-body energy %12.4f pre-pruned since energy is NaN.", i, ri, j, rj, get2Body(i, ri, j, rj)));
eliminateRotamerPair(residues, i, ri, j, rj, print);
}
}
}
}
}
// prune pairs
if (prunePairClashes) {
prunePairClashes(residues);
}
}
/**
* Remap to sequential integer keys.
*/
condenseEnergyMap(twoBodyEnergyMap);
if (loaded >= 3) {
if (twoBodyEnergyMap.size() > 0) {
if (master) {
logger.warning("Double-check that parameters match original run! Found trimers in restart file, but pairs job queue is non-empty.");
}
}
HashMap<String, Integer> reverseJobMapTrimers = new HashMap<>();
threeBodyEnergyMap.clear();
// allocate 3-Body energy array, fill in 3-Body energies from the restart file.
int trimerJobIndex = 0;
threeBodyEnergy = new double[nResidues][][][][][];
for (int i = 0; i < nResidues; i++) {
Residue resi = residues[i];
Rotamer[] roti = resi.getRotamers(library);
threeBodyEnergy[i] = new double[roti.length][][][][];
for (int ri = 0; ri < roti.length; ri++) {
if (check(i, ri)) {
continue;
}
threeBodyEnergy[i][ri] = new double[nResidues][][][];
for (int j = i + 1; j < nResidues; j++) {
Residue resj = residues[j];
Rotamer[] rotj = resj.getRotamers(library);
threeBodyEnergy[i][ri][j] = new double[rotj.length][][];
for (int rj = 0; rj < rotj.length; rj++) {
if (check(j, rj) || check(i, ri, j, rj)) {
continue;
}
threeBodyEnergy[i][ri][j][rj] = new double[nResidues][];
for (int k = j + 1; k < nResidues; k++) {
Residue resk = residues[k];
Rotamer[] rotk = resk.getRotamers(library);
threeBodyEnergy[i][ri][j][rj][k] = new double[rotk.length];
for (int rk = 0; rk < rotk.length; rk++) {
if (check(k, rk) || check(i, ri, k, rk) || check(j, rj, k, rk)) {
// Not implemented: check(i, ri, j, rj, k, rk).
continue;
}
Integer[] trimerJob = { i, ri, j, rj, k, rk };
if (decomposeOriginal && (ri != 0 || rj != 0 || rk != 0)) {
continue;
}
threeBodyEnergyMap.put(trimerJobIndex, trimerJob);
String revKey = format("%d %d %d %d %d %d", i, ri, j, rj, k, rk);
reverseJobMapTrimers.put(revKey, trimerJobIndex);
trimerJobIndex++;
}
}
}
}
}
}
// fill in 3-Body energies from file while removing the corresponding jobs from threeBodyEnergyMap
for (String line : tripleLines) {
try {
String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
int i;
if (tok[1].contains("-")) {
i = nameToNumber(tok[1], residues);
} else {
i = Integer.parseInt(tok[1]);
}
int ri = Integer.parseInt(tok[2]);
int j;
if (tok[3].contains("-")) {
j = nameToNumber(tok[3], residues);
} else {
j = Integer.parseInt(tok[3]);
}
int rj = Integer.parseInt(tok[4]);
int k;
if (tok[5].contains("-")) {
k = nameToNumber(tok[5], residues);
} else {
k = Integer.parseInt(tok[5]);
}
int rk = Integer.parseInt(tok[6]);
double energy = Double.parseDouble(tok[7]);
try {
threeBodyEnergy[i][ri][j][rj][k][rk] = energy;
} catch (ArrayIndexOutOfBoundsException ex) {
if (verbose) {
logIfMaster(format(" Restart file out-of-bounds index: %s", line));
}
} catch (NullPointerException npe) {
if (verbose) {
logIfMaster(format(" NPE in loading 3-body energies: pruning " + "likely changed! 3-body %s-%d %s-%d %s-%d", residues[i].toFormattedString(false, true), ri, residues[j], rj, residues[k], rk));
}
}
if (verbose) {
logIfMaster(format(" From restart file: Trimer energy %3d %-2d, %3d %-2d, %3d %-2d: %s", i, ri, j, rj, k, rk, formatEnergy(energy)));
}
// remove that job from the pool
String revKey = format("%d %d %d %d %d %d", i, ri, j, rj, k, rk);
Integer[] ret = threeBodyEnergyMap.remove(reverseJobMapTrimers.get(revKey));
if (ret == null) {
// logIfMaster(format("(sdl %d) Restart file contained unnecessary value for %s", BOXNUM, revKey));
}
} catch (NumberFormatException ex) {
logger.log(Level.WARNING, format("Unparsable line in energy restart file: \n%s", line), ex);
}
}
logIfMaster(" Loaded trimer energies from restart file.");
}
/**
* Remap to sequential integer keys.
*/
condenseEnergyMap(threeBodyEnergyMap);
return loaded;
} catch (IOException ex) {
logger.log(Level.WARNING, "Exception while loading energy restart file.", ex);
}
return 0;
}
Also used :
Path(java.nio.file.Path)
HashMap(java.util.HashMap)
ArrayList(java.util.ArrayList)
RotamerLibrary.applyRotamer(ffx.potential.bonded.RotamerLibrary.applyRotamer)
Rotamer(ffx.potential.bonded.Rotamer)
IOException(java.io.IOException)
IOException(java.io.IOException)
NACorrectionException(ffx.potential.bonded.NACorrectionException)
Residue(ffx.potential.bonded.Residue)
MultiResidue(ffx.potential.bonded.MultiResidue)
Aggregations
Rotamer (ffx.potential.bonded.Rotamer)56
Residue (ffx.potential.bonded.Residue)44
MultiResidue (ffx.potential.bonded.MultiResidue)42
RotamerLibrary.applyRotamer (ffx.potential.bonded.RotamerLibrary.applyRotamer)40
IOException (java.io.IOException)12
NACorrectionException (ffx.potential.bonded.NACorrectionException)10
Atom (ffx.potential.bonded.Atom)8
ResidueState (ffx.potential.bonded.ResidueState)8
ArrayList (java.util.ArrayList)7
File (java.io.File)6
PDBFilter (ffx.potential.parsers.PDBFilter)4
BufferedWriter (java.io.BufferedWriter)4
FileWriter (java.io.FileWriter)4
Polymer (ffx.potential.bonded.Polymer)3
TitrationUtils.inactivateResidue (ffx.potential.extended.TitrationUtils.inactivateResidue)3
AminoAcid3 (ffx.potential.bonded.ResidueEnumerations.AminoAcid3)2
RotamerLibrary (ffx.potential.bonded.RotamerLibrary)2
Torsion (ffx.potential.bonded.Torsion)2
Test (org.junit.Test)2
BooleanBuf (edu.rit.mp.BooleanBuf)1
