Examples with AminoAcid3 ffx.potential.bonded.ResidueEnumerations.AminoAcid3 used on opensource projects

Example 1 with AminoAcid3

use of ffx.potential.bonded.ResidueEnumerations.AminoAcid3 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));
        }
    }
}

Example 2 with AminoAcid3

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

the class PhMD method tryComboStep.

/**
 * Attempt a combination titration/rotamer MC move.
 *
 * @param targetMulti
 * @return accept/reject
 */
private boolean tryComboStep(MultiResidue targetMulti) {
    if (CAUTIOUS) {
        throw new UnsupportedOperationException();
    }
    // Record the pre-change total energy.
    double previousTotalEnergy = currentTotalEnergy();
    double previousElectrostaticEnergy = currentElectrostaticEnergy();
    // Write the pre-combo snapshot.
    writeSnapshot(true, StepType.COMBO, config.snapshots);
    String startString = targetMulti.toString();
    String startName = targetMulti.getActive().getName();
    // Choose from the list of available titrations for the active residue.
    List<Titration> avail = titrationMap.get(targetMulti.getActive());
    Titration titration = avail.get(rng.nextInt(avail.size()));
    // Perform the chosen titration.
    TitrationType titrationType = performTitration(targetMulti, titration, config.inactivateBackground);
    // Change rotamer state, but first save coordinates so we can return to them if rejected.
    Residue residue = targetMulti.getActive();
    ArrayList<Atom> atoms = residue.getAtomList();
    ResidueState origState = residue.storeState();
    double[] chi = new double[4];
    RotamerLibrary.measureAARotamer(residue, chi, false);
    AminoAcid3 aa = AminoAcid3.valueOf(residue.getName());
    Rotamer origCoordsRotamer = new Rotamer(aa, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
    // Swap to the new rotamer.
    Rotamer[] rotamers = residue.getRotamers(library);
    int rotaRand = rng.nextInt(rotamers.length);
    RotamerLibrary.applyRotamer(residue, rotamers[rotaRand]);
    // Write the post-combo snapshot.
    writeSnapshot(false, StepType.COMBO, config.snapshots);
    // Evaluate both MC criteria.
    String endName = targetMulti.getActive().getName();
    // Evaluate the titration probability of the step.
    double pKaref = titration.pKa;
    double dG_ref = titration.refEnergy;
    double temperature = thermostat.getCurrentTemperature();
    double kT = BOLTZMANN * temperature;
    double dG_elec = currentElectrostaticEnergy() - previousElectrostaticEnergy;
    if (config.zeroReferenceEnergies) {
        dG_ref = 0.0;
    }
    double prefix = Math.log(10) * kT * (pH - pKaref);
    if (titrationType == TitrationType.DEPROT) {
        prefix = -prefix;
    }
    double postfix = dG_elec - dG_ref;
    double dG_titr = prefix + postfix;
    double titrCriterion = exp(-dG_titr / kT);
    // Evaluate the rotamer probability of the step.
    double dG_rota = currentTotalEnergy() - previousTotalEnergy;
    double rotaCriterion = exp(-dG_rota / kT);
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Assessing possible MC combo step:\n"));
    sb.append(String.format("     dG_elec: %16.8f\n", dG_elec));
    sb.append(String.format("     dG_titr: %16.8f\n", dG_titr));
    sb.append(String.format("     dG_rota: %16.8f\n", dG_rota));
    sb.append(String.format("     -----\n"));
    // Automatic acceptance if both energy changes are favorable.
    if (dG_titr < 0 && dG_rota < 0 && config.mcOverride != MCOverride.REJECT) {
        sb.append(String.format("     Accepted!"));
        logger.info(sb.toString());
        numMovesAccepted++;
        propagateInactiveResidues(titratingMultis, false);
        return true;
    } else {
        // Conditionally accept based on combined probabilities.
        if (dG_titr < 0 || config.mcOverride == MCOverride.ACCEPT) {
            titrCriterion = 1.0;
        }
        if (dG_rota < 0) {
            rotaCriterion = 1.0;
        }
        if (config.mcOverride == MCOverride.REJECT) {
            titrCriterion = 0.0;
        }
        double metropolis = random();
        double comboCriterion = titrCriterion * rotaCriterion;
        sb.append(String.format("     titrCrit:   %9.4f\n", titrCriterion));
        sb.append(String.format("     rotaCrit:   %9.4f\n", rotaCriterion));
        sb.append(String.format("     criterion:  %9.4f\n", comboCriterion));
        sb.append(String.format("     rng:        %9.4f\n", metropolis));
        if (metropolis < comboCriterion) {
            sb.append(String.format("     Accepted!"));
            logger.info(sb.toString());
            numMovesAccepted++;
            propagateInactiveResidues(titratingMultis, false);
            return true;
        } else {
            // Move was denied.
            sb.append(String.format("     Denied."));
            logger.info(sb.toString());
            // Undo both pieces of the rejected move IN THE RIGHT ORDER.
            RotamerLibrary.applyRotamer(residue, origCoordsRotamer);
            performTitration(targetMulti, titration, config.inactivateBackground);
            ffe.reInit();
            molDyn.reInit();
            return false;
        }
    }
}

Example 3 with AminoAcid3

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

the class PhMD method tryRotamerStep.

/**
 * Attempt a rotamer MC move.
 *
 * @param targetMulti
 * @return accept/reject
 */
private boolean tryRotamerStep(MultiResidue targetMulti) {
    if (CAUTIOUS) {
        throw new UnsupportedOperationException();
    }
    // Record the pre-change total energy.
    double previousTotalEnergy = currentTotalEnergy();
    // Write the before-step snapshot.
    writeSnapshot(true, StepType.ROTAMER, config.snapshots);
    // Save coordinates so we can return to them if move is rejected.
    Residue residue = targetMulti.getActive();
    ArrayList<Atom> atoms = residue.getAtomList();
    ResidueState origState = residue.storeState();
    double[] chi = new double[4];
    RotamerLibrary.measureAARotamer(residue, chi, false);
    AminoAcid3 aa = AminoAcid3.valueOf(residue.getName());
    Rotamer origCoordsRotamer = new Rotamer(aa, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
    // Select a new rotamer and swap to it.
    Rotamer[] rotamers = residue.getRotamers(library);
    int rotaRand = rng.nextInt(rotamers.length);
    RotamerLibrary.applyRotamer(residue, rotamers[rotaRand]);
    // Write the post-rotamer change snapshot.
    writeSnapshot(false, StepType.ROTAMER, config.snapshots);
    // Check the MC criterion.
    double temperature = thermostat.getCurrentTemperature();
    double kT = BOLTZMANN * temperature;
    double postTotalEnergy = currentTotalEnergy();
    double dG_tot = postTotalEnergy - previousTotalEnergy;
    double criterion = exp(-dG_tot / kT);
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Assessing possible MC rotamer step:\n"));
    sb.append(String.format("     prev:   %16.8f\n", previousTotalEnergy));
    sb.append(String.format("     post:   %16.8f\n", postTotalEnergy));
    sb.append(String.format("     dG_tot: %16.8f\n", dG_tot));
    sb.append(String.format("     -----\n"));
    // Automatic acceptance if energy change is favorable.
    if (dG_tot < 0) {
        sb.append(String.format("     Accepted!"));
        logger.info(sb.toString());
        numMovesAccepted++;
        propagateInactiveResidues(titratingMultis, true);
        return true;
    } else {
        // Conditional acceptance if energy change is positive.
        double metropolis = random();
        sb.append(String.format("     criterion:  %9.4f\n", criterion));
        sb.append(String.format("     rng:        %9.4f\n", metropolis));
        if (metropolis < criterion) {
            sb.append(String.format("     Accepted!"));
            logger.info(sb.toString());
            numMovesAccepted++;
            propagateInactiveResidues(titratingMultis, true);
            return true;
        } else {
            // Move was denied.
            sb.append(String.format("     Denied."));
            logger.info(sb.toString());
            // Undo the rejected move.
            RotamerLibrary.applyRotamer(residue, origCoordsRotamer);
            return false;
        }
    }
}

Example 4 with AminoAcid3

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

the class RosenbluthChiAllMove method engage_cheap.

/**
 * So named, yet inadequate. Final stats on the ctrl vs bias algorithm for
 * chis2,3 of LYS monomer are: calc "4.23846e+07 / 22422" for the biased run
 * calc "4.21623e+06 / 10000" for the control run, where numerator = total
 * timer; demon = accepted Possible accelerations are predicated on the fact
 * that the test above was performed using unbinned ie fully-continuous
 * rotamers and sampling was done with replacement. Rectifying either of
 * these breaks balance, however, when used with MD...
 */
private boolean engage_cheap() {
    report.append(String.format(" Rosenbluth CBMC Move: %4d  (%s)\n", moveNumber, target));
    AminoAcid3 name = AminoAcid3.valueOf(target.getName());
    double[] chi = RotamerLibrary.measureRotamer(target, false);
    HashMap<Integer, BackBondedList> map = createBackBondedMap(name);
    // For each chi, create a test set from Boltzmann distr on torsion energy (Ubond).
    // Select from among this set based on Boltzmann weight of REMAINING energy (Uext).
    // The Uext partition function of each test set (wn) becomes a factor of the overall Rosenbluth (Wn).
    // ^ NOPE. Instead, Rosenbluth is going to be calculated just once for the whole combined set of chis.
    List<Torsion> allTors = new ArrayList<>();
    for (int i = 0; i < chi.length; i++) {
        Torsion tors = map.get(i).torsion;
        allTors.add(tors);
    }
    TrialSet newTrialSet = cheapTorsionSet(allTors, testSetSize, "bkn");
    // yields uExt(1) + uExt(2) + ...
    Wn = newTrialSet.sumExtBolt();
    if (Wn <= 0) {
        report.append("WARNING: Numerical instability in CMBC.");
        report.append("  Test set uExt values:  ");
        for (int i = 0; i < newTrialSet.uExt.length; i++) {
            report.append(String.format("%5.2f,  ", newTrialSet.uExt[i]));
        }
        report.append("  Discarding move.\n");
        target.revertState(origState);
        Wn = -1.0;
        Wo = 1000;
        logger.info(report.toString());
        return false;
    }
    // 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];
            proposedChis = newTrialSet.theta;
            double prob = uExtBolt / Wn * 100;
            if (printTestSets) {
                report.append(String.format("       Chose %d   %7.4f\t  %4.1f%%\n", j + 1, newTrialSet.uExt[j], 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).
    double ouDep = 0.0;
    for (Torsion tors : allTors) {
        // original-conf uDep
        ouDep += tors.energy(false);
    }
    // original-conf uExt
    double ouExt = totalEnergy() - ouDep;
    double ouExtBolt = FastMath.exp(-beta * ouExt);
    if (printTestSets) {
        report.append(String.format("       %3s %d:  %9.5g  %9.5g  %9.5g\n", "bko", 0, ouDep, ouExt, ouExtBolt));
    }
    writeSnapshot("bko", true);
    TrialSet oldTrialSet = cheapTorsionSet(allTors, 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));
    RotamerLibrary.applyRotamer(target, proposedMove);
    proposedChis = RotamerLibrary.measureRotamer(target, false);
    finalEnergy = totalEnergy();
    if (finalEnergy < CATASTROPHE_THRESHOLD) {
        report.append("\nWARNING: Multipole catastrophe in CMBC.\n");
        report.append("  Discarding move.\n");
        target.revertState(origState);
        updateAll();
        Wn = -1.0;
        Wo = 1000;
        logger.info(report.toString());
        return false;
    }
    double criterion = Math.min(1, Wn / Wo);
    rng = ThreadLocalRandom.current().nextDouble();
    report.append(String.format("    rng:    %5.2f\n", rng));
    if (rng < criterion) {
        accepted = true;
        numAccepted++;
        updateAll();
        write();
        report.append(String.format(" Accepted! %5d    NewEnergy: %.4f    Chi:", numAccepted, finalEnergy));
        for (int k = 0; k < proposedChis.length; k++) {
            report.append(String.format(" %7.2f", proposedChis[k]));
        }
        report.append(String.format("\n"));
    } else {
        accepted = false;
        target.revertState(origState);
        updateAll();
        report.append(String.format(" Denied.   %5d    OldEnergy: %.4f    Chi:", numAccepted, origEnergy));
        for (int k = 0; k < chi.length; k++) {
            report.append(String.format(" %7.2f", chi[k]));
        }
        report.append(String.format("\n"));
    }
    endTime = System.nanoTime();
    double took = (endTime - startTime) * NS_TO_MS;
    if (logTimings) {
        report.append(String.format("   Timing (ms): %.2f", took));
    }
    logger.info(report.toString());
    return accepted;
}

Example 5 with AminoAcid3

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

the class RosenbluthChiAllMove method engage_expensive.

/**
 * Follows Frenkel/Smit's derivation precisely, which for AMOEBA requires
 * SCF calls in the inner loops.
 */
private void engage_expensive() {
    report.append(String.format(" Rosenbluth CBMC Move: %4d  %s\n", moveNumber, target));
    AminoAcid3 name = AminoAcid3.valueOf(target.getName());
    double[] chi = RotamerLibrary.measureRotamer(target, false);
    HashMap<Integer, BackBondedList> map = createBackBondedMap(name);
    // For each chi, create a test set from Boltzmann distr on torsion energy (Ubond).
    // Select from among this set based on Boltzmann weight of REMAINING energy (Uext).
    // The Uext partition function of each test set (wn) becomes a factor of the overall Rosenbluth (Wn).
    // factors of Wn
    double[] wn = new double[chi.length];
    double[] finalChi = new double[chi.length];
    for (int i = 0; i < chi.length; i++) {
        Torsion tors = map.get(i).torsion;
        // TrialSet trialSet = boltzmannTorsionSet(tors, i, testSetSize, "bkn");
        TrialSet trialSet = expensiveTorsionSet(tors, i, testSetSize, "bkn");
        wn[i] = trialSet.sumExtBolt();
        if (i == 0) {
            Wn = wn[i];
        } else {
            Wn *= wn[i];
        }
        report.append(String.format(" wn,W running: %.2g %.2g", wn[i], Wn));
        logger.info(report.toString());
        // Choose a proposal move from amongst this trial set (bn).
        double rng = rand.nextDouble(wn[i]);
        double running = 0.0;
        for (int j = 0; j < trialSet.uExt.length; j++) {
            double uExtBolt = FastMath.exp(-beta * trialSet.uExt[j]);
            running += uExtBolt;
            if (rng < running) {
                // Yes, I mean i then j.
                finalChi[i] = trialSet.theta[j];
                double prob = uExtBolt / wn[i] * 100;
                report.append(String.format("       Chose %d   %7.4f\t%7.4f\t  %4.1f%%\n", j, trialSet.uExt[j], uExtBolt, prob));
                break;
            }
        }
    }
    report.append(String.format("    Wn Total:  %g\n", Wn));
    proposedMove = createRotamer(name, finalChi);
    // Reprise the above procedure for the old configuration.
    // The existing conformation forms first member of each test set (wo).
    // Overall Rosenbluth Wo is product of uExt partition functions.
    // factors of Wo
    double[] wo = new double[chi.length];
    for (int i = 0; i < chi.length; i++) {
        Torsion tors = map.get(i).torsion;
        // original-conf uDep
        double ouDep = tors.energy(false);
        // original-conf uExt
        double ouExt = totalEnergy() - ouDep;
        double ouExtBolt = FastMath.exp(-beta * ouExt);
        TrialSet trialSet = expensiveTorsionSet(tors, i, testSetSize - 1, "bko");
        wo[i] = ouExtBolt + trialSet.sumExtBolt();
        if (i == 0) {
            Wo = wo[i];
        } else {
            Wo *= wo[i];
        }
    }
    report.append(String.format("    Wo Total:  %g\n", Wo));
    report.append(String.format("    Wn/Wo:     %g\n", Wn / Wo));
    target.revertState(origState);
    updateAll();
    if (verbose) {
        logger.info(report.toString());
    }
}