Search in sources :

Example 36 with Atom

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

the class ForceFieldEnergyOpenMM method addTorsionTorsionForce.

private void addTorsionTorsionForce() {
    TorsionTorsion[] torsionTorsions = super.getTorsionTorsions();
    if (torsionTorsions == null || torsionTorsions.length < 1) {
        return;
    }
    /**
     * Load the torsion-torsions.
     */
    int nTypes = 0;
    LinkedHashMap<String, TorsionTorsionType> torTorTypes = new LinkedHashMap<>();
    int nTorsionTorsions = torsionTorsions.length;
    amoebaTorsionTorsionForce = OpenMM_AmoebaTorsionTorsionForce_create();
    for (int i = 0; i < nTorsionTorsions; i++) {
        TorsionTorsion torsionTorsion = torsionTorsions[i];
        int ia = torsionTorsion.getAtom(0).getXyzIndex() - 1;
        int ib = torsionTorsion.getAtom(1).getXyzIndex() - 1;
        int ic = torsionTorsion.getAtom(2).getXyzIndex() - 1;
        int id = torsionTorsion.getAtom(3).getXyzIndex() - 1;
        int ie = torsionTorsion.getAtom(4).getXyzIndex() - 1;
        TorsionTorsionType torsionTorsionType = torsionTorsion.torsionTorsionType;
        String key = torsionTorsionType.getKey();
        /**
         * Check if the TorTor parameters have already been added to the
         * Hash.
         */
        int gridIndex = 0;
        if (torTorTypes.containsKey(key)) {
            /**
             * If the TorTor has been added, get its (ordered) index in the
             * Hash.
             */
            int index = 0;
            for (String entry : torTorTypes.keySet()) {
                if (entry.equalsIgnoreCase(key)) {
                    gridIndex = index;
                    break;
                } else {
                    index++;
                }
            }
        } else {
            /**
             * Add the new TorTor.
             */
            torTorTypes.put(key, torsionTorsionType);
            gridIndex = nTypes;
            nTypes++;
        }
        Atom atom = torsionTorsion.getChiralAtom();
        int iChiral = -1;
        if (atom != null) {
            iChiral = atom.getXyzIndex() - 1;
        }
        OpenMM_AmoebaTorsionTorsionForce_addTorsionTorsion(amoebaTorsionTorsionForce, ia, ib, ic, id, ie, iChiral, gridIndex);
    }
    /**
     * Load the Torsion-Torsion parameters.
     */
    PointerByReference values = OpenMM_DoubleArray_create(6);
    int gridIndex = 0;
    for (String key : torTorTypes.keySet()) {
        TorsionTorsionType torTorType = torTorTypes.get(key);
        int nx = torTorType.nx;
        int ny = torTorType.ny;
        double[] tx = torTorType.tx;
        double[] ty = torTorType.ty;
        double[] f = torTorType.energy;
        double[] dx = torTorType.dx;
        double[] dy = torTorType.dy;
        double[] dxy = torTorType.dxy;
        /**
         * Create the 3D grid.
         */
        PointerByReference grid3D = OpenMM_3D_DoubleArray_create(nx, ny, 6);
        int xIndex = 0;
        int yIndex = 0;
        for (int j = 0; j < nx * ny; j++) {
            int addIndex = 0;
            OpenMM_DoubleArray_set(values, addIndex++, tx[xIndex]);
            OpenMM_DoubleArray_set(values, addIndex++, ty[yIndex]);
            OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * f[j]);
            OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dx[j]);
            OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dy[j]);
            OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dxy[j]);
            OpenMM_3D_DoubleArray_set(grid3D, yIndex, xIndex, values);
            xIndex++;
            if (xIndex == nx) {
                xIndex = 0;
                yIndex++;
            }
        }
        OpenMM_AmoebaTorsionTorsionForce_setTorsionTorsionGrid(amoebaTorsionTorsionForce, gridIndex++, grid3D);
        OpenMM_3D_DoubleArray_destroy(grid3D);
    }
    OpenMM_DoubleArray_destroy(values);
    OpenMM_System_addForce(system, amoebaTorsionTorsionForce);
    logger.log(Level.INFO, " Added Torsion-Torsions ({0})", nTorsionTorsions);
}
Also used : TorsionTorsionType(ffx.potential.parameters.TorsionTorsionType) PointerByReference(com.sun.jna.ptr.PointerByReference) OpenMM_AmoebaTorsionTorsionForce_addTorsionTorsion(simtk.openmm.AmoebaOpenMMLibrary.OpenMM_AmoebaTorsionTorsionForce_addTorsionTorsion) TorsionTorsion(ffx.potential.bonded.TorsionTorsion) OpenMM_System_addConstraint(simtk.openmm.OpenMMLibrary.OpenMM_System_addConstraint) CoordRestraint(ffx.potential.nonbonded.CoordRestraint) Atom(ffx.potential.bonded.Atom) LinkedHashMap(java.util.LinkedHashMap)

Example 37 with Atom

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

the class ForceFieldEnergyOpenMM method createVirtualHydrogenSites.

/**
 * Experimental. Virtual hydrogen sites require creation of new particles,
 * which then need to be handled (ignored?) for the multiple force.
 */
private void createVirtualHydrogenSites() {
    VanDerWaals vdW = super.getVdwNode();
    if (vdW == null) {
        return;
    }
    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();
        if (vdwType.reductionFactor < 1.0) {
            double factor = vdwType.reductionFactor;
            // Create the virtual site.
            PointerByReference virtualSite = OpenMM_TwoParticleAverageSite_create(i, ired[i], factor, 1.0 - factor);
            // Create a massless particle for the hydrogen vdW site.
            int id = OpenMM_System_addParticle(system, 0.0);
            // Denote the massless particle is a virtual site
            OpenMM_System_setVirtualSite(system, id, virtualSite);
        }
    }
}
Also used : VDWType(ffx.potential.parameters.VDWType) VanDerWaals(ffx.potential.nonbonded.VanDerWaals) PointerByReference(com.sun.jna.ptr.PointerByReference) Atom(ffx.potential.bonded.Atom) OpenMM_System_addConstraint(simtk.openmm.OpenMMLibrary.OpenMM_System_addConstraint) CoordRestraint(ffx.potential.nonbonded.CoordRestraint)

Example 38 with Atom

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

the class ForceFieldEnergyOpenMM method addCustomGBForce.

private void addCustomGBForce() {
    GeneralizedKirkwood gk = super.getGK();
    if (gk == null) {
        return;
    }
    customGBForce = OpenMM_CustomGBForce_create();
    OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "q");
    OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "radius");
    OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "scale");
    OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "solventDielectric", 78.3);
    OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "soluteDielectric", 1.0);
    // Factor of 0.1 for Ang to nm.
    OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "dOffset", gk.getDielecOffset() * OpenMM_NmPerAngstrom);
    OpenMM_CustomGBForce_addComputedValue(customGBForce, "I", // "step(r+sr2-or1)*0.5*((1/L^3-1/U^3)/3+(1/U^4-1/L^4)/8*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);"
    "0.5*((1/L^3-1/U^3)/3.0+(1/U^4-1/L^4)/8.0*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);" + "U=r+sr2;" + // + "C=2*(1/or1-1/L)*step(sr2-r-or1);"
    "C=2/3*(1/or1^3-1/L^3)*step(sr2-r-or1);" + // + "D=step(r-sr2)*(r-sr2) + (1-step(r-sr2))*(sr2-r);"
    "L = step(sr2 - r1r)*sr2mr + (1 - step(sr2 - r1r))*L;" + "sr2mr = sr2 - r;" + "r1r = radius1 + r;" + "L = step(r1sr2 - r)*radius1 + (1 - step(r1sr2 - r))*L;" + "r1sr2 = radius1 + sr2;" + "L = r - sr2;" + "sr2 = scale2 * radius2;" + "or1 = radius1; or2 = radius2", OpenMM_CustomGBForce_ParticlePairNoExclusions);
    OpenMM_CustomGBForce_addComputedValue(customGBForce, "B", // "psi=I*or; or=radius-0.009"
    "step(BB-radius)*BB + (1 - step(BB-radius))*radius;" + "BB = 1 / ( (3.0*III)^(1.0/3.0) );" + "III = step(II)*II + (1 - step(II))*1.0e-9/3.0;" + "II = maxI - I;" + "maxI = 1/(3.0*radius^3)", OpenMM_CustomGBForce_SingleParticle);
    double sTens = gk.getSurfaceTension();
    logger.info(String.format(" FFX surface tension: %9.5g kcal/mol/Ang^2", sTens));
    sTens *= OpenMM_KJPerKcal;
    // 100 square Angstroms per square nanometer.
    sTens *= 100.0;
    logger.info(String.format(" OpenMM surface tension: %9.5g kJ/mol/nm^2", sTens));
    String surfaceTension = Double.toString(sTens);
    OpenMM_CustomGBForce_addEnergyTerm(customGBForce, surfaceTension + "*(radius+0.14+dOffset)^2*((radius+dOffset)/B)^6/6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B", OpenMM_CustomGBForce_SingleParticle);
    /**
     * Particle pair term is the generalized Born cross term.
     */
    OpenMM_CustomGBForce_addEnergyTerm(customGBForce, "-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;" + "f=sqrt(r^2+B1*B2*exp(-r^2/(2.455*B1*B2)))", OpenMM_CustomGBForce_ParticlePair);
    double[] baseRadii = gk.getBaseRadii();
    double[] overlapScale = gk.getOverlapScale();
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    PointerByReference doubleArray = OpenMM_DoubleArray_create(0);
    for (int i = 0; i < nAtoms; i++) {
        MultipoleType multipoleType = atoms[i].getMultipoleType();
        OpenMM_DoubleArray_append(doubleArray, multipoleType.charge);
        OpenMM_DoubleArray_append(doubleArray, OpenMM_NmPerAngstrom * baseRadii[i]);
        OpenMM_DoubleArray_append(doubleArray, overlapScale[i]);
        OpenMM_CustomGBForce_addParticle(customGBForce, doubleArray);
        OpenMM_DoubleArray_resize(doubleArray, 0);
    }
    OpenMM_DoubleArray_destroy(doubleArray);
    double cut = gk.getCutoff();
    OpenMM_CustomGBForce_setCutoffDistance(customGBForce, cut);
    OpenMM_Force_setForceGroup(customGBForce, 1);
    OpenMM_System_addForce(system, customGBForce);
    logger.log(Level.INFO, " Added generalized Born force");
}
Also used : GeneralizedKirkwood(ffx.potential.nonbonded.GeneralizedKirkwood) PointerByReference(com.sun.jna.ptr.PointerByReference) MultipoleType(ffx.potential.parameters.MultipoleType) Atom(ffx.potential.bonded.Atom) OpenMM_System_addConstraint(simtk.openmm.OpenMMLibrary.OpenMM_System_addConstraint) CoordRestraint(ffx.potential.nonbonded.CoordRestraint)

Example 39 with Atom

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

the class MolecularAssembly method getMolecule.

private Atom getMolecule(Atom atom, boolean create) {
    String resName = atom.getResidueName();
    int resNum = atom.getResidueNumber();
    Character chainID = atom.getChainID();
    String segID = atom.getSegID();
    ArrayList<MSNode> list = ions.getChildList();
    for (MSNode node : list) {
        Molecule m = (Molecule) node;
        if (m.getSegID().equalsIgnoreCase(segID) && m.getResidueName().equalsIgnoreCase(resName) && m.getResidueNumber() == resNum) {
            return (Atom) m.addMSNode(atom);
        }
    }
    list = water.getChildList();
    for (MSNode node : list) {
        Molecule m = (Molecule) node;
        if (m.getSegID().equalsIgnoreCase(segID) && m.getResidueName().equalsIgnoreCase(resName) && m.getResidueNumber() == resNum) {
            return (Atom) m.addMSNode(atom);
        }
    }
    list = molecules.getChildList();
    for (MSNode node : list) {
        Molecule m = (Molecule) node;
        if (m.getSegID().equalsIgnoreCase(segID) && m.getResidueName().equalsIgnoreCase(resName) && m.getResidueNumber() == resNum) {
            return (Atom) m.addMSNode(atom);
        }
    }
    if (create) {
        Molecule m = new Molecule(resName, resNum, chainID, segID);
        m.addMSNode(atom);
        if (resName.equalsIgnoreCase("DOD") || resName.equalsIgnoreCase("HOH") || resName.equalsIgnoreCase("WAT")) {
            water.add(m);
        // NA, K, MG, MG2, CA, CA2, CL
        } else if (resName.equalsIgnoreCase("NA") || resName.equalsIgnoreCase("K") || resName.equalsIgnoreCase("MG") || resName.equalsIgnoreCase("MG2") || resName.equalsIgnoreCase("CA") || resName.equalsIgnoreCase("CA2") || resName.equalsIgnoreCase("CL") || resName.equalsIgnoreCase("BR") || resName.equalsIgnoreCase("ZN") || resName.equalsIgnoreCase("ZN2")) {
            ions.add(m);
        } else {
            molecules.add(m);
        }
        return atom;
    } else {
        return null;
    }
}
Also used : Molecule(ffx.potential.bonded.Molecule) MSNode(ffx.potential.bonded.MSNode) Atom(ffx.potential.bonded.Atom)

Example 40 with Atom

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

the class MolecularAssembly method findAtom.

/**
 * <p>
 * findAtom</p>
 *
 * @param atom a {@link ffx.potential.bonded.Atom} object.
 * @return a {@link ffx.potential.bonded.Atom} object.
 */
public Atom findAtom(Atom atom) {
    if (!atom.isHetero() || atom.isModRes()) {
        Polymer polymer = getPolymer(atom.getChainID(), atom.getSegID(), false);
        if (polymer != null) {
            Residue res = polymer.getResidue(atom.getResidueName(), atom.getResidueNumber(), false);
            if (res != null) {
                MSNode node = res.getAtomNode();
                Atom root = (Atom) node.contains(atom);
                return root;
            }
        }
        return null;
    } else {
        ArrayList<Molecule> list = getMolecules();
        for (MSNode node : list) {
            Molecule m = (Molecule) node;
            if (m.getSegID().equalsIgnoreCase(atom.getSegID()) && m.getResidueName().equalsIgnoreCase(atom.getResidueName()) && m.getResidueNumber() == atom.getResidueNumber()) {
                Atom root = (Atom) node.contains(atom);
                return root;
            }
        }
        ArrayList<MSNode> ionList = getIons();
        for (MSNode node : ionList) {
            Molecule m = (Molecule) node;
            if (m.getSegID().equalsIgnoreCase(atom.getSegID()) && m.getResidueName().equalsIgnoreCase(atom.getResidueName()) && m.getResidueNumber() == atom.getResidueNumber()) {
                Atom root = (Atom) node.contains(atom);
                return root;
            }
        }
        ArrayList<MSNode> waterList = getWaters();
        for (MSNode node : waterList) {
            Molecule m = (Molecule) node;
            if (m.getSegID().equalsIgnoreCase(atom.getSegID()) && m.getResidueName().equalsIgnoreCase(atom.getResidueName()) && m.getResidueNumber() == atom.getResidueNumber()) {
                Atom root = (Atom) node.contains(atom);
                return root;
            }
        }
        return null;
    }
}
Also used : Molecule(ffx.potential.bonded.Molecule) MSNode(ffx.potential.bonded.MSNode) Residue(ffx.potential.bonded.Residue) Polymer(ffx.potential.bonded.Polymer) Atom(ffx.potential.bonded.Atom)

Aggregations

Atom (ffx.potential.bonded.Atom)206 Residue (ffx.potential.bonded.Residue)42 Bond (ffx.potential.bonded.Bond)37 CoordRestraint (ffx.potential.nonbonded.CoordRestraint)34 ArrayList (java.util.ArrayList)33 OpenMM_System_addConstraint (simtk.openmm.OpenMMLibrary.OpenMM_System_addConstraint)23 Polymer (ffx.potential.bonded.Polymer)22 IOException (java.io.IOException)22 MSNode (ffx.potential.bonded.MSNode)21 Crystal (ffx.crystal.Crystal)20 Molecule (ffx.potential.bonded.Molecule)19 File (java.io.File)19 MultiResidue (ffx.potential.bonded.MultiResidue)17 MultipoleType (ffx.potential.parameters.MultipoleType)13 MissingHeavyAtomException (ffx.potential.bonded.BondedUtils.MissingHeavyAtomException)12 PointerByReference (com.sun.jna.ptr.PointerByReference)11 MissingAtomTypeException (ffx.potential.bonded.BondedUtils.MissingAtomTypeException)11 AtomType (ffx.potential.parameters.AtomType)11 List (java.util.List)11 MolecularAssembly (ffx.potential.MolecularAssembly)10