Examples with CoordRestraint ffx.potential.nonbonded.CoordRestraint used on opensource projects

Example 6 with CoordRestraint

use of ffx.potential.nonbonded.CoordRestraint in project ffx by mjschnie.

the class ForceFieldEnergy method getdEdXdL.

/**
 * {@inheritDoc}
 *
 * @param gradients
 */
@Override
public void getdEdXdL(double[] gradients) {
    if (!lambdaBondedTerms) {
        if (vanderWaalsTerm) {
            vanderWaals.getdEdXdL(gradients);
        }
        if (multipoleTerm) {
            particleMeshEwald.getdEdXdL(gradients);
        }
        if (restraintBondTerm) {
            for (int i = 0; i < nRestraintBonds; i++) {
                restraintBonds[i].getdEdXdL(gradients);
            }
        }
        if (ncsTerm && ncsRestraint != null) {
            ncsRestraint.getdEdXdL(gradients);
        }
        if (restrainTerm && !coordRestraints.isEmpty()) {
            for (CoordRestraint restraint : coordRestraints) {
                restraint.getdEdXdL(gradients);
            }
        // autoCoordRestraint.getdEdXdL(gradients);
        }
        if (comTerm && comRestraint != null) {
            comRestraint.getdEdXdL(gradients);
        }
        if (lambdaTorsions) {
            double[] grad = new double[3];
            int index = 0;
            for (int i = 0; i < nAtoms; i++) {
                Atom a = atoms[i];
                if (a.isActive()) {
                    a.getLambdaXYZGradient(grad);
                    gradients[index++] += grad[0];
                    gradients[index++] += grad[1];
                    gradients[index++] += grad[2];
                }
            }
        }
    }
}

Example 7 with CoordRestraint

use of ffx.potential.nonbonded.CoordRestraint 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 8 with CoordRestraint

use of ffx.potential.nonbonded.CoordRestraint in project ffx by mjschnie.

the class ConversionFilter method applyAtomProperties.

/**
 * Automatically sets atom-specific flags, particularly nouse and inactive,
 * and apply harmonic restraints. Intended to be called at the end of
 * convert() implementations.
 */
public void applyAtomProperties() {
    /**
     * What may be a more elegant implementation is to make convert() a
     * public concrete, but skeletal method, and then have convert()
     * call a protected abstract readFile method for each implementation.
     */
    Atom[] molaAtoms = activeMolecularAssembly.getAtomArray();
    int nmolaAtoms = molaAtoms.length;
    String[] nouseKeys = properties.getStringArray("nouse");
    for (String nouseKey : nouseKeys) {
        /*try {
                int[] nouseRange = parseAtNumArg("nouse", nouseKey, nmolaAtoms);
                logger.log(Level.INFO, String.format(" Atoms %d-%d set to be not "
                        + "used", nouseRange[0] + 1, nouseRange[1] + 1));
                for (int i = nouseRange[0]; i <= nouseRange[1]; i++) {
                    molaAtoms[i].setUse(false);
                }
            } catch (IllegalArgumentException ex) {
                logger.log(Level.INFO, ex.getLocalizedMessage());
            }*/
        String[] toks = nouseKey.split("\\s+");
        for (String tok : toks) {
            try {
                int[] nouseRange = parseAtNumArg("restraint", tok, nmolaAtoms);
                logger.info(String.format(" Setting atoms %d-%d to not be used", nouseRange[0] + 1, nouseRange[1] + 1));
                for (int j = nouseRange[0]; j <= nouseRange[1]; j++) {
                    molaAtoms[j].setUse(false);
                }
            } catch (IllegalArgumentException ex) {
                boolean atomFound = false;
                for (Atom atom : molaAtoms) {
                    if (atom.getName().equalsIgnoreCase(tok)) {
                        atomFound = true;
                        atom.setUse(false);
                    }
                }
                if (atomFound) {
                    logger.info(String.format(" Setting atoms with name %s to not be used", tok));
                } else {
                    logger.log(Level.INFO, ex.getLocalizedMessage());
                }
            }
        }
    }
    String[] inactiveKeys = properties.getStringArray("inactive");
    for (String inactiveKey : inactiveKeys) {
        try {
            int[] inactiveRange = parseAtNumArg("inactive", inactiveKey, nmolaAtoms);
            logger.log(Level.INFO, String.format(" Atoms %d-%d set to be not " + "used", inactiveRange[0] + 1, inactiveRange[1] + 1));
            for (int i = inactiveRange[0]; i <= inactiveRange[1]; i++) {
                molaAtoms[i].setActive(false);
            }
        } catch (IllegalArgumentException ex) {
            logger.log(Level.INFO, ex.getLocalizedMessage());
        }
    }
    coordRestraints = new ArrayList<>();
    String[] cRestraintStrings = properties.getStringArray("restraint");
    for (String coordRestraint : cRestraintStrings) {
        String[] toks = coordRestraint.split("\\s+");
        double forceconst = Double.parseDouble(toks[0]);
        logger.info(String.format(" Adding lambda-disabled coordinate restraint " + "with force constant %10.4f", forceconst));
        Set<Atom> restraintAtoms = new HashSet<>();
        for (int i = 1; i < toks.length; i++) {
            try {
                int[] nouseRange = parseAtNumArg("restraint", toks[i], nmolaAtoms);
                logger.info(String.format(" Adding atoms %d-%d to restraint", nouseRange[0], nouseRange[1]));
                for (int j = nouseRange[0]; j <= nouseRange[1]; j++) {
                    restraintAtoms.add(molaAtoms[j]);
                }
            } catch (IllegalArgumentException ex) {
                boolean atomFound = false;
                for (Atom atom : molaAtoms) {
                    if (atom.getName().equalsIgnoreCase(toks[i])) {
                        atomFound = true;
                        restraintAtoms.add(atom);
                    }
                }
                if (atomFound) {
                    logger.info(String.format(" Added atoms with name %s to restraint", toks[i]));
                } else {
                    logger.log(Level.INFO, ex.getLocalizedMessage());
                }
            }
        }
        if (!restraintAtoms.isEmpty()) {
            Atom[] ats = restraintAtoms.toArray(new Atom[restraintAtoms.size()]);
            coordRestraints.add(new CoordRestraint(ats, forceField, false, forceconst));
        } else {
            logger.warning(String.format(" Empty or unparseable restraint argument %s", coordRestraint));
        }
    }
    String[] lamRestraintStrings = properties.getStringArray("lamrestraint");
    for (String coordRestraint : lamRestraintStrings) {
        String[] toks = coordRestraint.split("\\s+");
        double forceconst = Double.parseDouble(toks[0]);
        logger.info(String.format(" Adding lambda-enabled coordinate restraint " + "with force constant %10.4f", forceconst));
        Set<Atom> restraintAtoms = new HashSet<>();
        for (int i = 1; i < toks.length; i++) {
            try {
                int[] nouseRange = parseAtNumArg("restraint", toks[i], nmolaAtoms);
                logger.info(String.format(" Adding atoms %d-%d to restraint", nouseRange[0], nouseRange[1]));
                for (int j = nouseRange[0]; j <= nouseRange[1]; j++) {
                    restraintAtoms.add(molaAtoms[j]);
                }
            } catch (IllegalArgumentException ex) {
                boolean atomFound = false;
                for (Atom atom : molaAtoms) {
                    if (atom.getName().equalsIgnoreCase(toks[i])) {
                        atomFound = true;
                        restraintAtoms.add(atom);
                    }
                }
                if (atomFound) {
                    logger.info(String.format(" Added atoms with name %s to restraint", toks[i]));
                } else {
                    logger.log(Level.INFO, ex.getLocalizedMessage());
                }
            }
        }
        if (!restraintAtoms.isEmpty()) {
            Atom[] ats = restraintAtoms.toArray(new Atom[restraintAtoms.size()]);
            coordRestraints.add(new CoordRestraint(ats, forceField, true, forceconst));
        } else {
            logger.warning(String.format(" Empty or unparseable restraint argument %s", coordRestraint));
        }
    }
    String[] xyzRestStrings = properties.getStringArray("xyzRestraint");
    // Variables to let sequential and otherwise identical xyzRestraint strings to be globbed into one restraint.
    List<Atom> restraintAts = new ArrayList<>();
    List<double[]> coords = new ArrayList<>();
    double lastForceConst = 0;
    boolean lastUseLam = false;
    for (String xR : xyzRestStrings) {
        String[] toks = xR.split("\\s+");
        int nToks = toks.length;
        if (nToks != 6) {
            logger.info(" XYZ restraint rejected: must have force constant, lambda boolean (true/false), 3 coordinates, and an atom number");
            logger.info(" For a coordinate restraint centered on original coordinates, use restraint or lamrestraint keys.");
            logger.info(String.format(" Rejected restraint string: %s", xR));
        } else {
            try {
                double forceConst = Double.parseDouble(toks[0]);
                boolean useLambda = Boolean.parseBoolean(toks[1]);
                if (forceConst != lastForceConst || useLambda != lastUseLam) {
                    int nAts = restraintAts.size();
                    if (nAts != 0) {
                        double[][] restXYZ = new double[3][nAts];
                        Atom[] atArr = restraintAts.toArray(new Atom[nAts]);
                        for (int i = 0; i < 3; i++) {
                            for (int j = 0; j < nAts; j++) {
                                restXYZ[i][j] = coords.get(j)[i];
                            }
                        }
                        CoordRestraint thePin = new CoordRestraint(atArr, forceField, lastUseLam, lastForceConst);
                        thePin.setCoordinatePin(restXYZ);
                        thePin.setIgnoreHydrogen(false);
                        coordRestraints.add(thePin);
                    }
                    restraintAts = new ArrayList<>();
                    coords = new ArrayList<>();
                    lastForceConst = forceConst;
                    lastUseLam = useLambda;
                }
                double[] atXYZ = new double[3];
                for (int i = 0; i < 3; i++) {
                    atXYZ[i] = Double.parseDouble(toks[i + 2]);
                }
                int atNum = Integer.parseInt(toks[5]) - 1;
                restraintAts.add(molaAtoms[atNum]);
                coords.add(atXYZ);
            } catch (Exception ex) {
                logger.info(String.format(" Exception parsing xyzRestraint %s: %s", xR, ex.toString()));
            }
        }
    }
    int nAts = restraintAts.size();
    if (nAts != 0) {
        double[][] restXYZ = new double[3][nAts];
        Atom[] atArr = restraintAts.toArray(new Atom[nAts]);
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < nAts; j++) {
                restXYZ[i][j] = coords.get(j)[i];
            }
        }
        CoordRestraint thePin = new CoordRestraint(atArr, forceField, lastUseLam, lastForceConst);
        thePin.setCoordinatePin(restXYZ);
        thePin.setIgnoreHydrogen(false);
        coordRestraints.add(thePin);
    }
    String[] noElStrings = properties.getStringArray("noElectro");
    for (String noE : noElStrings) {
        String[] toks = noE.split("\\s+");
        for (String tok : toks) {
            try {
                int[] noERange = parseAtNumArg("noElectro", tok, nmolaAtoms);
                for (int i = noERange[0]; i <= noERange[1]; i++) {
                    molaAtoms[i].setElectrostatics(false);
                }
                logger.log(Level.INFO, String.format(" Disabled electrostatics " + "for atoms %d-%d", noERange[0] + 1, noERange[1] + 1));
            } catch (IllegalArgumentException ex) {
                boolean atomFound = false;
                for (Atom atom : molaAtoms) {
                    if (atom.getName().equalsIgnoreCase(tok)) {
                        atomFound = true;
                        atom.setElectrostatics(false);
                    }
                }
                if (atomFound) {
                    logger.info(String.format(" Disabled electrostatics for atoms with name %s", tok));
                } else {
                    logger.log(Level.INFO, String.format(" No electrostatics " + "input %s could not be parsed as a numerical " + "range or atom type present in assembly", tok));
                }
            }
        }
    }
}