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

Example 11 with Polymer

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

the class PhMD method readyup.

private void readyup() {
    // Create MultiTerminus objects to wrap termini.
    if (config.titrateTermini) {
        for (Residue res : mola.getResidueList()) {
            if (res.getPreviousResidue() == null || res.getNextResidue() == null) {
                MultiTerminus multiTerminus = new MultiTerminus(res, ff, ffe, mola);
                Polymer polymer = findResiduePolymer(res, mola);
                polymer.addMultiTerminus(res, multiTerminus);
                reInitialize(true, false);
                titratingTermini.add(multiTerminus);
                logger.info(String.format(" Titrating: %s", multiTerminus));
            }
        }
    }
    /* Create containers for titratables: MultiResidues for discrete, ExtendedVariables for continuous. */
    if (distribution == Distribution.CONTINUOUS) {
        esvSystem = new ExtendedSystem(mola);
        esvSystem.setConstantPh(pH);
        for (Residue res : chosenResidues) {
            MultiResidue multi = TitrationUtils.titratingMultiresidueFactory(mola, res);
            TitrationESV esv = new TitrationESV(esvSystem, multi);
            titratingESVs.add(esv);
            for (Residue background : multi.getInactive()) {
                inactivateResidue(background);
            }
            esvSystem.addVariable(esv);
        }
        ffe.attachExtendedSystem(esvSystem);
        logger.info(format(" Continuous pHMD readied with %d residues.", titratingESVs.size()));
    } else {
        for (Residue res : chosenResidues) {
            // Create MultiResidue objects to wrap titratables.
            MultiResidue multiRes = new MultiResidue(res, ff, ffe);
            Polymer polymer = findResiduePolymer(res, mola);
            polymer.addMultiResidue(multiRes);
            recursiveMap(res, multiRes);
            // Switch back to the original form and ready the ForceFieldEnergy.
            multiRes.setActiveResidue(res);
            reInitialize(true, false);
            titratingMultis.add(multiRes);
            logger.info(String.format(" Titrating: %s", multiRes));
        }
        logger.info(format(" Discrete MCMD readied with %d residues.", titratingMultis.size()));
    }
    switch(distribution) {
        default:
        case DISCRETE:
            molDyn.setMonteCarloListener(this, MonteCarloNotification.EACH_STEP);
            break;
        case CONTINUOUS:
            ffe.attachExtendedSystem(esvSystem);
            molDyn.attachExtendedSystem(esvSystem, 100);
            break;
    }
}

Example 12 with Polymer

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

the class RotamerOptimizationTest method testSelfEnergyElimination.

@Test
public void testSelfEnergyElimination() {
    // Load the test system.
    load();
    // Initialize Parallel Java
    try {
        String[] args = new String[0];
        Comm.init(args);
    } catch (Exception e) {
        String message = String.format(" Exception starting up the Parallel Java communication layer.");
        logger.log(Level.WARNING, message, e.toString());
        message = String.format(" Skipping rotamer optimization test.");
        logger.log(Level.WARNING, message, e.toString());
        return;
    }
    // Run the optimization.
    RotamerLibrary rLib = RotamerLibrary.getDefaultLibrary();
    rLib.setLibrary(RotamerLibrary.ProteinLibrary.Richardson);
    rLib.setUseOrigCoordsRotamer(useOriginalRotamers);
    int counter = 1;
    ArrayList<Residue> residueList = new ArrayList<Residue>();
    Polymer[] polymers = molecularAssembly.getChains();
    int nPolymers = polymers.length;
    for (int p = 0; p < nPolymers; p++) {
        Polymer polymer = polymers[p];
        ArrayList<Residue> residues = polymer.getResidues();
        for (int i = 0; i < endResID; i++) {
            Residue residue = residues.get(i);
            Rotamer[] rotamers = residue.getRotamers(rLib);
            if (rotamers != null) {
                int nrot = rotamers.length;
                if (nrot == 1) {
                    RotamerLibrary.applyRotamer(residue, rotamers[0]);
                }
                if (counter >= startResID) {
                    residueList.add(residue);
                }
            }
            counter++;
        }
    }
    RotamerOptimization rotamerOptimization = new RotamerOptimization(molecularAssembly, forceFieldEnergy, null);
    rotamerOptimization.setThreeBodyEnergy(useThreeBody);
    rotamerOptimization.setUseGoldstein(useGoldstein);
    rotamerOptimization.setPruning(pruningLevel);
    rotamerOptimization.setEnergyRestartFile(restartFile);
    rotamerOptimization.setResidues(residueList);
    double energy;
    int nRes = residueList.size();
    if (doOverallOpt) {
        rotamerOptimization.turnRotamerPairEliminationOff();
        rotamerOptimization.setTestOverallOpt(true);
        energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
        // System.out.println("The expected overall energy is: " + energy);
        assertEquals(info + " Total Energy", expectedEnergy, energy, tolerance);
    }
    if (doSelfOpt) {
        rotamerOptimization.turnRotamerPairEliminationOff();
        rotamerOptimization.setTestSelfEnergyEliminations(true);
        energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
        // System.out.println("The expected self is: " + energy);
        assertEquals(info + " Self-Energy", expectedSelfEnergy, energy, tolerance);
        // Check that optimized rotamers are equivalent to the lowest self-energy of each residue.
        int[] optimum = rotamerOptimization.getOptimumRotamers();
        // Loop over all residues
        for (int i = 0; i < nRes; i++) {
            Residue res = residueList.get(i);
            Rotamer[] rotI = res.getRotamers(rLib);
            int nRot = rotI.length;
            int rotCounter = 0;
            while (rotCounter < nRot && rotamerOptimization.checkPrunedSingles(i, rotCounter)) {
                rotCounter++;
            }
            double lowEnergy = rotamerOptimization.getSelf(i, rotCounter);
            int bestRot = rotCounter;
            for (int ri = 1; ri < nRot; ri++) {
                if (rotamerOptimization.checkPrunedSingles(i, ri)) {
                    continue;
                } else {
                    double selfEnergy = rotamerOptimization.getSelf(i, ri);
                    if (selfEnergy < lowEnergy) {
                        lowEnergy = selfEnergy;
                        bestRot = ri;
                    }
                }
            }
            assertEquals(String.format(" %s Self-Energy of residue %d", info, i), optimum[i], bestRot);
        }
    }
    if (doPairOpt) {
        rotamerOptimization.turnRotamerPairEliminationOff();
        rotamerOptimization.setTestPairEnergyEliminations(pairResidue);
        energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
        assertEquals(info + " Pair-Energy", expectedPairEnergy, energy, tolerance);
        // Check that optimized rotamers are equivalent to the lowest 2-Body energy sum for the "pairResidue".
        int[] optimum = rotamerOptimization.getOptimumRotamers();
        Residue resI = residueList.get(pairResidue);
        Rotamer[] rotI = resI.getRotamers(rLib);
        int ni = rotI.length;
        double minEnergy = Double.POSITIVE_INFINITY;
        int bestRotI = -1;
        // Loop over the pairResidue rotamers to find its lowest energy rotamer.
        for (int ri = 0; ri < ni; ri++) {
            double energyForRi = 0.0;
            if (rotamerOptimization.checkPrunedSingles(pairResidue, ri)) {
                continue;
            }
            // Loop over residue J
            for (int j = 0; j < nRes; j++) {
                if (j == pairResidue) {
                    continue;
                }
                Residue resJ = residueList.get(j);
                Rotamer[] rotJ = resJ.getRotamers(rLib);
                int nRot = rotJ.length;
                int rj = 0;
                while (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, ri, j, rj)) {
                    if (++rj >= nRot) {
                        logger.warning("RJ is too large.");
                    }
                }
                double lowEnergy = rotamerOptimization.get2Body(pairResidue, ri, j, rj);
                for (rj = 1; rj < nRot; rj++) {
                    if (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, ri, j, rj)) {
                        continue;
                    } else {
                        double pairEnergy = rotamerOptimization.get2Body(pairResidue, ri, j, rj);
                        if (pairEnergy < lowEnergy) {
                            lowEnergy = pairEnergy;
                        }
                    }
                }
                energyForRi += lowEnergy;
            }
            if (energyForRi < minEnergy) {
                minEnergy = energyForRi;
                bestRotI = ri;
            }
        }
        assertEquals(String.format(" %s Best 2-body energy sum for residue %d is with rotamer %d at %10.4f.", info, pairResidue, bestRotI, minEnergy), optimum[pairResidue], bestRotI);
        // Given the minimum energy rotamer for "pairResidue" is "bestRotI", we can check selected rotamers for all other residues.
        for (int j = 0; j < nRes; j++) {
            if (j == pairResidue) {
                continue;
            }
            Residue resJ = residueList.get(j);
            Rotamer[] rotJ = resJ.getRotamers(rLib);
            int nRotJ = rotJ.length;
            int rotCounter = 0;
            while (rotamerOptimization.checkPrunedPairs(pairResidue, bestRotI, j, rotCounter) && rotCounter < nRotJ) {
                rotCounter++;
            }
            double lowEnergy = rotamerOptimization.get2Body(pairResidue, bestRotI, j, rotCounter);
            int bestRotJ = rotCounter;
            for (int rj = 1; rj < nRotJ; rj++) {
                if (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, bestRotI, j, rj)) {
                    continue;
                } else {
                    double pairEnergy = rotamerOptimization.get2Body(pairResidue, bestRotI, j, rj);
                    if (pairEnergy < lowEnergy) {
                        lowEnergy = pairEnergy;
                        bestRotJ = rj;
                    }
                }
            }
            assertEquals(String.format(" %s Pair-Energy of residue (%d,%d) with residue %d", info, pairResidue, bestRotI, j), optimum[j], bestRotJ);
        }
    }
    // Test 3-Body Energy Eliminations.
    if (doTripleOpt) {
        rotamerOptimization.turnRotamerPairEliminationOff();
        rotamerOptimization.setTestTripleEnergyEliminations(tripleResidue1, tripleResidue2);
        try {
            energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
            assertEquals(info + " Triple-Energy", expectedTripleEnergy, energy, tolerance);
        } catch (Exception e) {
            e.fillInStackTrace();
            e.printStackTrace();
            logger.log(java.util.logging.Level.INFO, "Error in doTripleOpt", e);
        }
        // Check that optimized rotamers are equivalent to the lowest 3-body energy of each residue with the tripleResidue1 and 2.
        int[] optimum = rotamerOptimization.getOptimumRotamers();
        // fix residue 1 and gets its rotamers
        Residue resI = residueList.get(tripleResidue1);
        Rotamer[] rotI = resI.getRotamers(rLib);
        int ni = rotI.length;
        // fix residue 2 and get its rotamers
        Residue resJ = residueList.get(tripleResidue2);
        Rotamer[] rotJ = resJ.getRotamers(rLib);
        int nj = rotJ.length;
        double minEnergyIJ = Double.POSITIVE_INFINITY;
        int bestRotI = -1;
        int bestRotJ = -1;
        for (int ri = 0; ri < ni; ri++) {
            // loop through rot I
            if (rotamerOptimization.check(tripleResidue1, ri)) {
                continue;
            }
            for (int rj = 0; rj < nj; rj++) {
                // loop through rot J
                if (rotamerOptimization.checkPrunedSingles(tripleResidue2, rj) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, tripleResidue2, rj)) {
                    continue;
                }
                double currentEnergy = 0.0;
                for (int k = 0; k < nRes; k++) {
                    // loop through all other residues
                    if (k == tripleResidue1 || k == tripleResidue2) {
                        continue;
                    }
                    Residue resK = residueList.get(k);
                    Rotamer[] rotK = resK.getRotamers(rLib);
                    int nk = rotK.length;
                    int rkStart = 0;
                    while (rotamerOptimization.checkPrunedSingles(k, rkStart) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, k, rkStart) || rotamerOptimization.checkPrunedPairs(tripleResidue2, rj, k, rkStart)) {
                        if (++rkStart >= nk) {
                            logger.warning("RJ is too large.");
                        }
                    }
                    double lowEnergy = rotamerOptimization.get3Body(tripleResidue1, ri, tripleResidue2, rj, k, rkStart);
                    for (int rk = rkStart; rk < nk; rk++) {
                        if (rotamerOptimization.checkPrunedSingles(k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue2, rj, k, rk)) {
                            continue;
                        } else {
                            double tripleEnergy = rotamerOptimization.get3Body(tripleResidue1, ri, tripleResidue2, rj, k, rk);
                            if (tripleEnergy < lowEnergy) {
                                lowEnergy = tripleEnergy;
                            }
                        }
                    }
                    // adds lowest energy conformation of residue k to that of the rotamer I
                    currentEnergy += lowEnergy;
                }
                if (currentEnergy < minEnergyIJ) {
                    minEnergyIJ = currentEnergy;
                    bestRotI = ri;
                    bestRotJ = rj;
                }
            }
        }
        assertEquals(String.format(" %s Best three-body energy sum for residue %d is with rotamer %d at %10.4f.", info, tripleResidue1, bestRotI, minEnergyIJ), optimum[tripleResidue1], bestRotI);
        assertEquals(String.format(" %s Best three-body energy sum for residue %d is with rotamer %d at %10.4f.", info, tripleResidue2, bestRotJ, minEnergyIJ), optimum[tripleResidue2], bestRotJ);
        // loop over the residues to find the best rotamer per residue given bestRotI and bestRotJ
        for (int k = 0; k < nRes; k++) {
            if (k == tripleResidue1 || k == tripleResidue2) {
                continue;
            }
            Residue resK = residueList.get(k);
            Rotamer[] rotK = resK.getRotamers(rLib);
            int nk = rotK.length;
            int rotCounter = 0;
            while (rotamerOptimization.checkPrunedPairs(tripleResidue1, bestRotI, k, rotCounter) && rotamerOptimization.checkPrunedPairs(tripleResidue2, bestRotJ, k, rotCounter) && rotCounter < nk) {
                rotCounter++;
            }
            double lowEnergy = rotamerOptimization.get3Body(tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k, rotCounter);
            int bestRotK = rotCounter;
            for (int rk = 1; rk < nk; rk++) {
                if (rotamerOptimization.checkPrunedSingles(k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue1, bestRotI, k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue2, bestRotJ, k, rk)) {
                    continue;
                } else {
                    double tripleEnergy = rotamerOptimization.get3Body(tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k, rk);
                    if (tripleEnergy < lowEnergy) {
                        lowEnergy = tripleEnergy;
                        bestRotK = rk;
                    }
                }
            }
            assertEquals(String.format(" %s Triple-Energy of residue (%d,%d) and residue (%d,%d) with residue %d", info, tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k), optimum[k], bestRotK);
        }
    }
}

Example 13 with Polymer

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

the class RotamerOptimizationTest method testPairEnergyElimination.

@Test
public void testPairEnergyElimination() {
    // Load the test system.
    load();
    // Initialize Parallel Java
    try {
        String[] args = new String[0];
        Comm.init(args);
    } catch (Exception e) {
        String message = String.format(" Exception starting up the Parallel Java communication layer.");
        logger.log(Level.WARNING, message, e.toString());
        message = String.format(" Skipping rotamer optimization test.");
        logger.log(Level.WARNING, message, e.toString());
        return;
    }
    // Run the optimization.
    RotamerLibrary rLib = RotamerLibrary.getDefaultLibrary();
    rLib.setLibrary(RotamerLibrary.ProteinLibrary.Richardson);
    rLib.setUseOrigCoordsRotamer(useOriginalRotamers);
    int counter = 1;
    ArrayList<Residue> residueList = new ArrayList<Residue>();
    Polymer[] polymers = molecularAssembly.getChains();
    int nPolymers = polymers.length;
    for (int p = 0; p < nPolymers; p++) {
        Polymer polymer = polymers[p];
        ArrayList<Residue> residues = polymer.getResidues();
        for (int i = 0; i < endResID; i++) {
            Residue residue = residues.get(i);
            Rotamer[] rotamers = residue.getRotamers(rLib);
            if (rotamers != null) {
                int nrot = rotamers.length;
                if (nrot == 1) {
                    RotamerLibrary.applyRotamer(residue, rotamers[0]);
                }
                if (counter >= startResID) {
                    residueList.add(residue);
                }
            }
            counter++;
        }
    }
    RotamerOptimization rotamerOptimization = new RotamerOptimization(molecularAssembly, forceFieldEnergy, null);
    rotamerOptimization.setThreeBodyEnergy(useThreeBody);
    rotamerOptimization.setUseGoldstein(useGoldstein);
    rotamerOptimization.setPruning(pruningLevel);
    rotamerOptimization.setEnergyRestartFile(restartFile);
    rotamerOptimization.setResidues(residueList);
    double energy;
    int nRes = residueList.size();
    if (doOverallOpt) {
        rotamerOptimization.turnRotamerSingleEliminationOff();
        energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
        // System.out.println("The expected overall energy is: " + energy);
        assertEquals(info + " Total Energy", expectedEnergy, energy, tolerance);
    }
    // ToDo: Test self-energy use for rotamer 2-body eliminations.
    if (doSelfOpt) {
        rotamerOptimization.turnRotamerSingleEliminationOff();
        rotamerOptimization.setTestSelfEnergyEliminations(true);
        energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
        // System.out.println("The expected self energy is: " + energy);
        assertEquals(info + " Self-Energy", expectedSelfEnergy, energy, tolerance);
        // Check that optimized rotamers are equivalent to the lowest self-energy of each residue.
        int[] optimum = rotamerOptimization.getOptimumRotamers();
        // Loop over all residues
        for (int i = 0; i < nRes; i++) {
            Residue res = residueList.get(i);
            Rotamer[] rotI = res.getRotamers(rLib);
            int nRot = rotI.length;
            int rotCounter = 0;
            while (rotCounter < nRot && rotamerOptimization.checkPrunedSingles(i, rotCounter)) {
                rotCounter++;
            }
            double lowEnergy = rotamerOptimization.getSelf(i, rotCounter);
            int bestRot = rotCounter;
            for (int ri = 1; ri < nRot; ri++) {
                if (rotamerOptimization.checkPrunedSingles(i, ri)) {
                    continue;
                } else {
                    double selfEnergy = rotamerOptimization.getSelf(i, ri);
                    if (selfEnergy < lowEnergy) {
                        lowEnergy = selfEnergy;
                        bestRot = ri;
                    }
                }
            }
            assertEquals(String.format(" %s Self-Energy of residue %d", info, i), optimum[i], bestRot);
        }
    }
    // ToDo: Test 2-body energy use for rotamer pair eliminations.
    if (doPairOpt) {
        rotamerOptimization.turnRotamerSingleEliminationOff();
        rotamerOptimization.setTestPairEnergyEliminations(pairResidue);
        energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
        // System.out.println("The expected 2-body energy is: " + energy);
        assertEquals(info + " Pair-Energy", expectedPairEnergy, energy, tolerance);
        // Check that optimized rotamers are equivalent to the lowest 2-body energy sum for the "pairResidue".
        int[] optimum = rotamerOptimization.getOptimumRotamers();
        Residue resI = residueList.get(pairResidue);
        Rotamer[] rotI = resI.getRotamers(rLib);
        int ni = rotI.length;
        double minEnergy = Double.POSITIVE_INFINITY;
        int bestRotI = -1;
        // Loop over the pairResidue rotamers to find its lowest energy rotamer.
        for (int ri = 0; ri < ni; ri++) {
            double energyForRi = 0.0;
            if (rotamerOptimization.checkPrunedSingles(pairResidue, ri)) {
                continue;
            }
            // Loop over residue J
            for (int j = 0; j < nRes; j++) {
                if (j == pairResidue) {
                    continue;
                }
                Residue resJ = residueList.get(j);
                Rotamer[] rotJ = resJ.getRotamers(rLib);
                int nRot = rotJ.length;
                int rj = 0;
                while (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, ri, j, rj)) {
                    if (++rj >= nRot) {
                        logger.warning("RJ is too large.");
                    }
                }
                double lowEnergy = rotamerOptimization.get2Body(pairResidue, ri, j, rj);
                for (rj = 1; rj < nRot; rj++) {
                    if (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, ri, j, rj)) {
                        continue;
                    } else {
                        double pairEnergy = rotamerOptimization.get2Body(pairResidue, ri, j, rj);
                        if (pairEnergy < lowEnergy) {
                            lowEnergy = pairEnergy;
                        }
                    }
                }
                energyForRi += lowEnergy;
            }
            if (energyForRi < minEnergy) {
                minEnergy = energyForRi;
                bestRotI = ri;
            }
        }
        assertEquals(String.format(" %s Best 2-body energy sum for residue %d is with rotamer %d at %10.4f.", info, pairResidue, bestRotI, minEnergy), optimum[pairResidue], bestRotI);
        // Given the minimum energy rotamer for "pairResidue" is "bestRotI", we can check selected rotamers for all other residues.
        for (int j = 0; j < nRes; j++) {
            if (j == pairResidue) {
                continue;
            }
            Residue resJ = residueList.get(j);
            Rotamer[] rotJ = resJ.getRotamers(rLib);
            int nRotJ = rotJ.length;
            int rotCounter = 0;
            while (rotamerOptimization.checkPrunedPairs(pairResidue, bestRotI, j, rotCounter) && rotCounter < nRotJ) {
                rotCounter++;
            }
            double lowEnergy = rotamerOptimization.get2Body(pairResidue, bestRotI, j, rotCounter);
            int bestRotJ = rotCounter;
            for (int rj = 1; rj < nRotJ; rj++) {
                if (rotamerOptimization.checkPrunedSingles(j, rj) || rotamerOptimization.checkPrunedPairs(pairResidue, bestRotI, j, rj)) {
                    continue;
                } else {
                    double pairEnergy = rotamerOptimization.get2Body(pairResidue, bestRotI, j, rj);
                    if (pairEnergy < lowEnergy) {
                        lowEnergy = pairEnergy;
                        bestRotJ = rj;
                    }
                }
            }
            if (bestRotJ != optimum[j]) {
                // Check if 2-body energies are equal.
                if (lowEnergy == rotamerOptimization.get2Body(pairResidue, bestRotI, j, optimum[j])) {
                    logger.warning(String.format(" Identical 2-body energies for %s: resi %d-%d, resj %d, best rotamer J %d, optimum J %d, 2-body energy (both) %10.6f", info, pairResidue, bestRotI, j, bestRotJ, optimum[j], lowEnergy));
                } else {
                    assertEquals(String.format(" %s Pair-Energy of residue (%d,%d) with residue %d", info, pairResidue, bestRotI, j), optimum[j], bestRotJ);
                }
            }
        }
    }
    // ToDo: Test 3-Body use for rotamer pair eliminations.
    if (doTripleOpt) {
        rotamerOptimization.turnRotamerSingleEliminationOff();
        rotamerOptimization.setTestTripleEnergyEliminations(tripleResidue1, tripleResidue2);
        try {
            energy = rotamerOptimization.optimize(RotamerOptimization.Algorithm.ALL);
            assertEquals(info + " Triple-Energy", expectedTripleEnergy, energy, tolerance);
        } catch (Exception e) {
            e.fillInStackTrace();
            e.printStackTrace();
            logger.log(java.util.logging.Level.INFO, "Error in doTripleOpt", e);
        }
        // Check that optimized rotamers are equivalent to the lowest 3-body energy of each residue with the tripleResidue1 and 2.
        int[] optimum = rotamerOptimization.getOptimumRotamers();
        // fix residue 1 and gets its rotamers
        Residue resI = residueList.get(tripleResidue1);
        Rotamer[] rotI = resI.getRotamers(rLib);
        int ni = rotI.length;
        // fix residue 2 and get its rotamers
        Residue resJ = residueList.get(tripleResidue2);
        Rotamer[] rotJ = resJ.getRotamers(rLib);
        int nj = rotJ.length;
        double minEnergyIJ = Double.POSITIVE_INFINITY;
        int bestRotI = -1;
        int bestRotJ = -1;
        for (int ri = 0; ri < ni; ri++) {
            // loop through rot I
            if (rotamerOptimization.check(tripleResidue1, ri)) {
                continue;
            }
            for (int rj = 0; rj < nj; rj++) {
                // loop through rot J
                if (rotamerOptimization.checkPrunedSingles(tripleResidue2, rj) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, tripleResidue2, rj)) {
                    continue;
                }
                double currentEnergy = 0.0;
                for (int k = 0; k < nRes; k++) {
                    // loop through all other residues
                    if (k == tripleResidue1 || k == tripleResidue2) {
                        continue;
                    }
                    Residue resK = residueList.get(k);
                    Rotamer[] rotK = resK.getRotamers(rLib);
                    int nk = rotK.length;
                    int rkStart = 0;
                    while (rotamerOptimization.checkPrunedSingles(k, rkStart) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, k, rkStart) || rotamerOptimization.checkPrunedPairs(tripleResidue2, rj, k, rkStart)) {
                        if (++rkStart >= nk) {
                            logger.warning("RJ is too large.");
                        }
                    }
                    double lowEnergy = rotamerOptimization.get3Body(tripleResidue1, ri, tripleResidue2, rj, k, rkStart);
                    for (int rk = rkStart; rk < nk; rk++) {
                        if (rotamerOptimization.checkPrunedSingles(k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue1, ri, k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue2, rj, k, rk)) {
                            continue;
                        } else {
                            double tripleEnergy = rotamerOptimization.get3Body(tripleResidue1, ri, tripleResidue2, rj, k, rk);
                            if (tripleEnergy < lowEnergy) {
                                lowEnergy = tripleEnergy;
                            }
                        }
                    }
                    // adds lowest energy conformation of residue k to that of the rotamer I
                    currentEnergy += lowEnergy;
                }
                if (currentEnergy < minEnergyIJ) {
                    minEnergyIJ = currentEnergy;
                    bestRotI = ri;
                    bestRotJ = rj;
                }
            }
        }
        assertEquals(String.format(" %s Best three-body energy sum for residue %d is with rotamer %d at %10.4f.", info, tripleResidue1, bestRotI, minEnergyIJ), optimum[tripleResidue1], bestRotI);
        assertEquals(String.format(" %s Best three-body energy sum for residue %d is with rotamer %d at %10.4f.", info, tripleResidue2, bestRotJ, minEnergyIJ), optimum[tripleResidue2], bestRotJ);
        // loop over the residues to find the best rotamer per residue given bestRotI and bestRotJ
        for (int k = 0; k < nRes; k++) {
            if (k == tripleResidue1 || k == tripleResidue2) {
                continue;
            }
            Residue resK = residueList.get(k);
            Rotamer[] rotK = resK.getRotamers(rLib);
            int nk = rotK.length;
            int rotCounter = 0;
            while (rotamerOptimization.checkPrunedPairs(tripleResidue1, bestRotI, k, rotCounter) && rotamerOptimization.checkPrunedPairs(tripleResidue2, bestRotJ, k, rotCounter) && rotCounter < nk) {
                rotCounter++;
            }
            double lowEnergy = rotamerOptimization.get3Body(tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k, rotCounter);
            int bestRotK = rotCounter;
            for (int rk = 1; rk < nk; rk++) {
                if (rotamerOptimization.checkPrunedSingles(k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue1, bestRotI, k, rk) || rotamerOptimization.checkPrunedPairs(tripleResidue2, bestRotJ, k, rk)) {
                    continue;
                } else {
                    double tripleEnergy = rotamerOptimization.get3Body(tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k, rk);
                    if (tripleEnergy < lowEnergy) {
                        lowEnergy = tripleEnergy;
                        bestRotK = rk;
                    }
                }
            }
            assertEquals(String.format(" %s Triple-Energy of residue (%d,%d) and residue (%d,%d) with residue %d", info, tripleResidue1, bestRotI, tripleResidue2, bestRotJ, k), optimum[k], bestRotK);
        }
    }
}

Example 14 with Polymer

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

the class ExtendedSystem method populate.

public void populate(List<String> residueIDs) {
    // Locate the Residue identified by the given resid.
    Polymer[] polymers = mola.getChains();
    for (String token : residueIDs) {
        char chainID = token.charAt(0);
        int resNum = Integer.parseInt(token.substring(1));
        Residue target = null;
        for (Polymer p : polymers) {
            char pid = p.getChainID().charValue();
            if (pid == chainID) {
                for (Residue res : p.getResidues()) {
                    if (res.getResidueNumber() == resNum) {
                        target = res;
                        break;
                    }
                }
                if (target != null) {
                    break;
                }
            }
        }
        if (target == null) {
            logger.severe("Couldn't find target residue " + token);
        }
        MultiResidue titrating = TitrationUtils.titratingMultiresidueFactory(mola, target);
        TitrationESV esv = new TitrationESV(this, titrating);
        this.addVariable(esv);
    }
}

Example 15 with Polymer

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

the class TitrationUtils method chooseTitratables.

/**
 * Identify titratable residues and choose them all.
 */
public static List<Residue> chooseTitratables(MolecularAssembly searchMe) {
    List<Residue> chosen = new ArrayList<>();
    Polymer[] polymers = searchMe.getChains();
    for (int i = 0; i < polymers.length; i++) {
        ArrayList<Residue> residues = polymers[i].getResidues();
        for (int j = 0; j < residues.size(); j++) {
            Residue res = residues.get(j);
            Titration[] avail = Titration.multiLookup(res);
            if (avail != null) {
                chosen.add(residues.get(j));
            }
        }
    }
    return chosen;
}