Search in sources :

Example 1 with AtomType

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

the class SimulationFilter method readFile.

/**
 * {@inheritDoc}
 */
@Override
public boolean readFile() {
    // information
    for (int i = 0; i < system.numatoms; i++) {
        AtomType atomType = atomTypes.get(system.types[i]);
        if (atomType == null) {
            atomType = new AtomType(system.types[i], -1, system.name[i], system.story[i], system.atomic[i], system.mass[i], 0);
            atomTypes.put(system.types[i], atomType);
        }
    }
    atomList = new ArrayList<Atom>();
    Vector<Integer> bonds1 = new Vector<Integer>();
    Vector<Integer> bonds2 = new Vector<Integer>();
    double[] d = new double[3];
    int[] b = new int[4];
    for (int i = 0; i < system.numatoms; i++) {
        d[0] = system.coordinates[0][i];
        d[1] = system.coordinates[1][i];
        d[2] = system.coordinates[2][i];
        String s = new String("" + system.types[i]);
        AtomType atomType = atomTypes.get(s);
        Atom a = new Atom(i + 1, new String("" + atomType.type), atomType, d);
        atomList.add(a);
        int b1 = i + 1;
        b[0] = system.connectivity[0][i];
        b[1] = system.connectivity[1][i];
        b[2] = system.connectivity[2][i];
        b[3] = system.connectivity[3][i];
        int j = 0;
        while (j < 4 && b[j] != 0) {
            int b2 = b[j];
            bonds1.add(b1);
            bonds2.add(b2);
            j++;
        }
    }
    bondList = new ArrayList<Bond>();
    for (int i = 0; i < bonds1.size(); i++) {
        int a1 = bonds1.get(i);
        int a2 = bonds2.get(i);
        if (a1 < a2) {
            Atom atom1 = atomList.get(a1 - 1);
            Atom atom2 = atomList.get(a2 - 1);
            bondList.add(new Bond(atom1, atom2));
        }
    }
    setFileRead(true);
    return true;
}
Also used : AtomType(ffx.potential.parameters.AtomType) Bond(ffx.potential.bonded.Bond) Vector(java.util.Vector) Atom(ffx.potential.bonded.Atom)

Example 2 with AtomType

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

the class AminoAcidUtils method assignAminoAcidSideChain.

/**
 * Assign atom types to a single amino acid side chain.
 *
 * @param position The position of this amino acid in the chain.
 * @param aminoAcid The amino acid to use.
 * @param residue The residue node.
 * @param CA The C-alpha carbon of this residue.
 * @param N The peptide nitrogen of this residue.
 * @param C The peptide carbonyl carbon.
 * @param forceField
 * @param bondList
 * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException this
 * exception is thrown if when heavy is atom is missing that cannot be
 * built.
 */
public static void assignAminoAcidSideChain(ResiduePosition position, AminoAcid3 aminoAcid, Residue residue, Atom CA, Atom N, Atom C, ForceField forceField, ArrayList<Bond> bondList) throws MissingHeavyAtomException {
    int k = AA_CB[aminoAcid.ordinal()];
    switch(aminoAcid) {
        case GLY:
            switch(position) {
                case FIRST_RESIDUE:
                    k = AA_HA[0][k];
                    break;
                case LAST_RESIDUE:
                    k = AA_HA[2][k];
                    break;
                default:
                    k = AA_HA[1][k];
            }
            buildHydrogen(residue, "HA3", CA, 1.10, N, 109.5, C, 109.5, 1, k, forceField, bondList);
            break;
        case ALA:
            buildAlanine(residue, CA, N, C, k, forceField, bondList);
            break;
        case VAL:
            buildValine(residue, CA, N, C, k, forceField, bondList);
            break;
        case LEU:
            buildLeucine(residue, CA, N, C, k, forceField, bondList);
            break;
        case ILE:
            buildIsoleucine(residue, CA, N, C, k, forceField, bondList);
            break;
        case SER:
            buildSerine(residue, CA, N, C, k, forceField, bondList);
            break;
        case THR:
            buildThreonine(residue, CA, N, C, k, forceField, bondList);
            break;
        case CYS:
            buildCysteine(residue, CA, N, C, k, forceField, bondList);
            break;
        case CYX:
            buildCystine(residue, CA, N, C, k, forceField, bondList);
            break;
        case CYD:
            buildDeprotonatedCysteine(residue, CA, N, C, k, forceField, bondList);
            break;
        case PRO:
            buildProline(residue, CA, N, C, k, position, forceField, bondList);
            break;
        case PHE:
            buildPhenylalanine(residue, CA, N, C, k, forceField, bondList);
            break;
        case TYR:
            buildTyrosine(residue, CA, N, C, k, forceField, bondList);
            break;
        case TYD:
            buildDeprotonatedTyrosine(residue, CA, N, C, k, forceField, bondList);
            break;
        case TRP:
            buildTryptophan(residue, CA, N, C, k, forceField, bondList);
            break;
        case HIS:
            buildHistidine(residue, CA, N, C, k, forceField, bondList);
            break;
        case HID:
            buildNeutralHistidineD(residue, CA, N, C, k, forceField, bondList);
            break;
        case HIE:
            buildNeutralHistidineE(residue, CA, N, C, k, forceField, bondList);
            break;
        case ASP:
            buildAspartate(residue, CA, N, C, k, forceField, bondList);
            break;
        case ASH:
            buildNeutralAsparticAcid(residue, CA, N, C, k, forceField, bondList);
            break;
        case ASN:
            buildAsparagine(residue, CA, N, C, k, forceField, bondList);
            break;
        case GLU:
            buildGlutamate(residue, CA, N, C, k, forceField, bondList);
            break;
        case GLH:
            buildNeutralGlutamicAcid(residue, CA, N, C, k, forceField, bondList);
            break;
        case GLN:
            buildGlutamine(residue, CA, N, C, k, forceField, bondList);
            break;
        case MET:
            buildMethionine(residue, CA, N, C, k, forceField, bondList);
            break;
        case LYS:
            buildLysine(residue, CA, N, C, k, forceField, bondList);
            break;
        case LYD:
            buildDeprotonatedLysine(residue, CA, N, C, k, forceField, bondList);
            break;
        case ARG:
            buildArginine(residue, CA, N, C, k, forceField, bondList);
            break;
        case ORN:
            buildOrnithine(residue, CA, N, C, k, forceField, bondList);
            break;
        case AIB:
            buildAIB(residue, CA, N, C, k, forceField, bondList);
            break;
        case PCA:
            buildPCA(residue, CA, N, C, k, forceField, bondList);
            break;
        case UNK:
            String residueName = residue.getName();
            logger.log(Level.INFO, " Patching side-chain {0}", residueName);
            HashMap<String, AtomType> types = forceField.getAtomTypes(residueName);
            if (!types.isEmpty()) {
                boolean patched = true;
                ArrayList<Atom> residueAtoms = residue.getAtomList();
                // Assign atom types for side-chain atoms.
                for (Atom atom : residueAtoms) {
                    String atomName = atom.getName().toUpperCase();
                    AtomType type = atom.getAtomType();
                    if (type == null) {
                        type = types.get(atomName);
                        atom.setAtomType(type);
                        types.remove(atomName);
                    }
                }
                // Create bonds between known atoms.
                if (patched) {
                    for (Atom atom : residueAtoms) {
                        String atomName = atom.getName();
                        String[] bonds = forceField.getBonds(residueName, atomName);
                        if (bonds != null) {
                            for (String name : bonds) {
                                Atom atom2 = (Atom) residue.getAtomNode(name);
                                if (atom2 != null && !atom.isBonded(atom2)) {
                                    buildBond(atom, atom2, forceField, bondList);
                                }
                            }
                        }
                    }
                }
                // Create missing hydrogen atoms.
                if (patched && !types.isEmpty()) {
                    // Create a hashmap of the molecule's atoms
                    HashMap<String, Atom> atomMap = new HashMap<>();
                    for (Atom atom : residueAtoms) {
                        atomMap.put(atom.getName().toUpperCase(), atom);
                    }
                    for (String atomName : types.keySet()) {
                        AtomType type = types.get(atomName);
                        String[] bonds = forceField.getBonds(residueName, atomName.toUpperCase());
                        if (bonds == null || bonds.length != 1) {
                            patched = false;
                            logger.log(Level.INFO, " Check biotype for hydrogen {0}.", 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.log(Level.FINE, " Created hydrogen {0}.", atomName);
                        hydrogen.setAtomType(type);
                        hydrogen.setHetero(true);
                        residue.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.log(Level.FINE, " Bonding {0} to {1} ({2} of {3}).", new Object[] { atomName, ia.getName(), numBonds, 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, 0);
                                        double[] e1 = hydrogen.getXYZ(null);
                                        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 = hydrogen.getXYZ(null);
                                        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.log(Level.INFO, " Check biotype for hydrogen {0}.", atomName);
                                        patched = false;
                                }
                                break;
                            case 3:
                                switch(numBonds) {
                                    case 2:
                                        intxyz(hydrogen, ia, 1.0, ib, 120.0, ic, 180.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.log(Level.INFO, " Check biotype for hydrogen {0}.", 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.log(Level.INFO, " Check biotype for hydrogen {0}.", 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.log(Level.INFO, " Check biotype for hydrogen {0}.", atomName);
                                        patched = false;
                                }
                                break;
                            default:
                                logger.log(Level.INFO, " Check biotype for hydrogen {0}.", atomName);
                                patched = false;
                        }
                        if (!patched) {
                            break;
                        } else {
                            buildBond(ia, hydrogen, forceField, bondList);
                        }
                    }
                }
                if (!patched) {
                    logger.log(Level.SEVERE, format(" Could not patch %s.", residueName));
                } else {
                    logger.log(Level.INFO, " Patch for {0} succeeded.", residueName);
                    residueAtoms = residue.getAtomList();
                    // Assign atom types for side-chain atoms.
                    double charge = 0.0;
                    for (Atom atom : residueAtoms) {
                        logger.info(atom.toString() + " -> " + atom.getAtomType().toString());
                    }
                }
            } else {
                switch(position) {
                    case FIRST_RESIDUE:
                        buildHydrogen(residue, "HA2", CA, 1.10e0, N, 109.5e0, C, 109.5e0, 1, 355, forceField, bondList);
                        break;
                    case LAST_RESIDUE:
                        buildHydrogen(residue, "HA2", CA, 1.10e0, N, 109.5e0, C, 109.5e0, 1, 506, forceField, bondList);
                        break;
                    default:
                        buildHydrogen(residue, "HA2", CA, 1.10e0, N, 109.5e0, C, 109.5e0, 1, 6, forceField, bondList);
                }
            }
            break;
    }
}
Also used : HashMap(java.util.HashMap) BondedUtils.buildHydrogenAtom(ffx.potential.bonded.BondedUtils.buildHydrogenAtom) AtomType(ffx.potential.parameters.AtomType) BondedUtils.findAtomType(ffx.potential.bonded.BondedUtils.findAtomType) BondedUtils.buildBond(ffx.potential.bonded.BondedUtils.buildBond)

Example 3 with AtomType

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

the class AminoAcidUtils method assignAminoAcidAtomTypes.

public static void assignAminoAcidAtomTypes(Residue residue, Residue previousResidue, Residue nextResidue, ForceField forceField, ArrayList<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
    String residueName = residue.getName().toUpperCase();
    int j = 1;
    ResiduePosition position = MIDDLE_RESIDUE;
    if (previousResidue == null) {
        j = 0;
        position = FIRST_RESIDUE;
    } else if (nextResidue == null) {
        j = 2;
        position = LAST_RESIDUE;
        /**
         * If the last residue only contains a nitrogen turn it into an NH2
         * group.
         */
        Atom N = (Atom) residue.getAtomNode("N");
        if (residue.getAtomNodeList().size() == 1 && N != null) {
            residueName = "NH2".intern();
            residue.setName(residueName);
        }
    }
    AminoAcid3 aminoAcid = getAminoAcid(residueName);
    int aminoAcidNumber = getAminoAcidNumber(residueName);
    /**
     * Non-standard Amino Acid; use ALA backbone types.
     */
    boolean nonStandard = false;
    if (aminoAcid == AminoAcid3.UNK) {
        aminoAcidNumber = getAminoAcidNumber("ALA");
        nonStandard = true;
    }
    /**
     * Only the last residue in a chain should have an OXT/OT2 atom.
     */
    if (nextResidue != null) {
        removeOXT_OT2(residue);
    }
    /**
     * Only the first nitrogen should have H1, H2 and H3 atoms, unless it's
     * an NME cap.
     */
    if (previousResidue != null) {
        removeH1_H2_H3(aminoAcid, residue);
    }
    /**
     * Check for missing heavy atoms. This check ignores special terminating
     * groups like FOR, NH2, etc.
     */
    if (!nonStandard) {
        try {
            checkForMissingHeavyAtoms(aminoAcidNumber, aminoAcid, position, residue);
        } catch (BondedUtils.MissingHeavyAtomException e) {
            logger.log(Level.INFO, " {0} could not be parsed.", residue.toString());
            logger.warning("MissingHeavyAtomException incoming from 194.");
            throw e;
        }
    }
    Atom pC = null;
    Atom pCA = null;
    if (previousResidue != null) {
        pC = (Atom) previousResidue.getAtomNode("C");
        pCA = (Atom) previousResidue.getAtomNode("CA");
    }
    /**
     * Backbone heavy atoms.
     */
    Atom N = (Atom) residue.getAtomNode("N");
    if (N != null) {
        N.setAtomType(BondedUtils.findAtomType(AA_N[j][aminoAcidNumber], forceField));
        if (position != FIRST_RESIDUE) {
            buildBond(pC, N, forceField, bondList);
        }
    }
    Atom CA = null;
    Atom C = null;
    Atom O = null;
    if (!(position == LAST_RESIDUE && aminoAcid == AminoAcid3.NH2)) {
        if (aminoAcid == AminoAcid3.ACE || aminoAcid == AminoAcid3.NME) {
            CA = buildHeavy(residue, "CH3", N, AA_CA[j][aminoAcidNumber], forceField, bondList);
        } else {
            CA = buildHeavy(residue, "CA", N, AA_CA[j][aminoAcidNumber], forceField, bondList);
        }
        if (!(position == LAST_RESIDUE && aminoAcid == AminoAcid3.NME)) {
            C = buildHeavy(residue, "C", CA, AA_C[j][aminoAcidNumber], forceField, bondList);
            O = (Atom) residue.getAtomNode("O");
            if (O == null) {
                O = (Atom) residue.getAtomNode("OT1");
            }
            AtomType atomType = findAtomType(AA_O[j][aminoAcidNumber], forceField);
            if (O == null) {
                MissingHeavyAtomException missingHeavyAtom = new MissingHeavyAtomException("O", atomType, C);
                logger.warning(" MissingHeavyAtomException incoming from 234.");
                throw missingHeavyAtom;
            }
            O.setAtomType(atomType);
            buildBond(C, O, forceField, bondList);
        }
    }
    /**
     * Nitrogen hydrogen atoms.
     */
    AtomType atomType = findAtomType(AA_HN[j][aminoAcidNumber], forceField);
    switch(position) {
        case FIRST_RESIDUE:
            switch(aminoAcid) {
                case PRO:
                    buildHydrogenAtom(residue, "H2", N, 1.02, CA, 109.5, C, 0.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H3", N, 1.02, CA, 109.5, C, -120.0, 0, atomType, forceField, bondList);
                    break;
                case PCA:
                    buildHydrogenAtom(residue, "H", N, 1.02, CA, 109.5, C, -60.0, 0, atomType, forceField, bondList);
                    break;
                case ACE:
                    break;
                default:
                    buildHydrogenAtom(residue, "H1", N, 1.02, CA, 109.5, C, 180.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H2", N, 1.02, CA, 109.5, C, 60.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H3", N, 1.02, CA, 109.5, C, -60.0, 0, atomType, forceField, bondList);
            }
            break;
        case LAST_RESIDUE:
            switch(aminoAcid) {
                case NH2:
                    buildHydrogenAtom(residue, "H1", N, 1.02, pC, 119.0, pCA, 0.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H2", N, 1.02, pC, 119.0, pCA, 180.0, 0, atomType, forceField, bondList);
                    break;
                case NME:
                    buildHydrogenAtom(residue, "H", N, 1.02, pC, 118.0, CA, 121.0, 1, atomType, forceField, bondList);
                    break;
                default:
                    buildHydrogenAtom(residue, "H", N, 1.02, pC, 119.0, CA, 119.0, 1, atomType, forceField, bondList);
            }
            break;
        default:
            // Mid-chain nitrogen hydrogen.
            buildHydrogenAtom(residue, "H", N, 1.02, pC, 119.0, CA, 119.0, 1, atomType, forceField, bondList);
    }
    /**
     * C-alpha hydrogen atoms.
     */
    String haName = "HA";
    if (aminoAcid == AminoAcid3.GLY) {
        haName = "HA2";
    }
    atomType = findAtomType(AA_HA[j][aminoAcidNumber], forceField);
    switch(position) {
        case FIRST_RESIDUE:
            switch(aminoAcid) {
                case FOR:
                    buildHydrogenAtom(residue, "H", C, 1.12, O, 0.0, null, 0.0, 0, atomType, forceField, bondList);
                    break;
                case ACE:
                    buildHydrogenAtom(residue, "H1", CA, 1.10, C, 109.5, O, 180.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H2", CA, 1.10, C, 109.5, O, 60.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H3", CA, 1.10, C, 109.5, O, -60.0, 0, atomType, forceField, bondList);
                    break;
                default:
                    buildHydrogenAtom(residue, haName, CA, 1.10, N, 109.5, C, 109.5, -1, atomType, forceField, bondList);
                    break;
            }
            break;
        case LAST_RESIDUE:
            switch(aminoAcid) {
                case NME:
                    buildHydrogenAtom(residue, "H1", CA, 1.10, N, 109.5, pC, 180.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H2", CA, 1.10, N, 109.5, pC, 60.0, 0, atomType, forceField, bondList);
                    buildHydrogenAtom(residue, "H3", CA, 1.10, N, 109.5, pC, -60.0, 0, atomType, forceField, bondList);
                    break;
                default:
                    buildHydrogenAtom(residue, haName, CA, 1.10, N, 109.5, C, 109.5, -1, atomType, forceField, bondList);
            }
            break;
        default:
            buildHydrogenAtom(residue, haName, CA, 1.10, N, 109.5, C, 109.0, -1, atomType, forceField, bondList);
    }
    /**
     * Build the amino acid side chain.
     */
    assignAminoAcidSideChain(position, aminoAcid, residue, CA, N, C, forceField, bondList);
    /**
     * Build the terminal oxygen if the residue is not NH2 or NME.
     */
    if (position == LAST_RESIDUE && !(aminoAcid == AminoAcid3.NH2 || aminoAcid == AminoAcid3.NME)) {
        atomType = findAtomType(AA_O[2][aminoAcidNumber], forceField);
        Atom OXT = (Atom) residue.getAtomNode("OXT");
        if (OXT == null) {
            OXT = (Atom) residue.getAtomNode("OT2");
            if (OXT != null) {
                OXT.setName("OXT");
            }
        }
        if (OXT == null) {
            String resName = C.getResidueName();
            int resSeq = C.getResidueNumber();
            Character chainID = C.getChainID();
            Character altLoc = C.getAltLoc();
            String segID = C.getSegID();
            double occupancy = C.getOccupancy();
            double tempFactor = C.getTempFactor();
            OXT = new Atom(0, "OXT", altLoc, new double[3], resName, resSeq, chainID, occupancy, tempFactor, segID);
            OXT.setAtomType(atomType);
            residue.addMSNode(OXT);
            intxyz(OXT, C, 1.25, CA, 117.0, O, 126.0, 1);
        } else {
            OXT.setAtomType(atomType);
        }
        buildBond(C, OXT, forceField, bondList);
    }
    /**
     * Do some checks on the current residue to make sure all atoms have
     * been assigned an atom type.
     */
    List<Atom> resAtoms = residue.getAtomList();
    for (Atom atom : resAtoms) {
        atomType = atom.getAtomType();
        if (atomType == null) {
            /**
             * Sometimes, with deuterons, a proton has been constructed in
             * its place, so we have a "dummy" deuteron still hanging
             * around.
             */
            String protonEq = atom.getName().replaceFirst("D", "H");
            Atom correspH = (Atom) residue.getAtomNode(protonEq);
            if (correspH == null || correspH.getAtomType() == null) {
                MissingAtomTypeException missingAtomTypeException = new MissingAtomTypeException(residue, atom);
                logger.warning("MissingAtomTypeException incoming from 393.");
                throw missingAtomTypeException;
            } else {
                correspH.setName(atom.getName());
                atom.removeFromParent();
                atom = correspH;
                atomType = atom.getAtomType();
            }
        }
        int numberOfBonds = atom.getNumBonds();
        if (numberOfBonds != atomType.valence) {
            if (atom == C && numberOfBonds == atomType.valence - 1 && position != LAST_RESIDUE) {
                continue;
            }
            logger.warning(format(" An atom for residue %s has the wrong number of bonds:\n %s", residueName, atom.toString()));
            logger.warning(format(" Expected: %d Actual: %d.", atomType.valence, numberOfBonds));
        }
    }
}
Also used : ResiduePosition(ffx.potential.bonded.Residue.ResiduePosition) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) BondedUtils.buildHydrogenAtom(ffx.potential.bonded.BondedUtils.buildHydrogenAtom) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) AtomType(ffx.potential.parameters.AtomType) BondedUtils.findAtomType(ffx.potential.bonded.BondedUtils.findAtomType)

Example 4 with AtomType

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

the class BondedUtils method buildHeavy.

public static Atom buildHeavy(MSGroup residue, String atomName, Atom bondedTo, int key, ForceField forceField, ArrayList<Bond> bondList) throws MissingHeavyAtomException {
    Atom atom = (Atom) residue.getAtomNode(atomName);
    AtomType atomType = findAtomType(key, forceField);
    if (atom == null) {
        MissingHeavyAtomException missingHeavyAtom = new MissingHeavyAtomException(atomName, atomType, bondedTo);
        throw missingHeavyAtom;
    }
    atom.setAtomType(atomType);
    if (bondedTo != null) {
        buildBond(atom, bondedTo, forceField, bondList);
    }
    return atom;
}
Also used : AtomType(ffx.potential.parameters.AtomType)

Example 5 with AtomType

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

the class RosenbluthChiAllMove method torsionSampler.

private void torsionSampler(List<Torsion> allTors) {
    // Collects data for plot of uTors vs. theta.
    // This is for finding the appropriate offset to add to each uTors such that the minimum is zero.
    // double origChi[] = RotamerLibrary.measureRotamer(target, false);
    updateAll();
    double[] origChi = measureLysine(target, false);
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Torsion Sampling for %s\n", target.getName()));
    sb.append(String.format("   origChi:"));
    for (int k = 0; k < origChi.length; k++) {
        sb.append(String.format(" %6.2f", origChi[k]));
    }
    sb.append("\n");
    for (int k = 0; k < origChi.length; k++) {
        Torsion tors = allTors.get(k);
        AtomType type1 = tors.getAtomArray()[0].getAtomType();
        AtomType type2 = tors.getAtomArray()[1].getAtomType();
        AtomType type3 = tors.getAtomArray()[2].getAtomType();
        AtomType type4 = tors.getAtomArray()[3].getAtomType();
        sb.append(String.format("   %d:    \"(%3d %3d %3s)  (%3d %3d %3s)  (%3d %3d %3s)  (%3d %3d %3s)\"\n", k, type1.type, type1.atomClass, type1.name, type2.type, type2.atomClass, type2.name, type3.type, type3.atomClass, type3.name, type4.type, type4.atomClass, type4.name));
    }
    logger.info(sb.toString());
    sb = new StringBuilder();
    String type = "combinatorial";
    if (type.equals("combinatorial")) {
        sb.append(String.format(" Resi chi0 chi1 chi2 chi3 List<uTors> uTorsSum uTotal\n"));
        logger.info(sb.toString());
        sb = new StringBuilder();
        boolean doChi0 = false, doChi1 = false;
        boolean doChi2 = true, doChi3 = true;
        double increment = 1.0;
        for (double chi0 = -180.0; chi0 < +180.0; chi0 += increment) {
            for (double chi1 = -180.0; chi1 <= +180.0; chi1 += increment) {
                for (double chi2 = -180.0; chi2 <= +180.0; chi2 += increment) {
                    for (double chi3 = -180.0; chi3 <= +180.0; chi3 += increment) {
                        sb.append(String.format(" %3s", target.getName()));
                        double[] newChi = new double[4];
                        if (doChi0) {
                            newChi[0] = chi0;
                        } else {
                            newChi[0] = origChi[0];
                        }
                        if (doChi1) {
                            newChi[1] = chi1;
                        } else {
                            newChi[1] = origChi[1];
                        }
                        if (doChi2) {
                            newChi[2] = chi2;
                        } else {
                            newChi[2] = origChi[2];
                        }
                        if (doChi3) {
                            newChi[3] = chi3;
                        } else {
                            newChi[3] = origChi[3];
                        }
                        Torsion chi0Tors = allTors.get(0);
                        Torsion chi1Tors = allTors.get(1);
                        Torsion chi2Tors = allTors.get(2);
                        Torsion chi3Tors = allTors.get(3);
                        Rotamer newState = createRotamer(target, newChi);
                        RotamerLibrary.applyRotamer(target, newState);
                        for (int wut = 0; wut < origChi.length; wut++) {
                            sb.append(String.format(" %3.0f", newChi[wut]));
                        }
                        writeSnapshot(String.format("%.0f", chi1), false);
                        double uTorsSum = 0.0;
                        for (int k = 0; k < allTors.size(); k++) {
                            double uTors = allTors.get(k).energy(false);
                            sb.append(String.format(" %5.2f", uTors));
                            uTorsSum += uTors;
                        }
                        sb.append(String.format(" %5.2f", uTorsSum));
                        double totalE = 1000.0;
                        try {
                            totalE = totalEnergy();
                        } catch (Exception ex) {
                        }
                        if (totalE > 1000.0) {
                            totalE = 1000.0;
                        }
                        sb.append(String.format(" %10.4f", totalE));
                        sb.append(String.format("\n"));
                        logger.info(sb.toString());
                        sb = new StringBuilder();
                        if (!doChi3) {
                            break;
                        }
                    }
                    if (!doChi2) {
                        break;
                    }
                }
                if (!doChi1) {
                    break;
                }
            }
            if (!doChi0) {
                break;
            }
        }
    } else {
        sb.append(String.format(" Resi Theta ( chi ) List<uTors,uTotal> \n"));
        logger.info(sb.toString());
        sb = new StringBuilder();
        for (double testTheta = -180.0; testTheta <= +180.0; testTheta += 0.01) {
            sb.append(String.format(" %3s %6.2f", target.getName(), testTheta));
            for (int k = 0; k < origChi.length; k++) {
                double[] newChi = new double[origChi.length];
                System.arraycopy(origChi, 0, newChi, 0, origChi.length);
                Torsion tors = allTors.get(k);
                newChi[k] = testTheta;
                Rotamer newState = createRotamer(target, newChi);
                RotamerLibrary.applyRotamer(target, newState);
                sb.append(" (");
                for (int wut = 0; wut < origChi.length; wut++) {
                    sb.append(String.format(" %6.2f", newChi[wut]));
                }
                sb.append(" )");
                writeSnapshot(String.format("%.0f", testTheta), false);
                double uTors = allTors.get(k).energy(false);
                double totalE = 0.0;
                try {
                    totalE = totalEnergy();
                } catch (Exception ex) {
                }
                sb.append(String.format(" %9.6f %9.6f", uTors, totalE));
            }
            sb.append(String.format("\n"));
            logger.info(sb.toString());
            sb = new StringBuilder();
        }
    }
    if (System.getProperty("cbmc-torsionSampler") != null) {
        try {
            File output = new File(System.getProperty("cbmc-torsionSampler"));
            BufferedWriter bw = new BufferedWriter(new FileWriter(output));
            bw.write(sb.toString());
            bw.close();
        } catch (IOException ex) {
        }
    }
    System.exit(0);
}
Also used : FileWriter(java.io.FileWriter) Rotamer(ffx.potential.bonded.Rotamer) IOException(java.io.IOException) IOException(java.io.IOException) BufferedWriter(java.io.BufferedWriter) AtomType(ffx.potential.parameters.AtomType) Torsion(ffx.potential.bonded.Torsion) File(java.io.File)

Aggregations

AtomType (ffx.potential.parameters.AtomType)19 Atom (ffx.potential.bonded.Atom)11 Bond (ffx.potential.bonded.Bond)7 MissingAtomTypeException (ffx.potential.bonded.BondedUtils.MissingAtomTypeException)6 MissingHeavyAtomException (ffx.potential.bonded.BondedUtils.MissingHeavyAtomException)6 IOException (java.io.IOException)5 AminoAcid3 (ffx.potential.bonded.ResidueEnumerations.AminoAcid3)4 NucleicAcid3 (ffx.potential.bonded.ResidueEnumerations.NucleicAcid3)4 HashMap (java.util.HashMap)4 Residue (ffx.potential.bonded.Residue)3 ResiduePosition (ffx.potential.bonded.Residue.ResiduePosition)3 ForceFieldString (ffx.potential.parameters.ForceField.ForceFieldString)3 File (java.io.File)3 Angle (ffx.potential.bonded.Angle)2 BondedUtils.buildHydrogenAtom (ffx.potential.bonded.BondedUtils.buildHydrogenAtom)2 BondedUtils.findAtomType (ffx.potential.bonded.BondedUtils.findAtomType)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