use of ffx.potential.bonded.ResidueState in project ffx by mjschnie.
the class RotamerOptimization method evaluateDistance.
/**
* Evaluates the pairwise distance between two residues' rotamers under any
* symmetry operator; does "lazy loading" for the distance matrix.
*
* @param i Residue i
* @param ri Rotamer for i
* @param j Residue j
* @param rj Rotamer for j
* @return Shortest distance
*/
private double evaluateDistance(int i, int ri, int j, int rj) {
Residue resi = allResiduesArray[i];
Rotamer[] rotamersI = resi.getRotamers(library);
Rotamer roti = rotamersI[ri];
double[][] xi;
if (roti.equals(resi.getRotamer())) {
xi = resi.storeCoordinateArray();
} else {
ResidueState origI = resi.storeState();
RotamerLibrary.applyRotamer(resi, roti);
xi = resi.storeCoordinateArray();
resi.revertState(origI);
}
Residue resj = allResiduesArray[j];
Rotamer[] rotamersJ = resj.getRotamers(library);
Rotamer rotj = rotamersJ[rj];
double[][] xj;
if (rotj.equals(resj.getRotamer())) {
xj = resj.storeCoordinateArray();
} else {
ResidueState origJ = resj.storeState();
RotamerLibrary.applyRotamer(resj, rotj);
xj = resj.storeCoordinateArray();
resj.revertState(origJ);
}
Crystal crystal = molecularAssembly.getCrystal();
int nSymm = crystal.spaceGroup.getNumberOfSymOps();
double minDist = Double.MAX_VALUE;
for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
double dist = interResidueDistance(xi, xj, symOp);
minDist = dist < minDist ? dist : minDist;
}
return minDist;
}
use of ffx.potential.bonded.ResidueState in project ffx by mjschnie.
the class RotamerOptimization method rotamerOptimizationDEE.
/**
* A global optimization over side-chain rotamers using a recursive
* algorithm and information about eliminated rotamers, rotamer pairs and
* rotamer triples
*
* @param molecularAssembly
* @param residues
* @param i
* @param currentRotamers
* @param lowEnergy
* @param optimum Optimum set of rotamers.
* @param permutationEnergies Energies of visited permutations or null.
* @return current energy.
*/
private double rotamerOptimizationDEE(MolecularAssembly molecularAssembly, Residue[] residues, int i, int[] currentRotamers, double lowEnergy, int[] optimum, double[] permutationEnergies) {
// This is the initialization condition.
if (i == 0) {
evaluatedPermutations = 0;
}
int nResidues = residues.length;
Residue residuei = residues[i];
Rotamer[] rotamersi = residuei.getRotamers(library);
int lenri = rotamersi.length;
double currentEnergy = Double.MAX_VALUE;
List<Residue> resList = Arrays.asList(residues);
/**
* As long as there are more residues, continue the recursion for each
* rotamer of the current residue.
*/
if (i < nResidues - 1) {
/**
* Loop over rotamers of residue i.
*/
for (int ri = 0; ri < lenri; ri++) {
/**
* Check if rotamer ri has been eliminated by DEE.
*/
if (check(i, ri)) {
continue;
}
/**
* Check if rotamer ri has been eliminated by an upstream
* rotamer (any residue's rotamer from j = 0 .. i-1).
*/
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;
}
applyRotamer(residuei, rotamersi[ri]);
currentRotamers[i] = ri;
double rotEnergy = rotamerOptimizationDEE(molecularAssembly, residues, i + 1, currentRotamers, lowEnergy, optimum, permutationEnergies);
if (rotEnergy < currentEnergy) {
currentEnergy = rotEnergy;
}
if (rotEnergy < lowEnergy) {
optimum[i] = ri;
lowEnergy = rotEnergy;
}
}
} else {
if (ensembleStates == null) {
ensembleStates = new ArrayList<>();
}
/**
* At the end of the recursion, compute the potential energy for
* each rotamer of the final residue. If a lower potential energy is
* discovered, the rotamers of each residue will be collected as the
* recursion returns up the chain.
*/
for (int ri = 0; ri < lenri; ri++) {
/**
* Check if rotamer ri has been eliminated by DEE.
*/
if (check(i, ri)) {
continue;
}
currentRotamers[i] = ri;
/**
* Check if rotamer ri has been eliminated by an upstream
* rotamer (any residue's rotamer from 0 .. i-1.
*/
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;
}
applyRotamer(residuei, rotamersi[ri]);
// Compute the energy based on a 3-body approximation
double approximateEnergy = computeEnergy(residues, currentRotamers, false);
double comparisonEnergy = approximateEnergy;
evaluatedPermutations++;
// Compute the AMOEBA energy
if (useFullAMOEBAEnergy) {
double amoebaEnergy = Double.NaN;
try {
amoebaEnergy = currentEnergy(resList);
} catch (ArithmeticException ex) {
logger.warning(String.format(" Exception %s in calculating full AMOEBA energy for permutation %d", ex.toString(), evaluatedPermutations));
}
comparisonEnergy = amoebaEnergy;
}
if (permutationEnergies != null) {
permutationEnergies[evaluatedPermutations - 1] = comparisonEnergy;
}
if (algorithmListener != null) {
algorithmListener.algorithmUpdate(molecularAssembly);
}
if (ensembleNumber > 1) {
if (master && printFiles) {
try {
FileWriter fw = new FileWriter(ensembleFile, true);
BufferedWriter bw = new BufferedWriter(fw);
bw.write(format("MODEL %d", evaluatedPermutations));
for (int j = 0; j < 75; j++) {
bw.write(" ");
}
bw.newLine();
bw.flush();
ensembleFilter.writeFile(ensembleFile, true);
bw.write(format("ENDMDL"));
for (int j = 0; j < 64; j++) {
bw.write(" ");
}
bw.newLine();
bw.close();
} catch (IOException e) {
logger.warning(format("Exception writing to file: %s", ensembleFile.getName()));
}
}
ResidueState[] states = ResidueState.storeAllCoordinates(residues);
ensembleStates.add(new ObjectPair<>(states, comparisonEnergy));
}
if (comparisonEnergy < currentEnergy) {
currentEnergy = comparisonEnergy;
}
if (comparisonEnergy < lowEnergy) {
lowEnergy = comparisonEnergy;
optimum[i] = ri;
}
if (useFullAMOEBAEnergy) {
// Log current results
logIfMaster(format(" %6e AMOEBA: %12.4f 3-Body: %12.4f Neglected: %12.4f (%12.4f)", (double) evaluatedPermutations, comparisonEnergy, approximateEnergy, comparisonEnergy - approximateEnergy, lowEnergy));
} else {
if (threeBodyTerm) {
;
logIfMaster(format(" %12s %5s %25f %5s %25f %5s", evaluatedPermutations, "|", approximateEnergy, "|", lowEnergy, "|"));
// logIfMaster(format(" %6e 3-Body: %12.4f (%12.4f)",
// (double) evaluatedPermutations, approximateEnergy, lowEnergy));
} else {
logIfMaster(format(" %12s %5s %25f %5s %25f %5s", evaluatedPermutations, "|", approximateEnergy, "|", lowEnergy, "|"));
// logIfMaster(format(" %6e 2-Body: %12.4f (%12.4f)",
// (double) evaluatedPermutations, approximateEnergy, lowEnergy));
}
}
}
ensembleStates.sort(null);
}
return currentEnergy;
}
use of ffx.potential.bonded.ResidueState in project ffx by mjschnie.
the class PhMD method tryComboStep.
/**
* Attempt a combination titration/rotamer MC move.
*
* @param targetMulti
* @return accept/reject
*/
private boolean tryComboStep(MultiResidue targetMulti) {
if (CAUTIOUS) {
throw new UnsupportedOperationException();
}
// Record the pre-change total energy.
double previousTotalEnergy = currentTotalEnergy();
double previousElectrostaticEnergy = currentElectrostaticEnergy();
// Write the pre-combo snapshot.
writeSnapshot(true, StepType.COMBO, config.snapshots);
String startString = targetMulti.toString();
String startName = targetMulti.getActive().getName();
// Choose from the list of available titrations for the active residue.
List<Titration> avail = titrationMap.get(targetMulti.getActive());
Titration titration = avail.get(rng.nextInt(avail.size()));
// Perform the chosen titration.
TitrationType titrationType = performTitration(targetMulti, titration, config.inactivateBackground);
// Change rotamer state, but first save coordinates so we can return to them if rejected.
Residue residue = targetMulti.getActive();
ArrayList<Atom> atoms = residue.getAtomList();
ResidueState origState = residue.storeState();
double[] chi = new double[4];
RotamerLibrary.measureAARotamer(residue, chi, false);
AminoAcid3 aa = AminoAcid3.valueOf(residue.getName());
Rotamer origCoordsRotamer = new Rotamer(aa, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
// Swap to the new rotamer.
Rotamer[] rotamers = residue.getRotamers(library);
int rotaRand = rng.nextInt(rotamers.length);
RotamerLibrary.applyRotamer(residue, rotamers[rotaRand]);
// Write the post-combo snapshot.
writeSnapshot(false, StepType.COMBO, config.snapshots);
// Evaluate both MC criteria.
String endName = targetMulti.getActive().getName();
// Evaluate the titration probability of the step.
double pKaref = titration.pKa;
double dG_ref = titration.refEnergy;
double temperature = thermostat.getCurrentTemperature();
double kT = BOLTZMANN * temperature;
double dG_elec = currentElectrostaticEnergy() - previousElectrostaticEnergy;
if (config.zeroReferenceEnergies) {
dG_ref = 0.0;
}
double prefix = Math.log(10) * kT * (pH - pKaref);
if (titrationType == TitrationType.DEPROT) {
prefix = -prefix;
}
double postfix = dG_elec - dG_ref;
double dG_titr = prefix + postfix;
double titrCriterion = exp(-dG_titr / kT);
// Evaluate the rotamer probability of the step.
double dG_rota = currentTotalEnergy() - previousTotalEnergy;
double rotaCriterion = exp(-dG_rota / kT);
StringBuilder sb = new StringBuilder();
sb.append(String.format(" Assessing possible MC combo step:\n"));
sb.append(String.format(" dG_elec: %16.8f\n", dG_elec));
sb.append(String.format(" dG_titr: %16.8f\n", dG_titr));
sb.append(String.format(" dG_rota: %16.8f\n", dG_rota));
sb.append(String.format(" -----\n"));
// Automatic acceptance if both energy changes are favorable.
if (dG_titr < 0 && dG_rota < 0 && config.mcOverride != MCOverride.REJECT) {
sb.append(String.format(" Accepted!"));
logger.info(sb.toString());
numMovesAccepted++;
propagateInactiveResidues(titratingMultis, false);
return true;
} else {
// Conditionally accept based on combined probabilities.
if (dG_titr < 0 || config.mcOverride == MCOverride.ACCEPT) {
titrCriterion = 1.0;
}
if (dG_rota < 0) {
rotaCriterion = 1.0;
}
if (config.mcOverride == MCOverride.REJECT) {
titrCriterion = 0.0;
}
double metropolis = random();
double comboCriterion = titrCriterion * rotaCriterion;
sb.append(String.format(" titrCrit: %9.4f\n", titrCriterion));
sb.append(String.format(" rotaCrit: %9.4f\n", rotaCriterion));
sb.append(String.format(" criterion: %9.4f\n", comboCriterion));
sb.append(String.format(" rng: %9.4f\n", metropolis));
if (metropolis < comboCriterion) {
sb.append(String.format(" Accepted!"));
logger.info(sb.toString());
numMovesAccepted++;
propagateInactiveResidues(titratingMultis, false);
return true;
} else {
// Move was denied.
sb.append(String.format(" Denied."));
logger.info(sb.toString());
// Undo both pieces of the rejected move IN THE RIGHT ORDER.
RotamerLibrary.applyRotamer(residue, origCoordsRotamer);
performTitration(targetMulti, titration, config.inactivateBackground);
ffe.reInit();
molDyn.reInit();
return false;
}
}
}
use of ffx.potential.bonded.ResidueState in project ffx by mjschnie.
the class PhMD method tryRotamerStep.
/**
* Attempt a rotamer MC move.
*
* @param targetMulti
* @return accept/reject
*/
private boolean tryRotamerStep(MultiResidue targetMulti) {
if (CAUTIOUS) {
throw new UnsupportedOperationException();
}
// Record the pre-change total energy.
double previousTotalEnergy = currentTotalEnergy();
// Write the before-step snapshot.
writeSnapshot(true, StepType.ROTAMER, config.snapshots);
// Save coordinates so we can return to them if move is rejected.
Residue residue = targetMulti.getActive();
ArrayList<Atom> atoms = residue.getAtomList();
ResidueState origState = residue.storeState();
double[] chi = new double[4];
RotamerLibrary.measureAARotamer(residue, chi, false);
AminoAcid3 aa = AminoAcid3.valueOf(residue.getName());
Rotamer origCoordsRotamer = new Rotamer(aa, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
// Select a new rotamer and swap to it.
Rotamer[] rotamers = residue.getRotamers(library);
int rotaRand = rng.nextInt(rotamers.length);
RotamerLibrary.applyRotamer(residue, rotamers[rotaRand]);
// Write the post-rotamer change snapshot.
writeSnapshot(false, StepType.ROTAMER, config.snapshots);
// Check the MC criterion.
double temperature = thermostat.getCurrentTemperature();
double kT = BOLTZMANN * temperature;
double postTotalEnergy = currentTotalEnergy();
double dG_tot = postTotalEnergy - previousTotalEnergy;
double criterion = exp(-dG_tot / kT);
StringBuilder sb = new StringBuilder();
sb.append(String.format(" Assessing possible MC rotamer step:\n"));
sb.append(String.format(" prev: %16.8f\n", previousTotalEnergy));
sb.append(String.format(" post: %16.8f\n", postTotalEnergy));
sb.append(String.format(" dG_tot: %16.8f\n", dG_tot));
sb.append(String.format(" -----\n"));
// Automatic acceptance if energy change is favorable.
if (dG_tot < 0) {
sb.append(String.format(" Accepted!"));
logger.info(sb.toString());
numMovesAccepted++;
propagateInactiveResidues(titratingMultis, true);
return true;
} else {
// Conditional acceptance if energy change is positive.
double metropolis = random();
sb.append(String.format(" criterion: %9.4f\n", criterion));
sb.append(String.format(" rng: %9.4f\n", metropolis));
if (metropolis < criterion) {
sb.append(String.format(" Accepted!"));
logger.info(sb.toString());
numMovesAccepted++;
propagateInactiveResidues(titratingMultis, true);
return true;
} else {
// Move was denied.
sb.append(String.format(" Denied."));
logger.info(sb.toString());
// Undo the rejected move.
RotamerLibrary.applyRotamer(residue, origCoordsRotamer);
return false;
}
}
}
use of ffx.potential.bonded.ResidueState 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;
}
Aggregations