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

Example 56 with Atom

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

the class RotamerOptimization method independent.

private double independent(List<Residue> residues) {
    if (x == null) {
        Atom[] atoms = molecularAssembly.getAtomArray();
        int nAtoms = atoms.length;
        x = new double[nAtoms * 3];
    }
    double e = Double.MAX_VALUE;
    List<Residue> rList = new ArrayList<>(Collections.nCopies(1, null));
    for (Residue residue : residues) {
        rList.set(0, residue);
        logger.info(format(" Optimizing %s side-chain.", residue));
        Rotamer[] rotamers = residue.getRotamers(library);
        potential.getCoordinates(x);
        e = Double.MAX_VALUE;
        int bestRotamer = -1;
        for (int j = 0; j < rotamers.length; j++) {
            Rotamer rotamer = rotamers[j];
            RotamerLibrary.applyRotamer(residue, rotamer);
            if (algorithmListener != null) {
                algorithmListener.algorithmUpdate(molecularAssembly);
            }
            double newE = Double.NaN;
            try {
                newE = currentEnergy(rList);
            } catch (ArithmeticException ex) {
                logger.fine(String.format(" Exception %s in energy calculations during independent for %s-%d", ex.toString(), residue, j));
            }
            if (newE < e) {
                bestRotamer = j;
            }
        }
        if (bestRotamer > -1) {
            Rotamer rotamer = rotamers[bestRotamer];
            RotamerLibrary.applyRotamer(residue, rotamer);
        }
        if (algorithmListener != null) {
            algorithmListener.algorithmUpdate(molecularAssembly);
        }
        if (terminate) {
            logger.info(format("\n Terminating after residue %s.\n", residue));
            break;
        }
    }
    return e;
}

Example 57 with Atom

use of ffx.potential.bonded.Atom 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 58 with Atom

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

the class PhMD method mcUpdate.

/**
 * The primary driver. Called by the MD engine at each dynamics step.
 */
@Override
public boolean mcUpdate(double temperature) {
    startTime = System.nanoTime();
    if (thermostat.getCurrentTemperature() > config.meltdownTemperature) {
        meltdown();
    }
    if (thermostat.getCurrentTemperature() > config.warningTemperature) {
        Atom[] atoms = mola.getAtomArray();
        logger.info(format(" System heating! Dumping atomic velocities for %d D.o.F.:", ffe.getNumberOfVariables()));
        double[] velocity = new double[3];
        for (Atom atom : atoms) {
            atom.getVelocity(velocity);
            logger.info(format(" %s: %s", atom.describe(Atom.Descriptions.Trim), Arrays.toString(velocity)));
        }
    }
    esvSystem.setTemperature(temperature);
    propagateInactiveResidues(titratingMultis);
    stepCount++;
    // Decide on the type of step to be taken.
    StepType stepType;
    if (stepCount % mcStepFrequency == 0 && stepCount % rotamerStepFrequency == 0) {
        stepType = StepType.COMBO;
    } else if (stepCount % mcStepFrequency == 0) {
        stepType = StepType.TITRATE;
    } else if (stepCount % rotamerStepFrequency == 0) {
        stepType = StepType.ROTAMER;
    } else {
        // Not yet time for an MC step, return to MD.
        if (config.logTimings) {
            long took = System.nanoTime() - startTime;
            logger.info(String.format(" CpHMD propagation time: %.6f", took * NS_TO_SEC));
        }
        return false;
    }
    logger.info(format("TitratingMultis: %d", titratingMultis.size()));
    // Randomly choose a target titratable residue to attempt protonation switch.
    int random = (config.titrateTermini) ? rng.nextInt(titratingMultis.size() + titratingTermini.size()) : rng.nextInt(titratingMultis.size());
    if (random >= titratingMultis.size()) {
        Residue target = titratingTermini.get(random - titratingMultis.size());
        boolean accepted = tryTerminusTitration((MultiTerminus) target);
        snapshotIndex++;
        if (accepted) {
            molDyn.reInit();
            previousTarget = target;
        }
        return accepted;
    }
    MultiResidue targetMulti = titratingMultis.get(random);
    // Check whether rotamer moves are possible for the selected residue.
    Residue targetMultiActive = targetMulti.getActive();
    Rotamer[] targetMultiRotamers = targetMultiActive.getRotamers(library);
    if (targetMultiRotamers == null || targetMultiRotamers.length <= 1) {
        if (stepType == StepType.ROTAMER) {
            return false;
        } else if (stepType == StepType.COMBO) {
            stepType = StepType.TITRATE;
        }
    }
    // Perform the MC move.
    boolean accepted;
    switch(stepType) {
        case TITRATE:
            accepted = tryTitrationStep(targetMulti);
            break;
        case ROTAMER:
            accepted = (config.useConformationalBias) ? tryCBMCStep(targetMulti) : tryRotamerStep(targetMulti);
            break;
        case COMBO:
            accepted = (config.useConformationalBias) ? tryCBMCStep(targetMulti) || tryTitrationStep(targetMulti) : tryComboStep(targetMulti);
            break;
        default:
            accepted = false;
            throw new IllegalStateException();
    }
    snapshotIndex++;
    if (accepted) {
        previousTarget = targetMulti;
    }
    if (config.logTimings) {
        long took = System.nanoTime() - startTime;
        logger.info(String.format(" CpHMD step time:        %.6f", took * NS_TO_SEC));
    }
    return accepted;
}

Example 59 with Atom

use of ffx.potential.bonded.Atom 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 60 with Atom

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