Examples with Bond ffx.potential.bonded.Bond used on opensource projects

Example 16 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class ForceFieldEnergy method reInit.

/**
 * Need to remove degree of freedom that are lost to prevent heating.
 */
public void reInit() {
    int[] molecule;
    if (esvTerm) {
        atoms = esvSystem.getExtendedAndBackgroundAtoms();
        molecule = esvSystem.getExtendedAndBackgroundMolecule();
    } else {
        atoms = molecularAssembly.getAtomArray();
        molecule = molecularAssembly.getMoleculeNumbers();
    }
    nAtoms = atoms.length;
    /* TODO Decide on only growing vs. always modifying xyz.
        if (xyz.length < 3 * nAtoms) {
            xyz = new double[nAtoms * 3];
        }   */
    xyz = new double[nAtoms * 3];
    getCoordinates(xyz);
    // Check that atom ordering is correct and count number of Active atoms.
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        int index = atom.getXyzIndex() - 1;
        if (index != i) {
            atom.setXyzIndex(i + 1);
        }
    }
    // Collect, count, pack and sort bonds.
    if (bondTerm) {
        ArrayList<ROLS> bond = molecularAssembly.getBondList();
        nBonds = 0;
        for (ROLS r : bond) {
            if (keep((Bond) r)) {
                nBonds++;
            }
        }
        if (nBonds > bonds.length) {
            bonds = new Bond[nBonds];
        }
        Arrays.fill(bonds, null);
        nBonds = 0;
        for (ROLS r : bond) {
            if (keep((Bond) r)) {
                bonds[nBonds++] = (Bond) r;
            }
        }
        Arrays.sort(bonds, 0, nBonds);
        if (nBonds > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Bonds:                             %10d", nBonds));
        }
    } else {
        nBonds = 0;
        bonds = null;
    }
    // Collect, count, pack and sort angles.
    if (angleTerm) {
        ArrayList<ROLS> angle = molecularAssembly.getAngleList();
        nAngles = 0;
        for (ROLS r : angle) {
            if (keep((Angle) r)) {
                nAngles++;
            }
        }
        if (nAngles > angles.length) {
            angles = new Angle[nAngles];
        }
        Arrays.fill(angles, null);
        nAngles = 0;
        for (ROLS r : angle) {
            if (keep((Angle) r)) {
                angles[nAngles++] = (Angle) r;
            }
        }
        Arrays.sort(angles, 0, nAngles);
        if (nAngles > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Angles:                            %10d", nAngles));
        }
    } else {
        nAngles = 0;
        angles = null;
    }
    // Collect, count, pack and sort stretch-bends.
    if (stretchBendTerm) {
        ArrayList<ROLS> stretchBend = molecularAssembly.getStretchBendList();
        nStretchBends = 0;
        for (ROLS r : stretchBend) {
            if (keep((StretchBend) r)) {
                nStretchBends++;
            }
        }
        if (nStretchBends > stretchBends.length) {
            stretchBends = new StretchBend[nStretchBends];
        }
        Arrays.fill(stretchBends, null);
        nStretchBends = 0;
        for (ROLS r : stretchBend) {
            if (keep((StretchBend) r)) {
                stretchBends[nStretchBends++] = (StretchBend) r;
            }
        }
        Arrays.sort(stretchBends, 0, nStretchBends);
        if (nStretchBends > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Stretch-Bends:                     %10d", nStretchBends));
        }
    } else {
        nStretchBends = 0;
        stretchBends = null;
    }
    // Collect, count, pack and sort Urey-Bradleys.
    if (ureyBradleyTerm) {
        ArrayList<ROLS> ureyBradley = molecularAssembly.getUreyBradleyList();
        nUreyBradleys = 0;
        for (ROLS r : ureyBradley) {
            if (keep((UreyBradley) r)) {
                nUreyBradleys++;
            }
        }
        if (nUreyBradleys > ureyBradleys.length) {
            ureyBradleys = new UreyBradley[nUreyBradleys];
        }
        fill(ureyBradleys, null);
        nUreyBradleys = 0;
        for (ROLS r : ureyBradley) {
            if (keep((UreyBradley) r)) {
                ureyBradleys[nUreyBradleys++] = (UreyBradley) r;
            }
        }
        Arrays.sort(ureyBradleys, 0, nUreyBradleys);
        if (nUreyBradleys > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Urey-Bradleys:                     %10d", nUreyBradleys));
        }
    } else {
        nUreyBradleys = 0;
        ureyBradleys = null;
    }
    /**
     * Set a multiplier on the force constants of bonded terms containing
     * hydrogens.
     */
    if (rigidHydrogens) {
        if (bonds != null) {
            for (Bond bond : bonds) {
                if (bond.containsHydrogen()) {
                    bond.setRigidScale(rigidScale);
                }
            }
        }
        if (angles != null) {
            for (Angle angle : angles) {
                if (angle.containsHydrogen()) {
                    angle.setRigidScale(rigidScale);
                }
            }
        }
        if (stretchBends != null) {
            for (StretchBend stretchBend : stretchBends) {
                if (stretchBend.containsHydrogen()) {
                    stretchBend.setRigidScale(rigidScale);
                }
            }
        }
        if (ureyBradleys != null) {
            for (UreyBradley ureyBradley : ureyBradleys) {
                if (ureyBradley.containsHydrogen()) {
                    ureyBradley.setRigidScale(rigidScale);
                }
            }
        }
    }
    // Collect, count, pack and sort out-of-plane bends.
    if (outOfPlaneBendTerm) {
        ArrayList<ROLS> outOfPlaneBend = molecularAssembly.getOutOfPlaneBendList();
        nOutOfPlaneBends = 0;
        for (ROLS r : outOfPlaneBend) {
            if (keep((OutOfPlaneBend) r)) {
                nOutOfPlaneBends++;
            }
        }
        if (nOutOfPlaneBends > outOfPlaneBends.length) {
            outOfPlaneBends = new OutOfPlaneBend[nOutOfPlaneBends];
        }
        fill(outOfPlaneBends, null);
        nOutOfPlaneBends = 0;
        for (ROLS r : outOfPlaneBend) {
            if (keep((OutOfPlaneBend) r)) {
                outOfPlaneBends[nOutOfPlaneBends++] = (OutOfPlaneBend) r;
            }
        }
        Arrays.sort(outOfPlaneBends, 0, nOutOfPlaneBends);
        if (nOutOfPlaneBends > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Out-of-Plane Bends:                %10d", nOutOfPlaneBends));
        }
    } else {
        nOutOfPlaneBends = 0;
        outOfPlaneBends = null;
    }
    // Collect, count, pack and sort torsions.
    if (torsionTerm) {
        ArrayList<ROLS> torsion = molecularAssembly.getTorsionList();
        nTorsions = 0;
        for (ROLS r : torsion) {
            if (keep((Torsion) r)) {
                nTorsions++;
            }
        }
        if (nTorsions >= torsions.length) {
            torsions = new Torsion[nTorsions];
        }
        fill(torsions, null);
        nTorsions = 0;
        for (ROLS r : torsion) {
            if (keep((Torsion) r)) {
                torsions[nTorsions++] = (Torsion) r;
            }
        }
        // Arrays.sort(torsions);
        if (nTorsions > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Torsions:                          %10d", nTorsions));
        }
    } else {
        nTorsions = 0;
        torsions = null;
    }
    // Collect, count, pack and sort pi-orbital torsions.
    if (piOrbitalTorsionTerm) {
        ArrayList<ROLS> piOrbitalTorsion = molecularAssembly.getPiOrbitalTorsionList();
        nPiOrbitalTorsions = 0;
        for (ROLS r : piOrbitalTorsion) {
            if (keep((PiOrbitalTorsion) r)) {
                nPiOrbitalTorsions++;
            }
        }
        if (nPiOrbitalTorsions >= piOrbitalTorsions.length) {
            piOrbitalTorsions = new PiOrbitalTorsion[nPiOrbitalTorsions];
        }
        fill(piOrbitalTorsions, null);
        nPiOrbitalTorsions = 0;
        for (ROLS r : piOrbitalTorsion) {
            if (keep((PiOrbitalTorsion) r)) {
                piOrbitalTorsions[nPiOrbitalTorsions++] = (PiOrbitalTorsion) r;
            }
        }
        if (nPiOrbitalTorsions > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Pi-Orbital Torsions:               %10d", nPiOrbitalTorsions));
        }
    } else {
        nPiOrbitalTorsions = 0;
        piOrbitalTorsions = null;
    }
    // Collect, count, pack and sort torsion-torsions.
    if (torsionTorsionTerm) {
        ArrayList<ROLS> torsionTorsion = molecularAssembly.getTorsionTorsionList();
        nTorsionTorsions = 0;
        for (ROLS r : torsionTorsion) {
            if (keep((TorsionTorsion) r)) {
                nTorsionTorsions++;
            }
        }
        if (nTorsionTorsions >= torsionTorsions.length) {
            torsionTorsions = new TorsionTorsion[nTorsionTorsions];
        }
        fill(torsionTorsions, null);
        nTorsionTorsions = 0;
        for (ROLS r : torsionTorsion) {
            if (keep((TorsionTorsion) r)) {
                torsionTorsions[nTorsionTorsions++] = (TorsionTorsion) r;
            }
        }
        if (nTorsionTorsions > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Torsion-Torsions:                  %10d", nTorsionTorsions));
        }
    } else {
        nTorsionTorsions = 0;
        torsionTorsions = null;
    }
    // Collect, count, pack and sort improper torsions.
    if (improperTorsionTerm) {
        ArrayList<ROLS> improperTorsion = molecularAssembly.getImproperTorsionList();
        nImproperTorsions = 0;
        for (ROLS r : improperTorsion) {
            if (keep((ImproperTorsion) r)) {
                nImproperTorsions++;
            }
        }
        if (nImproperTorsions >= improperTorsions.length) {
            improperTorsions = new ImproperTorsion[nImproperTorsions];
        }
        fill(improperTorsions, null);
        nImproperTorsions = 0;
        for (ROLS r : improperTorsion) {
            if (keep((ImproperTorsion) r)) {
                improperTorsions[nImproperTorsions++] = (ImproperTorsion) r;
            }
        }
        if (nImproperTorsions > 0 && logger.isLoggable(Level.FINEST)) {
            logger.finest(format("  Improper Torsions:                 %10d", nImproperTorsions));
        }
    } else {
        nImproperTorsions = 0;
        improperTorsions = null;
    }
    if (vanderWaalsTerm) {
        if (esvTerm) {
            vanderWaals.setAtoms(esvSystem.getExtendedAtoms(), esvSystem.getExtendedMolecule());
        } else {
            vanderWaals.setAtoms(atoms, molecule);
        }
    }
    if (multipoleTerm) {
        if (esvTerm) {
            particleMeshEwald.setAtoms(esvSystem.getExtendedAtoms(), esvSystem.getExtendedMolecule());
        } else {
            particleMeshEwald.setAtoms(atoms, molecule);
        }
    }
    if (ncsTerm) {
        logger.severe(" NCS energy term cannot be used with variable systems sizes.");
    }
    if (restrainTerm) {
        logger.severe(" Restrain energy term cannot be used with variable systems sizes.");
    }
    if (comTerm) {
        logger.severe(" COM restrain energy term cannot be used with variable systems sizes.");
    }
    bondedRegion = new BondedRegion();
}

Example 17 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addAmoebaMultipoleForce.

private void addAmoebaMultipoleForce() {
    ParticleMeshEwald pme = super.getPmeNode();
    if (pme == null) {
        return;
    }
    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);
    amoebaMultipoleForce = OpenMM_AmoebaMultipoleForce_create();
    double polarScale = 1.0;
    if (pme.getPolarizationType() != Polarization.MUTUAL) {
        OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_Direct);
        if (pme.getPolarizationType() == Polarization.NONE) {
            polarScale = 0.0;
        }
    } else {
        ForceField forceField = molecularAssembly.getForceField();
        String algorithm = forceField.getString(ForceField.ForceFieldString.SCF_ALGORITHM, "CG");
        ParticleMeshEwald.SCFAlgorithm scfAlgorithm;
        try {
            algorithm = algorithm.replaceAll("-", "_").toUpperCase();
            scfAlgorithm = ParticleMeshEwald.SCFAlgorithm.valueOf(algorithm);
        } catch (Exception e) {
            scfAlgorithm = ParticleMeshEwald.SCFAlgorithm.CG;
        }
        switch(scfAlgorithm) {
            case EPT:
                logger.info(" Using extrapolated perturbation theory approximation instead of full SCF calculations. Not supported in FFX reference implementation.");
                OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_Extrapolated);
                PointerByReference exptCoefficients = OpenMM_DoubleArray_create(4);
                OpenMM_DoubleArray_set(exptCoefficients, 0, -0.154);
                OpenMM_DoubleArray_set(exptCoefficients, 1, 0.017);
                OpenMM_DoubleArray_set(exptCoefficients, 2, 0.657);
                OpenMM_DoubleArray_set(exptCoefficients, 3, 0.475);
                OpenMM_AmoebaMultipoleForce_setExtrapolationCoefficients(amoebaMultipoleForce, exptCoefficients);
                OpenMM_DoubleArray_destroy(exptCoefficients);
                break;
            case CG:
            case SOR:
            default:
                OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_Mutual);
                break;
        }
    }
    PointerByReference dipoles = OpenMM_DoubleArray_create(3);
    PointerByReference quadrupoles = OpenMM_DoubleArray_create(9);
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        MultipoleType multipoleType = atom.getMultipoleType();
        PolarizeType polarType = atom.getPolarizeType();
        /**
         * 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;
        }
        double useFactor = 1.0;
        if (!atoms[i].getUse() || !atoms[i].getElectrostatics()) {
            // if (!atoms[i].getUse()) {
            useFactor = 0.0;
        }
        // Should be 1.0 at this point.
        double lambdaScale = lambda;
        if (!atom.applyLambda()) {
            lambdaScale = 1.0;
        }
        useFactor *= lambdaScale;
        /**
         * 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 * useFactor / 3.0);
            }
        }
        // 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;
            logger.info(String.format(" Atom type %s", atom.getAtomType().toString()));
        }
        double charge = multipoleType.charge * useFactor;
        /**
         * Add the multipole.
         */
        OpenMM_AmoebaMultipoleForce_addMultipole(amoebaMultipoleForce, charge, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis, polarType.thole, polarType.pdamp * dampingFactorConversion, polarType.polarizability * polarityConversion * polarScale);
    }
    OpenMM_DoubleArray_destroy(dipoles);
    OpenMM_DoubleArray_destroy(quadrupoles);
    Crystal crystal = super.getCrystal();
    if (!crystal.aperiodic()) {
        OpenMM_AmoebaMultipoleForce_setNonbondedMethod(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_PME);
        OpenMM_AmoebaMultipoleForce_setCutoffDistance(amoebaMultipoleForce, pme.getEwaldCutoff() * OpenMM_NmPerAngstrom);
        OpenMM_AmoebaMultipoleForce_setAEwald(amoebaMultipoleForce, pme.getEwaldCoefficient() / OpenMM_NmPerAngstrom);
        double ewaldTolerance = 1.0e-04;
        OpenMM_AmoebaMultipoleForce_setEwaldErrorTolerance(amoebaMultipoleForce, ewaldTolerance);
        PointerByReference gridDimensions = OpenMM_IntArray_create(3);
        ReciprocalSpace recip = pme.getReciprocalSpace();
        OpenMM_IntArray_set(gridDimensions, 0, recip.getXDim());
        OpenMM_IntArray_set(gridDimensions, 1, recip.getYDim());
        OpenMM_IntArray_set(gridDimensions, 2, recip.getZDim());
        OpenMM_AmoebaMultipoleForce_setPmeGridDimensions(amoebaMultipoleForce, gridDimensions);
        OpenMM_IntArray_destroy(gridDimensions);
    } else {
        OpenMM_AmoebaMultipoleForce_setNonbondedMethod(amoebaMultipoleForce, OpenMM_AmoebaMultipoleForce_NoCutoff);
    }
    OpenMM_AmoebaMultipoleForce_setMutualInducedMaxIterations(amoebaMultipoleForce, 500);
    OpenMM_AmoebaMultipoleForce_setMutualInducedTargetEpsilon(amoebaMultipoleForce, pme.getPolarEps());
    int[][] ip11 = pme.getPolarization11();
    int[][] ip12 = pme.getPolarization12();
    int[][] ip13 = pme.getPolarization13();
    ArrayList<Integer> list12 = new ArrayList<>();
    ArrayList<Integer> list13 = new ArrayList<>();
    ArrayList<Integer> list14 = new ArrayList<>();
    PointerByReference covalentMap = OpenMM_IntArray_create(0);
    for (int i = 0; i < nAtoms; i++) {
        Atom ai = atoms[i];
        list12.clear();
        list13.clear();
        list14.clear();
        for (Bond bond : ai.getBonds()) {
            int index = bond.get1_2(ai).getIndex() - 1;
            OpenMM_IntArray_append(covalentMap, index);
            list12.add(index);
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent12, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (Angle angle : ai.getAngles()) {
            Atom ak = angle.get1_3(ai);
            if (ak != null) {
                int index = ak.getIndex() - 1;
                if (!list12.contains(index)) {
                    list13.add(index);
                    OpenMM_IntArray_append(covalentMap, index);
                }
            }
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent13, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (Torsion torsion : ai.getTorsions()) {
            Atom ak = torsion.get1_4(ai);
            if (ak != null) {
                int index = ak.getIndex() - 1;
                if (!list12.contains(index) && !list13.contains(index)) {
                    list14.add(index);
                    OpenMM_IntArray_append(covalentMap, index);
                }
            }
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent14, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (Atom ak : ai.get1_5s()) {
            int index = ak.getIndex() - 1;
            if (!list12.contains(index) && !list13.contains(index) && !list14.contains(index)) {
                OpenMM_IntArray_append(covalentMap, index);
            }
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_Covalent15, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
        for (int j = 0; j < ip11[i].length; j++) {
            OpenMM_IntArray_append(covalentMap, ip11[i][j]);
        }
        OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i, OpenMM_AmoebaMultipoleForce_PolarizationCovalent11, covalentMap);
        OpenMM_IntArray_resize(covalentMap, 0);
    // for (int j = 0; j < ip12[i].length; j++) {
    // OpenMM_IntArray_append(covalentMap, ip12[i][j]);
    // }
    // OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
    // OpenMM_AmoebaMultipoleForce_PolarizationCovalent12, covalentMap);
    // OpenMM_IntArray_resize(covalentMap, 0);
    // 
    // for (int j = 0; j < ip13[i].length; j++) {
    // OpenMM_IntArray_append(covalentMap, ip13[i][j]);
    // }
    // OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
    // OpenMM_AmoebaMultipoleForce_PolarizationCovalent13, covalentMap);
    // OpenMM_IntArray_resize(covalentMap, 0);
    // 
    // OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
    // OpenMM_AmoebaMultipoleForce_PolarizationCovalent14, covalentMap);
    }
    OpenMM_IntArray_destroy(covalentMap);
    OpenMM_System_addForce(system, amoebaMultipoleForce);
    OpenMM_Force_setForceGroup(amoebaMultipoleForce, 1);
    logger.log(Level.INFO, " Added polarizable multipole force.");
    GeneralizedKirkwood gk = super.getGK();
    if (gk != null) {
        addGKForce();
    }
}

Example 18 with Bond

use of ffx.potential.bonded.Bond 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 19 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class MolecularAssembly method renderWire.

private Shape3D renderWire() {
    ArrayList<ROLS> bonds = getBondList();
    int numbonds = bonds.size();
    if (numbonds < 1) {
        return null;
    }
    Vector3d bondmidpoint = new Vector3d();
    double[] mid = { 0, 0, 0 };
    Vector3d v1 = new Vector3d();
    Vector3d v2 = new Vector3d();
    float[] a1 = { 0, 0, 0 };
    float[] a2 = { 0, 0, 0 };
    float[] col = new float[4];
    Bond bond;
    Atom atom1, atom2;
    LineArray la = new LineArray(4 * numbonds, GeometryArray.COORDINATES | GeometryArray.COLOR_4 | GeometryArray.NORMALS);
    la.setCapability(LineArray.ALLOW_COORDINATE_WRITE);
    la.setCapability(LineArray.ALLOW_COORDINATE_READ);
    la.setCapability(LineArray.ALLOW_COLOR_WRITE);
    la.setCapability(LineArray.ALLOW_COUNT_READ);
    la.setCapability(LineArray.ALLOW_INTERSECT);
    la.setCapability(LineArray.ALLOW_FORMAT_READ);
    atomLookUp = new Atom[4 * numbonds];
    int i = 0;
    col[3] = 0.9f;
    for (ListIterator<ROLS> li = bonds.listIterator(); li.hasNext(); ) {
        bond = (Bond) li.next();
        bond.setWire(la, i);
        atom1 = bond.getAtom(0);
        atom2 = bond.getAtom(1);
        atom1.getV3D(v1);
        atom2.getV3D(v2);
        a1[0] = (float) v1.x;
        a1[1] = (float) v1.y;
        a1[2] = (float) v1.z;
        a2[0] = (float) v2.x;
        a2[1] = (float) v2.y;
        a2[2] = (float) v2.z;
        // Find the bond center
        bondmidpoint.add(v1, v2);
        bondmidpoint.scale(0.5d);
        bondmidpoint.get(mid);
        // Atom #1
        Atom.AtomColor.get(atom1.getAtomicNumber()).get(col);
        atomLookUp[i] = atom1;
        la.setCoordinate(i, a1);
        la.setColor(i, col);
        la.setNormal(i, a2);
        i++;
        atomLookUp[i] = atom1;
        la.setCoordinate(i, mid);
        la.setColor(i, col);
        la.setNormal(i, a2);
        i++;
        // Atom #2
        Atom.AtomColor.get(atom2.getAtomicNumber()).get(col);
        atomLookUp[i] = atom2;
        la.setCoordinate(i, a2);
        la.setColor(i, col);
        la.setNormal(i, a1);
        i++;
        atomLookUp[i] = atom2;
        la.setCoordinate(i, mid);
        la.setColor(i, col);
        la.setNormal(i, a1);
        i++;
    }
    ColoringAttributes cola = new ColoringAttributes(new Color3f(), ColoringAttributes.SHADE_GOURAUD);
    Appearance app = new Appearance();
    lineAttributes = new LineAttributes();
    lineAttributes.setLineWidth(RendererCache.bondwidth);
    lineAttributes.setCapability(LineAttributes.ALLOW_WIDTH_WRITE);
    lineAttributes.setLineAntialiasingEnable(true);
    app.setLineAttributes(lineAttributes);
    app.setCapability(Appearance.ALLOW_LINE_ATTRIBUTES_READ);
    app.setCapability(Appearance.ALLOW_LINE_ATTRIBUTES_WRITE);
    RenderingAttributes ra = new RenderingAttributes();
    ra.setAlphaTestValue(0.1f);
    ra.setAlphaTestFunction(RenderingAttributes.GREATER);
    ra.setDepthBufferEnable(true);
    ra.setDepthBufferWriteEnable(true);
    app.setRenderingAttributes(ra);
    app.setColoringAttributes(cola);
    Shape3D wireframe = new Shape3D(la, app);
    // PickTool.setCapabilities(wire, PickTool.INTERSECT_COORD);
    wireframe.setUserData(this);
    wireframe.setBounds(new BoundingSphere(new Point3d(0, 0, 0), 1000.0));
    try {
        wireframe.setBoundsAutoCompute(false);
    } catch (Exception e) {
        e.printStackTrace();
    }
    wireframe.setCapability(Shape3D.ALLOW_GEOMETRY_READ);
    wireframe.setCapability(Shape3D.ALLOW_APPEARANCE_READ);
    wireframe.setCapability(Shape3D.ALLOW_LOCAL_TO_VWORLD_READ);
    return wireframe;
}

Example 20 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class Utilities method findPolymer.

/**
 * <p>
 * findPolymer</p>
 *
 * @param atoms List
 * @param currentAtom Atom
 * @param path List
 * @return List
 */
public static List<Atom> findPolymer(List<Atom> atoms, Atom currentAtom, List<Atom> path) {
    // Atom has no bonds to follow
    if (currentAtom.getBonds() == null) {
        path = getAtomListFromPool();
        path.add(currentAtom);
        return path;
    }
    // End of Recursion conditions
    if (currentAtom.getParent() != null) {
        return null;
    }
    int anum = currentAtom.getAtomicNumber();
    // Only C,N,O,P in a DNA/RNA/protein backbone
    if (anum != 6 && anum != 7 && anum != 8 && anum != 15) {
        return null;
    }
    // Allow the search to make it out of side chains, but not enter them...
    if (path != null && path.size() > 7) {
        // group
        if (anum == 8) {
            if (!formsBondsWith(currentAtom, 15)) {
                return null;
            }
        // Nitrogen is only in the backbone in peptide bonds
        } else if (anum == 7) {
            Atom carbonyl = findCarbonyl(currentAtom);
            if (carbonyl == null) {
                return null;
            }
            Atom alphaCarbon = findAlphaCarbon(currentAtom);
            if (alphaCarbon == null) {
                return null;
            }
        // Avoid more than 3 carbons in a row (phenyl groups, etc.)
        } else if (anum == 6) {
            Atom a;
            int anum2, anum3, anum4;
            int size = path.size();
            a = path.get(size - 1);
            anum2 = a.getAtomicNumber();
            if (anum2 == 6) {
                a = path.get(size - 2);
                anum3 = a.getAtomicNumber();
                if (anum3 == 6) {
                    a = path.get(size - 3);
                    anum4 = a.getAtomicNumber();
                    if (anum4 == 6) {
                        return null;
                    }
                }
            }
        }
    }
    // Atoms with only one bond are at the end of a Polymer
    Atom previousAtom = null;
    if (path != null) {
        previousAtom = path.get(path.size() - 1);
    }
    List<Bond> bonds = currentAtom.getBonds();
    if (bonds.size() == 1 && previousAtom != null) {
        return null;
    }
    // Initialization
    if (path == null) {
        path = getAtomListFromPool();
        previousAtom = null;
    // Or Continuation
    } else {
        List<Atom> pathclone = getAtomListFromPool();
        pathclone.addAll(path);
        path = pathclone;
    }
    // Add the currentAtom to the growing path
    path.add(currentAtom);
    // Continue search in each bond direction, but no backtracking over
    // previousAtom
    Atom nextAtom;
    List<Atom> newPolymer, maxPolymer = getAtomListFromPool();
    for (Bond b : bonds) {
        nextAtom = b.get1_2(currentAtom);
        // Check to avoid returning in the same direction and loops
        if (nextAtom != previousAtom && !path.contains(nextAtom)) {
            newPolymer = findPolymer(atoms, nextAtom, path);
            if (newPolymer != null) {
                // and if so, use the shorter Polymer (avoids loops)
                if (haveCommonAtom(newPolymer, maxPolymer)) {
                    if (newPolymer.size() < maxPolymer.size()) {
                        addAtomListToPool(maxPolymer);
                        maxPolymer = newPolymer;
                    }
                } else if (newPolymer.size() > maxPolymer.size()) {
                    addAtomListToPool(maxPolymer);
                    maxPolymer = newPolymer;
                }
            }
        }
    }
    // Add the currentAtom to the longest discovered chain and return
    maxPolymer.add(0, currentAtom);
    return maxPolymer;
}