use of ffx.potential.utils.PotentialsUtils in project ffx by mjschnie.
the class CrystalReciprocalSpaceTest method test1N7SPermanent.
/**
* Test of permanent method, of class CrystalReciprocalSpace.
*/
@Test
public void test1N7SPermanent() {
String filename = "ffx/xray/structures/1N7S.pdb";
int index = filename.lastIndexOf(".");
String name = filename.substring(0, index);
// load the structure
ClassLoader cl = this.getClass().getClassLoader();
File structure = new File(cl.getResource(filename).getPath());
PotentialsUtils potutil = new PotentialsUtils();
MolecularAssembly mola = potutil.open(structure);
CompositeConfiguration properties = mola.getProperties();
Crystal crystal = new Crystal(39.767, 51.750, 132.938, 90.00, 90.00, 90.00, "P212121");
Resolution resolution = new Resolution(1.45);
ReflectionList reflectionList = new ReflectionList(crystal, resolution);
DiffractionRefinementData refinementData = new DiffractionRefinementData(properties, reflectionList);
mola.finalize(true, mola.getForceField());
ForceFieldEnergy energy = mola.getPotentialEnergy();
List<Atom> atomList = mola.getAtomList();
Atom[] atomArray = atomList.toArray(new Atom[atomList.size()]);
// set up FFT and run it
ParallelTeam parallelTeam = new ParallelTeam();
CrystalReciprocalSpace crs = new CrystalReciprocalSpace(reflectionList, atomArray, parallelTeam, parallelTeam);
crs.computeAtomicDensity(refinementData.fc);
// tests
ComplexNumber b = new ComplexNumber(-828.584, -922.704);
HKL hkl = reflectionList.getHKL(1, 1, 4);
ComplexNumber a = refinementData.getFc(hkl.index());
System.out.println("1 1 4: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
assertEquals("1 1 4 reflection should be correct", -753.4722104328416, a.re(), 0.0001);
assertEquals("1 1 4 reflection should be correct", -1012.1341308707799, a.im(), 0.0001);
b.re(-70.4582);
b.im(-486.142);
hkl = reflectionList.getHKL(2, 1, 10);
a = refinementData.getFc(hkl.index());
System.out.println("2 1 10: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
assertEquals("2 1 10 reflection should be correct", -69.39660884054359, a.re(), 0.0001);
assertEquals("2 1 10 reflection should be correct", -412.0147625765328, a.im(), 0.0001);
}
use of ffx.potential.utils.PotentialsUtils in project ffx by mjschnie.
the class MolecularDynamics method writeStoredSnapshots.
/**
* Performs the inner loop of writing snapshots to disk; used by both detectAtypicalEnergy and a try-catch in dynamics.
*/
private void writeStoredSnapshots() {
int numSnaps = lastSnapshots.size();
File origFile = molecularAssembly.getFile();
String timeString = LocalDateTime.now().format(DateTimeFormatter.ofPattern("HH_mm_ss"));
PotentialsFunctions ef = new PotentialsUtils();
String filename = String.format("%s-%s-SNAP.pdb", FilenameUtils.removeExtension(molecularAssembly.getFile().getName()), timeString);
for (int is = 0; is < numSnaps; is++) {
DynamicsState oldState = lastSnapshots.poll();
oldState.revertState();
ef.saveAsPDB(molecularAssembly, new File(ef.versionFile(filename)));
}
molecularAssembly.setFile(origFile);
}
use of ffx.potential.utils.PotentialsUtils in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method energyAndGradient.
@Override
public double energyAndGradient(double[] x, double[] g, boolean verbose) {
if (lambdaBondedTerms) {
return 0.0;
}
/**
* Un-scale the coordinates.
*/
if (optimizationScaling != null) {
int len = x.length;
for (int i = 0; i < len; i++) {
x[i] /= optimizationScaling[i];
}
}
setCoordinates(x);
setOpenMMPositions(x, x.length);
int infoMask = OpenMM_State_Energy;
infoMask += OpenMM_State_Forces;
if (vdwLambdaTerm) {
infoMask += OpenMM_State_ParameterDerivatives;
}
state = OpenMM_Context_getState(context, infoMask, enforcePBC);
double e = OpenMM_State_getPotentialEnergy(state) / OpenMM_KJPerKcal;
if (vdwLambdaTerm) {
PointerByReference parameterArray = OpenMM_State_getEnergyParameterDerivatives(state);
int numDerives = OpenMM_ParameterArray_getSize(parameterArray);
if (numDerives > 0) {
vdwdUdL = OpenMM_ParameterArray_get(parameterArray, pointerForString("vdw_lambda")) / OpenMM_KJPerKcal;
}
}
if (maxDebugGradient < Double.MAX_VALUE) {
boolean extremeGrad = Arrays.stream(g).anyMatch((double gi) -> {
return (gi > maxDebugGradient || gi < -maxDebugGradient);
});
if (extremeGrad) {
File origFile = molecularAssembly.getFile();
String timeString = LocalDateTime.now().format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
String filename = String.format("%s-LARGEGRAD-%s.pdb", FilenameUtils.removeExtension(molecularAssembly.getFile().getName()), timeString);
PotentialsFunctions ef = new PotentialsUtils();
filename = ef.versionFile(filename);
logger.warning(String.format(" Excessively large gradients detected; printing snapshot to file %s", filename));
ef.saveAsPDB(molecularAssembly, new File(filename));
molecularAssembly.setFile(origFile);
}
}
if (verbose) {
logger.log(Level.INFO, String.format(" OpenMM Energy: %14.10g", e));
}
forces = OpenMM_State_getForces(state);
fillGradients(g);
/**
* Scale the coordinates and gradients.
*/
if (optimizationScaling != null) {
int len = x.length;
for (int i = 0; i < len; i++) {
x[i] *= optimizationScaling[i];
g[i] /= optimizationScaling[i];
}
}
OpenMM_State_destroy(state);
return e;
}
use of ffx.potential.utils.PotentialsUtils in project ffx by mjschnie.
the class ForceFieldEnergy method energyAndGradient.
/**
* {@inheritDoc}
*/
@Override
public double energyAndGradient(double[] x, double[] g, boolean verbose) {
/**
* Un-scale the coordinates.
*/
if (optimizationScaling != null) {
int len = x.length;
for (int i = 0; i < len; i++) {
x[i] /= optimizationScaling[i];
}
}
setCoordinates(x);
double e = energy(true, verbose);
// need to try-catch fillGradients.
try {
fillGradients(g);
/**
* Scale the coordinates and gradients.
*/
if (optimizationScaling != null) {
int len = x.length;
for (int i = 0; i < len; i++) {
x[i] *= optimizationScaling[i];
g[i] /= optimizationScaling[i];
}
}
if (maxDebugGradient < Double.MAX_VALUE) {
boolean extremeGrad = Arrays.stream(g).anyMatch((double gi) -> {
return (gi > maxDebugGradient || gi < -maxDebugGradient);
});
if (extremeGrad) {
File origFile = molecularAssembly.getFile();
String timeString = LocalDateTime.now().format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
String filename = String.format("%s-LARGEGRAD-%s.pdb", FilenameUtils.removeExtension(molecularAssembly.getFile().getName()), timeString);
PotentialsFunctions ef = new PotentialsUtils();
filename = ef.versionFile(filename);
logger.warning(String.format(" Excessively large gradients detected; printing snapshot to file %s", filename));
ef.saveAsPDB(molecularAssembly, new File(filename));
molecularAssembly.setFile(origFile);
}
}
return e;
} catch (EnergyException ex) {
ex.printStackTrace();
if (printOnFailure) {
String timeString = LocalDateTime.now().format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
String filename = String.format("%s-ERROR-%s.pdb", FilenameUtils.removeExtension(molecularAssembly.getFile().getName()), timeString);
PotentialsFunctions ef = new PotentialsUtils();
filename = ef.versionFile(filename);
logger.info(String.format(" Writing on-error snapshot to file %s", filename));
ef.saveAsPDB(molecularAssembly, new File(filename));
}
if (ex.doCauseSevere()) {
Utilities.printStackTrace(ex);
logger.log(Level.SEVERE, " Error in calculating energies or gradients", ex);
} else {
ex.printStackTrace();
// Rethrow exception
throw ex;
}
// Should ordinarily be unreachable.
throw ex;
}
}
use of ffx.potential.utils.PotentialsUtils in project ffx by mjschnie.
the class ForceFieldEnergy method energy.
/**
* <p>
* energy</p>
*
* @param gradient a boolean.
* @param print a boolean.
* @return a double.
*/
public double energy(boolean gradient, boolean print) {
try {
bondTime = 0;
angleTime = 0;
stretchBendTime = 0;
ureyBradleyTime = 0;
outOfPlaneBendTime = 0;
torsionTime = 0;
piOrbitalTorsionTime = 0;
torsionTorsionTime = 0;
improperTorsionTime = 0;
vanDerWaalsTime = 0;
electrostaticTime = 0;
restraintBondTime = 0;
ncsTime = 0;
coordRestraintTime = 0;
totalTime = System.nanoTime();
// Zero out the potential energy of each bonded term.
bondEnergy = 0.0;
angleEnergy = 0.0;
stretchBendEnergy = 0.0;
ureyBradleyEnergy = 0.0;
outOfPlaneBendEnergy = 0.0;
torsionEnergy = 0.0;
piOrbitalTorsionEnergy = 0.0;
torsionTorsionEnergy = 0.0;
improperTorsionEnergy = 0.0;
totalBondedEnergy = 0.0;
// Zero out potential energy of restraint terms
restraintBondEnergy = 0.0;
ncsEnergy = 0.0;
restrainEnergy = 0.0;
// Zero out bond and angle RMSDs.
bondRMSD = 0.0;
angleRMSD = 0.0;
// Zero out the potential energy of each non-bonded term.
vanDerWaalsEnergy = 0.0;
permanentMultipoleEnergy = 0.0;
permanentRealSpaceEnergy = 0.0;
permanentSelfEnergy = 0.0;
permanentReciprocalEnergy = 0.0;
polarizationEnergy = 0.0;
inducedRealSpaceEnergy = 0.0;
inducedSelfEnergy = 0.0;
inducedReciprocalEnergy = 0.0;
totalMultipoleEnergy = 0.0;
totalNonBondedEnergy = 0.0;
// Zero out the solvation energy.
solvationEnergy = 0.0;
// Zero out the relative solvation energy (sequence optimization)
relativeSolvationEnergy = 0.0;
nRelativeSolvations = 0;
esvBias = 0.0;
// Zero out the total potential energy.
totalEnergy = 0.0;
// Zero out the Cartesian coordinate gradient for each atom.
if (gradient) {
for (int i = 0; i < nAtoms; i++) {
atoms[i].setXYZGradient(0.0, 0.0, 0.0);
atoms[i].setLambdaXYZGradient(0.0, 0.0, 0.0);
}
}
/**
* Computed the bonded energy terms in parallel.
*/
try {
bondedRegion.setGradient(gradient);
parallelTeam.execute(bondedRegion);
} catch (RuntimeException ex) {
logger.warning("Runtime exception during bonded term calculation.");
throw ex;
} catch (Exception ex) {
Utilities.printStackTrace(ex);
logger.severe(ex.toString());
}
if (!lambdaBondedTerms) {
/**
* Compute restraint terms.
*/
if (ncsTerm) {
ncsTime = -System.nanoTime();
ncsEnergy = ncsRestraint.residual(gradient, print);
ncsTime += System.nanoTime();
}
if (restrainTerm && !coordRestraints.isEmpty()) {
coordRestraintTime = -System.nanoTime();
for (CoordRestraint restraint : coordRestraints) {
restrainEnergy += restraint.residual(gradient, print);
}
coordRestraintTime += System.nanoTime();
}
if (comTerm) {
comRestraintTime = -System.nanoTime();
comRestraintEnergy = comRestraint.residual(gradient, print);
comRestraintTime += System.nanoTime();
}
/**
* Compute non-bonded terms.
*/
if (vanderWaalsTerm) {
vanDerWaalsTime = -System.nanoTime();
vanDerWaalsEnergy = vanderWaals.energy(gradient, print);
nVanDerWaalInteractions = this.vanderWaals.getInteractions();
vanDerWaalsTime += System.nanoTime();
}
if (multipoleTerm) {
electrostaticTime = -System.nanoTime();
totalMultipoleEnergy = particleMeshEwald.energy(gradient, print);
permanentMultipoleEnergy = particleMeshEwald.getPermanentEnergy();
permanentRealSpaceEnergy = particleMeshEwald.getPermRealEnergy();
permanentSelfEnergy = particleMeshEwald.getPermSelfEnergy();
permanentReciprocalEnergy = particleMeshEwald.getPermRecipEnergy();
polarizationEnergy = particleMeshEwald.getPolarizationEnergy();
inducedRealSpaceEnergy = particleMeshEwald.getIndRealEnergy();
inducedSelfEnergy = particleMeshEwald.getIndSelfEnergy();
inducedReciprocalEnergy = particleMeshEwald.getIndRecipEnergy();
nPermanentInteractions = particleMeshEwald.getInteractions();
solvationEnergy = particleMeshEwald.getGKEnergy();
nGKInteractions = particleMeshEwald.getGKInteractions();
electrostaticTime += System.nanoTime();
}
}
if (relativeSolvationTerm) {
List<Residue> residuesList = molecularAssembly.getResidueList();
for (Residue residue : residuesList) {
if (residue instanceof MultiResidue) {
Atom refAtom = residue.getSideChainAtoms().get(0);
if (refAtom != null && refAtom.getUse()) {
/**
* Reasonably confident that it should be -=, as we
* are trying to penalize residues with strong
* solvation energy.
*/
double thisSolvation = relativeSolvation.getSolvationEnergy(residue, false);
relativeSolvationEnergy -= thisSolvation;
if (thisSolvation != 0) {
nRelativeSolvations++;
}
}
}
}
}
totalTime = System.nanoTime() - totalTime;
totalBondedEnergy = bondEnergy + restraintBondEnergy + angleEnergy + stretchBendEnergy + ureyBradleyEnergy + outOfPlaneBendEnergy + torsionEnergy + piOrbitalTorsionEnergy + improperTorsionEnergy + torsionTorsionEnergy + ncsEnergy + restrainEnergy;
totalNonBondedEnergy = vanDerWaalsEnergy + totalMultipoleEnergy + relativeSolvationEnergy;
totalEnergy = totalBondedEnergy + totalNonBondedEnergy + solvationEnergy;
if (esvTerm) {
esvBias = esvSystem.getBiasEnergy();
totalEnergy += esvBias;
}
} catch (EnergyException ex) {
if (printOnFailure) {
File origFile = molecularAssembly.getFile();
String timeString = LocalDateTime.now().format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
String filename = String.format("%s-ERROR-%s.pdb", FilenameUtils.removeExtension(molecularAssembly.getFile().getName()), timeString);
PotentialsFunctions ef = new PotentialsUtils();
filename = ef.versionFile(filename);
logger.info(String.format(" Writing on-error snapshot to file %s", filename));
ef.saveAsPDB(molecularAssembly, new File(filename));
molecularAssembly.setFile(origFile);
}
if (ex.doCauseSevere()) {
logger.log(Level.SEVERE, " Error in calculating energies or gradients", ex);
return 0.0;
} else {
// Rethrow exception
throw ex;
}
}
if (print || printOverride) {
if (printCompact) {
logger.info(this.toString());
} else {
StringBuilder sb = new StringBuilder();
if (gradient) {
sb.append("\n Computed Potential Energy and Atomic Coordinate Gradients\n");
} else {
sb.append("\n Computed Potential Energy\n");
}
sb.append(this);
logger.info(sb.toString());
}
}
return totalEnergy;
}
Aggregations