Examples with VDWType ffx.potential.parameters.VDWType used on opensource projects

Example 6 with VDWType

use of ffx.potential.parameters.VDWType in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addWCAForce.

private void addWCAForce() {
    double epso = 0.1100;
    double epsh = 0.0135;
    double rmino = 1.7025;
    double rminh = 1.3275;
    double awater = 0.033428;
    double slevy = 1.0;
    double dispoff = 0.26;
    double shctd = 0.81;
    VanDerWaals vdW = super.getVdwNode();
    VanDerWaalsForm vdwForm = vdW.getVDWForm();
    double radScale = 1.0;
    if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
        radScale = 0.5;
    }
    amoebaWcaDispersionForce = OpenMM_AmoebaWcaDispersionForce_create();
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        // cdispTotal += nonpol__.cdisp[ii];
        Atom atom = atoms[i];
        VDWType vdwType = atom.getVDWType();
        double radius = vdwType.radius;
        double eps = vdwType.wellDepth;
        OpenMM_AmoebaWcaDispersionForce_addParticle(amoebaWcaDispersionForce, OpenMM_NmPerAngstrom * radius * radScale, OpenMM_KJPerKcal * eps);
    }
    OpenMM_AmoebaWcaDispersionForce_setEpso(amoebaWcaDispersionForce, epso * OpenMM_KJPerKcal);
    OpenMM_AmoebaWcaDispersionForce_setEpsh(amoebaWcaDispersionForce, epsh * OpenMM_KJPerKcal);
    OpenMM_AmoebaWcaDispersionForce_setRmino(amoebaWcaDispersionForce, rmino * OpenMM_NmPerAngstrom);
    OpenMM_AmoebaWcaDispersionForce_setRminh(amoebaWcaDispersionForce, rminh * OpenMM_NmPerAngstrom);
    OpenMM_AmoebaWcaDispersionForce_setDispoff(amoebaWcaDispersionForce, dispoff * OpenMM_NmPerAngstrom);
    OpenMM_AmoebaWcaDispersionForce_setAwater(amoebaWcaDispersionForce, awater / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
    OpenMM_AmoebaWcaDispersionForce_setSlevy(amoebaWcaDispersionForce, slevy);
    OpenMM_AmoebaWcaDispersionForce_setShctd(amoebaWcaDispersionForce, shctd);
    OpenMM_System_addForce(system, amoebaWcaDispersionForce);
    logger.log(Level.INFO, " Added WCA dispersion force.");
}

Example 7 with VDWType

use of ffx.potential.parameters.VDWType in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method updateWCAForce.

/**
 * Updates the WCA force for change in Use flags.
 *
 * @param atoms Array of all Atoms in the system
 */
private void updateWCAForce(Atom[] atoms) {
    VanDerWaals vdW = super.getVdwNode();
    VanDerWaalsForm vdwForm = vdW.getVDWForm();
    double radScale = 1.0;
    if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
        radScale = 0.5;
    }
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        double useFactor = 1.0;
        if (!atoms[i].getUse()) {
            useFactor = 0.0;
        }
        double lambdaScale = lambda;
        if (!atoms[i].applyLambda()) {
            lambdaScale = 1.0;
        }
        useFactor *= lambdaScale;
        Atom atom = atoms[i];
        VDWType vdwType = atom.getVDWType();
        double radius = vdwType.radius;
        double eps = vdwType.wellDepth;
        OpenMM_AmoebaWcaDispersionForce_setParticleParameters(amoebaWcaDispersionForce, i, OpenMM_NmPerAngstrom * radius * radScale, OpenMM_KJPerKcal * eps * useFactor);
    }
    OpenMM_AmoebaWcaDispersionForce_updateParametersInContext(amoebaWcaDispersionForce, context);
}

Example 8 with VDWType

use of ffx.potential.parameters.VDWType 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();
    }
}

Example 9 with VDWType

use of ffx.potential.parameters.VDWType in project ffx by mjschnie.

the class VanDerWaals method initAtomArrays.

/**
 * Allocate coordinate arrays and set up reduction indices and values.
 */
private void initAtomArrays() {
    if (esvTerm) {
        atoms = esvSystem.getExtendedAtoms();
        nAtoms = atoms.length;
    }
    if (atomClass == null || nAtoms > atomClass.length || lambdaTerm || esvTerm) {
        atomClass = new int[nAtoms];
        coordinates = new double[nAtoms * 3];
        reduced = new double[nSymm][nAtoms * 3];
        reducedXYZ = reduced[0];
        reductionIndex = new int[nAtoms];
        reductionValue = new double[nAtoms];
        bondMask = new int[nAtoms][];
        angleMask = new int[nAtoms][];
        if (vdwForm.vdwType == VanDerWaalsForm.VDW_TYPE.LENNARD_JONES) {
            torsionMask = new int[nAtoms][];
        } else {
            torsionMask = null;
        }
        use = new boolean[nAtoms];
        isSoft = new boolean[nAtoms];
        softCore = new boolean[2][nAtoms];
        lambdaGradX = null;
        lambdaGradY = null;
        lambdaGradZ = null;
        switch(atomicDoubleArrayImpl) {
            case MULTI:
                gradX = new MultiDoubleArray(threadCount, nAtoms);
                gradY = new MultiDoubleArray(threadCount, nAtoms);
                gradZ = new MultiDoubleArray(threadCount, nAtoms);
                if (lambdaTerm) {
                    lambdaGradX = new MultiDoubleArray(threadCount, nAtoms);
                    lambdaGradY = new MultiDoubleArray(threadCount, nAtoms);
                    lambdaGradZ = new MultiDoubleArray(threadCount, nAtoms);
                }
                break;
            case PJ:
                gradX = new PJDoubleArray(threadCount, nAtoms);
                gradY = new PJDoubleArray(threadCount, nAtoms);
                gradZ = new PJDoubleArray(threadCount, nAtoms);
                if (lambdaTerm) {
                    lambdaGradX = new PJDoubleArray(threadCount, nAtoms);
                    lambdaGradY = new PJDoubleArray(threadCount, nAtoms);
                    lambdaGradZ = new PJDoubleArray(threadCount, nAtoms);
                }
                break;
            case ADDER:
            default:
                gradX = new AdderDoubleArray(threadCount, nAtoms);
                gradY = new AdderDoubleArray(threadCount, nAtoms);
                gradZ = new AdderDoubleArray(threadCount, nAtoms);
                if (lambdaTerm) {
                    lambdaGradX = new AdderDoubleArray(threadCount, nAtoms);
                    lambdaGradY = new AdderDoubleArray(threadCount, nAtoms);
                    lambdaGradZ = new AdderDoubleArray(threadCount, nAtoms);
                }
                break;
        }
    }
    /**
     * Initialize all atoms to be used in the energy.
     */
    fill(use, true);
    fill(isSoft, false);
    fill(softCore[HARD], false);
    fill(softCore[SOFT], false);
    softCoreInit = false;
    // Needs initialized regardless of esvTerm.
    esvAtoms = new boolean[nAtoms];
    esvLambda = new double[nAtoms];
    atomEsvID = new int[nAtoms];
    fill(esvAtoms, false);
    fill(esvLambda, 1.0);
    fill(atomEsvID, -1);
    if (esvTerm) {
        updateEsvLambda();
    }
    lambdaFactors = new LambdaFactors[threadCount];
    for (int i = 0; i < threadCount; i++) {
        if (esvTerm) {
            lambdaFactors[i] = new LambdaFactorsESV();
        } else if (lambdaTerm) {
            lambdaFactors[i] = new LambdaFactorsOSRW();
        } else {
            lambdaFactors[i] = new LambdaFactors();
        }
    }
    for (int i = 0; i < nAtoms; i++) {
        Atom ai = atoms[i];
        assert (i == ai.getXyzIndex() - 1);
        double[] xyz = ai.getXYZ(null);
        int i3 = i * 3;
        coordinates[i3 + XX] = xyz[XX];
        coordinates[i3 + YY] = xyz[YY];
        coordinates[i3 + ZZ] = xyz[ZZ];
        AtomType atomType = ai.getAtomType();
        if (atomType == null) {
            logger.severe(ai.toString());
            // Severe no longer guarantees program crash.
            continue;
        }
        String vdwIndex = forceField.getString(ForceField.ForceFieldString.VDWINDEX, "Class");
        if (vdwIndex.equalsIgnoreCase("Type")) {
            atomClass[i] = atomType.type;
        } else {
            atomClass[i] = atomType.atomClass;
        }
        VDWType type = forceField.getVDWType(Integer.toString(atomClass[i]));
        if (type == null) {
            logger.info(" No VdW type for atom class " + atomClass[i]);
            logger.severe(" No VdW type for atom " + ai.toString());
            return;
        }
        ai.setVDWType(type);
        ArrayList<Bond> bonds = ai.getBonds();
        int numBonds = bonds.size();
        if (type.reductionFactor > 0.0 && numBonds == 1) {
            Bond bond = bonds.get(0);
            Atom heavyAtom = bond.get1_2(ai);
            // Atom indexes start at 1
            reductionIndex[i] = heavyAtom.getIndex() - 1;
            reductionValue[i] = type.reductionFactor;
        } else {
            reductionIndex[i] = i;
            reductionValue[i] = 0.0;
        }
        bondMask[i] = new int[numBonds];
        for (int j = 0; j < numBonds; j++) {
            Bond bond = bonds.get(j);
            bondMask[i][j] = bond.get1_2(ai).getIndex() - 1;
        }
        ArrayList<Angle> angles = ai.getAngles();
        int numAngles = 0;
        for (Angle angle : angles) {
            Atom ak = angle.get1_3(ai);
            if (ak != null) {
                numAngles++;
            }
        }
        angleMask[i] = new int[numAngles];
        int j = 0;
        for (Angle angle : angles) {
            Atom ak = angle.get1_3(ai);
            if (ak != null) {
                angleMask[i][j++] = ak.getIndex() - 1;
            }
        }
        if (vdwForm.scale14 != 1.0) {
            ArrayList<Torsion> torsions = ai.getTorsions();
            int numTorsions = 0;
            for (Torsion torsion : torsions) {
                Atom ak = torsion.get1_4(ai);
                if (ak != null) {
                    numTorsions++;
                }
            }
            torsionMask[i] = new int[numTorsions];
            j = 0;
            for (Torsion torsion : torsions) {
                Atom ak = torsion.get1_4(ai);
                if (ak != null) {
                    torsionMask[i][j++] = ak.getIndex() - 1;
                }
            }
        }
    }
}