Examples with Crystal ffx.crystal.Crystal used on opensource projects

Example 6 with Crystal

use of ffx.crystal.Crystal in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addAmoebaVDWForce.

private void addAmoebaVDWForce() {
    VanDerWaals vdW = super.getVdwNode();
    if (vdW == null) {
        return;
    }
    amoebaVDWForce = OpenMM_AmoebaVdwForce_create();
    OpenMM_System_addForce(system, amoebaVDWForce);
    OpenMM_Force_setForceGroup(amoebaVDWForce, 1);
    VanDerWaalsForm vdwForm = vdW.getVDWForm();
    NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
    Crystal crystal = super.getCrystal();
    double radScale = 1.0;
    if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
        radScale = 0.5;
    }
    /**
     * Note that the API says it wants a SIGMA value.
     */
    if (vdwForm.radiusType == VanDerWaalsForm.RADIUS_TYPE.R_MIN) {
    // radScale *= 1.122462048309372981;
    }
    int[] ired = vdW.getReductionIndex();
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        VDWType vdwType = atom.getVDWType();
        OpenMM_AmoebaVdwForce_addParticle(amoebaVDWForce, ired[i], OpenMM_NmPerAngstrom * vdwType.radius * radScale, OpenMM_KJPerKcal * vdwType.wellDepth, vdwType.reductionFactor);
    }
    // OpenMM_AmoebaVdwForce_setSigmaCombiningRule(amoebaVdwForce, toPropertyForm(vdwForm.radiusRule.name()));
    // OpenMM_AmoebaVdwForce_setEpsilonCombiningRule(amoebaVdwForce, toPropertyForm(vdwForm.epsilonRule.name()));
    OpenMM_AmoebaVdwForce_setCutoffDistance(amoebaVDWForce, nonbondedCutoff.off * OpenMM_NmPerAngstrom);
    OpenMM_AmoebaVdwForce_setUseDispersionCorrection(amoebaVDWForce, OpenMM_Boolean.OpenMM_False);
    if (crystal.aperiodic()) {
        OpenMM_AmoebaVdwForce_setNonbondedMethod(amoebaVDWForce, OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_NoCutoff);
    } else {
        OpenMM_AmoebaVdwForce_setNonbondedMethod(amoebaVDWForce, OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_CutoffPeriodic);
    }
    /**
     * Create exclusion lists.
     */
    PointerByReference exclusions = OpenMM_IntArray_create(0);
    double[] mask = new double[nAtoms];
    Arrays.fill(mask, 1.0);
    for (int i = 0; i < nAtoms; i++) {
        OpenMM_IntArray_append(exclusions, i);
        vdW.applyMask(mask, i);
        for (int j = 0; j < nAtoms; j++) {
            if (mask[j] == 0.0) {
                OpenMM_IntArray_append(exclusions, j);
            }
        }
        vdW.removeMask(mask, i);
        OpenMM_AmoebaVdwForce_setParticleExclusions(amoebaVDWForce, i, exclusions);
        OpenMM_IntArray_resize(exclusions, 0);
    }
    OpenMM_IntArray_destroy(exclusions);
    logger.log(Level.INFO, " Added van der Waals force.");
}

Example 7 with Crystal

use of ffx.crystal.Crystal in project ffx by mjschnie.

the class ForceFieldEnergyOpenMM method addCustomNonbondedSoftcoreForce.

/**
 * 1.) Handle interactions between non-alchemical atoms with our default
 * OpenMM NonBondedForce. Note that alchemical atoms must have eps=0 to turn
 * them off in this force.
 * <p>
 * 2.) Handle interactions between alchemical atoms and mixed non-alchemical
 * <-> alchemical interactions with an OpenMM CustomNonBondedForce.
 */
private void addCustomNonbondedSoftcoreForce() {
    VanDerWaals vdW = super.getVdwNode();
    if (vdW == null) {
        return;
    }
    /**
     * Only 6-12 LJ with arithmetic mean to define sigma and geometric mean
     * for epsilon is supported.
     */
    VanDerWaalsForm vdwForm = vdW.getVDWForm();
    if (vdwForm.vdwType != LENNARD_JONES || vdwForm.radiusRule != ARITHMETIC || vdwForm.epsilonRule != GEOMETRIC) {
        logger.info(format(" VDW Type:         %s", vdwForm.vdwType));
        logger.info(format(" VDW Radius Rule:  %s", vdwForm.radiusRule));
        logger.info(format(" VDW Epsilon Rule: %s", vdwForm.epsilonRule));
        logger.log(Level.SEVERE, String.format(" Unsuppporterd van der Waals functional form."));
        return;
    }
    // Sterics mixing rules.
    String stericsMixingRules = " epsilon = sqrt(epsilon1*epsilon2);";
    stericsMixingRules += " rmin = 0.5 * (sigma1 + sigma2) * 1.122462048309372981;";
    // Softcore Lennard-Jones, with a form equivalent to that used in FFX VanDerWaals class.
    String stericsEnergyExpression = "(vdw_lambda^beta)*epsilon*x*(x-2.0);";
    // Effective softcore distance for sterics.
    stericsEnergyExpression += " x = 1.0 / (alpha*(1.0-vdw_lambda)^2.0 + (r/rmin)^6.0);";
    // Define energy expression for sterics.
    String energyExpression = stericsEnergyExpression + stericsMixingRules;
    fixedChargeSoftcore = OpenMM_CustomNonbondedForce_create(energyExpression);
    // Get the Alpha and Beta constants from the VanDerWaals instance.
    double alpha = vdW.getAlpha();
    double beta = vdW.getBeta();
    logger.info(format(" Custom non-bonded force with alpha = %8.6f and beta = %8.6f", alpha, beta));
    OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "vdw_lambda", 1.0);
    OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "alpha", alpha);
    OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "beta", beta);
    OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "sigma");
    OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "epsilon");
    /**
     * Add particles.
     */
    PointerByReference alchemicalGroup = OpenMM_IntSet_create();
    PointerByReference nonAlchemicalGroup = OpenMM_IntSet_create();
    DoubleByReference charge = new DoubleByReference();
    DoubleByReference sigma = new DoubleByReference();
    DoubleByReference eps = new DoubleByReference();
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    for (int i = 0; i < nAtoms; i++) {
        Atom atom = atoms[i];
        OpenMM_NonbondedForce_getParticleParameters(fixedChargeNonBondedForce, i, charge, sigma, eps);
        if (atom.applyLambda()) {
            OpenMM_IntSet_insert(alchemicalGroup, i);
            logger.info(format(" Adding alchemical atom %s.", atom));
        } else {
            OpenMM_IntSet_insert(nonAlchemicalGroup, i);
        }
        PointerByReference particleParameters = OpenMM_DoubleArray_create(0);
        OpenMM_DoubleArray_append(particleParameters, sigma.getValue());
        OpenMM_DoubleArray_append(particleParameters, eps.getValue());
        OpenMM_CustomNonbondedForce_addParticle(fixedChargeSoftcore, particleParameters);
        OpenMM_DoubleArray_destroy(particleParameters);
    }
    OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup, alchemicalGroup);
    OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup, nonAlchemicalGroup);
    OpenMM_IntSet_destroy(alchemicalGroup);
    OpenMM_IntSet_destroy(nonAlchemicalGroup);
    Crystal crystal = super.getCrystal();
    if (crystal.aperiodic()) {
        OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore, OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_NoCutoff);
    } else {
        OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore, OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_CutoffPeriodic);
    }
    NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
    double off = nonbondedCutoff.off;
    double cut = nonbondedCutoff.cut;
    OpenMM_CustomNonbondedForce_setCutoffDistance(fixedChargeSoftcore, OpenMM_NmPerAngstrom * off);
    OpenMM_CustomNonbondedForce_setUseSwitchingFunction(fixedChargeSoftcore, OpenMM_True);
    OpenMM_CustomNonbondedForce_setSwitchingDistance(fixedChargeSoftcore, OpenMM_NmPerAngstrom * cut);
    if (cut == off) {
        logger.warning(" OpenMM does not properly handle cutoffs where cut == off!");
        if (cut == Double.MAX_VALUE || cut == Double.POSITIVE_INFINITY) {
            logger.info(" Detected infinite or max-value cutoff; setting cut to 1E+40 for OpenMM.");
            cut = 1E40;
        } else {
            logger.info(String.format(" Detected cut %8.4g == off %8.4g; scaling cut to 0.99 of off for OpenMM.", cut, off));
            cut *= 0.99;
        }
    }
    // Add energy parameter derivative
    OpenMM_CustomNonbondedForce_addEnergyParameterDerivative(fixedChargeSoftcore, "vdw_lambda");
    OpenMM_System_addForce(system, fixedChargeSoftcore);
    logger.log(Level.INFO, String.format(" Added fixed charge softcore sterics force."));
    GeneralizedKirkwood gk = super.getGK();
    if (gk != null) {
        logger.severe(" OpenMM alchemical methods are not supported for GB.");
        addCustomGBForce();
    }
    // Not entirely sure how to initialize this portion
    alchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(stericsEnergyExpression);
    nonAlchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(stericsEnergyExpression);
    // allStericsForce = (alchemicalAlchemicalStericsForce + nonAlchemicalAlchemicalStericsForce);
    // Can be reduced to two lines if I can figure out how to combine the two custom bonded sterics forces
    OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "rmin");
    OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "epsilon");
    OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "vdw_lambda", 1.0);
    OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "alpha", alpha);
    OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "beta", beta);
    OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "rmin");
    OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "epsilon");
    OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "vdw_lambda", 1.0);
    OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "alpha", alpha);
    OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "beta", beta);
    int range = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
    IntByReference atomi = new IntByReference();
    IntByReference atomj = new IntByReference();
    int[][] torsionMask = vdW.getTorsionMask();
    for (int i = 0; i < range; i++) {
        OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, atomi, atomj, charge, sigma, eps);
        OpenMM_CustomNonbondedForce_addExclusion(fixedChargeSoftcore, atomi.getValue(), atomj.getValue());
        int[] maskI = torsionMask[atomi.getValue()];
        int jID = atomj.getValue();
        boolean epsException = false;
        for (int j = 0; j < maskI.length; j++) {
            if (maskI[j] == jID) {
                epsException = true;
                break;
            }
        }
        Atom atom1 = atoms[atomi.getValue()];
        Atom atom2 = atoms[atomj.getValue()];
        boolean bothAlchemical = false;
        boolean oneAlchemical = false;
        if (atom1.applyLambda() && atom2.applyLambda()) {
            bothAlchemical = true;
        } else if ((atom1.applyLambda() && !atom2.applyLambda()) || (!atom1.applyLambda() && atom2.applyLambda())) {
            oneAlchemical = true;
        }
        if (bothAlchemical) {
            if (epsException) {
                PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
                OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
                OpenMM_DoubleArray_append(bondParameters, eps.getValue());
                OpenMM_CustomBondForce_addBond(alchemicalAlchemicalStericsForce, atomi.getValue(), atomj.getValue(), bondParameters);
                OpenMM_DoubleArray_destroy(bondParameters);
            }
        } else if (oneAlchemical) {
            if (epsException) {
                PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
                OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
                OpenMM_DoubleArray_append(bondParameters, eps.getValue());
                OpenMM_CustomBondForce_addBond(nonAlchemicalAlchemicalStericsForce, atomi.getValue(), atomj.getValue(), bondParameters);
                OpenMM_DoubleArray_destroy(bondParameters);
            }
        }
    }
    /**
     *        for (int i = 0; i < range; i++){
     *            OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, atomi, atomj, charge, sigma, eps);
     *
     *            Atom atom1 = atoms[atomi.getValue()];
     *            Atom atom2 = atoms[atomj.getValue()];
     *
     *            if (atom1.applyLambda() || atom2.applyLambda()){
     *                OpenMM_NonbondedForce_setExceptionParameters(fixedChargeNonBondedForce, i, atomi.getValue(), atomj.getValue(), abs(0.0*charge.getValue()), sigma.getValue(), abs(0.0*eps.getValue()));
     *            }
     *        }
     */
    OpenMM_CustomBondForce_addEnergyParameterDerivative(alchemicalAlchemicalStericsForce, "vdw_lambda");
    OpenMM_CustomBondForce_addEnergyParameterDerivative(nonAlchemicalAlchemicalStericsForce, "vdw_lambda");
    OpenMM_System_addForce(system, alchemicalAlchemicalStericsForce);
    OpenMM_System_addForce(system, nonAlchemicalAlchemicalStericsForce);
}

Example 8 with Crystal

use of ffx.crystal.Crystal in project ffx by mjschnie.

the class RotamerOptimization method generateSuperbox.

/**
 * Returns the superbox used to generate the boxes for sliding box. If
 * superbox coordinates manually set, uses them plus the defined buffer.
 * Else, if an aperiodic system, uses maximum and minimum C alpha (or N1/9)
 * coordinates plus superboxBuffer (by default, 8A, longer than a lysine
 * side chain or N1/N9 distance to any other atom). Else, it just uses the
 * ordinary crystal.
 *
 * @param residueList List of residues to incorporate.
 * @return Superbox crystal.
 */
private Crystal generateSuperbox(List<Residue> residueList) {
    double[] maxXYZ = new double[3];
    double[] minXYZ = new double[3];
    Crystal originalCrystal = molecularAssembly.getCrystal();
    if (manualSuperbox) {
        for (int i = 0; i < maxXYZ.length; i++) {
            int ii = 2 * i;
            minXYZ[i] = boxDimensions[ii] - superboxBuffer;
            maxXYZ[i] = boxDimensions[ii + 1] + superboxBuffer;
        }
    } else if (originalCrystal.aperiodic()) {
        if (residueList == null || residueList.isEmpty()) {
            throw new IllegalArgumentException(" Null or empty residue list when generating superbox.");
        }
        Atom initializerAtom = residueList.get(0).getReferenceAtom();
        initializerAtom.getXYZ(minXYZ);
        initializerAtom.getXYZ(maxXYZ);
        for (Residue residue : residueList) {
            Atom refAtom = residue.getReferenceAtom();
            double[] refAtomCoords = new double[3];
            refAtom.getXYZ(refAtomCoords);
            for (int i = 0; i < 3; i++) {
                maxXYZ[i] = (refAtomCoords[i] > maxXYZ[i] ? refAtomCoords[i] : maxXYZ[i]);
                minXYZ[i] = (refAtomCoords[i] < minXYZ[i] ? refAtomCoords[i] : minXYZ[i]);
            }
        }
        for (int i = 0; i < 3; i++) {
            minXYZ[i] -= superboxBuffer;
            maxXYZ[i] += superboxBuffer;
        }
    } else {
        return originalCrystal;
    }
    double newA = maxXYZ[0] - minXYZ[0];
    double newB = maxXYZ[1] - minXYZ[1];
    double newC = maxXYZ[2] - minXYZ[2];
    if (manualSuperbox) {
        logger.info(format(" Manual superbox set over (minX, maxX, minY, " + "maxY, minZ, maxZ): %f, %f, %f, %f, %f, %f", minXYZ[0], maxXYZ[0], minXYZ[1], maxXYZ[1], minXYZ[2], maxXYZ[2]));
        logger.info(format(" Buffer size (included in dimensions): %f\n", superboxBuffer));
    } else {
        // Crystal systems will have already returned.
        logger.info(" System is aperiodic: protein box generated over these coordinates (minX, maxX, minY, maxY, minZ, maxZ):");
        String message = " Aperiodic box dimensions: ";
        for (int i = 0; i < minXYZ.length; i++) {
            message = message.concat(format("%f,%f,", minXYZ[i], maxXYZ[i]));
        }
        message = message.substring(0, message.length() - 1);
        logger.info(message);
        logger.info(format(" Buffer size (included in dimensions): %f\n", superboxBuffer));
    }
    return new Crystal(newA, newB, newC, 90.0, 90.0, 90.0, "P1");
}

Example 9 with Crystal

use of ffx.crystal.Crystal in project ffx by mjschnie.

the class RotamerOptimization method evaluateDistance.

/**
 * Evaluates the pairwise distance between two residues' rotamers under any
 * symmetry operator; does "lazy loading" for the distance matrix.
 *
 * @param i Residue i
 * @param ri Rotamer for i
 * @param j Residue j
 * @param rj Rotamer for j
 * @return Shortest distance
 */
private double evaluateDistance(int i, int ri, int j, int rj) {
    Residue resi = allResiduesArray[i];
    Rotamer[] rotamersI = resi.getRotamers(library);
    Rotamer roti = rotamersI[ri];
    double[][] xi;
    if (roti.equals(resi.getRotamer())) {
        xi = resi.storeCoordinateArray();
    } else {
        ResidueState origI = resi.storeState();
        RotamerLibrary.applyRotamer(resi, roti);
        xi = resi.storeCoordinateArray();
        resi.revertState(origI);
    }
    Residue resj = allResiduesArray[j];
    Rotamer[] rotamersJ = resj.getRotamers(library);
    Rotamer rotj = rotamersJ[rj];
    double[][] xj;
    if (rotj.equals(resj.getRotamer())) {
        xj = resj.storeCoordinateArray();
    } else {
        ResidueState origJ = resj.storeState();
        RotamerLibrary.applyRotamer(resj, rotj);
        xj = resj.storeCoordinateArray();
        resj.revertState(origJ);
    }
    Crystal crystal = molecularAssembly.getCrystal();
    int nSymm = crystal.spaceGroup.getNumberOfSymOps();
    double minDist = Double.MAX_VALUE;
    for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
        SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
        double dist = interResidueDistance(xi, xj, symOp);
        minDist = dist < minDist ? dist : minDist;
    }
    return minDist;
}

Example 10 with Crystal

use of ffx.crystal.Crystal in project ffx by mjschnie.

the class RotamerOptimization method interResidueDistance.

/**
 * Calculates the minimum distance between two sets of coordinates in a
 * given symmetry operator.
 *
 * @param resi Coordinates of i by [atom][xyz]
 * @param resj Coordinates of j by [atom][xyz]
 * @param symOp Symmetry operator to apply
 * @return Minimum distance
 */
private double interResidueDistance(double[][] resi, double[][] resj, SymOp symOp) {
    double dist = Double.MAX_VALUE;
    Crystal crystal = molecularAssembly.getCrystal();
    int ni = resi.length;
    for (int i = 0; i < ni; i++) {
        double[] xi = resi[i];
        int nj = resj.length;
        for (int j = 0; j < nj; j++) {
            double[] xj = resj[j];
            if (symOp != null) {
                crystal.applySymOp(xj, xj, symOp);
            }
            // Generally: compare on square-of-distance, and square root only at return.
            // double r = Math.sqrt(crystal.image(xi[0] - xj[0], xi[1] - xj[1], xi[2] - xj[2]));
            double r = crystal.image(xi[0] - xj[0], xi[1] - xj[1], xi[2] - xj[2]);
            if (r < dist) {
                dist = r;
            }
        }
    }
    return sqrt(dist);
}