use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class RotamerOptimization method pruneSingleClashes.
/**
* Uses calculated energies to prune rotamers based on a threshold distance
* from that residue's minimum energy rotamer (by default 20 kcal/mol). The
* threshold can be modulated by presence of nucleic acids or MultiResidues,
* which require more generous pruning criteria.
*
* @param residues Residues to prune rotamers over.
*/
private void pruneSingleClashes(Residue[] residues) {
if (!pruneClashes) {
return;
}
for (int i = 0; i < residues.length; i++) {
Residue residue = residues[i];
Rotamer[] rotamers = residue.getRotamers(library);
int nrot = rotamers.length;
double minEnergy = Double.MAX_VALUE;
int minRot = -1;
for (int ri = 0; ri < nrot; ri++) {
if (!check(i, ri) && getSelf(i, ri) < minEnergy) {
minEnergy = getSelf(i, ri);
minRot = ri;
}
}
/**
* Regular: ep = minEnergy + clashThreshold Nucleic acids: ep =
* minEnergy + (clashThreshold * factor * factor) MultiResidues: ep
* = minEnergy + multiResClashThreshold
*
* Nucleic acids are bigger than amino acids, and MultiResidues can
* have wild swings in energy on account of chemical perturbation.
*/
double energyToPrune = (residue instanceof MultiResidue) ? multiResClashThreshold : clashThreshold;
energyToPrune = (residue.getResidueType() == NA) ? energyToPrune * singletonNAPruningFactor * pruningFactor : energyToPrune;
energyToPrune += minEnergy;
for (int ri = 0; ri < nrot; ri++) {
if (!check(i, ri) && (getSelf(i, ri) > energyToPrune)) {
if (eliminateRotamer(residues, i, ri, print)) {
logIfMaster(format(" Rotamer (%7s,%2d) self-energy %s pruned by (%7s,%2d) %s.", residue, ri, formatEnergy(getSelf(i, ri)), residue, minRot, formatEnergy(minEnergy)));
}
}
}
}
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class RotamerOptimization method rotamerEnergies.
/**
* Turn off non-bonded contributions from all residues except for one.
* Compute the self-energy for each residue relative to the backbone
* contribution.
*
* @param residues A list of residues that we undergo rotamer optimization.
* @return template energy
*/
private double rotamerEnergies(Residue[] residues) {
if (residues == null) {
logger.warning(" Attempt to compute rotamer energies for an empty array of residues.");
return 0.0;
}
int nResidues = residues.length;
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
allocateEliminationMemory(residues);
if (decomposeOriginal) {
assert library.getUsingOrigCoordsRotamer();
for (int i = 0; i < nResidues; i++) {
Residue resi = residues[i];
Rotamer[] rotsi = resi.getRotamers(library);
int lenri = rotsi.length;
// Leave the 0'th original-coordinates rotamer alone.
for (int ri = 1; ri < lenri; ri++) {
eliminateRotamer(residues, i, ri, false);
}
}
}
/**
* Initialize all atoms to be used.
*/
for (int i = 0; i < nAtoms; i++) {
atoms[i].setUse(true);
}
if (!usingBoxOptimization) {
energiesToWrite = Collections.synchronizedList(new ArrayList<String>());
receiveThread = new ReceiveThread(residues);
receiveThread.start();
if (master && writeEnergyRestart && printFiles) {
energyWriterThread = new EnergyWriterThread(receiveThread);
energyWriterThread.start();
}
}
int loaded = 0;
if (loadEnergyRestart) {
if (usingBoxOptimization) {
loaded = loadEnergyRestart(energyRestartFile, residues, boxLoadIndex, boxLoadCellIndices);
} else {
loaded = loadEnergyRestart(energyRestartFile, residues);
}
}
long energyStartTime = System.nanoTime();
SelfEnergyRegion singlesRegion = new SelfEnergyRegion(residues);
TwoBodyEnergyRegion pairsRegion = new TwoBodyEnergyRegion(residues);
ThreeBodyEnergyRegion triplesRegion = new ThreeBodyEnergyRegion(residues);
FourBodyEnergyRegion quadsRegion = new FourBodyEnergyRegion(residues);
try {
if (loaded < 1) {
selfEnergyMap.clear();
// allocate selfEnergy
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++) {
if (!check(i, ri)) {
Integer[] selfJob = { i, ri };
if (decomposeOriginal && ri != 0) {
continue;
}
selfEnergyMap.put(singleJobIndex++, selfJob);
}
}
}
}
// broadcast that this proc is done with startup and allocation; ready for singles
boolean thisProcReady = true;
BooleanBuf thisProcReadyBuf = BooleanBuf.buffer(thisProcReady);
multicastBuf(thisProcReadyBuf);
// launch parallel singles calculation
while (!readyForSingles) {
Thread.sleep(POLLING_FREQUENCY);
}
logger.info(String.format(" Number of self energies to calculate: %d", selfEnergyMap.size()));
energyWorkerTeam.execute(singlesRegion);
long singlesTime = System.nanoTime() - energyStartTime;
logIfMaster(format(" Time for single energies: %12.4g", (singlesTime * 1.0E-9)));
if (loaded < 2) {
twoBodyEnergyMap.clear();
// allocate twoBodyEnergy and create jobs
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 (checkToJ(i, ri, j, rj)) {
continue;
}
Integer[] pairJob = { i, ri, j, rj };
if (decomposeOriginal && (ri != 0 || rj != 0)) {
continue;
}
twoBodyEnergyMap.put(pairJobIndex++, pairJob);
}
}
}
}
}
// broadcast that this proc is done with pruning and allocation; ready for pairs
multicastBuf(thisProcReadyBuf);
// launch parallel pair calculation
while (!readyFor2Body) {
Thread.sleep(POLLING_FREQUENCY);
}
logger.info(String.format(" Number of 2-body energies to calculate: %d", twoBodyEnergyMap.size()));
energyWorkerTeam.execute(pairsRegion);
long pairsTime = System.nanoTime() - (singlesTime + energyStartTime);
long triplesTime = 0;
long quadsTime = 0;
logIfMaster(format(" Time for 2-body energies: %12.4g", (pairsTime * 1.0E-9)));
if (threeBodyTerm) {
if (loaded < 3) {
threeBodyEnergyMap.clear();
// allocate threeBodyEnergy and create jobs
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 (checkToJ(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 (checkToK(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);
}
}
}
}
}
}
}
// broadcast that this proc is done with pruning and allocation; ready for trimers
multicastBuf(thisProcReadyBuf);
// launch parallel 3-Body calculation
while (!readyFor3Body) {
Thread.sleep(POLLING_FREQUENCY);
}
logger.info(String.format(" Number of 3-Body energies to calculate: %d", threeBodyEnergyMap.size()));
energyWorkerTeam.execute(triplesRegion);
triplesTime = System.nanoTime() - (pairsTime + singlesTime + energyStartTime);
logIfMaster(format(" Time for 3-Body energies: %12.4g", (triplesTime * 1.0E-9)));
}
if (compute4BodyEnergy) {
logger.info(" Creating 4-Body jobs...");
fourBodyEnergyMap.clear();
boolean maxedOut = false;
// create 4-Body jobs (no memory allocation)
int quadJobIndex = 0;
for (int i = 0; i < nResidues; i++) {
Residue resi = residues[i];
Rotamer[] roti = resi.getRotamers(library);
for (int ri = 0; ri < roti.length; ri++) {
if (check(i, ri)) {
continue;
}
for (int j = i + 1; j < nResidues; j++) {
Residue resj = residues[j];
Rotamer[] rotj = resj.getRotamers(library);
for (int rj = 0; rj < rotj.length; rj++) {
/*if (check(j, rj) || check(i, ri, j, rj)) {
continue;
}*/
if (checkToJ(i, ri, j, rj)) {
continue;
}
for (int k = j + 1; k < nResidues; k++) {
Residue resk = residues[k];
Rotamer[] rotk = resk.getRotamers(library);
for (int rk = 0; rk < rotk.length; rk++) {
/*if (check(k, rk) || check(i, ri, k, rk) || check(j, rj, k, rk) || check(i, ri, j, rj, k, rk)) {
continue;
}*/
if (checkToK(i, ri, j, rj, k, rk)) {
continue;
}
for (int l = k + 1; l < nResidues; l++) {
Residue resl = residues[l];
Rotamer[] rotl = resl.getRotamers(library);
for (int rl = 0; rl < rotl.length; rl++) {
/*if (check(l, rl) || check(i, ri, l, rl) ||
check(j, rj, l, rl) || check(k, rk, l, rl) ||
check(i, ri, j, rj, l, rl) || check(i, ri, k, rk, l, rl) ||
check(j, rj, k, rk, l, rl)) {
continue;
}*/
if (checkToL(i, ri, j, rj, k, rk, l, rl)) {
continue;
}
Integer[] quadJob = { i, ri, j, rj, k, rk, l, rl };
if (decomposeOriginal && (ri != 0 || rj != 0 || rk != 0 || rl != 0)) {
continue;
}
fourBodyEnergyMap.put(quadJobIndex++, quadJob);
if (fourBodyEnergyMap.size() > max4BodyCount) {
maxedOut = true;
break;
}
}
if (maxedOut) {
break;
}
}
if (maxedOut) {
break;
}
}
if (maxedOut) {
break;
}
}
if (maxedOut) {
break;
}
}
if (maxedOut) {
break;
}
}
if (maxedOut) {
break;
}
}
if (maxedOut) {
break;
}
}
// broadcast that this proc is done with pruning and allocation; ready for quads
// logger.info(format(" Proc %d broadcasting ready for quads.", world.rank()));
multicastBuf(thisProcReadyBuf);
// launch parallel 3-Body calculation
int waiting = 0;
while (!readyFor4Body) {
Thread.sleep(POLLING_FREQUENCY);
}
logger.info(format(" Running quads: %d jobs.", fourBodyEnergyMap.size()));
energyWorkerTeam.execute(quadsRegion);
quadsTime = System.nanoTime() - (triplesTime + pairsTime + singlesTime + energyStartTime);
logIfMaster(format(" Time for 4-Body energies: %12.4g", quadsTime * 1.0E-9));
}
long allTime = singlesTime + pairsTime + triplesTime + quadsTime;
logIfMaster(format(" Time for all energies: %12.4g", allTime * 1.0E-9));
} catch (Exception ex) {
String message = " Exception computing rotamer energies in parallel.";
logger.log(Level.SEVERE, message, ex);
}
// Turn on all atoms.
for (int i = 0; i < atoms.length; i++) {
atoms[i].setUse(true);
}
// Print the energy with all rotamers in their default conformation.
if (verboseEnergies && master) {
try {
double defaultEnergy = currentEnergy(residues);
logger.info(format(" Energy of the system with rotamers in their default conformation: %s", formatEnergy(defaultEnergy)));
} catch (ArithmeticException ex) {
logger.severe(String.format(" Exception %s in calculating default energy; FFX shutting down", ex.toString()));
}
}
return backboneEnergy;
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class RotamerOptimization method reconcileNARotamersWithSubsequentResidues.
/**
* For NA residues inside some optimization window, prune any rotamers which
* would be incompatible with the established rotamers of downstream NA
* residues. Could in theory be done by self energies, but every rotamer
* which can be eliminated without calculating a self energy makes the
* optimization much faster. Technically, this works by pinning it to its
* current pucker, but if the input structure is any good, it would be the
* current pucker anyways.
*
* @param residues Residues to check for incompatible rotamers.
* @return Number of rotamers eliminated for each Residue.
*/
private void reconcileNARotamersWithSubsequentResidues(Residue[] residues, int[] eliminatedRotamers) {
for (int i = 0; i < residues.length; i++) {
Residue residuei = residues[i];
switch(residuei.getResidueType()) {
case NA:
Rotamer[] rotamers = residues[i].getRotamers(library);
Residue nextResidue = residuei.getNextResidue();
if (rotamers == null || nextResidue == null) {
break;
}
boolean isInList = false;
for (Residue residue : residues) {
if (residue.equals(nextResidue)) {
isInList = true;
break;
}
}
if (isInList) {
break;
}
double delta = RotamerLibrary.measureDelta(residuei);
if (RotamerLibrary.checkPucker(delta) == 1) {
for (int j = 0; j < rotamers.length; j++) {
if (RotamerLibrary.checkPucker(rotamers[j].chi7) != 1) {
if (print) {
logIfMaster(format(" Rotamer %d of residue %s eliminated " + "for incompatibility with the sugar pucker of previous " + "residue %s outside the window.", j, residuei.toFormattedString(false, true), nextResidue.toString()));
}
eliminateRotamer(residues, i, j, print);
eliminatedRotamers[i]++;
}
}
} else {
for (int j = 0; j < rotamers.length; j++) {
if (RotamerLibrary.checkPucker(rotamers[j].chi7) != 2) {
if (print) {
logIfMaster(format(" Rotamer %d of residue %s eliminated " + "for incompatibility with the sugar pucker of previous " + "residue %s outside the window.", j, residuei.toFormattedString(false, true), nextResidue.toString()));
}
eliminateRotamer(residues, i, j, print);
eliminatedRotamers[i]++;
}
}
}
break;
case AA:
default:
break;
}
}
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class RotamerOptimization method dryRun.
/**
* A global optimization over side-chain rotamers using a recursive
* algorithm and information about eliminated rotamers, rotamer pairs and
* rotamer triples.
*
* @param residues
* @param i
* @param currentRotamers
* @return 0.
*/
private double dryRun(Residue[] residues, int i, int[] currentRotamers) {
// This is the initialization condition.
if (i == 0) {
evaluatedPermutations = 0;
evaluatedPermutationsPrint = 1000;
}
if (evaluatedPermutations >= evaluatedPermutationsPrint) {
if (evaluatedPermutations % evaluatedPermutationsPrint == 0) {
logIfMaster(format(" The permutations have reached %10.4e.", (double) evaluatedPermutationsPrint));
evaluatedPermutationsPrint *= 10;
}
}
int nResidues = residues.length;
Residue residuei = residues[i];
Rotamer[] rotamersi = residuei.getRotamers(library);
int lenri = rotamersi.length;
if (i < nResidues - 1) {
for (int ri = 0; ri < lenri; ri++) {
if (check(i, ri)) {
continue;
}
boolean deadEnd = false;
for (int j = 0; j < i; j++) {
int rj = currentRotamers[j];
deadEnd = check(j, rj, i, ri);
if (deadEnd) {
break;
}
}
if (deadEnd) {
continue;
}
currentRotamers[i] = ri;
dryRun(residues, i + 1, currentRotamers);
}
} else {
/**
* At the end of the recursion, check each rotamer of the final
* residue.
*/
for (int ri = 0; ri < lenri; ri++) {
if (check(i, ri)) {
continue;
}
currentRotamers[i] = ri;
boolean deadEnd = false;
for (int j = 0; j < i; j++) {
int rj = currentRotamers[j];
deadEnd = check(j, rj, i, ri);
if (deadEnd) {
break;
}
}
if (!deadEnd) {
evaluatedPermutations++;
}
}
}
return 0.0;
}
use of ffx.potential.bonded.Residue in project ffx by mjschnie.
the class RotamerOptimization method validateDEE.
private boolean validateDEE(Residue[] residues) {
int nres = eliminatedSingles.length;
// Validate residues
for (int i = 0; i < nres; i++) {
Residue residuei = residues[i];
int ni = eliminatedSingles[i].length;
boolean valid = false;
for (int ri = 0; ri < ni; ri++) {
if (!check(i, ri)) {
valid = true;
}
}
if (!valid) {
logger.severe(format(" Coding error: all %d rotamers for residue %s eliminated.", ni, residuei));
}
}
// Validate pairs
for (int i = 0; i < nres; i++) {
Residue residuei = residues[i];
Rotamer[] rotamersi = residuei.getRotamers(library);
int ni = rotamersi.length;
for (int j = i + 1; j < nres; j++) {
Residue residuej = residues[j];
Rotamer[] rotamersj = residuej.getRotamers(library);
int nj = rotamersj.length;
boolean valid = false;
for (int ri = 0; ri < ni; ri++) {
for (int rj = 0; rj < nj; rj++) {
if (!check(i, ri, j, rj)) {
valid = true;
}
}
}
if (!valid) {
logger.severe(format(" Coding error: all pairs for %s with residue %s eliminated.", residuei.toFormattedString(false, true), residuej));
}
}
}
return true;
}
Aggregations