use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method updateAmoebaMultipoleForce.
/**
* Updates the Amoeba electrostatic multipolar force for change in Use
* flags.
*
* @param atoms Array of all Atoms in the system
*/
private void updateAmoebaMultipoleForce(Atom[] atoms) {
ParticleMeshEwald pme = super.getPmeNode();
int[][] axisAtom = pme.getAxisAtoms();
double dipoleConversion = OpenMM_NmPerAngstrom;
double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
double polarScale = 1.0;
if (pme.getPolarizationType() == Polarization.NONE) {
polarScale = 0.0;
}
PointerByReference dipoles = OpenMM_DoubleArray_create(3);
PointerByReference quadrupoles = OpenMM_DoubleArray_create(9);
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom atom = atoms[i];
MultipoleType multipoleType = atom.getMultipoleType();
PolarizeType polarType = atom.getPolarizeType();
double useFactor = 1.0;
if (!atoms[i].getUse() || !atoms[i].getElectrostatics()) {
// if (!atoms[i].getUse()) {
useFactor = 0.0;
}
double lambdaScale = lambda;
if (!atom.applyLambda()) {
lambdaScale = 1.0;
}
useFactor *= lambdaScale;
/**
* Define the frame definition.
*/
int axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
switch(multipoleType.frameDefinition) {
case ZONLY:
axisType = OpenMM_AmoebaMultipoleForce_ZOnly;
break;
case ZTHENX:
axisType = OpenMM_AmoebaMultipoleForce_ZThenX;
break;
case BISECTOR:
axisType = OpenMM_AmoebaMultipoleForce_Bisector;
break;
case ZTHENBISECTOR:
axisType = OpenMM_AmoebaMultipoleForce_ZBisect;
break;
case TRISECTOR:
axisType = OpenMM_AmoebaMultipoleForce_ThreeFold;
break;
default:
break;
}
/**
* Load local multipole coefficients.
*/
for (int j = 0; j < 3; j++) {
OpenMM_DoubleArray_set(dipoles, j, multipoleType.dipole[j] * dipoleConversion * useFactor);
}
int l = 0;
for (int j = 0; j < 3; j++) {
for (int k = 0; k < 3; k++) {
OpenMM_DoubleArray_set(quadrupoles, l++, multipoleType.quadrupole[j][k] * quadrupoleConversion / 3.0 * useFactor);
}
}
// int zaxis = 0;
int zaxis = 1;
// int xaxis = 0;
int xaxis = 1;
// int yaxis = 0;
int yaxis = 1;
int[] refAtoms = axisAtom[i];
if (refAtoms != null) {
zaxis = refAtoms[0];
if (refAtoms.length > 1) {
xaxis = refAtoms[1];
if (refAtoms.length > 2) {
yaxis = refAtoms[2];
}
}
} else {
axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
}
/**
* Add the multipole.
*/
OpenMM_AmoebaMultipoleForce_setMultipoleParameters(amoebaMultipoleForce, i, multipoleType.charge * useFactor, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis, polarType.thole, polarType.pdamp * dampingFactorConversion, polarType.polarizability * polarityConversion * polarScale * useFactor);
}
OpenMM_DoubleArray_destroy(dipoles);
OpenMM_DoubleArray_destroy(quadrupoles);
OpenMM_AmoebaMultipoleForce_updateParametersInContext(amoebaMultipoleForce, context);
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method fillGradients.
/**
* Private method for internal use, so we don't have subclasses calling
* super.energy, and this class delegating to the subclass's getGradients
* method.
*
* @param g Gradient array to fill.
* @return Gradient array.
*/
public double[] fillGradients(double[] g) {
assert (g != null);
int n = getNumberOfVariables();
if (g.length < n) {
g = new double[n];
}
int index = 0;
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom a = atoms[i];
if (a.isActive()) {
OpenMM_Vec3 posInNm = OpenMM_Vec3Array_get(forces, i);
/**
* Convert OpenMM Forces in KJ/Nm into an FFX gradient in
* Kcal/A.
*/
double gx = -posInNm.x * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
double gy = -posInNm.y * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
double gz = -posInNm.z * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
if (Double.isNaN(gx) || Double.isInfinite(gx) || Double.isNaN(gy) || Double.isInfinite(gy) || Double.isNaN(gz) || Double.isInfinite(gz)) {
/*String message = format("The gradient of atom %s is (%8.3f,%8.3f,%8.3f).",
a.toString(), gx, gy, gz);*/
StringBuilder sb = new StringBuilder(format("The gradient of atom %s is (%8.3f,%8.3f,%8.3f).", a.toString(), gx, gy, gz));
double[] vals = new double[3];
a.getVelocity(vals);
sb.append(format("\n Velocities: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
a.getAcceleration(vals);
sb.append(format("\n Accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
a.getPreviousAcceleration(vals);
sb.append(format("\n Previous accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
// logger.severe(message);
throw new EnergyException(sb.toString());
}
a.setXYZGradient(gx, gy, gz);
g[index++] = gx;
g[index++] = gy;
g[index++] = gz;
}
}
return g;
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method getOpenMMAccelerations.
/**
* getOpenMMAccelerations takes in a PointerByReference containing the force
* information of the context. This method creates a Vec3Array that contains
* the three dimensional information of the forces on the atoms. This method
* then adds these values (divided by mass, effectively turning them into
* accelerations) to a new double array a and returns it to the method call
*
* @param forces
* @param numberOfVariables
* @param mass
* @return
*/
public double[] getOpenMMAccelerations(PointerByReference forces, int numberOfVariables, double[] mass, double[] a) {
assert numberOfVariables == getNumberOfVariables();
if (a == null || a.length < numberOfVariables) {
a = new double[numberOfVariables];
}
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
int offset = i * 3;
OpenMM_Vec3 acc = OpenMM_Vec3Array_get(forces, i);
a[offset] = (acc.x * 10.0) / mass[i];
a[offset + 1] = (acc.y * 10.0) / mass[i + 1];
a[offset + 2] = (acc.z * 10.0) / mass[i + 2];
Atom atom = atoms[i];
double[] acceleration = { a[offset], a[offset + 1], a[offset + 2] };
atom.setAcceleration(acceleration);
}
return a;
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method addAmoebaVDWForce.
private void addAmoebaVDWForce() {
VanDerWaals vdW = super.getVdwNode();
if (vdW == null) {
return;
}
amoebaVDWForce = OpenMM_AmoebaVdwForce_create();
OpenMM_System_addForce(system, amoebaVDWForce);
OpenMM_Force_setForceGroup(amoebaVDWForce, 1);
VanDerWaalsForm vdwForm = vdW.getVDWForm();
NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
Crystal crystal = super.getCrystal();
double radScale = 1.0;
if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
radScale = 0.5;
}
/**
* Note that the API says it wants a SIGMA value.
*/
if (vdwForm.radiusType == VanDerWaalsForm.RADIUS_TYPE.R_MIN) {
// radScale *= 1.122462048309372981;
}
int[] ired = vdW.getReductionIndex();
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom atom = atoms[i];
VDWType vdwType = atom.getVDWType();
OpenMM_AmoebaVdwForce_addParticle(amoebaVDWForce, ired[i], OpenMM_NmPerAngstrom * vdwType.radius * radScale, OpenMM_KJPerKcal * vdwType.wellDepth, vdwType.reductionFactor);
}
// OpenMM_AmoebaVdwForce_setSigmaCombiningRule(amoebaVdwForce, toPropertyForm(vdwForm.radiusRule.name()));
// OpenMM_AmoebaVdwForce_setEpsilonCombiningRule(amoebaVdwForce, toPropertyForm(vdwForm.epsilonRule.name()));
OpenMM_AmoebaVdwForce_setCutoffDistance(amoebaVDWForce, nonbondedCutoff.off * OpenMM_NmPerAngstrom);
OpenMM_AmoebaVdwForce_setUseDispersionCorrection(amoebaVDWForce, OpenMM_Boolean.OpenMM_False);
if (crystal.aperiodic()) {
OpenMM_AmoebaVdwForce_setNonbondedMethod(amoebaVDWForce, OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_NoCutoff);
} else {
OpenMM_AmoebaVdwForce_setNonbondedMethod(amoebaVDWForce, OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_CutoffPeriodic);
}
/**
* Create exclusion lists.
*/
PointerByReference exclusions = OpenMM_IntArray_create(0);
double[] mask = new double[nAtoms];
Arrays.fill(mask, 1.0);
for (int i = 0; i < nAtoms; i++) {
OpenMM_IntArray_append(exclusions, i);
vdW.applyMask(mask, i);
for (int j = 0; j < nAtoms; j++) {
if (mask[j] == 0.0) {
OpenMM_IntArray_append(exclusions, j);
}
}
vdW.removeMask(mask, i);
OpenMM_AmoebaVdwForce_setParticleExclusions(amoebaVDWForce, i, exclusions);
OpenMM_IntArray_resize(exclusions, 0);
}
OpenMM_IntArray_destroy(exclusions);
logger.log(Level.INFO, " Added van der Waals force.");
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method addCustomNonbondedSoftcoreForce.
/**
* 1.) Handle interactions between non-alchemical atoms with our default
* OpenMM NonBondedForce. Note that alchemical atoms must have eps=0 to turn
* them off in this force.
* <p>
* 2.) Handle interactions between alchemical atoms and mixed non-alchemical
* <-> alchemical interactions with an OpenMM CustomNonBondedForce.
*/
private void addCustomNonbondedSoftcoreForce() {
VanDerWaals vdW = super.getVdwNode();
if (vdW == null) {
return;
}
/**
* Only 6-12 LJ with arithmetic mean to define sigma and geometric mean
* for epsilon is supported.
*/
VanDerWaalsForm vdwForm = vdW.getVDWForm();
if (vdwForm.vdwType != LENNARD_JONES || vdwForm.radiusRule != ARITHMETIC || vdwForm.epsilonRule != GEOMETRIC) {
logger.info(format(" VDW Type: %s", vdwForm.vdwType));
logger.info(format(" VDW Radius Rule: %s", vdwForm.radiusRule));
logger.info(format(" VDW Epsilon Rule: %s", vdwForm.epsilonRule));
logger.log(Level.SEVERE, String.format(" Unsuppporterd van der Waals functional form."));
return;
}
// Sterics mixing rules.
String stericsMixingRules = " epsilon = sqrt(epsilon1*epsilon2);";
stericsMixingRules += " rmin = 0.5 * (sigma1 + sigma2) * 1.122462048309372981;";
// Softcore Lennard-Jones, with a form equivalent to that used in FFX VanDerWaals class.
String stericsEnergyExpression = "(vdw_lambda^beta)*epsilon*x*(x-2.0);";
// Effective softcore distance for sterics.
stericsEnergyExpression += " x = 1.0 / (alpha*(1.0-vdw_lambda)^2.0 + (r/rmin)^6.0);";
// Define energy expression for sterics.
String energyExpression = stericsEnergyExpression + stericsMixingRules;
fixedChargeSoftcore = OpenMM_CustomNonbondedForce_create(energyExpression);
// Get the Alpha and Beta constants from the VanDerWaals instance.
double alpha = vdW.getAlpha();
double beta = vdW.getBeta();
logger.info(format(" Custom non-bonded force with alpha = %8.6f and beta = %8.6f", alpha, beta));
OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "vdw_lambda", 1.0);
OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "alpha", alpha);
OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "beta", beta);
OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "sigma");
OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "epsilon");
/**
* Add particles.
*/
PointerByReference alchemicalGroup = OpenMM_IntSet_create();
PointerByReference nonAlchemicalGroup = OpenMM_IntSet_create();
DoubleByReference charge = new DoubleByReference();
DoubleByReference sigma = new DoubleByReference();
DoubleByReference eps = new DoubleByReference();
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom atom = atoms[i];
OpenMM_NonbondedForce_getParticleParameters(fixedChargeNonBondedForce, i, charge, sigma, eps);
if (atom.applyLambda()) {
OpenMM_IntSet_insert(alchemicalGroup, i);
logger.info(format(" Adding alchemical atom %s.", atom));
} else {
OpenMM_IntSet_insert(nonAlchemicalGroup, i);
}
PointerByReference particleParameters = OpenMM_DoubleArray_create(0);
OpenMM_DoubleArray_append(particleParameters, sigma.getValue());
OpenMM_DoubleArray_append(particleParameters, eps.getValue());
OpenMM_CustomNonbondedForce_addParticle(fixedChargeSoftcore, particleParameters);
OpenMM_DoubleArray_destroy(particleParameters);
}
OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup, alchemicalGroup);
OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup, nonAlchemicalGroup);
OpenMM_IntSet_destroy(alchemicalGroup);
OpenMM_IntSet_destroy(nonAlchemicalGroup);
Crystal crystal = super.getCrystal();
if (crystal.aperiodic()) {
OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore, OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_NoCutoff);
} else {
OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore, OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_CutoffPeriodic);
}
NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
double off = nonbondedCutoff.off;
double cut = nonbondedCutoff.cut;
OpenMM_CustomNonbondedForce_setCutoffDistance(fixedChargeSoftcore, OpenMM_NmPerAngstrom * off);
OpenMM_CustomNonbondedForce_setUseSwitchingFunction(fixedChargeSoftcore, OpenMM_True);
OpenMM_CustomNonbondedForce_setSwitchingDistance(fixedChargeSoftcore, OpenMM_NmPerAngstrom * cut);
if (cut == off) {
logger.warning(" OpenMM does not properly handle cutoffs where cut == off!");
if (cut == Double.MAX_VALUE || cut == Double.POSITIVE_INFINITY) {
logger.info(" Detected infinite or max-value cutoff; setting cut to 1E+40 for OpenMM.");
cut = 1E40;
} else {
logger.info(String.format(" Detected cut %8.4g == off %8.4g; scaling cut to 0.99 of off for OpenMM.", cut, off));
cut *= 0.99;
}
}
// Add energy parameter derivative
OpenMM_CustomNonbondedForce_addEnergyParameterDerivative(fixedChargeSoftcore, "vdw_lambda");
OpenMM_System_addForce(system, fixedChargeSoftcore);
logger.log(Level.INFO, String.format(" Added fixed charge softcore sterics force."));
GeneralizedKirkwood gk = super.getGK();
if (gk != null) {
logger.severe(" OpenMM alchemical methods are not supported for GB.");
addCustomGBForce();
}
// Not entirely sure how to initialize this portion
alchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(stericsEnergyExpression);
nonAlchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(stericsEnergyExpression);
// allStericsForce = (alchemicalAlchemicalStericsForce + nonAlchemicalAlchemicalStericsForce);
// Can be reduced to two lines if I can figure out how to combine the two custom bonded sterics forces
OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "rmin");
OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "epsilon");
OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "vdw_lambda", 1.0);
OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "alpha", alpha);
OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "beta", beta);
OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "rmin");
OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "epsilon");
OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "vdw_lambda", 1.0);
OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "alpha", alpha);
OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "beta", beta);
int range = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
IntByReference atomi = new IntByReference();
IntByReference atomj = new IntByReference();
int[][] torsionMask = vdW.getTorsionMask();
for (int i = 0; i < range; i++) {
OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, atomi, atomj, charge, sigma, eps);
OpenMM_CustomNonbondedForce_addExclusion(fixedChargeSoftcore, atomi.getValue(), atomj.getValue());
int[] maskI = torsionMask[atomi.getValue()];
int jID = atomj.getValue();
boolean epsException = false;
for (int j = 0; j < maskI.length; j++) {
if (maskI[j] == jID) {
epsException = true;
break;
}
}
Atom atom1 = atoms[atomi.getValue()];
Atom atom2 = atoms[atomj.getValue()];
boolean bothAlchemical = false;
boolean oneAlchemical = false;
if (atom1.applyLambda() && atom2.applyLambda()) {
bothAlchemical = true;
} else if ((atom1.applyLambda() && !atom2.applyLambda()) || (!atom1.applyLambda() && atom2.applyLambda())) {
oneAlchemical = true;
}
if (bothAlchemical) {
if (epsException) {
PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
OpenMM_DoubleArray_append(bondParameters, eps.getValue());
OpenMM_CustomBondForce_addBond(alchemicalAlchemicalStericsForce, atomi.getValue(), atomj.getValue(), bondParameters);
OpenMM_DoubleArray_destroy(bondParameters);
}
} else if (oneAlchemical) {
if (epsException) {
PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
OpenMM_DoubleArray_append(bondParameters, eps.getValue());
OpenMM_CustomBondForce_addBond(nonAlchemicalAlchemicalStericsForce, atomi.getValue(), atomj.getValue(), bondParameters);
OpenMM_DoubleArray_destroy(bondParameters);
}
}
}
/**
* for (int i = 0; i < range; i++){
* OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, atomi, atomj, charge, sigma, eps);
*
* Atom atom1 = atoms[atomi.getValue()];
* Atom atom2 = atoms[atomj.getValue()];
*
* if (atom1.applyLambda() || atom2.applyLambda()){
* OpenMM_NonbondedForce_setExceptionParameters(fixedChargeNonBondedForce, i, atomi.getValue(), atomj.getValue(), abs(0.0*charge.getValue()), sigma.getValue(), abs(0.0*eps.getValue()));
* }
* }
*/
OpenMM_CustomBondForce_addEnergyParameterDerivative(alchemicalAlchemicalStericsForce, "vdw_lambda");
OpenMM_CustomBondForce_addEnergyParameterDerivative(nonAlchemicalAlchemicalStericsForce, "vdw_lambda");
OpenMM_System_addForce(system, alchemicalAlchemicalStericsForce);
OpenMM_System_addForce(system, nonAlchemicalAlchemicalStericsForce);
}
Aggregations