Examples with NonbondedCutoff ffx.potential.nonbonded.NonbondedCutoff used on opensource projects

Example 1 with NonbondedCutoff

use of ffx.potential.nonbonded.NonbondedCutoff 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.");
}

Example 2 with NonbondedCutoff

use of ffx.potential.nonbonded.NonbondedCutoff 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);
}

Example 3 with NonbondedCutoff

use of ffx.potential.nonbonded.NonbondedCutoff in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addFixedChargeNonBondedForce.

/**
 * Uses arithmetic mean to define sigma and geometric mean for epsilon.
 */
private void addFixedChargeNonBondedForce() {
    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;
    }
    fixedChargeNonBondedForce = OpenMM_NonbondedForce_create();
    /**
     * 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;
    }
    /**
     * Add particles.
     */
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        VDWType vdwType = atom.getVDWType();
        double sigma = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
        double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
        double charge = 0.0;
        MultipoleType multipoleType = atom.getMultipoleType();
        if (multipoleType != null && atoms[i].getElectrostatics()) {
            charge = multipoleType.charge;
        }
        OpenMM_NonbondedForce_addParticle(fixedChargeNonBondedForce, charge, sigma, eps);
    }
    /**
     * Define 1-4 scale factors.
     */
    double lj14Scale = vdwForm.getScale14();
    double coulomb14Scale = 1.0 / 1.2;
    ParticleMeshEwald pme = super.getPmeNode();
    Bond[] bonds = super.getBonds();
    if (bonds != null && bonds.length > 0) {
        int nBonds = bonds.length;
        PointerByReference bondArray;
        bondArray = OpenMM_BondArray_create(0);
        for (int i = 0; i < nBonds; i++) {
            Bond bond = bonds[i];
            int i1 = bond.getAtom(0).getXyzIndex() - 1;
            int i2 = bond.getAtom(1).getXyzIndex() - 1;
            OpenMM_BondArray_append(bondArray, i1, i2);
        }
        if (pme != null) {
            coulomb14Scale = pme.getScale14();
        }
        OpenMM_NonbondedForce_createExceptionsFromBonds(fixedChargeNonBondedForce, bondArray, coulomb14Scale, lj14Scale);
        OpenMM_BondArray_destroy(bondArray);
        int num = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
        chargeExclusion = new boolean[num];
        vdWExclusion = new boolean[num];
        exceptionChargeProd = new double[num];
        exceptionEps = new double[num];
        IntByReference particle1 = new IntByReference();
        IntByReference particle2 = new IntByReference();
        DoubleByReference chargeProd = new DoubleByReference();
        DoubleByReference sigma = new DoubleByReference();
        DoubleByReference eps = new DoubleByReference();
        for (int i = 0; i < num; i++) {
            OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, particle1, particle2, chargeProd, sigma, eps);
            if (abs(chargeProd.getValue()) > 0.0) {
                chargeExclusion[i] = false;
                exceptionChargeProd[i] = chargeProd.getValue();
            } else {
                exceptionChargeProd[i] = 0.0;
                chargeExclusion[i] = true;
            }
            if (abs(eps.getValue()) > 0.0) {
                vdWExclusion[i] = false;
                exceptionEps[i] = eps.getValue();
            } else {
                vdWExclusion[i] = true;
                exceptionEps[i] = 0.0;
            }
        }
    }
    Crystal crystal = super.getCrystal();
    if (crystal.aperiodic()) {
        OpenMM_NonbondedForce_setNonbondedMethod(fixedChargeNonBondedForce, OpenMM_NonbondedForce_NonbondedMethod.OpenMM_NonbondedForce_NoCutoff);
    } else {
        OpenMM_NonbondedForce_setNonbondedMethod(fixedChargeNonBondedForce, OpenMM_NonbondedForce_NonbondedMethod.OpenMM_NonbondedForce_PME);
        if (pme != null) {
            // Units of the Ewald coefficient are A^-1; Multiply by AngstromsPerNM to convert to (Nm^-1).
            double aEwald = OpenMM_AngstromsPerNm * pme.getEwaldCoefficient();
            int nx = pme.getReciprocalSpace().getXDim();
            int ny = pme.getReciprocalSpace().getYDim();
            int nz = pme.getReciprocalSpace().getZDim();
            OpenMM_NonbondedForce_setPMEParameters(fixedChargeNonBondedForce, aEwald, nx, ny, nz);
        }
        NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
        double off = nonbondedCutoff.off;
        double cut = nonbondedCutoff.cut;
        OpenMM_NonbondedForce_setCutoffDistance(fixedChargeNonBondedForce, OpenMM_NmPerAngstrom * off);
        OpenMM_NonbondedForce_setUseSwitchingFunction(fixedChargeNonBondedForce, OpenMM_True);
        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;
            }
        }
        OpenMM_NonbondedForce_setSwitchingDistance(fixedChargeNonBondedForce, OpenMM_NmPerAngstrom * cut);
    }
    OpenMM_NonbondedForce_setUseDispersionCorrection(fixedChargeNonBondedForce, OpenMM_False);
    // OpenMM_Force_setForceGroup(fixedChargeNonBondedForce, 1);
    OpenMM_System_addForce(system, fixedChargeNonBondedForce);
    logger.log(Level.INFO, String.format(" Added fixed charge non-bonded force."));
    GeneralizedKirkwood gk = super.getGK();
    if (gk != null) {
        addCustomGBForce();
    }
}