Search in sources :

Example 16 with AminoAcid3

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

the class PhMD method recursiveMap.

/**
 * Finds Titration definitions for the given Residue and adds them to the given MultiResidue.
 * For three-state transitions, simply populate the enumeration with multiple titrations
 * from a shared state and this will include them in MultiResidue construction.
 */
private void recursiveMap(Residue member, MultiResidue multiRes) {
    // Map titrations for this member.
    Titration[] titrations = Titration.multiLookup(member);
    titrationMap.put(member, Arrays.asList(titrations));
    // For each titration, check whether it needs added as a MultiResidue option.
    for (Titration titration : titrations) {
        // Allow manual override of Histidine treatment.
        if ((titration.deprotForm == AminoAcid3.HID && config.histidineMode == HistidineMode.HIE_ONLY) || (titration.deprotForm == AminoAcid3.HIE && config.histidineMode == HistidineMode.HID_ONLY)) {
            continue;
        }
        // Find all the choices currently available to this MultiResidue.
        List<AminoAcid3> choices = new ArrayList<>();
        for (Residue choice : multiRes.getConsideredResidues()) {
            choices.add(choice.getAminoAcid3());
        }
        // If this Titration target is not a choice for the MultiResidue, then add it.
        if (!choices.contains(titration.protForm) || !(choices.contains(titration.deprotForm))) {
            String targetName = (member.getAminoAcid3() == titration.protForm) ? titration.deprotForm.toString() : titration.protForm.toString();
            int resNumber = member.getResidueNumber();
            ResidueType resType = member.getResidueType();
            Residue newChoice = new Residue(targetName, resNumber, resType);
            multiRes.addResidue(newChoice);
            titrationMap.put(newChoice, Arrays.asList(Titration.multiLookup(newChoice)));
        }
    }
}
Also used : ResidueType(ffx.potential.bonded.Residue.ResidueType) TitrationUtils.inactivateResidue(ffx.potential.extended.TitrationUtils.inactivateResidue) MultiResidue(ffx.potential.bonded.MultiResidue) Residue(ffx.potential.bonded.Residue) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) ArrayList(java.util.ArrayList) Titration(ffx.potential.extended.TitrationUtils.Titration) TitrationUtils.performTitration(ffx.potential.extended.TitrationUtils.performTitration)

Example 17 with AminoAcid3

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

the class RosenbluthChiAllMove method engage_diffs.

private void engage_diffs() {
    report.append(String.format(" Rosenbluth CBMC Move: %4d\n", moveNumber));
    report.append(String.format("    residue:   %s\n", target.toString()));
    double origFastEnergy = fastEnergy();
    double origSlowEnergy = slowEnergy();
    AminoAcid3 name = AminoAcid3.valueOf(target.getName());
    double[] chi = RotamerLibrary.measureRotamer(target, false);
    // Now we're really just wingin' it: uDep is now "fast energy" (like RESPA does),
    // uExt is the slow part, and we operate solely on diffs... so the bko_zero term -> unity.
    TrialSet newTrialSet = cheapTorsionSet_diffs(testSetSize, "bkn");
    // <-- explore difference between prod W(n) and sum wi
    Wn = newTrialSet.sumExtBolt();
    if (Wn <= 0) {
        StringBuilder sb = new StringBuilder();
        sb.append("Numerical instability in CMBC:");
        sb.append("  Test set uExt values:  ");
        for (int i = 0; i < newTrialSet.uExt.length; i++) {
            sb.append(String.format("%5.2f,  ", newTrialSet.uExt[i]));
        }
        logger.warning(sb.toString());
    }
    // Choose a proposal move from amongst this trial set (bn).
    double rng = rand.nextDouble(Wn);
    double running = 0.0;
    for (int j = 0; j < newTrialSet.uExt.length; j++) {
        double uExtBolt = FastMath.exp(-beta * newTrialSet.uExt[j]);
        running += uExtBolt;
        if (rng < running) {
            proposedMove = newTrialSet.rotamer[j];
            double prob = uExtBolt / Wn * 100;
            report.append(String.format("       Chose %d   %7.4f\t%7.4f\t  %4.1f%%\n", j + 1, newTrialSet.uExt[j], uExtBolt, prob));
            break;
        }
    }
    report.append(String.format("    Wn Total:  %g\n", Wn));
    // Reprise the above procedure for the old configuration.
    // The existing conformation forms first member of each test set (wo).
    target.revertState(origState);
    // original-conf uDep
    double ouDep = slowEnergy() - origSlowEnergy;
    // original-conf uExt
    double ouExt = fastEnergy() - origFastEnergy;
    if (Math.abs(ouDep) > tolerance || Math.abs(ouExt) > tolerance) {
        report.append(" WAIT WHAT (v2).  ouDep/ouExt remainder: " + ouDep + "," + ouExt + "\n");
    }
    double ouExtBolt = FastMath.exp(-beta * ouExt);
    report.append(String.format("       %3s %d:  %9.5g  %9.5g  %9.5g\n", "bko", 0, ouDep, ouExt, ouExtBolt));
    writeSnapshot("bko", true);
    TrialSet oldTrialSet = cheapTorsionSet_diffs(testSetSize - 1, "bko");
    Wo = ouExtBolt + oldTrialSet.sumExtBolt();
    report.append(String.format("    Wo Total:  %11.4g\n", Wo));
    report.append(String.format("    Wn/Wo:     %11.4g", Wn / Wo));
    target.revertState(origState);
    updateAll();
    if (verbose) {
        logger.info(report.toString());
    }
}
Also used : AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3)

Example 18 with AminoAcid3

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

the class PDBFilter method renameAminoAcidToPDBStandard.

/**
 * Renames the Atoms in an amino acid to PDB standard.
 *
 * @param residue Residue to fix atom names of.
 */
private static void renameAminoAcidToPDBStandard(Residue residue) {
    final Atom N = findNitrogenAtom(residue);
    AminoAcid3 aa3 = residue.getAminoAcid3();
    if (N != null) {
        N.setName("N");
        Atom CA = getAlphaCarbon(residue, N);
        CA.setName("CA");
        List<Atom> has = findBondedAtoms(CA, 1);
        switch(aa3) {
            case NME:
                // Do all renaming here then return out of the method.
                findBondedAtoms(N, 1).get(0).setName("H");
                CA.setName("CH3");
                for (int i = 1; i <= 3; i++) {
                    has.get(i - 1).setName(String.format("H%d", i));
                }
                return;
            case GLY:
                has.get(0).setName("HA2");
                has.get(1).setName("HA3");
                break;
            default:
                has.get(0).setName("HA");
                break;
        }
        CA.setName("CA");
        Atom C = null;
        Atom CB = null;
        for (Atom carb : findBondedAtoms(CA, 6)) {
            // Second check is because of serine/threonine OG bonded straight to CB.
            if (hasAttachedAtom(carb, 8) && !hasAttachedAtom(carb, 1)) {
                C = carb;
                C.setName("C");
            } else {
                CB = carb;
                CB.setName("CB");
            }
        }
        if (C == null) {
            throw new IllegalArgumentException(String.format(" Could not find carbonyl carbon for residue %s!", residue));
        }
        if (CB == null && aa3 != AminoAcid3.GLY) {
            throw new IllegalArgumentException(String.format(" Could not find beta carbon for residue %s!", residue));
        }
        List<Atom> ctermOxygens = findBondedAtoms(C, 8);
        switch(ctermOxygens.size()) {
            case 1:
                ctermOxygens.get(0).setName("O");
                break;
            case 2:
                Atom O = null;
                for (Atom oxy : ctermOxygens) {
                    if (oxy.getBonds().size() == 2) {
                        O = oxy;
                        O.setName("O");
                        findBondedAtoms(O, 1).get(0).setName("HO");
                    }
                }
                if (O == null) {
                    ctermOxygens.get(0).setName("O");
                    ctermOxygens.get(1).setName("OXT");
                }
        }
        List<Atom> amideProtons = findBondedAtoms(N, 1);
        switch(amideProtons.size()) {
            case 1:
                amideProtons.get(0).setName("H");
                break;
            default:
                // Should catch both N-termini and proline.
                for (int i = 1; i <= amideProtons.size(); i++) {
                    amideProtons.get(i - 1).setName(String.format("H%d", i));
                }
                break;
        }
        /**
         * All common atoms are now named: N, H[1-3], CA, HA[2-3], CB, C, O[XT], [HO]
         */
        renameCommonAminoAcids(residue, aa3, CA, CB);
    } else {
        switch(aa3) {
            case ACE:
                {
                    Atom O = findAtomsOfElement(residue, 8).get(0);
                    O.setName("O");
                    Atom C = findBondedAtoms(O, 6).get(0);
                    C.setName("C");
                    Atom CH3 = findBondedAtoms(C, 6).get(0);
                    CH3.setName("CH3");
                    List<Atom> hydrogens = findBondedAtoms(CH3, 1);
                    for (int i = 1; i <= 3; i++) {
                        hydrogens.get(i - 1).setName(String.format("H%d", i));
                    }
                }
                break;
            default:
                throw new IllegalArgumentException(String.format(" Could not find nitrogen atom for residue %s!", residue));
        }
    }
}
Also used : AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) ResidueEnumerations.nucleicAcidList(ffx.potential.bonded.ResidueEnumerations.nucleicAcidList) ResidueEnumerations.aminoAcidList(ffx.potential.bonded.ResidueEnumerations.aminoAcidList) List(java.util.List) ArrayList(java.util.ArrayList) Atom(ffx.potential.bonded.Atom)

Example 19 with AminoAcid3

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

the class BiojavaFilter method convert.

/**
 * {@inheritDoc}
 *
 * Parse the Biojava Structure
 *
 * @return true if the structure is successfully converted
 */
@Override
public boolean convert() {
    // First atom is #1, to match xyz file format
    int xyzIndex = 1;
    setConverted(false);
    systems.add(activeMolecularAssembly);
    if (mutate) {
        List<Character> chainIDs = new ArrayList<>();
        for (Chain chain : structure.getChains()) {
            chainIDs.add(chain.getChainID().charAt(0));
        }
        if (!chainIDs.contains(mutateChainID)) {
            if (chainIDs.size() == 1) {
                logger.warning(String.format(" Chain ID %c for " + "mutation not found: only one chain %c " + "found.", mutateChainID, chainIDs.get(0)));
                mutateChainID = chainIDs.get(0);
            } else {
                logger.warning(String.format(" Chain ID %c for " + "mutation not found: mutation will not " + "proceed.", mutateChainID));
            }
        }
    }
    /**
     * Echo the alternate location being parsed.
     */
    if (currentAltLoc == 'A') {
        logger.info(String.format(" Reading %s", structure.getName()));
    } else {
        logger.info(String.format(" Reading %s alternate location %s", structure.getName(), currentAltLoc));
    }
    org.biojava.bio.structure.Atom[] bjAtoms = StructureTools.getAllAtomArray(structure);
    int nAtoms = bjAtoms.length;
    Atom[] ffxAtoms = new Atom[nAtoms];
    /**
     * Reset the current chain and segID.
     */
    currentChainID = null;
    currentSegID = null;
    PDBCrystallographicInfo cInfo = structure.getCrystallographicInfo();
    if (cInfo.isCrystallographic()) {
        // I do not think we need to check if it already has these properties,
        // but it can be done.
        properties.addProperty("a-axis", cInfo.getA());
        properties.addProperty("b-axis", cInfo.getB());
        properties.addProperty("c-axis", cInfo.getC());
        properties.addProperty("alpha", cInfo.getAlpha());
        properties.addProperty("beta", cInfo.getBeta());
        properties.addProperty("gamma", cInfo.getGamma());
        properties.addProperty("spacegroup", SpaceGroup.pdb2ShortName(cInfo.getSpaceGroup()));
    }
    for (org.biojava.bio.structure.Atom atom : bjAtoms) {
        String name = atom.getName().toUpperCase().trim();
        double[] xyz = new double[3];
        xyz[0] = atom.getX();
        xyz[1] = atom.getY();
        xyz[2] = atom.getZ();
        char altLoc = atom.getAltLoc();
        if (!altLocs.contains(altLoc)) {
            altLocs.add(altLoc);
        }
        if (altLoc != ' ' && altLoc != 'A' && altLoc != currentAltLoc) {
            break;
        }
        if (name.contains("1H") || name.toUpperCase().contains("2H") || name.toUpperCase().contains("3H")) {
            // VERSION3_2 is presently just a placeholder for "anything non-standard".
            fileStandard = VERSION3_2;
        }
        Group group = atom.getGroup();
        ResidueNumber resnum = group.getResidueNumber();
        int resSeq = resnum.getSeqNum();
        String resName = group.getPDBName().trim().toUpperCase();
        Chain chain = group.getChain();
        char chainID = chain.getChainID().charAt(0);
        String segID = getSegID(chainID);
        boolean printAtom = false;
        if (mutate && chainID == mutateChainID && mutateResID == resSeq) {
            if (name.equals("N") || name.equals("C") || name.equals("O") || name.equals("CA")) {
                printAtom = true;
                name = mutateToResname;
            } else {
                logger.info(String.format(" Deleting atom %s of %s %d", name, resName, resSeq));
                break;
            }
        }
        Atom newAtom = new Atom(0, name, altLoc, xyz, resName, resSeq, chainID, atom.getOccupancy(), atom.getTempFactor(), segID);
        /* Biojava sets at least some capping groups, and possibly nonstandard
             amino acids to be heteroatoms. */
        boolean hetatm = true;
        for (AminoAcid3 aa3Name : AminoAcid3.values()) {
            if (aa3Name.name().equals(resName)) {
                hetatm = false;
                break;
            }
        }
        newAtom.setHetero(hetatm);
        // Look Ma, Biojava doesn't care about anisou!
        Atom returnedAtom = (Atom) activeMolecularAssembly.addMSNode(newAtom);
        if (returnedAtom != newAtom) {
            atoms.put(atom.getPDBserial(), returnedAtom);
            if (logger.isLoggable(Level.FINE)) {
                logger.fine(String.format("%s has been retained over\n%s", returnedAtom.toString(), newAtom.toString()));
            }
        } else {
            atoms.put(atom.getPDBserial(), newAtom);
            if (newAtom.getIndex() == 0) {
                newAtom.setXyzIndex(xyzIndex++);
            }
            if (printAtom) {
                logger.info(newAtom.toString());
            }
        }
    }
    List<Bond> ssBondList = new ArrayList<>();
    for (SSBond ssBond : structure.getSSBonds()) {
        Polymer c1 = activeMolecularAssembly.getChain(ssBond.getChainID1());
        Polymer c2 = activeMolecularAssembly.getChain(ssBond.getChainID2());
        int rn1;
        int rn2;
        try {
            rn1 = Integer.parseInt(ssBond.getResnum1());
            rn2 = Integer.parseInt(ssBond.getResnum1());
        } catch (NumberFormatException ex) {
            logger.warning(String.format(" Could not parse SSbond %d", ssBond.getSerNum()));
            continue;
        }
        Residue r1 = c1.getResidue(rn1);
        Residue r2 = c2.getResidue(rn2);
        List<Atom> atoms1 = r1.getAtomList();
        List<Atom> atoms2 = r2.getAtomList();
        Atom SG1 = null;
        Atom SG2 = null;
        for (Atom atom : atoms1) {
            if (atom.getName().equalsIgnoreCase("SG")) {
                SG1 = atom;
                break;
            }
        }
        for (Atom atom : atoms2) {
            if (atom.getName().equalsIgnoreCase("SG")) {
                SG2 = atom;
                break;
            }
        }
        if (SG1 == null) {
            logger.warning(String.format(" SG atom 1 of SS-bond %s is null", ssBond.toString()));
        }
        if (SG2 == null) {
            logger.warning(String.format(" SG atom 2 of SS-bond %s is null", ssBond.toString()));
        }
        if (SG1 == null || SG2 == null) {
            continue;
        }
        double d = VectorMath.dist(SG1.getXYZ(null), SG2.getXYZ(null));
        if (d < 3.0) {
            r1.setName("CYX");
            r2.setName("CYX");
            for (Atom atom : atoms1) {
                atom.setResName("CYX");
            }
            for (Atom atom : atoms2) {
                atom.setResName("CYX");
            }
            Bond bond = new Bond(SG1, SG2);
            ssBondList.add(bond);
        } else {
            String message = String.format("Ignoring [%s]\n due to distance %8.3f A.", ssBond.toString(), d);
            logger.log(Level.WARNING, message);
        }
    }
    int pdbAtoms = activeMolecularAssembly.getAtomArray().length;
    assignAtomTypes();
    StringBuilder sb = new StringBuilder(" Disulfide Bonds:");
    for (Bond bond : ssBondList) {
        Atom a1 = bond.getAtom(0);
        Atom a2 = bond.getAtom(1);
        int[] c = new int[2];
        c[0] = a1.getAtomType().atomClass;
        c[1] = a2.getAtomType().atomClass;
        String key = ffx.potential.parameters.BondType.sortKey(c);
        ffx.potential.parameters.BondType bondType = forceField.getBondType(key);
        if (bondType == null) {
            logger.severe(String.format("No BondType for key: %s\n %s\n %s", key, a1, a2));
        } else {
            bond.setBondType(bondType);
        }
        double d = VectorMath.dist(a1.getXYZ(null), a2.getXYZ(null));
        Polymer c1 = activeMolecularAssembly.getChain(a1.getSegID());
        Polymer c2 = activeMolecularAssembly.getChain(a2.getSegID());
        Residue r1 = c1.getResidue(a1.getResidueNumber());
        Residue r2 = c2.getResidue(a2.getResidueNumber());
        sb.append(String.format("\n S-S distance of %6.2f for %s and %s.", d, r1.toString(), r2.toString()));
        bondList.add(bond);
    }
    if (ssBondList.size() > 0) {
        logger.info(sb.toString());
    }
    /**
     * Finally, re-number the atoms if missing atoms were created.
     */
    if (pdbAtoms != activeMolecularAssembly.getAtomArray().length) {
        numberAtoms();
    }
    return true;
}
Also used : Chain(org.biojava.bio.structure.Chain) Group(org.biojava.bio.structure.Group) MSGroup(ffx.potential.bonded.MSGroup) SpaceGroup(ffx.crystal.SpaceGroup) ArrayList(java.util.ArrayList) PDBCrystallographicInfo(org.biojava.bio.structure.PDBCrystallographicInfo) SSBond(org.biojava.bio.structure.SSBond) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) Polymer(ffx.potential.bonded.Polymer) Atom(ffx.potential.bonded.Atom) Residue(ffx.potential.bonded.Residue) ResidueNumber(org.biojava.bio.structure.ResidueNumber) Bond(ffx.potential.bonded.Bond) SSBond(org.biojava.bio.structure.SSBond)

Aggregations

AminoAcid3 (ffx.potential.bonded.ResidueEnumerations.AminoAcid3)19 Atom (ffx.potential.bonded.Atom)9 ArrayList (java.util.ArrayList)9 Residue (ffx.potential.bonded.Residue)7 MissingAtomTypeException (ffx.potential.bonded.BondedUtils.MissingAtomTypeException)4 MissingHeavyAtomException (ffx.potential.bonded.BondedUtils.MissingHeavyAtomException)4 MultiResidue (ffx.potential.bonded.MultiResidue)4 NucleicAcid3 (ffx.potential.bonded.ResidueEnumerations.NucleicAcid3)4 AtomType (ffx.potential.parameters.AtomType)4 Bond (ffx.potential.bonded.Bond)3 Polymer (ffx.potential.bonded.Polymer)3 ResidueEnumerations.aminoAcidList (ffx.potential.bonded.ResidueEnumerations.aminoAcidList)3 ResidueEnumerations.nucleicAcidList (ffx.potential.bonded.ResidueEnumerations.nucleicAcidList)3 Torsion (ffx.potential.bonded.Torsion)3 TitrationUtils.inactivateResidue (ffx.potential.extended.TitrationUtils.inactivateResidue)3 List (java.util.List)3 MSNode (ffx.potential.bonded.MSNode)2 Molecule (ffx.potential.bonded.Molecule)2 ResiduePosition (ffx.potential.bonded.Residue.ResiduePosition)2 ResidueState (ffx.potential.bonded.ResidueState)2