Examples with PDBFilter ffx.potential.parsers.PDBFilter used on opensource projects

Example 16 with PDBFilter

use of ffx.potential.parsers.PDBFilter in project ffx by mjschnie.

the class RotamerOptimization method globalOptimization.

/**
 * The main driver for optimizing a block of residues using DEE.
 *
 * @param residueList Residues to optimize.
 * @return Final energy.
 */
private double globalOptimization(List<Residue> residueList) {
    int currentEnsemble = Integer.MAX_VALUE;
    Residue[] residues = residueList.toArray(new Residue[residueList.size()]);
    int nResidues = residues.length;
    int[] currentRotamers = new int[nResidues];
    int iterations = 0;
    boolean finalTry = false;
    int bestEnsembleTargetDiffThusFar = Integer.MAX_VALUE;
    double bestBufferThusFar = ensembleBuffer;
    double startingBuffer = ensembleBuffer;
    optimum = new int[nResidues];
    if (ensembleEnergy > 0.0) {
        ensembleBuffer = ensembleEnergy;
        applyEliminationCriteria(residues, true, true);
        if (x == null) {
            Atom[] atoms = molecularAssembly.getAtomArray();
            int nAtoms = atoms.length;
            x = new double[nAtoms * 3];
        }
        /**
         * Compute the number of permutations without eliminating dead-ends
         * and compute the number of permutations using singleton
         * elimination.
         */
        double permutations = 1;
        double singletonPermutations = 1;
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int nr = rotamers.length;
            if (nr > 1) {
                int nrot = 0;
                for (int ri = 0; ri < nr; ri++) {
                    if (!eliminatedSingles[i][ri]) {
                        nrot++;
                    }
                }
                permutations *= rotamers.length;
                if (nrot > 1) {
                    singletonPermutations *= nrot;
                }
            }
        }
        dryRun(residues, 0, currentRotamers);
        double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
        double afterPairElim = singletonPermutations - pairTotalElimination;
        if (evaluatedPermutations == 0) {
            logger.severe(" No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
        }
        if (master && printFiles && ensembleFile == null) {
            File file = molecularAssembly.getFile();
            String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
            ensembleFile = new File(filename + ".ens");
            if (ensembleFile.exists()) {
                for (int i = 2; i < 1000; i++) {
                    ensembleFile = new File(filename + ".ens_" + i);
                    if (!ensembleFile.exists()) {
                        break;
                    }
                }
                if (ensembleFile.exists()) {
                    logger.warning(format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
                }
            }
            ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null, null);
            logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
        }
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Condition", "|", "Number of Permutations Left", "|", "Removed", "|"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "No Eliminations", "|", permutations, "|", "", "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single Eliminations", "|", singletonPermutations, "|", permutations - singletonPermutations, "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Pair Eliminations", "|", afterPairElim, "|", pairTotalElimination, "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single and Pair Eliminations", "|", (double) evaluatedPermutations, "|", pairTotalElimination + (permutations - singletonPermutations), "|"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------\n"));
        logIfMaster(format(" Energy of permutations:"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Permutation", "|", "Energy", "|", "Lowest Possible Energy", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        double e;
        if (useMonteCarlo()) {
            firstValidPerm(residues, 0, currentRotamers);
            System.arraycopy(currentRotamers, 0, optimum, 0, nResidues);
            rotamerOptimizationMC(residues, optimum, currentRotamers, nMCsteps, false, mcUseAll);
            logIfMaster(" Ensembles not currently compatible with Monte Carlo search");
        /**
         * Not currently compatible with ensembles.
         */
        } else {
            double[] permutationEnergies = new double[evaluatedPermutations];
            ensembleStates = new ArrayList<>();
            e = rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE, optimum, permutationEnergies);
            int[][] acceptedPermutations = new int[evaluatedPermutations][];
            for (int i = 0; i < acceptedPermutations.length; i++) {
                acceptedPermutations[i] = null;
            }
            logIfMaster(format("\n Checking permutations for distance < %5.3f kcal/mol from GMEC energy %10.8f kcal/mol", ensembleEnergy, e));
            dryRunForEnsemble(residues, 0, currentRotamers, e, permutationEnergies, acceptedPermutations);
            int numAcceptedPermutations = 0;
            for (int i = 0; i < acceptedPermutations.length; i++) {
                if (acceptedPermutations[i] != null) {
                    ++numAcceptedPermutations;
                    logIfMaster(format(" Accepting permutation %d at %8.6f < %8.6f", i, permutationEnergies[i] - e, ensembleEnergy));
                    for (int j = 0; j < nResidues; j++) {
                        Residue residuej = residues[j];
                        Rotamer[] rotamersj = residuej.getRotamers(library);
                        RotamerLibrary.applyRotamer(residuej, rotamersj[acceptedPermutations[i][j]]);
                    }
                    ResidueState[] states = ResidueState.storeAllCoordinates(residues);
                    ensembleStates.add(new ObjectPair<>(states, permutationEnergies[i]));
                    if (printFiles && master) {
                        try {
                            FileWriter fw = new FileWriter(ensembleFile, true);
                            BufferedWriter bw = new BufferedWriter(fw);
                            bw.write(format("MODEL        %d", numAcceptedPermutations));
                            for (int j = 0; j < 75; j++) {
                                bw.write(" ");
                            }
                            bw.newLine();
                            bw.flush();
                            ensembleFilter.writeFile(ensembleFile, true);
                            bw.write(format("ENDMDL"));
                            for (int j = 0; j < 64; j++) {
                                bw.write(" ");
                            }
                            bw.newLine();
                            bw.close();
                        } catch (IOException ex) {
                            logger.warning(format(" Exception writing to file: %s", ensembleFile.getName()));
                        }
                    }
                }
            }
            logIfMaster(format(" Number of permutations within %5.3f kcal/mol of GMEC energy: %6.4e", ensembleEnergy, (double) numAcceptedPermutations));
            ensembleStates.sort(null);
        }
        logIfMaster(format(" Final rotamers:"));
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
        logIfMaster(format("%14s %3s %10s %3s %9s %3s %9s %3s %9s %3s", "Residue", "|", "Chi 1", "|", "Chi 2", "|", "Chi 3", "|", "Chi 4", "|"));
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int ri = optimum[i];
            Rotamer rotamer = rotamers[ri];
            logIfMaster(format(" %c (%7s,2d) | %s", residue.getChainID(), residue, ri, rotamer.toAngleString()));
            RotamerLibrary.applyRotamer(residue, rotamer);
        }
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------\n"));
        double sumSelfEnergy = 0;
        double sumPairEnergy = 0;
        double sumTrimerEnergy = 0;
        for (int i = 0; i < nResidues; i++) {
            int ri = optimum[i];
            sumSelfEnergy += getSelf(i, ri);
            logIfMaster(format(" Final self Energy (%8s,%2d): %12.4f", residues[i].toFormattedString(false, true), ri, getSelf(i, ri)));
        }
        for (int i = 0; i < nResidues - 1; i++) {
            int ri = optimum[i];
            for (int j = i + 1; j < nResidues; j++) {
                int rj = optimum[j];
                sumPairEnergy += get2Body(i, ri, j, rj);
                if (get2Body(i, ri, j, rj) > 10.0) {
                    logIfMaster(format(" Large Final Pair Energy (%8s,%2d) (%8s,%2d): %12.4f", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, get2Body(i, ri, j, rj)));
                }
            }
        }
        try {
            e = currentEnergy(residueList);
        } catch (ArithmeticException ex) {
            e = Double.NaN;
            logger.severe(String.format(" Exception %s in calculating current energy at the end of triples", ex.toString()));
        }
        logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Type", "|", "Energy", "|", "Lowest Possible Energy", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Self:", "|", sumSelfEnergy, "|", "", "|"));
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Pair:", "|", sumPairEnergy, "|", "", "|"));
        double approximateEnergy = backboneEnergy + sumSelfEnergy + sumPairEnergy;
        if (threeBodyTerm) {
            for (int i = 0; i < nResidues - 2; i++) {
                int ri = optimum[i];
                for (int j = i + 1; j < nResidues - 1; j++) {
                    int rj = optimum[j];
                    for (int k = j + 1; k < nResidues; k++) {
                        int rk = optimum[k];
                        try {
                            sumTrimerEnergy += get3Body(i, ri, j, rj, k, rk);
                        } catch (Exception ex) {
                            logger.warning(ex.toString());
                        }
                    }
                }
            }
            approximateEnergy += sumTrimerEnergy;
            double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - backboneEnergy;
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Trimer:", "|", sumTrimerEnergy, "|", "", "|"));
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
        } else {
            double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - backboneEnergy;
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
        }
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Approximate:", "|", approximateEnergy, "|", "", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------\n"));
        return e;
    }
    /**
     * Permutations used only to set maximum bound on ensembleNumber, thus
     * it is safe here to put that value in a 32-bit int.
     */
    int nPerms = 1;
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int nr = rotamers.length;
        if (nr > 1) {
            nPerms *= rotamers.length;
        }
        if (nPerms > ensembleNumber) {
            break;
        }
    }
    if (nPerms < ensembleNumber) {
        logger.info(format(" Requested an ensemble of %d, but only %d permutations exist; returning full ensemble", ensembleNumber, nPerms));
        ensembleNumber = nPerms;
    }
    while (currentEnsemble != ensembleNumber) {
        if (monteCarlo) {
            logIfMaster(" Ensemble search not currently compatible with Monte Carlo");
            ensembleNumber = 1;
        }
        if (iterations == 0) {
            applyEliminationCriteria(residues, true, true);
        } else {
            applyEliminationCriteria(residues, false, false);
        }
        if (x == null) {
            Atom[] atoms = molecularAssembly.getAtomArray();
            int nAtoms = atoms.length;
            x = new double[nAtoms * 3];
        }
        /**
         * Compute the number of permutations without eliminating dead-ends
         * and compute the number of permutations using singleton
         * elimination.
         */
        double permutations = 1;
        double singletonPermutations = 1;
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int nr = rotamers.length;
            if (nr > 1) {
                int nrot = 0;
                for (int ri = 0; ri < nr; ri++) {
                    if (!eliminatedSingles[i][ri]) {
                        nrot++;
                    }
                }
                permutations *= rotamers.length;
                if (nrot > 1) {
                    singletonPermutations *= nrot;
                }
            }
        }
        logIfMaster(format(" Collecting Permutations:"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
        dryRun(residues, 0, currentRotamers);
        double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
        double afterPairElim = singletonPermutations - pairTotalElimination;
        currentEnsemble = (int) evaluatedPermutations;
        if (ensembleNumber == 1 && currentEnsemble == 0) {
            logger.severe(" No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
        }
        if (ensembleNumber > 1) {
            if (master && printFiles && ensembleFile == null) {
                File file = molecularAssembly.getFile();
                String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
                ensembleFile = new File(filename + ".ens");
                if (ensembleFile.exists()) {
                    for (int i = 2; i < 1000; i++) {
                        ensembleFile = new File(filename + ".ens_" + i);
                        if (!ensembleFile.exists()) {
                            break;
                        }
                    }
                    if (ensembleFile.exists()) {
                        logger.warning(format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
                    }
                }
                ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null, null);
                logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
            }
            logIfMaster(format(" Ensemble Search Stats: (buffer: %5.3f, current: %d, target: %d)", ensembleBuffer, currentEnsemble, ensembleNumber));
        }
        if (ensembleNumber == 1 || finalTry) {
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Condition", "|", "Number of Permutations Left", "|", "Number of Permutations Removed", "|"));
            logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "No Eliminations", "|", permutations, "|", "", "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single Eliminations", "|", singletonPermutations, "|", permutations - singletonPermutations, "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Pair Eliminations", "|", afterPairElim, "|", pairTotalElimination, "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single and Pair Eliminations", "|", (double) evaluatedPermutations, "|", pairTotalElimination + (permutations - singletonPermutations), "|"));
            logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------\n"));
            logIfMaster(format(" Energy of permutations:"));
            logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
            logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Permutation", "|", "Energy", "|", "Lowest Possible Energy", "|"));
            logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
            break;
        }
        if (Math.abs(currentEnsemble - ensembleNumber) < bestEnsembleTargetDiffThusFar) {
            bestEnsembleTargetDiffThusFar = Math.abs(currentEnsemble - ensembleNumber);
            bestBufferThusFar = ensembleBuffer;
        }
        if (currentEnsemble > ensembleNumber) {
            ensembleBuffer -= ensembleBufferStep;
            ensembleBufferStep -= (ensembleBufferStep * 0.01);
            iterations++;
        } else if (currentEnsemble < ensembleNumber) {
            ensembleBuffer += ensembleBufferStep;
            ensembleBufferStep -= (ensembleBufferStep * 0.01);
            iterations++;
        }
        if (iterations > 100) {
            if (currentEnsemble == 0) {
                // TODO: Decide whether we like these next four lines.  Has the potential to produce a crazy amount of permutations.
                logIfMaster(" Ensemble still empty; increasing buffer energy.");
                startingBuffer = 3 * startingBuffer;
                setEnsemble(10, startingBuffer);
                iterations = 0;
            } else {
                ensembleBuffer = bestBufferThusFar;
                finalTry = true;
            }
        }
    }
    if (currentEnsemble == 0) {
        logger.warning(" No valid rotamer permutations found; results will be unreliable.  Try increasing the starting ensemble buffer.");
    }
    double[] permutationEnergyStub = null;
    if (useMonteCarlo()) {
        firstValidPerm(residues, 0, currentRotamers);
        rotamerOptimizationMC(residues, optimum, currentRotamers, nMCsteps, false, mcUseAll);
    } else {
        rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE, optimum, permutationEnergyStub);
    }
    double[] residueEnergy = new double[nResidues];
    double sumSelfEnergy = 0;
    double sumLowSelfEnergy = 0;
    logIfMaster(format("%s", " ----------------------------------------------------------------------------------\n"));
    logIfMaster(format(" Energy contributions:"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------"));
    logIfMaster(format(" %15s %5s %25s %5s %25s %5s", "Type", "|", "Energy", "|", "Lowest Possible Energy", "|"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------"));
    for (int i = 0; i < nResidues; i++) {
        int ri = optimum[i];
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        turnOnAtoms(residue);
        RotamerLibrary.applyRotamer(residue, rotamers[ri]);
        double self = getSelf(i, ri);
        residueEnergy[i] = self;
        sumSelfEnergy += self;
        double lowest = lowestSelfEnergy(residues, i);
        sumLowSelfEnergy += lowest;
        if (self - lowest > 10.0) {
            logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self (" + residues[i] + "," + ri + "):", "|", self, "|", lowest, "|"));
        }
    }
    double sumPairEnergy = 0.0;
    double sumLowPairEnergy = 0.0;
    double[] resPairEnergy = new double[nResidues];
    double[] lowPairEnergy = new double[nResidues];
    for (int i = 0; i < nResidues - 1; i++) {
        StringBuilder sb = new StringBuilder();
        int ri = optimum[i];
        double sumPairEnergyI = 0;
        double sumLowPairEnergyI = 0;
        for (int j = i + 1; j < nResidues; j++) {
            int rj = optimum[j];
            double pair = get2Body(i, ri, j, rj);
            residueEnergy[i] += 0.5 * pair;
            residueEnergy[j] += 0.5 * pair;
            sumPairEnergy += pair;
            sumPairEnergyI += pair;
            double lowest = lowestPairEnergy(residues, i, ri, j);
            sumLowPairEnergy += lowest;
            sumLowPairEnergyI += lowest;
            resPairEnergy[i] = 0.5 * pair;
            resPairEnergy[j] = 0.5 * pair;
            lowPairEnergy[i] = 0.5 * lowest;
            lowPairEnergy[j] = 0.5 * lowest;
            if (resPairEnergy[i] - lowPairEnergy[i] > 10.0) {
                sb.append(format("  Pair Energy (%8s,%2d) (%8s,%2d): %12.4f (Lowest: %12.4f).\n", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, pair, lowest));
            }
        }
        if (sumPairEnergyI - sumLowPairEnergyI > 10.0) {
            logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self (" + residues[i] + "," + ri + "):", "|", sumPairEnergyI, "|", sumLowPairEnergyI, "|"));
            sb.trimToSize();
            if (!sb.toString().isEmpty()) {
                logIfMaster(sb.toString());
            }
        }
    }
    double e = Double.NaN;
    try {
        e = currentEnergy(residueList);
    } catch (ArithmeticException ex) {
        logger.severe(String.format(" Exception %s in calculating current energy at the end of self and pairs", ex.toString()));
    }
    logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Backbone:", "|", backboneEnergy, "|", "", "|"));
    logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self:", "|", sumSelfEnergy, "|", sumLowSelfEnergy, "|"));
    logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Pair:", "|", sumPairEnergy, "|", sumLowPairEnergy, "|"));
    double approximateEnergy = backboneEnergy + sumSelfEnergy + sumPairEnergy;
    double sumTrimerEnergy = 0;
    if (threeBodyTerm) {
        for (int i = 0; i < nResidues - 2; i++) {
            int ri = optimum[i];
            for (int j = i + 1; j < nResidues - 1; j++) {
                int rj = optimum[j];
                for (int k = j + 1; k < nResidues; k++) {
                    int rk = optimum[k];
                    try {
                        double triple = get3Body(i, ri, j, rj, k, rk);
                        double thirdTrip = triple / 3.0;
                        residueEnergy[i] += thirdTrip;
                        residueEnergy[j] += thirdTrip;
                        residueEnergy[k] += thirdTrip;
                        sumTrimerEnergy += triple;
                    } catch (Exception ex) {
                        logger.warning(ex.toString());
                    }
                }
            }
        }
        approximateEnergy += sumTrimerEnergy;
        double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - backboneEnergy;
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Trimer:", "|", sumTrimerEnergy, "|", "", "|"));
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
    } else {
        double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - backboneEnergy;
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
    }
    logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Approximate:", "|", approximateEnergy, "|", "", "|"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------\n"));
    logIfMaster(format(" Final rotamers:"));
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
    logIfMaster(format("%17s %3s %10s %3s %9s %3s %9s %3s %9s %3s %10s %3s", "Residue", "|", "Chi 1", "|", "Chi 2", "|", "Chi 3", "|", "Chi 4", "|", "Energy", "|"));
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int ri = optimum[i];
        Rotamer rotamer = rotamers[ri];
        logIfMaster(format(" %3d %c (%7s,%2d) | %s %12.4f |", i + 1, residue.getChainID(), residue, ri, rotamer.toAngleString(), residueEnergy[i]));
        RotamerLibrary.applyRotamer(residue, rotamer);
    }
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------\n"));
    return e;
}

Example 17 with PDBFilter

use of ffx.potential.parsers.PDBFilter in project ffx by mjschnie.

the class PhMD method writeSnapshot.

private void writeSnapshot(String extension) {
    String filename = FilenameUtils.removeExtension(mola.getFile().toString()) + extension + snapshotIndex;
    if (config.snapshots == Snapshots.INTERLEAVED) {
        filename = mola.getFile().getAbsolutePath();
        if (!filename.contains("dyn")) {
            filename = FilenameUtils.removeExtension(filename) + "_dyn.pdb";
        }
    }
    File file = new File(filename);
    PDBFilter writer = new PDBFilter(file, mola, null, null);
    writer.writeFile(file, false);
}

Example 18 with PDBFilter

use of ffx.potential.parsers.PDBFilter in project ffx by mjschnie.

the class PhMD method writeSnapshot.

private void writeSnapshot(boolean beforeChange, StepType stepType, Snapshots snapshotsType) {
    // Write the after-step snapshot.
    if (snapshotsType != Snapshots.NONE) {
        String postfixA = ".";
        switch(stepType) {
            case TITRATE:
                postfixA = ".pro";
                break;
            case ROTAMER:
                postfixA = ".rot";
                break;
            case COMBO:
                postfixA = ".cbo";
                break;
        }
        String postfixB = (beforeChange) ? "S-" : "F-";
        String filename = FilenameUtils.removeExtension(mola.getFile().toString()) + postfixA + postfixB + snapshotIndex;
        if (snapshotsType == Snapshots.INTERLEAVED) {
            filename = mola.getFile().getAbsolutePath();
            if (!filename.contains("dyn")) {
                filename = FilenameUtils.removeExtension(filename) + "_dyn.pdb";
            }
        }
        File afterFile = new File(filename);
        PDBFilter afterWriter = new PDBFilter(afterFile, mola, null, null);
        afterWriter.writeFile(afterFile, false);
    }
}

Example 19 with PDBFilter

use of ffx.potential.parsers.PDBFilter in project ffx by mjschnie.

the class TransitionTemperedOSRW method writeTraversal.

private void writeTraversal() {
    double heldTraversalLambda = 0.5;
    if (!traversalInHand.isEmpty()) {
        heldTraversalLambda = Double.parseDouble(traversalInHand.get(0).split(",")[0]);
        if ((lambda > 0.2 && traversalSnapshotTarget == 0) || (lambda < 0.8 && traversalSnapshotTarget == 1)) {
            int snapshotCounts = Integer.parseInt(traversalInHand.get(0).split(",")[1]);
            traversalInHand.remove(0);
            File fileToWrite;
            int numStructures;
            if (traversalSnapshotTarget == 0) {
                fileToWrite = lambdaZeroFile;
                numStructures = ++lambdaZeroStructures;
            } else {
                fileToWrite = lambdaOneFile;
                numStructures = ++lambdaOneStructures;
            }
            try {
                FileWriter fw = new FileWriter(fileToWrite, true);
                BufferedWriter bw = new BufferedWriter(fw);
                bw.write(String.format("MODEL        %d          L=%.4f  counts=%d", numStructures, heldTraversalLambda, snapshotCounts));
                for (int i = 0; i < 50; i++) {
                    bw.write(" ");
                }
                bw.newLine();
                for (int i = 0; i < traversalInHand.size(); i++) {
                    bw.write(traversalInHand.get(i));
                    bw.newLine();
                }
                bw.write(String.format("ENDMDL"));
                for (int i = 0; i < 75; i++) {
                    bw.write(" ");
                }
                bw.newLine();
                bw.close();
                logger.info(String.format(" Wrote traversal structure L=%.4f", heldTraversalLambda));
            } catch (Exception exception) {
                logger.warning(String.format("Exception writing to file: %s", fileToWrite.getName()));
            }
            heldTraversalLambda = 0.5;
            traversalInHand.clear();
            traversalSnapshotTarget = 1 - traversalSnapshotTarget;
        }
    }
    if (((lambda < 0.1 && traversalInHand.isEmpty()) || (lambda < heldTraversalLambda - 0.025 && !traversalInHand.isEmpty())) && (traversalSnapshotTarget == 0 || traversalSnapshotTarget == -1)) {
        if (lambdaZeroFilter == null) {
            lambdaZeroFilter = new PDBFilter(lambdaZeroFile, lambdaZeroAssembly, null, null);
            lambdaZeroFilter.setListMode(true);
        }
        lambdaZeroFilter.clearListOutput();
        lambdaZeroFilter.writeFileWithHeader(lambdaFile, format("%.4f,%d", lambda, totalWeight));
        traversalInHand = lambdaZeroFilter.getListOutput();
        traversalSnapshotTarget = 0;
    } else if (((lambda > 0.9 && traversalInHand.isEmpty()) || (lambda > heldTraversalLambda + 0.025 && !traversalInHand.isEmpty())) && (traversalSnapshotTarget == 1 || traversalSnapshotTarget == -1)) {
        if (lambdaOneFilter == null) {
            lambdaOneFilter = new PDBFilter(lambdaOneFile, lambdaOneAssembly, null, null);
            lambdaOneFilter.setListMode(true);
        }
        lambdaOneFilter.clearListOutput();
        lambdaOneFilter.writeFileWithHeader(lambdaFile, format("%.4f,%d", lambda, totalWeight));
        traversalInHand = lambdaOneFilter.getListOutput();
        traversalSnapshotTarget = 1;
    }
}

Example 20 with PDBFilter

use of ffx.potential.parsers.PDBFilter in project ffx by mjschnie.

the class RosenbluthChiAllMove method engage_control.

/**
 * For validation. Performs Monte Carlo chi moves WITHOUT biasing. Randomly
 * select one chi, give it a random theta. Accept on the vanilla Metropolis
 * criterion.
 */
private boolean engage_control() {
    report.append(String.format(" Rosenbluth Control Move: %4d  %s\n", moveNumber, target));
    double origEnergy = totalEnergy();
    double[] origChi = RotamerLibrary.measureRotamer(target, false);
    int chiIndex = rand.nextInt(origChi.length);
    double theta = rand.nextDouble(360.0) - 180;
    double[] newChi = new double[origChi.length];
    System.arraycopy(origChi, 0, newChi, 0, origChi.length);
    newChi[chiIndex] = theta;
    Rotamer newState = createRotamer(target, newChi);
    RotamerLibrary.applyRotamer(target, newState);
    double finalEnergy = totalEnergy();
    double dU = finalEnergy - origEnergy;
    double criterion = FastMath.exp(-beta * dU);
    double rng = rand.nextDouble();
    report.append(String.format("    move (chi,theta): %d %5.1f\n", chiIndex, theta));
    report.append(String.format("    orig, final, dU:  %.2g %.2g %.2g\n", origEnergy, finalEnergy, dU));
    report.append(String.format("    crit, rng:        %.2g %.2g\n", criterion, rng));
    if (rng < criterion) {
        accepted = true;
        report.append(String.format("    Accepted!\n"));
        PDBFilter writer = new PDBFilter(mola.getFile(), mola, null, null);
        String filename = FilenameUtils.removeExtension(mola.getFile().toString());
        filename = mola.getFile().getAbsolutePath();
        if (!filename.contains("_mc")) {
            filename = FilenameUtils.removeExtension(filename) + "_mc.pdb";
        }
        File file = new File(filename);
        writer.writeFile(file, false);
    } else {
        accepted = false;
        report.append(String.format("    Denied.\n"));
        target.revertState(origState);
    }
    updateAll();
    if (verbose) {
        logger.info(report.toString());
    }
    return (rng < criterion);
}