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

Example 31 with Crystal

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

the class ParticleMeshEwaldCart method ligandElec.

/**
 * 3.) Aperiodic ligand electrostatics.
 *
 * A.) Real space with an Ewald coefficient of 0.0 (no reciprocal space).
 *
 * B.) Polarization scaled as in Step 2 by (1 - lambda).
 */
private double ligandElec() {
    for (int i = 0; i < nAtoms; i++) {
        use[i] = atoms[i].applyLambda();
    }
    /**
     * Scale the permanent vacuum electrostatics. The softcore alpha is not
     * necessary (nothing in vacuum to collide with).
     */
    doPermanentRealSpace = true;
    permanentScale = 1.0 - lPowPerm;
    dEdLSign = -1.0;
    double lAlphaBack = lAlpha;
    double dlAlphaBack = dlAlpha;
    double d2lAlphaBack = d2lAlpha;
    lAlpha = 0.0;
    dlAlpha = 0.0;
    d2lAlpha = 0.0;
    /**
     * If we are past the end of the polarization lambda window, then only
     * the condensed phase is necessary.
     */
    if (lambda <= polLambdaEnd) {
        doPolarization = true;
        polarizationScale = 1.0 - lPowPol;
    } else {
        doPolarization = false;
        polarizationScale = 0.0;
    }
    /**
     * Save the current real space PME parameters.
     */
    double offBack = off;
    double aewaldBack = aewald;
    off = Double.MAX_VALUE;
    aewald = 0.0;
    setEwaldParameters(off, aewald);
    /**
     * Save the current parallelization schedule.
     */
    IntegerSchedule permanentScheduleBack = permanentSchedule;
    IntegerSchedule ewaldScheduleBack = realSpaceSchedule;
    Range[] rangesBack = realSpaceRanges;
    permanentSchedule = vaporPermanentSchedule;
    realSpaceSchedule = vaporEwaldSchedule;
    realSpaceRanges = vacuumRanges;
    /**
     * Use vacuum crystal / vacuum neighborLists.
     */
    Crystal crystalBack = crystal;
    int nSymmBack = nSymm;
    int[][][] listsBack = neighborLists;
    neighborLists = vaporLists;
    crystal = vaporCrystal;
    nSymm = 1;
    /**
     * Turn off GK if in use.
     */
    boolean gkBack = generalizedKirkwoodTerm;
    /**
     * Turn off Pre-conditioned conjugate gradient SCF solver.
     */
    SCFAlgorithm scfBack = scfAlgorithm;
    scfAlgorithm = SCFAlgorithm.SOR;
    if (doLigandGKElec) {
        generalizedKirkwoodTerm = true;
        generalizedKirkwood.setNeighborList(vaporLists);
        generalizedKirkwood.setLambda(lambda);
        generalizedKirkwood.setCutoff(off);
        generalizedKirkwood.setCrystal(vaporCrystal);
        generalizedKirkwood.setLambdaFunction(polarizationScale, dEdLSign * dlPowPol, dEdLSign * d2lPowPol);
    } else {
        generalizedKirkwoodTerm = false;
    }
    double energy = computeEnergy(false);
    /**
     * Revert to the saved parameters.
     */
    off = offBack;
    aewald = aewaldBack;
    setEwaldParameters(off, aewald);
    neighborLists = listsBack;
    crystal = crystalBack;
    nSymm = nSymmBack;
    permanentSchedule = permanentScheduleBack;
    realSpaceSchedule = ewaldScheduleBack;
    realSpaceRanges = rangesBack;
    lAlpha = lAlphaBack;
    dlAlpha = dlAlphaBack;
    d2lAlpha = d2lAlphaBack;
    generalizedKirkwoodTerm = gkBack;
    scfAlgorithm = scfBack;
    fill(use, true);
    return energy;
}

Example 32 with Crystal

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

the class ParticleMeshEwaldCart method initSoftCoreInit.

/**
 * Initialize a boolean array of soft atoms and, if requested, ligand vapor
 * electrostatics.
 */
private void initSoftCoreInit(boolean print) {
    if (initSoftCore) {
        return;
    }
    /**
     * Initialize a boolean array that marks soft atoms.
     */
    StringBuilder sb = new StringBuilder("\n Softcore Atoms:\n");
    int count = 0;
    for (int i = 0; i < nAtoms; i++) {
        Atom ai = atoms[i];
        if (ai.applyLambda()) {
            isSoft[i] = true;
            if (print) {
                sb.append(ai.toString()).append("\n");
            }
            count++;
        }
    }
    if (count > 0 && print) {
        logger.info(sb.toString());
    }
    /**
     * Initialize boundary conditions, an n^2 neighbor list and parallel
     * scheduling for ligand vapor electrostatics.
     */
    if (doLigandVaporElec) {
        double maxr = 10.0;
        for (int i = 0; i < nAtoms; i++) {
            Atom ai = atoms[i];
            if (ai.applyLambda()) {
                /**
                 * Determine ligand size.
                 */
                for (int j = i + 1; j < nAtoms; j++) {
                    Atom aj = atoms[j];
                    if (aj.applyLambda()) {
                        double dx = ai.getX() - aj.getX();
                        double dy = ai.getY() - aj.getY();
                        double dz = ai.getZ() - aj.getZ();
                        double r = sqrt(dx * dx + dy * dy + dz * dz);
                        maxr = max(r, maxr);
                    }
                }
            }
        }
        double vacuumOff = 2 * maxr;
        vaporCrystal = new Crystal(3 * vacuumOff, 3 * vacuumOff, 3 * vacuumOff, 90.0, 90.0, 90.0, "P1");
        vaporCrystal.setAperiodic(true);
        NeighborList vacuumNeighborList = new NeighborList(null, vaporCrystal, atoms, vacuumOff, 2.0, parallelTeam);
        vacuumNeighborList.setIntermolecular(false, molecule);
        vaporLists = new int[1][nAtoms][];
        double[][] coords = new double[1][nAtoms * 3];
        for (int i = 0; i < nAtoms; i++) {
            coords[0][i * 3] = atoms[i].getX();
            coords[0][i * 3 + 1] = atoms[i].getY();
            coords[0][i * 3 + 2] = atoms[i].getZ();
        }
        vacuumNeighborList.buildList(coords, vaporLists, isSoft, true, true);
        vaporPermanentSchedule = vacuumNeighborList.getPairwiseSchedule();
        vaporEwaldSchedule = vaporPermanentSchedule;
        vacuumRanges = new Range[maxThreads];
        vacuumNeighborList.setDisableUpdates(forceField.getBoolean(ForceField.ForceFieldBoolean.DISABLE_NEIGHBOR_UPDATES, false));
    } else {
        vaporCrystal = null;
        vaporLists = null;
        vaporPermanentSchedule = null;
        vaporEwaldSchedule = null;
        vacuumRanges = null;
    }
    /**
     * Set this flag to true to avoid re-initialization.
     */
    initSoftCore = true;
}

Example 33 with Crystal

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

the class ParticleMeshEwaldQI method ligandElec.

/**
 * 3.) Ligand in vapor A.) Real space with an Ewald coefficient of 0.0 (no
 * reciprocal space). B.) Polarization scaled as in Step 2 by (1 - lambda).
 */
private double ligandElec() {
    for (int i = 0; i < nAtoms; i++) {
        use[i] = atoms[i].applyLambda();
    }
    /**
     * Scale the permanent vacuum electrostatics. The softcore alpha is not
     * necessary (nothing in vacuum to collide with).
     */
    doPermanentRealSpace = true;
    permanentScale = 1.0 - lPowPerm;
    dEdLSign = -1.0;
    double lAlphaBack = lAlpha;
    double dlAlphaBack = dlAlpha;
    double d2lAlphaBack = d2lAlpha;
    lAlpha = 0.0;
    dlAlpha = 0.0;
    d2lAlpha = 0.0;
    /**
     * If we are past the end of the polarization lambda window, then only
     * the condensed phase is necessary.
     */
    if (lambda <= polLambdaEnd) {
        doPolarization = true;
        polarizationScale = 1.0 - lPowPol;
    } else {
        doPolarization = false;
        polarizationScale = 0.0;
    }
    /**
     * Save the current real space PME parameters.
     */
    double offBack = off;
    double aewaldBack = aewald;
    off = Double.MAX_VALUE;
    aewald = 0.0;
    setEwaldParameters(off, aewald);
    /**
     * Save the current parallelization schedule.
     */
    IntegerSchedule permanentScheduleBack = permanentSchedule;
    IntegerSchedule ewaldScheduleBack = realSpaceSchedule;
    Range[] rangesBack = realSpaceRanges;
    permanentSchedule = vaporPermanentSchedule;
    realSpaceSchedule = vaporEwaldSchedule;
    realSpaceRanges = vacuumRanges;
    /**
     * Use vacuum crystal / vacuum neighborLists.
     */
    Crystal crystalBack = crystal;
    int nSymmBack = nSymm;
    int[][][] listsBack = neighborLists;
    neighborLists = vaporLists;
    crystal = vaporCrystal;
    nSymm = 1;
    for (LambdaFactors lf : lambdaFactors) {
        lf.setFactors();
    }
    /**
     * Turn off GK if it is in use, unless it's being used as the decoupling
     * target. If so, set its parameters and lambda (derivative) factors.
     */
    boolean gkBack = generalizedKirkwoodTerm;
    if (doLigandGKElec) {
        generalizedKirkwoodTerm = true;
        generalizedKirkwood.setNeighborList(vaporLists);
        generalizedKirkwood.setLambda(lambda);
        generalizedKirkwood.setCutoff(off);
        generalizedKirkwood.setCrystal(vaporCrystal);
        // TODO Decide whether to send LambdaFactors into generalizedKirkwood.
        generalizedKirkwood.setLambdaFunction(polarizationScale, dEdLSign * dlPowPol, dEdLSign * d2lPowPol);
    } else {
        generalizedKirkwoodTerm = false;
    }
    double energy = computeEnergy(false);
    /**
     * Revert to the saved parameters.
     */
    off = offBack;
    aewald = aewaldBack;
    setEwaldParameters(off, aewald);
    neighborLists = listsBack;
    crystal = crystalBack;
    nSymm = nSymmBack;
    permanentSchedule = permanentScheduleBack;
    realSpaceSchedule = ewaldScheduleBack;
    realSpaceRanges = rangesBack;
    lAlpha = lAlphaBack;
    dlAlpha = dlAlphaBack;
    d2lAlpha = d2lAlphaBack;
    generalizedKirkwoodTerm = gkBack;
    for (LambdaFactors lf : lambdaFactors) {
        lf.setFactors();
    }
    fill(use, true);
    return energy;
}

Example 34 with Crystal

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

the class DYNFilter method writeDYN.

/**
 * <p>
 * writeDYN</p>
 *
 * @param dynFile a {@link java.io.File} object.
 * @param unitCell a {@link ffx.crystal.Crystal} object.
 * @param x an array of double.
 * @param v an array of double.
 * @param a an array of double.
 * @param ap an array of double.
 * @return a boolean.
 */
public boolean writeDYN(File dynFile, Crystal crystal, double[] x, double[] v, double[] a, double[] ap) {
    FileWriter fw = null;
    BufferedWriter bw = null;
    try {
        fw = new FileWriter(dynFile);
        bw = new BufferedWriter(fw);
        bw.write(" Number of Atoms and Title :\n");
        assert (x.length % 3 == 0);
        int numberOfAtoms = x.length / 3;
        String output = format("%7d  %s\n", numberOfAtoms, label);
        bw.write(output);
        bw.write(" Periodic Box Dimensions :\n");
        Crystal unitCell = crystal.getUnitCell();
        bw.write(format("%26.16E%26.16E%26.16E\n", unitCell.a, unitCell.b, unitCell.c));
        bw.write(format("%26.16E%26.16E%26.16E\n", unitCell.alpha, unitCell.beta, unitCell.gamma));
        bw.write(" Current Atomic Positions :\n");
        for (int i = 0; i < numberOfAtoms; i++) {
            int k = i * 3;
            bw.write(format("%26.16E%26.16E%26.16E\n", x[k], x[k + 1], x[k + 2]));
        }
        bw.write(" Current Atomic Velocities :\n");
        for (int i = 0; i < numberOfAtoms; i++) {
            int k = i * 3;
            bw.write(format("%26.16E%26.16E%26.16E\n", v[k], v[k + 1], v[k + 2]));
        }
        bw.write(" Current Atomic Accelerations :\n");
        for (int i = 0; i < numberOfAtoms; i++) {
            int k = i * 3;
            bw.write(format("%26.16E%26.16E%26.16E\n", a[k], a[k + 1], a[k + 2]));
        }
        bw.write(" Previous Atomic Accelerations :\n");
        for (int i = 0; i < numberOfAtoms; i++) {
            int k = i * 3;
            bw.write(format("%26.16E%26.16E%26.16E\n", ap[k], ap[k + 1], ap[k + 2]));
        }
    } catch (IOException e) {
        String message = " Exception writing dynamic restart file " + dynFile;
        logger.log(Level.SEVERE, message, e);
        return false;
    } finally {
        try {
            bw.close();
            fw.close();
            return true;
        } catch (IOException e) {
            String message = " Exception closing dynamic restart file " + dynFile;
            logger.log(Level.WARNING, message, e);
            return false;
        }
    }
}

Example 35 with Crystal

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

the class CNSFilter method getReflectionList.

/**
 * {@inheritDoc}
 */
@Override
public ReflectionList getReflectionList(File cnsFile, CompositeConfiguration properties) {
    try {
        BufferedReader br = new BufferedReader(new FileReader(cnsFile));
        String string;
        while ((string = br.readLine()) != null) {
            String[] strArray = string.split("\\s+");
            if (strArray[0].equalsIgnoreCase("{")) {
                if (strArray[1].toLowerCase().startsWith("sg=")) {
                    spaceGroupName = strArray[1].substring(3);
                    cell[0] = parseDouble(strArray[2].substring(2));
                    cell[1] = parseDouble(strArray[3].substring(2));
                    cell[2] = parseDouble(strArray[4].substring(2));
                    cell[3] = parseDouble(strArray[5].substring(6));
                    cell[4] = parseDouble(strArray[6].substring(5));
                    cell[5] = parseDouble(strArray[7].substring(6));
                }
            } else if (strArray[0].equalsIgnoreCase("CRYST1")) {
                cell[0] = parseDouble(strArray[1]);
                cell[1] = parseDouble(strArray[2]);
                cell[2] = parseDouble(strArray[3]);
                cell[3] = parseDouble(strArray[4]);
                cell[4] = parseDouble(strArray[5]);
                cell[5] = parseDouble(strArray[6]);
                spaceGroupName = SpaceGroup.pdb2ShortName(string.substring(55, 65));
            } else if (strArray[0].toLowerCase().startsWith("inde")) {
                break;
            }
        }
        br.close();
    } catch (IOException e) {
        String message = " CNS IO Exception.";
        logger.log(Level.WARNING, message, e);
        return null;
    }
    Resolution resolution = null;
    if (properties != null) {
        resolution = Resolution.checkProperties(properties);
        resHigh = resolution.resolution;
    }
    if (spaceGroupName != null) {
        spaceGroupNum = SpaceGroup.spaceGroupNumber(spaceGroupName);
    }
    if (spaceGroupNum < 0 || cell[0] < 0 || resolution == null) {
        logger.info(" The CNS file contains insufficient information to generate the reflection list.");
        return null;
    }
    if (logger.isLoggable(Level.INFO)) {
        StringBuilder sb = new StringBuilder();
        sb.append(format("\n Opening %s\n", cnsFile.getName()));
        sb.append(format(" Setting up Reflection List based on CNS:\n"));
        sb.append(format("  Spacegroup #: %d (name: %s)\n", spaceGroupNum, SpaceGroup.spaceGroupNames[spaceGroupNum - 1]));
        sb.append(format("  Resolution:   %8.3f\n", resHigh));
        sb.append(format("  Cell:         %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n", cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]));
        logger.info(sb.toString());
    }
    Crystal crystal = new Crystal(cell[0], cell[1], cell[2], cell[3], cell[4], cell[5], SpaceGroup.spaceGroupNames[spaceGroupNum - 1]);
    ReflectionList reflectionlist = new ReflectionList(crystal, resolution, properties);
    return reflectionlist;
}