Search in sources :

Example 1 with NucleicAcid3

use of ffx.potential.bonded.ResidueEnumerations.NucleicAcid3 in project ffx by mjschnie.

the class NucleicAcidUtils method assignNucleicAcidAtomTypes.

/**
 * Assign atom types for a nucleic acid polymer.
 *
 * @param residues A list of residues that form the nucleic acid polymer.
 * @param forceField The ForceField in use.
 * @param bondList A list of created bonds.
 * @throws MissingHeavyAtomException
 * @throws MissingAtomTypeException
 */
public static void assignNucleicAcidAtomTypes(List<Residue> residues, ForceField forceField, ArrayList<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
    /**
     * A reference to the O3* atom of the previous base.
     */
    Atom pSugarO3 = null;
    /**
     * Loop over residues.
     */
    int numberOfResidues = residues.size();
    for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
        /**
         * Match the residue name to a known nucleic acid residue.
         */
        Residue residue = residues.get(residueNumber);
        String residueName = residue.getName().toUpperCase();
        NucleicAcid3 nucleicAcid = null;
        int naNumber = -1;
        for (NucleicAcid3 nucleic : nucleicAcidList) {
            naNumber++;
            String nuc3 = nucleic.toString();
            nuc3 = nuc3.substring(nuc3.length() - 3);
            if (nuc3.equalsIgnoreCase(residueName)) {
                nucleicAcid = nucleic;
                break;
            }
        }
        /**
         * Do atom name conversions.
         */
        List<Atom> resAtoms = residue.getAtomList();
        int natoms = resAtoms.size();
        for (int i = 0; i < natoms; i++) {
            Atom atom = resAtoms.get(i);
            String name = atom.getName();
            name = name.replace('*', '\'');
            // name = name.replace('D', 'H');
            atom.setName(name);
        }
        /**
         * Check if the sugar is deoxyribose and change the residue name if
         * necessary.
         */
        boolean isDNA = false;
        Atom sugarO2 = (Atom) residue.getAtomNode("O2\'");
        if (sugarO2 == null) {
            /**
             * Assume deoxyribose (DNA) since there is an O2* atom.
             */
            isDNA = true;
            if (!residueName.startsWith("D")) {
                switch(nucleicAcid) {
                    case ADE:
                        nucleicAcid = NucleicAcid3.DAD;
                        residueName = "DAD";
                        residue.setName(residueName);
                        break;
                    case CYT:
                        nucleicAcid = NucleicAcid3.DCY;
                        residueName = "DCY";
                        residue.setName(residueName);
                        break;
                    case GUA:
                        nucleicAcid = NucleicAcid3.DGU;
                        residueName = "DGU";
                        residue.setName(residueName);
                        break;
                    case THY:
                        nucleicAcid = NucleicAcid3.DTY;
                        residueName = "DTY";
                        residue.setName(residueName);
                        break;
                    default:
                }
            }
        } else /**
         * Assume ribose (RNA) since there is an O2* atom.
         */
        {
            if (residueName.startsWith("D")) {
                switch(nucleicAcid) {
                    case DAD:
                        nucleicAcid = NucleicAcid3.ADE;
                        residueName = "ADE";
                        residue.setName(residueName);
                        break;
                    case DCY:
                        nucleicAcid = NucleicAcid3.CYT;
                        residueName = "CYT";
                        residue.setName(residueName);
                        break;
                    case DGU:
                        nucleicAcid = NucleicAcid3.GUA;
                        residueName = "GUA";
                        residue.setName(residueName);
                        break;
                    case DTY:
                        nucleicAcid = NucleicAcid3.THY;
                        residueName = "THY";
                        residue.setName(residueName);
                        break;
                    default:
                }
            }
        }
        /**
         * Set a position flag.
         */
        Residue.ResiduePosition position = MIDDLE_RESIDUE;
        if (residueNumber == 0) {
            position = FIRST_RESIDUE;
        } else if (residueNumber == numberOfResidues - 1) {
            position = LAST_RESIDUE;
        }
        /**
         * Build the phosphate atoms of the current residue.
         */
        Atom phosphate = null;
        Atom sugarO5 = null;
        if (position == FIRST_RESIDUE) {
            /**
             * The 5' O5' oxygen of the nucleic acid is generally terminated
             * by 1.) A phosphate group PO3 (-3). 2.) A hydrogen.
             *
             * If the base has phosphate atom we will assume a PO3 group.
             */
            phosphate = (Atom) residue.getAtomNode("P");
            if (phosphate != null) {
                if (isDNA) {
                    phosphate = buildHeavy(residue, "P", null, 1247, forceField, bondList);
                    buildHeavy(residue, "OP1", phosphate, 1248, forceField, bondList);
                    buildHeavy(residue, "OP2", phosphate, 1248, forceField, bondList);
                    buildHeavy(residue, "OP3", phosphate, 1248, forceField, bondList);
                    sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1246, forceField, bondList);
                } else {
                    phosphate = buildHeavy(residue, "P", null, 1235, forceField, bondList);
                    buildHeavy(residue, "OP1", phosphate, 1236, forceField, bondList);
                    buildHeavy(residue, "OP2", phosphate, 1236, forceField, bondList);
                    buildHeavy(residue, "OP3", phosphate, 1236, forceField, bondList);
                    sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1234, forceField, bondList);
                }
            } else if (isDNA) {
                sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1244, forceField, bondList);
            } else {
                sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1232, forceField, bondList);
            }
        } else {
            phosphate = buildHeavy(residue, "P", pSugarO3, NA_P[naNumber], forceField, bondList);
            buildHeavy(residue, "OP1", phosphate, NA_OP[naNumber], forceField, bondList);
            buildHeavy(residue, "OP2", phosphate, NA_OP[naNumber], forceField, bondList);
            sugarO5 = buildHeavy(residue, "O5\'", phosphate, NA_O5[naNumber], forceField, bondList);
        }
        /**
         * Build the ribose sugar atoms of the current base.
         */
        Atom sugarC5 = buildHeavy(residue, "C5\'", sugarO5, NA_C5[naNumber], forceField, bondList);
        Atom sugarC4 = buildHeavy(residue, "C4\'", sugarC5, NA_C4[naNumber], forceField, bondList);
        Atom sugarO4 = buildHeavy(residue, "O4\'", sugarC4, NA_O4[naNumber], forceField, bondList);
        Atom sugarC1 = buildHeavy(residue, "C1\'", sugarO4, NA_C1[naNumber], forceField, bondList);
        Atom sugarC3 = buildHeavy(residue, "C3\'", sugarC4, NA_C3[naNumber], forceField, bondList);
        Atom sugarC2 = buildHeavy(residue, "C2\'", sugarC3, NA_C2[naNumber], forceField, bondList);
        buildBond(sugarC2, sugarC1, forceField, bondList);
        Atom sugarO3 = null;
        if (position == LAST_RESIDUE || numberOfResidues == 1) {
            if (isDNA) {
                sugarO3 = buildHeavy(residue, "O3\'", sugarC3, 1249, forceField, bondList);
            } else {
                sugarO3 = buildHeavy(residue, "O3\'", sugarC3, 1237, forceField, bondList);
            }
        } else {
            sugarO3 = buildHeavy(residue, "O3\'", sugarC3, NA_O3[naNumber], forceField, bondList);
        }
        if (!isDNA) {
            sugarO2 = buildHeavy(residue, "O2\'", sugarC2, NA_O2[naNumber], forceField, bondList);
        }
        /**
         * Build the backbone hydrogen atoms.
         */
        if (position == FIRST_RESIDUE && NA_P == null) {
            buildHydrogen(residue, "H5T", sugarO5, 1.00e0, sugarC5, 109.5e0, sugarC4, 180.0e0, 0, NA_H5T[naNumber], forceField, bondList);
        }
        buildHydrogen(residue, "H5\'1", sugarC5, 1.09e0, sugarO5, 109.5e0, sugarC4, 109.5e0, 1, NA_H51[naNumber], forceField, bondList);
        buildHydrogen(residue, "H5\'2", sugarC5, 1.09e0, sugarO5, 109.5e0, sugarC4, 109.5e0, -1, NA_H52[naNumber], forceField, bondList);
        buildHydrogen(residue, "H4\'", sugarC4, 1.09e0, sugarC5, 109.5e0, sugarC3, 109.5e0, -1, NA_H4[naNumber], forceField, bondList);
        buildHydrogen(residue, "H3\'", sugarC3, 1.09e0, sugarC4, 109.5e0, sugarC2, 109.5e0, -1, NA_H3[naNumber], forceField, bondList);
        if (isDNA) {
            buildHydrogen(residue, "H2\'1", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, -1, NA_H21[naNumber], forceField, bondList);
            buildHydrogen(residue, "H2\'2", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, 1, NA_H22[naNumber], forceField, bondList);
        } else {
            buildHydrogen(residue, "H2\'", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, -1, NA_H21[naNumber], forceField, bondList);
            // Add the NA_O2' Methyl for OMC and OMG
            if (nucleicAcid == NucleicAcid3.OMC || nucleicAcid == NucleicAcid3.OMG) {
                Atom CM2 = buildHeavy(residue, "CM2", sugarO2, 1427, forceField, bondList);
                Atom HM21 = buildHydrogen(residue, "HM21", CM2, 1.08e0, sugarO2, 109.5e0, sugarC2, 0.0e0, 0, 1428, forceField, bondList);
                buildHydrogen(residue, "HM22", CM2, 1.08e0, sugarO2, 109.5e0, HM21, 109.5e0, 1, 1429, forceField, bondList);
                buildHydrogen(residue, "HM23", CM2, 1.08e0, sugarO2, 109.5e0, HM21, 109.5e0, -1, 1430, forceField, bondList);
            } else {
                buildHydrogen(residue, "HO\'", sugarO2, 1.00e0, sugarC2, 109.5e0, sugarC3, 180.0e0, 0, NA_H22[naNumber], forceField, bondList);
            }
        }
        buildHydrogen(residue, "H1\'", sugarC1, 1.09e0, sugarO4, 109.5e0, sugarC2, 109.5e0, -1, NA_H1[naNumber], forceField, bondList);
        if (position == LAST_RESIDUE || numberOfResidues == 1) {
            buildHydrogen(residue, "H3T", sugarO3, 1.00e0, sugarC3, 109.5e0, sugarC4, 180.0e0, 0, NA_H3T[naNumber], forceField, bondList);
        // Else, if it is terminated by a 3' phosphate cap:
        // Will need to see how PDB would label a 3' phosphate cap.
        }
        /**
         * Build the nucleic acid base.
         */
        try {
            assignNucleicAcidBaseAtomTypes(nucleicAcid, residue, sugarC1, sugarO4, sugarC2, forceField, bondList);
        } catch (MissingHeavyAtomException missingHeavyAtomException) {
            throw missingHeavyAtomException;
        }
        /**
         * Do some checks on the current base to make sure all atoms have
         * been assigned an atom type.
         */
        resAtoms = residue.getAtomList();
        for (Atom atom : resAtoms) {
            AtomType atomType = atom.getAtomType();
            if (atomType == null) {
                MissingAtomTypeException missingAtomTypeException = new MissingAtomTypeException(residue, atom);
                throw missingAtomTypeException;
            }
            int numberOfBonds = atom.getNumBonds();
            if (numberOfBonds != atomType.valence) {
                if (atom == sugarO3 && numberOfBonds == atomType.valence - 1 && position != LAST_RESIDUE && numberOfResidues != 1) {
                    continue;
                }
                logger.log(Level.WARNING, format(" An atom for residue %s has the wrong number of bonds:\n %s", residueName, atom.toString()));
                logger.log(Level.WARNING, format(" Expected: %d Actual: %d.", atomType.valence, numberOfBonds));
            }
        }
        /**
         * Save a reference to the current O3* oxygen.
         */
        pSugarO3 = sugarO3;
    }
}
Also used : MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) NucleicAcid3(ffx.potential.bonded.ResidueEnumerations.NucleicAcid3) AtomType(ffx.potential.parameters.AtomType) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException)

Example 2 with NucleicAcid3

use of ffx.potential.bonded.ResidueEnumerations.NucleicAcid3 in project ffx by mjschnie.

the class PDBFilter method assignAtomTypes.

/**
 * Assign force field atoms types to common chemistries using "biotype"
 * records.
 */
private void assignAtomTypes() {
    /**
     * Create a list to store bonds defined by 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;
            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;
                        checkHydrogenAtomNames(residue);
                        break;
                    }
                }
                // Check for a patch.
                if (!aa) {
                    logger.info(" Checking for non-standard amino acid patch " + name);
                    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, forceField, bondList);
                } catch (MissingHeavyAtomException | MissingAtomTypeException e) {
                    logger.severe(e.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 {
                logger.fine(" " + atom.toString() + " -> " + atomType.toString());
                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.");
        }
    }
    resolvePolymerLinks(molecules);
}
Also used : NucleicAcid3(ffx.potential.bonded.ResidueEnumerations.NucleicAcid3) HashMap(java.util.HashMap) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) MSNode(ffx.potential.bonded.MSNode) AtomType(ffx.potential.parameters.AtomType) ResidueEnumerations.nucleicAcidList(ffx.potential.bonded.ResidueEnumerations.nucleicAcidList) ResidueEnumerations.aminoAcidList(ffx.potential.bonded.ResidueEnumerations.aminoAcidList) List(java.util.List) ArrayList(java.util.ArrayList) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) Polymer(ffx.potential.bonded.Polymer) Atom(ffx.potential.bonded.Atom) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) IOException(java.io.IOException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) Molecule(ffx.potential.bonded.Molecule) Residue(ffx.potential.bonded.Residue) Bond(ffx.potential.bonded.Bond)

Example 3 with NucleicAcid3

use of ffx.potential.bonded.ResidueEnumerations.NucleicAcid3 in project ffx by mjschnie.

the class BiojavaFilter method assignNucleicAcidAtomTypes.

/**
 * Assign atom types for a nucleic acid polymer.
 *
 * @param residues
 * @throws ffx.potential.parsers.PDBFilter.MissingHeavyAtomException
 */
private void assignNucleicAcidAtomTypes(List<Residue> residues) throws MissingHeavyAtomException, MissingAtomTypeException {
    /**
     * A reference to the O3* atom of the previous base.
     */
    Atom pO3s = null;
    /**
     * Loop over residues.
     */
    int numberOfResidues = residues.size();
    for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
        /**
         * Match the residue name to a known nucleic acid residue.
         */
        Residue residue = residues.get(residueNumber);
        String residueName = residue.getName().toUpperCase();
        NucleicAcid3 nucleicAcid = null;
        int naNumber = -1;
        for (NucleicAcid3 nucleic : nucleicAcidList) {
            naNumber++;
            String nuc3 = nucleic.toString();
            nuc3 = nuc3.substring(nuc3.length() - 3);
            if (nuc3.equalsIgnoreCase(residueName)) {
                nucleicAcid = nucleic;
                break;
            }
        }
        /**
         * Do atom name conversions.
         */
        List<Atom> resAtoms = residue.getAtomList();
        int natoms = resAtoms.size();
        for (int i = 0; i < natoms; i++) {
            Atom atom = resAtoms.get(i);
            String name = atom.getName();
            name = name.replace('*', '\'');
            // name = name.replace('D', 'H');
            atom.setName(name);
        }
        /**
         * Check if this is a 3' phosphate being listed as its own residue.
         */
        /*if (residue.getAtomList().size() == 1) {
             Atom P3s = (Atom) residue.getAtomNode("NA_P");
             if (P3s != null) {
             Residue prevResidue = residue.getPreviousResidue();
             if (prevResidue != null) {
             Atom O2sPrev = (Atom) prevResidue.getAtomNode("NA_O2\'");
             if (O2sPrev == null) {
             P3s = buildHeavy(prevResidue, "P3s", null, 1247);
             } else {
             P3s = buildHeavy(prevResidue, "P3s", null, 1235);
             }
             } else {
             return;
             }
             } else {
             return;
             }
             }*/
        /**
         * Check if the sugar is deoxyribose and change the residue name if
         * necessary.
         */
        boolean isDNA = false;
        Atom O2s = (Atom) residue.getAtomNode("O2\'");
        if (O2s == null) {
            /**
             * Assume deoxyribose (DNA) since there is an O2* atom.
             */
            isDNA = true;
            if (!residueName.startsWith("D")) {
                switch(nucleicAcid) {
                    case ADE:
                        nucleicAcid = NucleicAcid3.DAD;
                        residueName = "DAD";
                        residue.setName(residueName);
                        break;
                    case CYT:
                        nucleicAcid = NucleicAcid3.DCY;
                        residueName = "DCY";
                        residue.setName(residueName);
                        break;
                    case GUA:
                        nucleicAcid = NucleicAcid3.DGU;
                        residueName = "DGU";
                        residue.setName(residueName);
                        break;
                    case THY:
                        nucleicAcid = NucleicAcid3.DTY;
                        residueName = "DTY";
                        residue.setName(residueName);
                        break;
                    default:
                }
            }
        } else /**
         * Assume ribose (RNA) since there is an O2* atom.
         */
        if (residueName.startsWith("D")) {
            switch(nucleicAcid) {
                case DAD:
                    nucleicAcid = NucleicAcid3.ADE;
                    residueName = "ADE";
                    residue.setName(residueName);
                    break;
                case DCY:
                    nucleicAcid = NucleicAcid3.CYT;
                    residueName = "CYT";
                    residue.setName(residueName);
                    break;
                case DGU:
                    nucleicAcid = NucleicAcid3.GUA;
                    residueName = "GUA";
                    residue.setName(residueName);
                    break;
                case DTY:
                    nucleicAcid = NucleicAcid3.THY;
                    residueName = "THY";
                    residue.setName(residueName);
                    break;
                default:
            }
        }
        /**
         * Set a position flag.
         */
        ResiduePosition position = MIDDLE_RESIDUE;
        if (residueNumber == 0) {
            position = FIRST_RESIDUE;
        } else if (residueNumber == numberOfResidues - 1) {
            position = LAST_RESIDUE;
        }
        /**
         * Build the phosphate atoms of the current residue.
         */
        Atom phosphate = null;
        Atom O5s = null;
        if (position == FIRST_RESIDUE) {
            /**
             * The 5' O5' oxygen of the nucleic acid is generally terminated
             * by 1.) A phosphate group PO3 (-3). 2.) A hydrogen.
             *
             * If the base has phosphate atom we will assume a PO3 group.
             */
            phosphate = (Atom) residue.getAtomNode("P");
            if (phosphate != null) {
                if (isDNA) {
                    phosphate = buildHeavy(residue, "P", null, 1247);
                    buildHeavy(residue, "OP1", phosphate, 1248);
                    buildHeavy(residue, "OP2", phosphate, 1248);
                    buildHeavy(residue, "OP3", phosphate, 1248);
                    O5s = buildHeavy(residue, "O5\'", phosphate, 1246);
                } else {
                    phosphate = buildHeavy(residue, "P", null, 1235);
                    buildHeavy(residue, "OP1", phosphate, 1236);
                    buildHeavy(residue, "OP2", phosphate, 1236);
                    buildHeavy(residue, "OP3", phosphate, 1236);
                    O5s = buildHeavy(residue, "O5\'", phosphate, 1234);
                }
            } else if (isDNA) {
                O5s = buildHeavy(residue, "O5\'", phosphate, 1244);
            } else {
                O5s = buildHeavy(residue, "O5\'", phosphate, 1232);
            }
        } else {
            phosphate = buildHeavy(residue, "P", pO3s, NA_P[naNumber]);
            buildHeavy(residue, "OP1", phosphate, NA_OP[naNumber]);
            buildHeavy(residue, "OP2", phosphate, NA_OP[naNumber]);
            O5s = buildHeavy(residue, "O5\'", phosphate, NA_O5[naNumber]);
        }
        /**
         * Build the ribose sugar atoms of the current base.
         */
        Atom C5s = buildHeavy(residue, "C5\'", O5s, NA_C5[naNumber]);
        Atom C4s = buildHeavy(residue, "C4\'", C5s, NA_C4[naNumber]);
        Atom O4s = buildHeavy(residue, "O4\'", C4s, NA_O4[naNumber]);
        Atom C1s = buildHeavy(residue, "C1\'", O4s, NA_C1[naNumber]);
        Atom C3s = buildHeavy(residue, "C3\'", C4s, NA_C3[naNumber]);
        Atom C2s = buildHeavy(residue, "C2\'", C3s, NA_C2[naNumber]);
        buildBond(C2s, C1s);
        Atom O3s = null;
        if (position == LAST_RESIDUE || numberOfResidues == 1) {
            if (isDNA) {
                O3s = buildHeavy(residue, "O3\'", C3s, 1249);
            } else {
                O3s = buildHeavy(residue, "O3\'", C3s, 1237);
            }
        } else {
            O3s = buildHeavy(residue, "O3\'", C3s, NA_O3[naNumber]);
        }
        if (!isDNA) {
            O2s = buildHeavy(residue, "O2\'", C2s, NA_O2[naNumber]);
        }
        /**
         * Build the backbone hydrogen atoms.
         */
        if (position == FIRST_RESIDUE && phosphate == null) {
            buildHydrogen(residue, "H5T", O5s, 1.00e0, C5s, 109.5e0, C4s, 180.0e0, 0, NA_H5T[naNumber]);
        }
        buildHydrogen(residue, "H5\'1", C5s, 1.09e0, O5s, 109.5e0, C4s, 109.5e0, 1, NA_H51[naNumber]);
        buildHydrogen(residue, "H5\'2", C5s, 1.09e0, O5s, 109.5e0, C4s, 109.5e0, -1, NA_H52[naNumber]);
        buildHydrogen(residue, "H4\'", C4s, 1.09e0, C5s, 109.5e0, C3s, 109.5e0, -1, NA_H4[naNumber]);
        buildHydrogen(residue, "H3\'", C3s, 1.09e0, C4s, 109.5e0, C2s, 109.5e0, -1, NA_H3[naNumber]);
        if (isDNA) {
            buildHydrogen(residue, "H2\'1", C2s, 1.09e0, C3s, 109.5e0, C1s, 109.5e0, -1, NA_H21[naNumber]);
            buildHydrogen(residue, "H2\'2", C2s, 1.09e0, C3s, 109.5e0, C1s, 109.5e0, 1, NA_H22[naNumber]);
        } else {
            buildHydrogen(residue, "H2\'", C2s, 1.09e0, C3s, 109.5e0, C1s, 109.5e0, -1, NA_H21[naNumber]);
            // Add the NA_O2' Methyl for OMC and OMG
            if (nucleicAcid == NucleicAcid3.OMC || nucleicAcid == NucleicAcid3.OMG) {
                Atom CM2 = buildHeavy(residue, "CM2", O2s, 1427);
                Atom HM21 = buildHydrogen(residue, "HM21", CM2, 1.08e0, O2s, 109.5e0, C2s, 0.0e0, 0, 1428);
                buildHydrogen(residue, "HM22", CM2, 1.08e0, O2s, 109.5e0, HM21, 109.5e0, 1, 1429);
                buildHydrogen(residue, "HM23", CM2, 1.08e0, O2s, 109.5e0, HM21, 109.5e0, -1, 1430);
            } else {
                buildHydrogen(residue, "HO\'", O2s, 1.00e0, C2s, 109.5e0, C3s, 180.0e0, 0, NA_H22[naNumber]);
            }
        }
        buildHydrogen(residue, "H1\'", C1s, 1.09e0, O4s, 109.5e0, C2s, 109.5e0, -1, NA_H1[naNumber]);
        if (position == LAST_RESIDUE || numberOfResidues == 1) {
            buildHydrogen(residue, "H3T", O3s, 1.00e0, C3s, 109.5e0, C4s, 180.0e0, 0, NA_H3T[naNumber]);
        // Else, if it is terminated by a 3' phosphate cap:
        // Will need to see how PDB would label a 3' phosphate cap.
        }
        /**
         * Build the nucleic acid base.
         */
        try {
            assignNucleicAcidBaseAtomTypes(nucleicAcid, residue, C1s, O4s, C2s);
        } catch (MissingHeavyAtomException missingHeavyAtomException) {
            logger.throwing(PDBFilter.class.getName(), "assignNucleicAcidAtomTypes", missingHeavyAtomException);
            throw missingHeavyAtomException;
        }
        /**
         * Do some checks on the current base to make sure all atoms have
         * been assigned an atom type.
         */
        resAtoms = residue.getAtomList();
        for (Atom atom : resAtoms) {
            AtomType atomType = atom.getAtomType();
            if (atomType == null) {
                MissingAtomTypeException missingAtomTypeException = new MissingAtomTypeException(residue, atom);
                logger.throwing(PDBFilter.class.getName(), "assignNucleicAcidAtomTypes", missingAtomTypeException);
                throw missingAtomTypeException;
            }
            int numberOfBonds = atom.getNumBonds();
            if (numberOfBonds != atomType.valence) {
                if (atom == O3s && numberOfBonds == atomType.valence - 1 && position != LAST_RESIDUE && numberOfResidues != 1) {
                    continue;
                }
                logger.log(Level.WARNING, format(" An atom for residue %s has the wrong number of bonds:\n %s", residueName, atom.toString()));
                logger.log(Level.WARNING, format(" Expected: %d Actual: %d.", atomType.valence, numberOfBonds));
            }
        }
        /**
         * Save a reference to the current O3* oxygen.
         */
        pO3s = O3s;
    }
}
Also used : ResiduePosition(ffx.potential.bonded.Residue.ResiduePosition) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) Residue(ffx.potential.bonded.Residue) NucleicAcid3(ffx.potential.bonded.ResidueEnumerations.NucleicAcid3) AtomType(ffx.potential.parameters.AtomType) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) Atom(ffx.potential.bonded.Atom)

Example 4 with NucleicAcid3

use of ffx.potential.bonded.ResidueEnumerations.NucleicAcid3 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.");
        }
    }
}
Also used : NucleicAcid3(ffx.potential.bonded.ResidueEnumerations.NucleicAcid3) HashMap(java.util.HashMap) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) MSNode(ffx.potential.bonded.MSNode) AtomType(ffx.potential.parameters.AtomType) ResidueEnumerations.nucleicAcidList(ffx.potential.bonded.ResidueEnumerations.nucleicAcidList) ResidueEnumerations.aminoAcidList(ffx.potential.bonded.ResidueEnumerations.aminoAcidList) ArrayList(java.util.ArrayList) List(java.util.List) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) Polymer(ffx.potential.bonded.Polymer) HetAtoms(ffx.potential.parsers.PDBFilter.HetAtoms) Atom(ffx.potential.bonded.Atom) IOException(java.io.IOException) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) Molecule(ffx.potential.bonded.Molecule) Residue(ffx.potential.bonded.Residue) Bond(ffx.potential.bonded.Bond) SSBond(org.biojava.bio.structure.SSBond)

Example 5 with NucleicAcid3

use of ffx.potential.bonded.ResidueEnumerations.NucleicAcid3 in project ffx by mjschnie.

the class Residue method getRotamers.

/**
 * Gets the Rotamers for this residue, potentially incorporating the
 * original coordinates if RotamerLibrary's original coordinates rotamer
 * flag has been set.
 *
 * @param library Rotamer library to use
 * @return An array of Rotamer.
 */
public Rotamer[] getRotamers(RotamerLibrary library) {
    /**
     * If the rotamers for this residue have been cached, return them.
     */
    if (rotamers != null) {
        return rotamers;
    }
    /**
     * Return rotamers for this residue from the RotamerLibrary.
     */
    Rotamer[] libRotamers = library.getRotamers(this);
    /**
     * If there are no rotamers, and addOrigRot is true, return an array with
     * only an original-coordinates rotamer. Else if there are no rotamers,
     * return (null) library rotamers. If there are rotamers, and original
     * coordinates are turned off, return (filled) library rotamers. Else,
     * continue generating the rotamers array.
     */
    if (libRotamers == null) {
        if (addOrigRot) {
            rotamers = new Rotamer[1];
            ResidueState origState = this.storeState();
            double[] chi = RotamerLibrary.measureRotamer(this, false);
            switch(residueType) {
                case AA:
                    AminoAcid3 aa3 = AminoAcid3.UNK;
                    try {
                        aa3 = AminoAcid3.valueOf(getName());
                    } catch (Exception e) {
                    }
                    originalRotamer = new Rotamer(aa3, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
                    rotamers[0] = originalRotamer;
                    break;
                case NA:
                    NucleicAcid3 na3 = NucleicAcid3.UNK;
                    try {
                        na3 = NucleicAcid3.valueOf(getName());
                    } catch (Exception e) {
                    }
                    originalRotamer = new Rotamer(na3, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0, chi[4], 0, chi[5], 0);
                    rotamers[0] = originalRotamer;
                    break;
                default:
                    originalRotamer = null;
                    // Resets to null.
                    rotamers = libRotamers;
            }
            return rotamers;
        } else {
            return libRotamers;
        }
    } else if (!library.getUsingOrigCoordsRotamer()) {
        return libRotamers;
    }
    /**
     * Define the current coordinates as a new rotamer.
     */
    ResidueState origState = storeState();
    double[] chi = RotamerLibrary.measureRotamer(this, false);
    switch(residueType) {
        case AA:
            AminoAcid3 aa3 = this.getAminoAcid3();
            originalRotamer = new Rotamer(aa3, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
            break;
        case NA:
            NucleicAcid3 na3 = this.getNucleicAcid3();
            originalRotamer = new Rotamer(na3, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0, chi[4], 0, chi[5], 0);
            break;
        default:
            double[] rotaValues = new double[chi.length * 2];
            for (int i = 0; i < chi.length; i++) {
                int ii = i * 2;
                rotaValues[ii] = chi[i];
                rotaValues[ii + 1] = 0.0;
            }
            originalRotamer = new Rotamer(origState, rotaValues);
            break;
    }
    /**
     * Add the new rotamer to those from the library and cache the result.
     */
    int nRots = libRotamers.length;
    rotamers = new Rotamer[nRots + 1];
    if (origAtEnd) {
        System.arraycopy(libRotamers, 0, rotamers, 0, nRots);
        rotamers[rotamers.length - 1] = originalRotamer;
    } else {
        System.arraycopy(libRotamers, 0, rotamers, 1, nRots);
        rotamers[0] = originalRotamer;
    }
    /**
     * Return the rotamer array.
     */
    return rotamers;
}
Also used : NucleicAcid3(ffx.potential.bonded.ResidueEnumerations.NucleicAcid3) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3)

Aggregations

NucleicAcid3 (ffx.potential.bonded.ResidueEnumerations.NucleicAcid3)8 MissingAtomTypeException (ffx.potential.bonded.BondedUtils.MissingAtomTypeException)4 MissingHeavyAtomException (ffx.potential.bonded.BondedUtils.MissingHeavyAtomException)4 AminoAcid3 (ffx.potential.bonded.ResidueEnumerations.AminoAcid3)4 AtomType (ffx.potential.parameters.AtomType)4 Atom (ffx.potential.bonded.Atom)3 Residue (ffx.potential.bonded.Residue)3 ArrayList (java.util.ArrayList)3 Bond (ffx.potential.bonded.Bond)2 MSNode (ffx.potential.bonded.MSNode)2 Molecule (ffx.potential.bonded.Molecule)2 Polymer (ffx.potential.bonded.Polymer)2 ResidueEnumerations.aminoAcidList (ffx.potential.bonded.ResidueEnumerations.aminoAcidList)2 ResidueEnumerations.nucleicAcidList (ffx.potential.bonded.ResidueEnumerations.nucleicAcidList)2 IOException (java.io.IOException)2 HashMap (java.util.HashMap)2 List (java.util.List)2 ResiduePosition (ffx.potential.bonded.Residue.ResiduePosition)1 HetAtoms (ffx.potential.parsers.PDBFilter.HetAtoms)1 SSBond (org.biojava.bio.structure.SSBond)1