use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergy method getdEdXdL.
/**
* {@inheritDoc}
*
* @param gradients
*/
@Override
public void getdEdXdL(double[] gradients) {
if (!lambdaBondedTerms) {
if (vanderWaalsTerm) {
vanderWaals.getdEdXdL(gradients);
}
if (multipoleTerm) {
particleMeshEwald.getdEdXdL(gradients);
}
if (restraintBondTerm) {
for (int i = 0; i < nRestraintBonds; i++) {
restraintBonds[i].getdEdXdL(gradients);
}
}
if (ncsTerm && ncsRestraint != null) {
ncsRestraint.getdEdXdL(gradients);
}
if (restrainTerm && !coordRestraints.isEmpty()) {
for (CoordRestraint restraint : coordRestraints) {
restraint.getdEdXdL(gradients);
}
// autoCoordRestraint.getdEdXdL(gradients);
}
if (comTerm && comRestraint != null) {
comRestraint.getdEdXdL(gradients);
}
if (lambdaTorsions) {
double[] grad = new double[3];
int index = 0;
for (int i = 0; i < nAtoms; i++) {
Atom a = atoms[i];
if (a.isActive()) {
a.getLambdaXYZGradient(grad);
gradients[index++] += grad[0];
gradients[index++] += grad[1];
gradients[index++] += grad[2];
}
}
}
}
}
use of ffx.potential.bonded.Atom 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;
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergy method getCoordinates.
/**
* {@inheritDoc}
*
* @param x
*/
@Override
public double[] getCoordinates(double[] x) {
int n = getNumberOfVariables();
if (x == null || x.length < n) {
x = new double[n];
}
int index = 0;
for (int i = 0; i < nAtoms; i++) {
Atom a = atoms[i];
if (a.isActive() && !a.isBackground()) {
x[index++] = a.getX();
x[index++] = a.getY();
x[index++] = a.getZ();
}
}
return x;
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method addGKForce.
private void addGKForce() {
GeneralizedKirkwood gk = super.getGK();
amoebaGeneralizedKirkwoodForce = OpenMM_AmoebaGeneralizedKirkwoodForce_create();
OpenMM_AmoebaGeneralizedKirkwoodForce_setSolventDielectric(amoebaGeneralizedKirkwoodForce, 78.3);
OpenMM_AmoebaGeneralizedKirkwoodForce_setSoluteDielectric(amoebaGeneralizedKirkwoodForce, 1.0);
double[] overlapScale = gk.getOverlapScale();
double[] baseRadii = gk.getBaseRadii();
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
MultipoleType multipoleType = atoms[i].getMultipoleType();
OpenMM_AmoebaGeneralizedKirkwoodForce_addParticle(amoebaGeneralizedKirkwoodForce, multipoleType.charge, OpenMM_NmPerAngstrom * baseRadii[i], overlapScale[i]);
}
OpenMM_AmoebaGeneralizedKirkwoodForce_setProbeRadius(amoebaGeneralizedKirkwoodForce, 1.4 * OpenMM_NmPerAngstrom);
NonPolar nonpolar = gk.getNonPolarModel();
switch(nonpolar) {
case BORN_SOLV:
case BORN_CAV_DISP:
default:
// Configure a Born Radii based surface area term.
double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal * OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm;
OpenMM_AmoebaGeneralizedKirkwoodForce_setIncludeCavityTerm(amoebaGeneralizedKirkwoodForce, OpenMM_True);
OpenMM_AmoebaGeneralizedKirkwoodForce_setSurfaceAreaFactor(amoebaGeneralizedKirkwoodForce, -surfaceTension);
break;
case CAV:
case CAV_DISP:
case HYDROPHOBIC_PMF:
case NONE:
// This NonPolar model does not use a Born Radii based surface area term.
OpenMM_AmoebaGeneralizedKirkwoodForce_setIncludeCavityTerm(amoebaGeneralizedKirkwoodForce, OpenMM_False);
break;
}
OpenMM_System_addForce(system, amoebaGeneralizedKirkwoodForce);
switch(nonpolar) {
case CAV_DISP:
case BORN_CAV_DISP:
addWCAForce();
break;
case CAV:
case HYDROPHOBIC_PMF:
case BORN_SOLV:
case NONE:
default:
}
logger.log(Level.INFO, " Added generalized Kirkwood force.");
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method addRestraintBonds.
/**
* Adds restraint bonds, if any.
*/
private void addRestraintBonds() {
List<RestraintBond> restraintBonds = super.getRestraintBonds();
if (restraintBonds != null && !restraintBonds.isEmpty()) {
// OpenMM's HarmonicBondForce class uses k, not 1/2*k as does FFX.
double kParameterConversion = BondType.units * 2.0 * OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
for (RestraintBond rbond : super.getRestraintBonds()) {
// Is not a valid substitute for a targeting computer.
PointerByReference theForce;
BondType bType = rbond.bondType;
BondType.BondFunction funct = bType.bondFunction;
if (restraintForces.containsKey(funct)) {
theForce = restraintForces.get(funct);
} else {
theForce = OpenMM_CustomBondForce_create(funct.toMathematicalForm());
OpenMM_CustomBondForce_addPerBondParameter(theForce, "k");
OpenMM_CustomBondForce_addPerBondParameter(theForce, "r0");
if (funct.hasFlatBottom()) {
OpenMM_CustomBondForce_addPerBondParameter(theForce, "fb");
}
// Wholly untested code.
switch(funct) {
case QUARTIC:
case FLAT_BOTTOM_QUARTIC:
OpenMM_CustomBondForce_addGlobalParameter(theForce, "cubic", BondType.cubic / OpenMM_NmPerAngstrom);
OpenMM_CustomBondForce_addGlobalParameter(theForce, "quartic", BondType.quartic / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
break;
default:
break;
}
OpenMM_System_addForce(system, theForce);
}
double forceConst = bType.forceConstant * kParameterConversion;
double equilDist = bType.distance * OpenMM_NmPerAngstrom;
Atom[] ats = rbond.getAtomArray();
int at1 = ats[0].getXyzIndex() - 1;
int at2 = ats[1].getXyzIndex() - 1;
PointerByReference bondParams = OpenMM_DoubleArray_create(0);
OpenMM_DoubleArray_append(bondParams, forceConst);
OpenMM_DoubleArray_append(bondParams, equilDist);
if (funct.hasFlatBottom()) {
OpenMM_DoubleArray_append(bondParams, bType.flatBottomRadius * OpenMM_NmPerAngstrom);
}
OpenMM_CustomBondForce_addBond(theForce, at1, at2, bondParams);
OpenMM_DoubleArray_destroy(bondParams);
}
}
}
Aggregations