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

Example 6 with MultiResidue

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

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

the class TitrationUtils method propagateInactiveResidues.

/**
 * Copies atomic coordinates from each active residue to its inactive
 * counterparts. Inactive hydrogen coordinates are updated by geometry with
 * the propagated heavies.
 */
public static void propagateInactiveResidues(MultiResidue multiRes, boolean propagateDynamics) {
    // Propagate all atom coordinates from active residues to their inactive counterparts.
    Residue active = multiRes.getActive();
    String activeResName = active.getName();
    List<Residue> inactives = multiRes.getInactive();
    for (Atom activeAtom : active.getAtomList()) {
        String activeName = activeAtom.getName();
        for (Residue inactive : inactives) {
            Atom inactiveAtom = (Atom) inactive.getAtomNode(activeName);
            if (inactiveAtom != null) {
                // Propagate position and gradient.
                double[] activeXYZ = activeAtom.getXYZ(null);
                inactiveAtom.setXYZ(activeXYZ);
                double[] grad = new double[3];
                activeAtom.getXYZGradient(grad);
                inactiveAtom.setXYZGradient(grad[0], grad[1], grad[2]);
                if (propagateDynamics) {
                    // Propagate velocity, acceleration, and previous acceleration.
                    double[] activeVelocity = new double[3];
                    activeAtom.getVelocity(activeVelocity);
                    inactiveAtom.setVelocity(activeVelocity);
                    double[] activeAccel = new double[3];
                    activeAtom.getAcceleration(activeAccel);
                    inactiveAtom.setAcceleration(activeAccel);
                    double[] activePrevAcc = new double[3];
                    activeAtom.getPreviousAcceleration(activePrevAcc);
                    inactiveAtom.setPreviousAcceleration(activePrevAcc);
                }
            } else {
                if (activeName.equals("C") || activeName.equals("O") || activeName.equals("N") || activeName.equals("CA") || activeName.equals("H") || activeName.equals("HA")) {
                // Backbone atoms aren't supposed to exist in inactive multiResidue components; so no problem.
                } else if (isTitratableHydrogen(activeAtom)) {
                /**
                 * i.e. ((activeResName.equals("LYS") &&
                 * activeName.equals("HZ3")) ||
                 * (activeResName.equals("TYR") &&
                 * activeName.equals("HH")) ||
                 * (activeResName.equals("CYS") &&
                 * activeName.equals("HG")) ||
                 * (activeResName.equals("HIS") &&
                 * (activeName.equals("HD1") ||
                 * activeName.equals("HE2"))) ||
                 * (activeResName.equals("HID") &&
                 * activeName.equals("HD1")) ||
                 * (activeResName.equals("HIE") &&
                 * activeName.equals("HE2")) ||
                 * (activeResName.equals("ASH") &&
                 * activeName.equals("HD2")) ||
                 * (activeResName.equals("GLH") &&
                 * activeName.equals("HE2")))
                 */
                // These titratable protons are handled below; so no problem.
                } else {
                    // Now we have a problem.
                    logger.warning(format("Couldn't propagate inactive MultiResidue atom: %s: %s, %s", multiRes, activeName, activeAtom));
                }
            }
        }
    }
    rebuildStrandedProtons(multiRes);
}

Example 8 with MultiResidue

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

the class TitrationUtils method rebuildStrandedProtons.

/**
 * Rebuild stranded titratable protons from ideal geometry. "Stranded
 * protons" are titrating H+ atoms on inactive MultiRes members; when
 * propagating new coordinates to inactive residues, no coords/velocity
 * exist for them.
 */
private static void rebuildStrandedProtons(MultiResidue multiRes) {
    // If inactive residue is a protonated form, move the stranded hydrogen to new coords (based on propagated heavies).
    // Also give the stranded hydrogen a maxwell velocity and remove its accelerations.
    List<Residue> inactives = multiRes.getInactive();
    for (Residue inactive : inactives) {
        List<Atom> resetMe = new ArrayList<>();
        switch(inactive.getName()) {
            case "LYS":
                {
                    Atom HZ3 = (Atom) inactive.getAtomNode("HZ3");
                    Atom NZ = (Atom) inactive.getAtomNode("NZ");
                    Atom CE = (Atom) inactive.getAtomNode("CE");
                    Atom HZ1 = (Atom) inactive.getAtomNode("HZ1");
                    BondedUtils.intxyz(HZ3, NZ, 1.02, CE, 109.5, HZ1, 109.5, -1);
                    resetMe.add(HZ3);
                    break;
                }
            case "ASH":
                {
                    Atom HD2 = (Atom) inactive.getAtomNode("HD2");
                    Atom OD2 = (Atom) inactive.getAtomNode("OD2");
                    Atom CG = (Atom) inactive.getAtomNode("CG");
                    Atom OD1 = (Atom) inactive.getAtomNode("OD1");
                    BondedUtils.intxyz(HD2, OD2, 0.98, CG, 108.7, OD1, 0.0, 0);
                    resetMe.add(HD2);
                    break;
                }
            case "GLH":
                {
                    Atom HE2 = (Atom) inactive.getAtomNode("HE2");
                    Atom OE2 = (Atom) inactive.getAtomNode("OE2");
                    Atom CD = (Atom) inactive.getAtomNode("CD");
                    Atom OE1 = (Atom) inactive.getAtomNode("OE1");
                    BondedUtils.intxyz(HE2, OE2, 0.98, CD, 108.7, OE1, 0.0, 0);
                    resetMe.add(HE2);
                    break;
                }
            case "HIS":
                {
                    Atom HE2 = (Atom) inactive.getAtomNode("HE2");
                    Atom NE2 = (Atom) inactive.getAtomNode("NE2");
                    Atom CD2 = (Atom) inactive.getAtomNode("CD2");
                    Atom CE1 = (Atom) inactive.getAtomNode("CE1");
                    Atom HD1 = (Atom) inactive.getAtomNode("HD1");
                    Atom ND1 = (Atom) inactive.getAtomNode("ND1");
                    Atom CG = (Atom) inactive.getAtomNode("CG");
                    Atom CB = (Atom) inactive.getAtomNode("CB");
                    BondedUtils.intxyz(HE2, NE2, 1.02, CD2, 126.0, CE1, 126.0, 1);
                    BondedUtils.intxyz(HD1, ND1, 1.02, CG, 126.0, CB, 0.0, 0);
                    resetMe.add(HE2);
                    resetMe.add(HD1);
                    break;
                }
            case "HID":
                {
                    Atom HD1 = (Atom) inactive.getAtomNode("HD1");
                    Atom ND1 = (Atom) inactive.getAtomNode("ND1");
                    Atom CG = (Atom) inactive.getAtomNode("CG");
                    Atom CB = (Atom) inactive.getAtomNode("CB");
                    BondedUtils.intxyz(HD1, ND1, 1.02, CG, 126.0, CB, 0.0, 0);
                    resetMe.add(HD1);
                    break;
                }
            case "HIE":
                {
                    Atom HE2 = (Atom) inactive.getAtomNode("HE2");
                    Atom NE2 = (Atom) inactive.getAtomNode("NE2");
                    Atom CD2 = (Atom) inactive.getAtomNode("CD2");
                    Atom CE1 = (Atom) inactive.getAtomNode("CE1");
                    BondedUtils.intxyz(HE2, NE2, 1.02, CD2, 126.0, CE1, 126.0, 1);
                    resetMe.add(HE2);
                    break;
                }
            case "CYS":
                {
                    Atom HG = (Atom) inactive.getAtomNode("HG");
                    Atom SG = (Atom) inactive.getAtomNode("SG");
                    Atom CB = (Atom) inactive.getAtomNode("CB");
                    Atom CA = (Atom) inactive.getAtomNode("CA");
                    BondedUtils.intxyz(HG, SG, 1.34, CB, 96.0, CA, 180.0, 0);
                    resetMe.add(HG);
                    break;
                }
            case "TYR":
                {
                    Atom HH = (Atom) inactive.getAtomNode("HH");
                    Atom OH = (Atom) inactive.getAtomNode("OH");
                    Atom CZ = (Atom) inactive.getAtomNode("CZ");
                    Atom CE2 = (Atom) inactive.getAtomNode("CE2");
                    BondedUtils.intxyz(HH, OH, 0.97, CZ, 108.0, CE2, 0.0, 0);
                    resetMe.add(HH);
                    break;
                }
            default:
        }
        for (Atom a : resetMe) {
            if (heavyStrandedDynamics) {
                // Use of heavy atom dynamics properties is in testing.
                a.setXYZGradient(0, 0, 0);
                double[] heavyVelocity = new double[3];
                double[] heavyAccel = new double[3];
                double[] heavyPrevAccel = new double[3];
                Atom heavy = a.getBonds().get(0).get1_2(a);
                heavy.getVelocity(heavyVelocity);
                heavy.getAcceleration(heavyAccel);
                heavy.getPreviousAcceleration(heavyPrevAccel);
                a.setVelocity(heavyVelocity);
                a.setAcceleration(heavyAccel);
                a.setPreviousAcceleration(heavyPrevAccel);
            } else {
                // PREVIOUSLY: draw vel from maxwell and set accel to zero
                a.setXYZGradient(0, 0, 0);
                a.setVelocity(ExtUtils.maxwellVelocity(a.getMass(), ExtConstants.roomTemperature));
                a.setAcceleration(new double[] { 0, 0, 0 });
                a.setPreviousAcceleration(new double[] { 0, 0, 0 });
            }
        }
    }
}

Example 9 with MultiResidue

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

the class ForceFieldEnergy method energy.

/**
 * <p>
 * energy</p>
 *
 * @param gradient a boolean.
 * @param print a boolean.
 * @return a double.
 */
public double energy(boolean gradient, boolean print) {
    try {
        bondTime = 0;
        angleTime = 0;
        stretchBendTime = 0;
        ureyBradleyTime = 0;
        outOfPlaneBendTime = 0;
        torsionTime = 0;
        piOrbitalTorsionTime = 0;
        torsionTorsionTime = 0;
        improperTorsionTime = 0;
        vanDerWaalsTime = 0;
        electrostaticTime = 0;
        restraintBondTime = 0;
        ncsTime = 0;
        coordRestraintTime = 0;
        totalTime = System.nanoTime();
        // Zero out the potential energy of each bonded term.
        bondEnergy = 0.0;
        angleEnergy = 0.0;
        stretchBendEnergy = 0.0;
        ureyBradleyEnergy = 0.0;
        outOfPlaneBendEnergy = 0.0;
        torsionEnergy = 0.0;
        piOrbitalTorsionEnergy = 0.0;
        torsionTorsionEnergy = 0.0;
        improperTorsionEnergy = 0.0;
        totalBondedEnergy = 0.0;
        // Zero out potential energy of restraint terms
        restraintBondEnergy = 0.0;
        ncsEnergy = 0.0;
        restrainEnergy = 0.0;
        // Zero out bond and angle RMSDs.
        bondRMSD = 0.0;
        angleRMSD = 0.0;
        // Zero out the potential energy of each non-bonded term.
        vanDerWaalsEnergy = 0.0;
        permanentMultipoleEnergy = 0.0;
        permanentRealSpaceEnergy = 0.0;
        permanentSelfEnergy = 0.0;
        permanentReciprocalEnergy = 0.0;
        polarizationEnergy = 0.0;
        inducedRealSpaceEnergy = 0.0;
        inducedSelfEnergy = 0.0;
        inducedReciprocalEnergy = 0.0;
        totalMultipoleEnergy = 0.0;
        totalNonBondedEnergy = 0.0;
        // Zero out the solvation energy.
        solvationEnergy = 0.0;
        // Zero out the relative solvation energy (sequence optimization)
        relativeSolvationEnergy = 0.0;
        nRelativeSolvations = 0;
        esvBias = 0.0;
        // Zero out the total potential energy.
        totalEnergy = 0.0;
        // Zero out the Cartesian coordinate gradient for each atom.
        if (gradient) {
            for (int i = 0; i < nAtoms; i++) {
                atoms[i].setXYZGradient(0.0, 0.0, 0.0);
                atoms[i].setLambdaXYZGradient(0.0, 0.0, 0.0);
            }
        }
        /**
         * Computed the bonded energy terms in parallel.
         */
        try {
            bondedRegion.setGradient(gradient);
            parallelTeam.execute(bondedRegion);
        } catch (RuntimeException ex) {
            logger.warning("Runtime exception during bonded term calculation.");
            throw ex;
        } catch (Exception ex) {
            Utilities.printStackTrace(ex);
            logger.severe(ex.toString());
        }
        if (!lambdaBondedTerms) {
            /**
             * Compute restraint terms.
             */
            if (ncsTerm) {
                ncsTime = -System.nanoTime();
                ncsEnergy = ncsRestraint.residual(gradient, print);
                ncsTime += System.nanoTime();
            }
            if (restrainTerm && !coordRestraints.isEmpty()) {
                coordRestraintTime = -System.nanoTime();
                for (CoordRestraint restraint : coordRestraints) {
                    restrainEnergy += restraint.residual(gradient, print);
                }
                coordRestraintTime += System.nanoTime();
            }
            if (comTerm) {
                comRestraintTime = -System.nanoTime();
                comRestraintEnergy = comRestraint.residual(gradient, print);
                comRestraintTime += System.nanoTime();
            }
            /**
             * Compute non-bonded terms.
             */
            if (vanderWaalsTerm) {
                vanDerWaalsTime = -System.nanoTime();
                vanDerWaalsEnergy = vanderWaals.energy(gradient, print);
                nVanDerWaalInteractions = this.vanderWaals.getInteractions();
                vanDerWaalsTime += System.nanoTime();
            }
            if (multipoleTerm) {
                electrostaticTime = -System.nanoTime();
                totalMultipoleEnergy = particleMeshEwald.energy(gradient, print);
                permanentMultipoleEnergy = particleMeshEwald.getPermanentEnergy();
                permanentRealSpaceEnergy = particleMeshEwald.getPermRealEnergy();
                permanentSelfEnergy = particleMeshEwald.getPermSelfEnergy();
                permanentReciprocalEnergy = particleMeshEwald.getPermRecipEnergy();
                polarizationEnergy = particleMeshEwald.getPolarizationEnergy();
                inducedRealSpaceEnergy = particleMeshEwald.getIndRealEnergy();
                inducedSelfEnergy = particleMeshEwald.getIndSelfEnergy();
                inducedReciprocalEnergy = particleMeshEwald.getIndRecipEnergy();
                nPermanentInteractions = particleMeshEwald.getInteractions();
                solvationEnergy = particleMeshEwald.getGKEnergy();
                nGKInteractions = particleMeshEwald.getGKInteractions();
                electrostaticTime += System.nanoTime();
            }
        }
        if (relativeSolvationTerm) {
            List<Residue> residuesList = molecularAssembly.getResidueList();
            for (Residue residue : residuesList) {
                if (residue instanceof MultiResidue) {
                    Atom refAtom = residue.getSideChainAtoms().get(0);
                    if (refAtom != null && refAtom.getUse()) {
                        /**
                         * Reasonably confident that it should be -=, as we
                         * are trying to penalize residues with strong
                         * solvation energy.
                         */
                        double thisSolvation = relativeSolvation.getSolvationEnergy(residue, false);
                        relativeSolvationEnergy -= thisSolvation;
                        if (thisSolvation != 0) {
                            nRelativeSolvations++;
                        }
                    }
                }
            }
        }
        totalTime = System.nanoTime() - totalTime;
        totalBondedEnergy = bondEnergy + restraintBondEnergy + angleEnergy + stretchBendEnergy + ureyBradleyEnergy + outOfPlaneBendEnergy + torsionEnergy + piOrbitalTorsionEnergy + improperTorsionEnergy + torsionTorsionEnergy + ncsEnergy + restrainEnergy;
        totalNonBondedEnergy = vanDerWaalsEnergy + totalMultipoleEnergy + relativeSolvationEnergy;
        totalEnergy = totalBondedEnergy + totalNonBondedEnergy + solvationEnergy;
        if (esvTerm) {
            esvBias = esvSystem.getBiasEnergy();
            totalEnergy += esvBias;
        }
    } catch (EnergyException ex) {
        if (printOnFailure) {
            File origFile = molecularAssembly.getFile();
            String timeString = LocalDateTime.now().format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
            String filename = String.format("%s-ERROR-%s.pdb", FilenameUtils.removeExtension(molecularAssembly.getFile().getName()), timeString);
            PotentialsFunctions ef = new PotentialsUtils();
            filename = ef.versionFile(filename);
            logger.info(String.format(" Writing on-error snapshot to file %s", filename));
            ef.saveAsPDB(molecularAssembly, new File(filename));
            molecularAssembly.setFile(origFile);
        }
        if (ex.doCauseSevere()) {
            logger.log(Level.SEVERE, " Error in calculating energies or gradients", ex);
            return 0.0;
        } else {
            // Rethrow exception
            throw ex;
        }
    }
    if (print || printOverride) {
        if (printCompact) {
            logger.info(this.toString());
        } else {
            StringBuilder sb = new StringBuilder();
            if (gradient) {
                sb.append("\n Computed Potential Energy and Atomic Coordinate Gradients\n");
            } else {
                sb.append("\n Computed Potential Energy\n");
            }
            sb.append(this);
            logger.info(sb.toString());
        }
    }
    return totalEnergy;
}

Example 10 with MultiResidue

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

the class TitrationUtils method performTitration.

/**
 * Perform the requested titration on the given MultiResidue. Remember to
 * call reInit() on affected ForceFieldEnergy and MolecularDynamics objects!
 *
 * @param multiRes
 * @param titration
 */
public static TitrationType performTitration(MultiResidue multiRes, Titration titration, boolean inactivateBackground) {
    AminoAcid3 current = multiRes.getActive().getAminoAcid3();
    final TitrationType type;
    final AminoAcid3 target;
    if (current == titration.protForm) {
        type = TitrationType.DEPROT;
        target = titration.deprotForm;
    } else if (current == titration.deprotForm) {
        type = TitrationType.PROT;
        target = titration.protForm;
    } else {
        throw new IllegalStateException();
    }
    boolean success = multiRes.requestSetActiveResidue(target);
    if (!success) {
        logger.severe(String.format("Couldn't perform requested titration for MultiRes: %s", multiRes.toString()));
    }
    List<Atom> oldAtoms = multiRes.getActive().getAtomList();
    List<Atom> newAtoms = multiRes.getActive().getAtomList();
    // identify which atoms were actually inserted/removed
    List<Atom> removedAtoms = new ArrayList<>();
    List<Atom> insertedAtoms = new ArrayList<>();
    for (Atom oldAtom : oldAtoms) {
        boolean found = false;
        for (Atom newAtom : newAtoms) {
            if (newAtom == oldAtom || (newAtom.getResidueNumber() == oldAtom.getResidueNumber() && newAtom.getName().equals(oldAtom.getName()))) {
                found = true;
                break;
            }
        }
        if (!found) {
            removedAtoms.add(oldAtom);
        }
    }
    for (Atom newAtom : newAtoms) {
        boolean found = false;
        for (Atom oldAtom : oldAtoms) {
            if (newAtom == oldAtom || (newAtom.getResidueNumber() == oldAtom.getResidueNumber() && newAtom.getName().equals(oldAtom.getName()))) {
                found = true;
                break;
            }
        }
        if (!found) {
            insertedAtoms.add(newAtom);
        }
    }
    if (insertedAtoms.size() + removedAtoms.size() > 1) {
        logger.warning("Protonate: removed + inserted atom count > 1.");
    }
    if (inactivateBackground) {
        activateResidue(multiRes.getActive());
        for (Residue res : multiRes.getInactive()) {
            inactivateResidue(res);
        }
    }
    return type;
}