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

Example 1 with Torsion

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

the class ForceFieldEnergyOpenMM method addTorsionForce.

private void addTorsionForce() {
    Torsion[] torsions = super.getTorsions();
    if (torsions == null || torsions.length < 1) {
        return;
    }
    int nTorsions = torsions.length;
    amoebaTorsionForce = OpenMM_PeriodicTorsionForce_create();
    for (int i = 0; i < nTorsions; i++) {
        Torsion torsion = torsions[i];
        int a1 = torsion.getAtom(0).getXyzIndex() - 1;
        int a2 = torsion.getAtom(1).getXyzIndex() - 1;
        int a3 = torsion.getAtom(2).getXyzIndex() - 1;
        int a4 = torsion.getAtom(3).getXyzIndex() - 1;
        TorsionType torsionType = torsion.torsionType;
        int nTerms = torsionType.phase.length;
        for (int j = 0; j < nTerms; j++) {
            OpenMM_PeriodicTorsionForce_addTorsion(amoebaTorsionForce, a1, a2, a3, a4, j + 1, torsionType.phase[j] * OpenMM_RadiansPerDegree, OpenMM_KJPerKcal * torsion.units * torsionType.amplitude[j]);
        }
    }
    OpenMM_System_addForce(system, amoebaTorsionForce);
    logger.log(Level.INFO, " Added Torsions ({0})", nTorsions);
}

Example 2 with Torsion

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

the class PhDiscount method crashDump.

/**
 * Attempt to print sources of catastrophic system heating.
 */
private void crashDump(RuntimeException error) {
    writeSnapshot(".meltdown-");
    ffe.energy(false, true);
    mola.getDescendants(BondedTerm.class).stream().filter(BondedTerm::isExtendedSystemMember).forEach(term -> {
        try {
            ((Bond) term).log();
        } catch (Exception ex) {
        }
        try {
            ((Angle) term).log();
        } catch (Exception ex) {
        }
        try {
            ((Torsion) term).log();
        } catch (Exception ex) {
        }
    });
    if (ffe.getVanDerWaalsEnergy() > 1000) {
        extendedAtoms = esvSystem.getExtendedAtoms();
        nAtomsExt = esvSystem.getExtendedAtoms().length;
        for (int i = 0; i < nAtomsExt; i++) {
            Atom ai = extendedAtoms[i];
            for (int j = 0; j < nAtomsExt; j++) {
                Atom aj = extendedAtoms[j];
                if (!esvSystem.isExtended(i) && !esvSystem.isExtended(j)) {
                    continue;
                }
                if (ai == aj || ai.getBond(aj) != null) {
                    continue;
                }
                double dist = FastMath.sqrt(FastMath.pow((aj.getX() - ai.getX()), 2) + FastMath.pow((aj.getY() - ai.getY()), 2) + FastMath.pow((aj.getZ() - ai.getZ()), 2));
                if (dist < 0.8 * (aj.getVDWR() + ai.getVDWR())) {
                    logger.warning(String.format("Close vdW contact for atoms: \n   %s\n   %s", aj, ai));
                }
            }
        }
    }
    throw error;
}

Example 3 with Torsion

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

Example 4 with Torsion

use of ffx.potential.bonded.Torsion 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 Torsion

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