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

Example 91 with Atom

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

the class BiojavaFilter method assignAtomTypes.

/**
 * Assign force field atoms types to common chemistries using "biotype"
 * records.
 */
public void assignAtomTypes() {
    /**
     * Create a new List to store bonds determined based on PDB atom names.
     */
    bondList = new ArrayList<>();
    /**
     * To Do: Look for cyclic peptides and disulfides.
     */
    Polymer[] polymers = activeMolecularAssembly.getChains();
    /**
     * Loop over chains.
     */
    if (polymers != null) {
        logger.info(format("\n Assigning atom types for %d chains.", polymers.length));
        for (Polymer polymer : polymers) {
            List<Residue> residues = polymer.getResidues();
            int numberOfResidues = residues.size();
            /**
             * Check if all residues are known amino acids.
             */
            boolean isProtein = true;
            if (!residues.isEmpty()) {
            // renameNTerminusHydrogens(residues.get(0)); Not safe to use until it distinguishes between true N-termini and N-terminal residues in general.
            }
            for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
                Residue residue = residues.get(residueNumber);
                String name = residue.getName().toUpperCase();
                boolean aa = false;
                for (AminoAcid3 amino : aminoAcidList) {
                    if (amino.toString().equalsIgnoreCase(name)) {
                        aa = true;
                        renameNonstandardHydrogens(residue);
                        break;
                    }
                }
                // Check for a patch.
                if (!aa) {
                    HashMap<String, AtomType> types = forceField.getAtomTypes(name);
                    if (types.isEmpty()) {
                        isProtein = false;
                        break;
                    } else {
                        logger.info(" Patch found for non-standard amino acid " + name);
                    }
                }
            }
            /**
             * If all the residues in this chain have known amino acids
             * names, then attempt to assign atom types.
             */
            if (isProtein) {
                try {
                    logger.info(format(" Amino acid chain %s", polymer.getName()));
                    double dist = properties.getDouble("chainbreak", 3.0);
                    // Detect main chain breaks!
                    List<List<Residue>> subChains = findChainBreaks(residues, dist);
                    for (List<Residue> subChain : subChains) {
                        assignAminoAcidAtomTypes(subChain);
                    }
                } catch (MissingHeavyAtomException missingHeavyAtomException) {
                    logger.severe(missingHeavyAtomException.toString());
                } catch (MissingAtomTypeException missingAtomTypeException) {
                    logger.severe(missingAtomTypeException.toString());
                }
                continue;
            }
            /**
             * Check if all residues have known nucleic acids names.
             */
            boolean isNucleicAcid = true;
            for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
                Residue residue = residues.get(residueNumber);
                String name = residue.getName().toUpperCase();
                /**
                 * Convert 1 and 2-character nucleic acid names to
                 * 3-character names.
                 */
                if (name.length() == 1) {
                    if (name.equals("A")) {
                        name = NucleicAcid3.ADE.toString();
                    } else if (name.equals("C")) {
                        name = NucleicAcid3.CYT.toString();
                    } else if (name.equals("G")) {
                        name = NucleicAcid3.GUA.toString();
                    } else if (name.equals("T")) {
                        name = NucleicAcid3.THY.toString();
                    } else if (name.equals("U")) {
                        name = NucleicAcid3.URI.toString();
                    }
                } else if (name.length() == 2) {
                    if (name.equals("YG")) {
                        name = NucleicAcid3.YYG.toString();
                    }
                }
                residue.setName(name);
                NucleicAcid3 nucleicAcid = null;
                for (NucleicAcid3 nucleic : nucleicAcidList) {
                    String nuc3 = nucleic.toString();
                    nuc3 = nuc3.substring(nuc3.length() - 3);
                    if (nuc3.equalsIgnoreCase(name)) {
                        nucleicAcid = nucleic;
                        break;
                    }
                }
                if (nucleicAcid == null) {
                    logger.info(format("Nucleic acid was not recognized %s.", name));
                    isNucleicAcid = false;
                    break;
                }
            }
            /**
             * If all the residues in this chain have known nucleic acids
             * names, then attempt to assign atom types.
             */
            if (isNucleicAcid) {
                try {
                    logger.info(format(" Nucleic acid chain %s", polymer.getName()));
                    assignNucleicAcidAtomTypes(residues);
                } catch (MissingHeavyAtomException missingHeavyAtomException) {
                    logger.severe(missingHeavyAtomException.toString());
                } catch (MissingAtomTypeException missingAtomTypeException) {
                    logger.severe(missingAtomTypeException.toString());
                }
            }
        }
    }
    // Assign ion atom types.
    ArrayList<MSNode> ions = activeMolecularAssembly.getIons();
    if (ions != null && ions.size() > 0) {
        logger.info(format(" Assigning atom types for %d ions.", ions.size()));
        for (MSNode m : ions) {
            Molecule ion = (Molecule) m;
            String name = ion.getResidueName().toUpperCase();
            HetAtoms hetatm = HetAtoms.valueOf(name);
            Atom atom = ion.getAtomList().get(0);
            if (ion.getAtomList().size() != 1) {
                logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
            }
            try {
                switch(hetatm) {
                    case NA:
                        atom.setAtomType(findAtomType(2003));
                        break;
                    case K:
                        atom.setAtomType(findAtomType(2004));
                        break;
                    case MG:
                    case MG2:
                        atom.setAtomType(findAtomType(2005));
                        break;
                    case CA:
                    case CA2:
                        atom.setAtomType(findAtomType(2006));
                        break;
                    case CL:
                        atom.setAtomType(findAtomType(2007));
                        break;
                    case ZN:
                    case ZN2:
                        atom.setAtomType(findAtomType(2008));
                        break;
                    case BR:
                        atom.setAtomType(findAtomType(2009));
                        break;
                    default:
                        logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
                }
            } catch (Exception e) {
                String message = "Error assigning atom types.";
                logger.log(Level.SEVERE, message, e);
            }
        }
    }
    // Assign water atom types.
    ArrayList<MSNode> water = activeMolecularAssembly.getWaters();
    if (water != null && water.size() > 0) {
        logger.info(format(" Assigning atom types for %d waters.", water.size()));
        for (MSNode m : water) {
            Molecule wat = (Molecule) m;
            try {
                Atom O = buildHeavy(wat, "O", null, 2001);
                Atom H1 = buildHydrogen(wat, "H1", O, 0.96e0, null, 109.5e0, null, 120.0e0, 0, 2002);
                H1.setHetero(true);
                Atom H2 = buildHydrogen(wat, "H2", O, 0.96e0, H1, 109.5e0, null, 120.0e0, 0, 2002);
                H2.setHetero(true);
            } catch (Exception e) {
                String message = "Error assigning atom types to a water.";
                logger.log(Level.SEVERE, message, e);
            }
        }
    }
    // Assign small molecule atom types.
    ArrayList<Molecule> molecules = activeMolecularAssembly.getMolecules();
    for (MSNode m : molecules) {
        Molecule molecule = (Molecule) m;
        String moleculeName = molecule.getResidueName();
        logger.info(" Attempting to patch " + moleculeName);
        ArrayList<Atom> moleculeAtoms = molecule.getAtomList();
        boolean patched = true;
        HashMap<String, AtomType> types = forceField.getAtomTypes(moleculeName);
        /**
         * Assign atom types for all known atoms.
         */
        for (Atom atom : moleculeAtoms) {
            String atomName = atom.getName().toUpperCase();
            AtomType atomType = types.get(atomName);
            if (atomType == null) {
                logger.info(" No atom type was found for " + atomName + " of " + moleculeName + ".");
                patched = false;
                break;
            } else {
                atom.setAtomType(atomType);
                types.remove(atomName);
            }
        }
        /**
         * Create missing hydrogen atoms. Check for missing heavy atoms.
         */
        if (patched && !types.isEmpty()) {
            for (AtomType type : types.values()) {
                if (type.atomicNumber != 1) {
                    logger.info(" Missing heavy atom " + type.name);
                    patched = false;
                    break;
                }
            }
        }
        // Create bonds between known atoms.
        if (patched) {
            for (Atom atom : moleculeAtoms) {
                String atomName = atom.getName();
                String[] bonds = forceField.getBonds(moleculeName, atomName);
                if (bonds != null) {
                    for (String name : bonds) {
                        Atom atom2 = molecule.getAtom(name);
                        if (atom2 != null && !atom.isBonded(atom2)) {
                            buildBond(atom, atom2);
                        }
                    }
                }
            }
        }
        // Create missing hydrogen atoms.
        if (patched && !types.isEmpty()) {
            // Create a hashmap of the molecule's atoms
            HashMap<String, Atom> atomMap = new HashMap<String, Atom>();
            for (Atom atom : moleculeAtoms) {
                atomMap.put(atom.getName().toUpperCase(), atom);
            }
            for (String atomName : types.keySet()) {
                AtomType type = types.get(atomName);
                String[] bonds = forceField.getBonds(moleculeName, atomName.toUpperCase());
                if (bonds == null || bonds.length != 1) {
                    patched = false;
                    logger.info(" Check biotype for hydrogen " + type.name + ".");
                    break;
                }
                // Get the heavy atom the hydrogen is bonded to.
                Atom ia = atomMap.get(bonds[0].toUpperCase());
                Atom hydrogen = new Atom(0, atomName, ia.getAltLoc(), new double[3], ia.getResidueName(), ia.getResidueNumber(), ia.getChainID(), ia.getOccupancy(), ia.getTempFactor(), ia.getSegID());
                logger.fine(" Created hydrogen " + atomName + ".");
                hydrogen.setAtomType(type);
                hydrogen.setHetero(true);
                molecule.addMSNode(hydrogen);
                int valence = ia.getAtomType().valence;
                List<Bond> aBonds = ia.getBonds();
                int numBonds = aBonds.size();
                /**
                 * Try to find the following configuration: ib-ia-ic
                 */
                Atom ib = null;
                Atom ic = null;
                Atom id = null;
                if (numBonds > 0) {
                    Bond bond = aBonds.get(0);
                    ib = bond.get1_2(ia);
                }
                if (numBonds > 1) {
                    Bond bond = aBonds.get(1);
                    ic = bond.get1_2(ia);
                }
                if (numBonds > 2) {
                    Bond bond = aBonds.get(2);
                    id = bond.get1_2(ia);
                }
                /**
                 * Building the hydrogens depends on hybridization and the
                 * locations of other bonded atoms.
                 */
                logger.fine(" Bonding " + atomName + " to " + ia.getName() + " (" + numBonds + " of " + valence + ").");
                switch(valence) {
                    case 4:
                        switch(numBonds) {
                            case 3:
                                // Find the average coordinates of atoms ib, ic and id.
                                double[] b = ib.getXYZ(null);
                                double[] c = ib.getXYZ(null);
                                double[] d = ib.getXYZ(null);
                                double[] a = new double[3];
                                a[0] = (b[0] + c[0] + d[0]) / 3.0;
                                a[1] = (b[1] + c[1] + d[1]) / 3.0;
                                a[2] = (b[2] + c[2] + d[2]) / 3.0;
                                // Place the hydrogen at chiral position #1.
                                intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 1);
                                double[] e1 = new double[3];
                                hydrogen.getXYZ(e1);
                                double[] ret = new double[3];
                                diff(a, e1, ret);
                                double l1 = r(ret);
                                // Place the hydrogen at chiral position #2.
                                intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, -1);
                                double[] e2 = new double[3];
                                hydrogen.getXYZ(e2);
                                diff(a, e2, ret);
                                double l2 = r(ret);
                                // Revert to #1 if it is farther from the average.
                                if (l1 > l2) {
                                    hydrogen.setXYZ(e1);
                                }
                                break;
                            case 2:
                                intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 0);
                                break;
                            case 1:
                                intxyz(hydrogen, ia, 1.0, ib, 109.5, null, 0.0, 0);
                                break;
                            case 0:
                                intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
                                break;
                            default:
                                logger.info(" Check biotype for hydrogen " + atomName + ".");
                                patched = false;
                        }
                        break;
                    case 3:
                        switch(numBonds) {
                            case 2:
                                intxyz(hydrogen, ia, 1.0, ib, 120.0, ic, 0.0, 0);
                                break;
                            case 1:
                                intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
                                break;
                            case 0:
                                intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
                                break;
                            default:
                                logger.info(" Check biotype for hydrogen " + atomName + ".");
                                patched = false;
                        }
                        break;
                    case 2:
                        switch(numBonds) {
                            case 1:
                                intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
                                break;
                            case 0:
                                intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
                                break;
                            default:
                                logger.info(" Check biotype for hydrogen " + atomName + ".");
                                patched = false;
                        }
                        break;
                    case 1:
                        switch(numBonds) {
                            case 0:
                                intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
                                break;
                            default:
                                logger.info(" Check biotype for hydrogen " + atomName + ".");
                                patched = false;
                        }
                        break;
                    default:
                        logger.info(" Check biotype for hydrogen " + atomName + ".");
                        patched = false;
                }
                if (!patched) {
                    break;
                } else {
                    buildBond(ia, hydrogen);
                }
            }
        }
        if (!patched) {
            logger.log(Level.WARNING, format(" Deleting unrecognized molecule %s.", m.toString()));
            activeMolecularAssembly.deleteMolecule((Molecule) m);
        } else {
            logger.info(" Patch for " + moleculeName + " succeeded.");
        }
    }
}

Example 92 with Atom

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

the class PotentialEnergyTest method testSoftCore.

public void testSoftCore() {
    logger.info(format(" %s softcore unity test on %s", pmeName, structure.getName()));
    boolean gradient = false;
    boolean print = false;
    double e = forceFieldEnergy.energy(gradient, print);
    Atom[] atoms = molecularAssembly.getAtomArray();
    int n = atoms.length;
    // Make the 3 atoms soft
    for (int i = n; i > n - 3; i--) {
        Atom ai = atoms[i - 1];
        ai.setApplyLambda(true);
    }
    // Compute the energy with Lambda = 1.0;
    double lambda = 1.0;
    forceFieldEnergy.setLambda(lambda);
    double e2 = forceFieldEnergy.energy(gradient, print);
    String msg = String.format("Lambda Unity Test: %g == %g", e, e2);
    assertEquals(msg, e, e2, tolerance);
}

Example 93 with Atom

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

the class MultipoleType method multipoleTypeFactory.

public static MultipoleType multipoleTypeFactory(Atom atom, ForceField forceField) {
    AtomType atomType = atom.getAtomType();
    if (atomType == null) {
        String message = " Multipoles can only be assigned to atoms that have been typed.";
        logger.severe(message);
        return null;
    }
    PolarizeType polarizeType = forceField.getPolarizeType(atomType.getKey());
    if (polarizeType != null) {
        atom.setPolarizeType(polarizeType);
    } else {
        String message = " No polarization type was found for " + atom.toString();
        logger.fine(message);
        double polarizability = 0.0;
        double thole = 0.0;
        int[] polarizationGroup = null;
        polarizeType = new PolarizeType(atomType.type, polarizability, thole, polarizationGroup);
        forceField.addForceFieldType(polarizeType);
        atom.setPolarizeType(polarizeType);
    }
    String key;
    // No reference atoms.
    key = atomType.getKey() + " 0 0";
    MultipoleType multipoleType = forceField.getMultipoleType(key);
    if (multipoleType != null) {
        atom.setMultipoleType(multipoleType);
        atom.setAxisAtoms(null);
        return multipoleType;
    }
    // No bonds.
    List<Bond> bonds = atom.getBonds();
    if (bonds == null || bonds.size() < 1) {
        String message = "Multipoles can only be assigned after bonded relationships are defined.\n";
        logger.severe(message);
    }
    // 1 reference atom.
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        key = atomType.getKey() + " " + atom2.getAtomType().getKey() + " 0";
        multipoleType = forceField.getMultipoleType(key);
        if (multipoleType != null) {
            atom.setMultipoleType(multipoleType);
            atom.setAxisAtoms(atom2);
            return multipoleType;
        }
    }
    // 2 reference atoms.
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        String key2 = atom2.getAtomType().getKey();
        for (Bond b2 : bonds) {
            if (b == b2) {
                continue;
            }
            Atom atom3 = b2.get1_2(atom);
            String key3 = atom3.getAtomType().getKey();
            key = atomType.getKey() + " " + key2 + " " + key3;
            multipoleType = forceField.getMultipoleType(key);
            if (multipoleType != null) {
                atom.setMultipoleType(multipoleType);
                atom.setAxisAtoms(atom2, atom3);
                return multipoleType;
            }
        }
    }
    /**
     * 3 reference atoms.
     */
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        String key2 = atom2.getAtomType().getKey();
        for (Bond b2 : bonds) {
            if (b == b2) {
                continue;
            }
            Atom atom3 = b2.get1_2(atom);
            String key3 = atom3.getAtomType().getKey();
            for (Bond b3 : bonds) {
                if (b == b3 || b2 == b3) {
                    continue;
                }
                Atom atom4 = b3.get1_2(atom);
                String key4 = atom4.getAtomType().getKey();
                key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
                multipoleType = forceField.getMultipoleType(key);
                if (multipoleType != null) {
                    atom.setMultipoleType(multipoleType);
                    atom.setAxisAtoms(atom2, atom3, atom4);
                    return multipoleType;
                }
            }
            List<Angle> angles = atom.getAngles();
            for (Angle angle : angles) {
                Atom atom4 = angle.get1_3(atom);
                if (atom4 != null) {
                    String key4 = atom4.getAtomType().getKey();
                    key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
                    multipoleType = forceField.getMultipoleType(key);
                    if (multipoleType != null) {
                        atom.setMultipoleType(multipoleType);
                        atom.setAxisAtoms(atom2, atom3, atom4);
                        return multipoleType;
                    }
                }
            }
        }
    }
    /**
     * Revert to a 2 reference atom definition that may include a 1-3 site.
     * For example a hydrogen on water.
     */
    for (Bond b : bonds) {
        Atom atom2 = b.get1_2(atom);
        String key2 = atom2.getAtomType().getKey();
        List<Angle> angles = atom.getAngles();
        for (Angle angle : angles) {
            Atom atom3 = angle.get1_3(atom);
            if (atom3 != null) {
                String key3 = atom3.getAtomType().getKey();
                key = atomType.getKey() + " " + key2 + " " + key3;
                multipoleType = forceField.getMultipoleType(key);
                if (multipoleType != null) {
                    atom.setMultipoleType(multipoleType);
                    atom.setAxisAtoms(atom2, atom3);
                    return multipoleType;
                }
                for (Angle angle2 : angles) {
                    Atom atom4 = angle2.get1_3(atom);
                    if (atom4 != null && atom4 != atom3) {
                        String key4 = atom4.getAtomType().getKey();
                        key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
                        multipoleType = forceField.getMultipoleType(key);
                        if (multipoleType != null) {
                            atom.setMultipoleType(multipoleType);
                            atom.setAxisAtoms(atom2, atom3, atom4);
                            return multipoleType;
                        }
                    }
                }
            }
        }
    }
    return null;
}

Example 94 with Atom

use of ffx.potential.bonded.Atom 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 95 with Atom

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

the class ForceFieldEnergy method getVelocity.

/**
 * Returns an array of velocity values for active atoms.
 *
 * @param velocity if the velocity array is null, it will be allocated.
 * @return the velocity array is returned.
 */
@Override
public double[] getVelocity(double[] velocity) {
    int n = getNumberOfVariables();
    if (velocity == null || velocity.length < n) {
        velocity = new double[n];
    }
    int index = 0;
    double[] v = new double[3];
    for (int i = 0; i < nAtoms; i++) {
        Atom a = atoms[i];
        if (a.isActive()) {
            a.getVelocity(v);
            velocity[index++] = v[0];
            velocity[index++] = v[1];
            velocity[index++] = v[2];
        }
    }
    return velocity;
}