use of ffx.potential.nonbonded.VanDerWaalsForm in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method updateFixedChargeNonBondedForce.
/**
* Updates the fixed-charge non-bonded force for change in Use flags.
*
* @param atoms Array of all Atoms in the system
* @param vdwLambdaTerm If true, set charges and eps values to zero for
* Lambda atoms
*/
private void updateFixedChargeNonBondedForce(Atom[] atoms, boolean vdwLambdaTerm) {
VanDerWaals vdW = super.getVdwNode();
/**
* 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.log(Level.SEVERE, String.format(" Unsuppporterd van der Waals functional form."));
return;
}
/**
* OpenMM vdW force requires a diameter (i.e. not radius).
*/
double radScale = 1.0;
if (vdwForm.radiusSize == RADIUS) {
radScale = 2.0;
}
/**
* OpenMM vdw force requires atomic sigma values (i.e. not r-min).
*/
if (vdwForm.radiusType == R_MIN) {
radScale /= 1.122462048309372981;
}
/**
* Update parameters.
*/
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom atom = atoms[i];
boolean applyLambda = atom.applyLambda();
double charge = Double.MIN_VALUE;
MultipoleType multipoleType = atom.getMultipoleType();
if (multipoleType != null && atoms[i].getElectrostatics()) {
charge = multipoleType.charge;
if (lambdaTerm && applyLambda) {
charge *= lambda;
}
}
VDWType vdwType = atom.getVDWType();
double sigma = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
if ((vdwLambdaTerm && applyLambda) || !atoms[i].getUse()) {
eps = 0.0;
charge = 0.0;
}
OpenMM_NonbondedForce_setParticleParameters(fixedChargeNonBondedForce, i, charge, sigma, eps);
}
// OpenMM_NonbondedForce_updateParametersInContext(fixedChargeNonBondedForce, context);
/**
* Update Exceptions.
*/
IntByReference particle1 = new IntByReference();
IntByReference particle2 = new IntByReference();
DoubleByReference chargeProd = new DoubleByReference();
DoubleByReference sigma = new DoubleByReference();
DoubleByReference eps = new DoubleByReference();
int numExceptions = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
for (int i = 0; i < numExceptions; i++) {
/**
* Only update exceptions.
*/
if (chargeExclusion[i] && vdWExclusion[i]) {
continue;
}
OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, particle1, particle2, chargeProd, sigma, eps);
int i1 = particle1.getValue();
int i2 = particle2.getValue();
double qq = exceptionChargeProd[i];
double epsilon = exceptionEps[i];
Atom atom1 = atoms[i1];
Atom atom2 = atoms[i2];
double lambdaValue = lambda;
if (lambda == 0.0) {
lambdaValue = 1.0e-6;
}
if (atom1.applyLambda()) {
qq *= lambdaValue;
if (vdwLambdaTerm) {
epsilon = 1.0e-6;
qq = 1.0e-6;
}
}
if (atom2.applyLambda()) {
qq *= lambdaValue;
if (vdwLambdaTerm) {
epsilon = 1.0e-6;
qq = 1.0e-6;
}
}
if (!atom1.getUse() || !atom2.getUse()) {
qq = 1.0e-6;
epsilon = 1.0e-6;
}
OpenMM_NonbondedForce_setExceptionParameters(fixedChargeNonBondedForce, i, i1, i2, qq, sigma.getValue(), epsilon);
/**
* logger.info(format(" B Exception %d %d %d q=%10.8f s=%10.8f
* e=%10.8f.", i, i1, i2, chargeProd.getValue(), sigma.getValue(),
* eps.getValue()));
*
* logger.info(format(" E Exception %d %d %d q=%10.8f s=%10.8f
* e=%10.8f.", i, i1, i2, qq, sigma.getValue(), epsilon));
*/
}
OpenMM_NonbondedForce_updateParametersInContext(fixedChargeNonBondedForce, context);
}
use of ffx.potential.nonbonded.VanDerWaalsForm in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method updateAmoebaVDWForce.
/**
* Updates the AMOEBA van der Waals force for change in Use flags.
*
* @param atoms Array of all Atoms in the system
*/
private void updateAmoebaVDWForce(Atom[] atoms) {
VanDerWaals vdW = super.getVdwNode();
VanDerWaalsForm vdwForm = vdW.getVDWForm();
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();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom atom = atoms[i];
VDWType vdwType = atom.getVDWType();
double useFactor = 1.0;
if (!atoms[i].getUse()) {
useFactor = 0.0;
}
double eps = OpenMM_KJPerKcal * vdwType.wellDepth * useFactor;
OpenMM_AmoebaVdwForce_setParticleParameters(amoebaVDWForce, i, ired[i], OpenMM_NmPerAngstrom * vdwType.radius * radScale, eps, vdwType.reductionFactor);
}
OpenMM_AmoebaVdwForce_updateParametersInContext(amoebaVDWForce, context);
}
use of ffx.potential.nonbonded.VanDerWaalsForm in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method getVDWRadius.
public double getVDWRadius() {
VanDerWaals vdW = super.getVdwNode();
if (vdW == null) {
return -1.0;
}
/**
* 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.log(Level.SEVERE, String.format(" Unsuppporterd van der Waals functional form."));
return -1.0;
}
/**
* OpenMM vdW force requires a diameter (i.e. not radius).
*/
double radScale = 1.0;
if (vdwForm.radiusSize == RADIUS) {
radScale = 2.0;
}
/**
* OpenMM vdw force requires atomic sigma values (i.e. not r-min).
*/
if (vdwForm.radiusType == R_MIN) {
radScale /= 1.122462048309372981;
}
return radScale;
}
use of ffx.potential.nonbonded.VanDerWaalsForm 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.nonbonded.VanDerWaalsForm 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