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

Example 1 with Crystal

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

the class CCP4MapFilter method getCrystal.

/**
 * {@inheritDoc}
 */
@Override
public Crystal getCrystal(String fileName, CompositeConfiguration properties) {
    int imapData;
    int spaceGroup = -1;
    double cellA = -1.0;
    double cellB = -1.0;
    double cellC = -1.0;
    double cellAlpha = -1.0;
    double cellBeta = -1.0;
    double cellGamma = -1.0;
    ByteOrder byteOrder = ByteOrder.nativeOrder();
    FileInputStream fileInputStream;
    DataInputStream dataInputStream;
    // first determine byte order of file versus system
    try {
        fileInputStream = new FileInputStream(fileName);
        dataInputStream = new DataInputStream(fileInputStream);
        dataInputStream.skipBytes(212);
        byte[] bytes = new byte[4];
        dataInputStream.read(bytes, 0, 4);
        ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
        imapData = byteBuffer.order(ByteOrder.BIG_ENDIAN).getInt();
        String stampString = Integer.toHexString(imapData);
        switch(stampString.charAt(0)) {
            case '1':
            case '3':
                if (byteOrder.equals(ByteOrder.LITTLE_ENDIAN)) {
                    byteOrder = ByteOrder.BIG_ENDIAN;
                }
                break;
            case '4':
                if (byteOrder.equals(ByteOrder.BIG_ENDIAN)) {
                    byteOrder = ByteOrder.LITTLE_ENDIAN;
                }
                break;
        }
        fileInputStream.close();
    } catch (Exception e) {
        String message = " Fatal exception reading CCP4 map.\n";
        logger.log(Level.SEVERE, message, e);
    }
    try {
        fileInputStream = new FileInputStream(fileName);
        dataInputStream = new DataInputStream(fileInputStream);
        dataInputStream.skipBytes(40);
        byte[] bytes = new byte[80];
        dataInputStream.read(bytes, 0, 80);
        ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
        cellA = byteBuffer.order(byteOrder).getFloat();
        cellB = byteBuffer.order(byteOrder).getFloat();
        cellC = byteBuffer.order(byteOrder).getFloat();
        cellAlpha = byteBuffer.order(byteOrder).getFloat();
        cellBeta = byteBuffer.order(byteOrder).getFloat();
        cellGamma = byteBuffer.order(byteOrder).getFloat();
        for (int i = 0; i < 3; i++) {
            byteBuffer.order(byteOrder).getInt();
        }
        for (int i = 0; i < 3; i++) {
            byteBuffer.order(byteOrder).getFloat();
        }
        spaceGroup = byteBuffer.order(byteOrder).getInt();
        fileInputStream.close();
    } catch (Exception e) {
        String message = " Fatal exception reading CCP4 map.\n";
        logger.log(Level.SEVERE, message, e);
    }
    return new Crystal(cellA, cellB, cellC, cellAlpha, cellBeta, cellGamma, SpaceGroup.spaceGroupNames[spaceGroup - 1]);
}

Example 2 with Crystal

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

the class MTZFilter method getReflectionList.

/**
 * {@inheritDoc}
 */
@Override
public ReflectionList getReflectionList(File mtzFile, CompositeConfiguration properties) {
    ByteOrder byteOrder = ByteOrder.nativeOrder();
    FileInputStream fileInputStream;
    DataInputStream dataInputStream;
    try {
        fileInputStream = new FileInputStream(mtzFile);
        dataInputStream = new DataInputStream(fileInputStream);
        byte[] headerOffset = new byte[4];
        byte[] bytes = new byte[80];
        int offset = 0;
        // Eat "MTZ" title.
        dataInputStream.read(bytes, offset, 4);
        String mtzstr = new String(bytes);
        // Header offset.
        dataInputStream.read(headerOffset, offset, 4);
        // Machine stamp.
        dataInputStream.read(bytes, offset, 4);
        ByteBuffer byteBuffer = ByteBuffer.wrap(bytes);
        int stamp = byteBuffer.order(ByteOrder.BIG_ENDIAN).getInt();
        String stampstr = Integer.toHexString(stamp);
        switch(stampstr.charAt(0)) {
            case '1':
            case '3':
                if (byteOrder.equals(ByteOrder.LITTLE_ENDIAN)) {
                    byteOrder = ByteOrder.BIG_ENDIAN;
                }
                break;
            case '4':
                if (byteOrder.equals(ByteOrder.BIG_ENDIAN)) {
                    byteOrder = ByteOrder.LITTLE_ENDIAN;
                }
                break;
        }
        byteBuffer = ByteBuffer.wrap(headerOffset);
        int headerOffsetI = byteBuffer.order(byteOrder).getInt();
        // skip to header and parse
        dataInputStream.skipBytes((headerOffsetI - 4) * 4);
        for (Boolean parsing = true; parsing; dataInputStream.read(bytes, offset, 80)) {
            mtzstr = new String(bytes);
            parsing = parseHeader(mtzstr);
        }
    } catch (EOFException e) {
        String message = " MTZ end of file reached.";
        logger.log(Level.WARNING, message, e);
        return null;
    } catch (IOException e) {
        String message = " MTZ IO exception.";
        logger.log(Level.WARNING, message, e);
        return null;
    }
    // column identifiers
    foString = sigFoString = rFreeString = null;
    if (properties != null) {
        foString = properties.getString("fostring", null);
        sigFoString = properties.getString("sigfostring", null);
        rFreeString = properties.getString("rfreestring", null);
    }
    h = k = l = fo = sigFo = rFree = -1;
    fPlus = sigFPlus = fMinus = sigFMinus = rFreePlus = rFreeMinus = -1;
    fc = phiC = -1;
    boolean print = false;
    parseColumns(print);
    parseFcColumns(print);
    if (fo < 0 && fPlus < 0 && sigFo < 0 && sigFPlus < 0 && fc < 0 && phiC < 0) {
        logger.info(" The MTZ header contains insufficient information to generate the reflection list.");
        logger.info(" For non-default column labels set fostring/sigfostring in the properties file.");
        return null;
    }
    Column column;
    if (fo > 0) {
        column = (Column) columns.get(fo);
    } else if (fPlus > 0) {
        column = (Column) columns.get(fPlus);
    } else {
        column = (Column) columns.get(fc);
    }
    Dataset dataSet = (Dataset) dataSets.get(column.id - dsetOffset);
    if (logger.isLoggable(Level.INFO)) {
        StringBuilder sb = new StringBuilder();
        sb.append(format("\n Reading %s\n\n", mtzFile.getName()));
        sb.append(format(" Setting up reflection list based on MTZ file.\n"));
        sb.append(format("  Space group number: %d (name: %s)\n", spaceGroupNum, SpaceGroup.spaceGroupNames[spaceGroupNum - 1]));
        sb.append(format("  Resolution:         %8.3f\n", 0.999999 * resHigh));
        sb.append(format("  Cell:               %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n", dataSet.cell[0], dataSet.cell[1], dataSet.cell[2], dataSet.cell[3], dataSet.cell[4], dataSet.cell[5]));
        logger.info(sb.toString());
    }
    Crystal crystal = new Crystal(dataSet.cell[0], dataSet.cell[1], dataSet.cell[2], dataSet.cell[3], dataSet.cell[4], dataSet.cell[5], SpaceGroup.spaceGroupNames[spaceGroupNum - 1]);
    double sampling = 0.6;
    if (properties != null) {
        sampling = properties.getDouble("sampling", 0.6);
    }
    Resolution resolution = new Resolution(0.999999 * resHigh, sampling);
    return new ReflectionList(crystal, resolution, properties);
}

Example 3 with Crystal

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

the class CrystalReciprocalSpaceTest method test1N7SPermanent.

/**
 * Test of permanent method, of class CrystalReciprocalSpace.
 */
@Test
public void test1N7SPermanent() {
    String filename = "ffx/xray/structures/1N7S.pdb";
    int index = filename.lastIndexOf(".");
    String name = filename.substring(0, index);
    // load the structure
    ClassLoader cl = this.getClass().getClassLoader();
    File structure = new File(cl.getResource(filename).getPath());
    PotentialsUtils potutil = new PotentialsUtils();
    MolecularAssembly mola = potutil.open(structure);
    CompositeConfiguration properties = mola.getProperties();
    Crystal crystal = new Crystal(39.767, 51.750, 132.938, 90.00, 90.00, 90.00, "P212121");
    Resolution resolution = new Resolution(1.45);
    ReflectionList reflectionList = new ReflectionList(crystal, resolution);
    DiffractionRefinementData refinementData = new DiffractionRefinementData(properties, reflectionList);
    mola.finalize(true, mola.getForceField());
    ForceFieldEnergy energy = mola.getPotentialEnergy();
    List<Atom> atomList = mola.getAtomList();
    Atom[] atomArray = atomList.toArray(new Atom[atomList.size()]);
    // set up FFT and run it
    ParallelTeam parallelTeam = new ParallelTeam();
    CrystalReciprocalSpace crs = new CrystalReciprocalSpace(reflectionList, atomArray, parallelTeam, parallelTeam);
    crs.computeAtomicDensity(refinementData.fc);
    // tests
    ComplexNumber b = new ComplexNumber(-828.584, -922.704);
    HKL hkl = reflectionList.getHKL(1, 1, 4);
    ComplexNumber a = refinementData.getFc(hkl.index());
    System.out.println("1 1 4: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
    assertEquals("1 1 4 reflection should be correct", -753.4722104328416, a.re(), 0.0001);
    assertEquals("1 1 4 reflection should be correct", -1012.1341308707799, a.im(), 0.0001);
    b.re(-70.4582);
    b.im(-486.142);
    hkl = reflectionList.getHKL(2, 1, 10);
    a = refinementData.getFc(hkl.index());
    System.out.println("2 1 10: " + a.toString() + " | " + b.toString() + " | " + a.divides(b).toString());
    assertEquals("2 1 10 reflection should be correct", -69.39660884054359, a.re(), 0.0001);
    assertEquals("2 1 10 reflection should be correct", -412.0147625765328, a.im(), 0.0001);
}

Example 4 with Crystal

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

the class MolecularDynamics method run.

/**
 * {@inheritDoc}
 */
@Override
public void run() {
    done = false;
    terminate = false;
    /**
     * Set the target temperature.
     */
    thermostat.setTargetTemperature(targetTemperature);
    thermostat.setQuiet(quiet);
    if (integrator instanceof Stochastic) {
        Stochastic stochastic = (Stochastic) integrator;
        stochastic.setTemperature(targetTemperature);
    }
    /**
     * Set the step size.
     */
    integrator.setTimeStep(dt);
    if (!initialized) {
        /**
         * Initialize from a restart file.
         */
        if (loadRestart) {
            Crystal crystal = molecularAssembly.getCrystal();
            if (!dynFilter.readDYN(restartFile, crystal, x, v, a, aPrevious)) {
                String message = " Could not load the restart file - dynamics terminated.";
                logger.log(Level.WARNING, message);
                done = true;
                return;
            } else {
                molecularAssembly.getPotentialEnergy().setCrystal(crystal);
            }
        } else {
            /**
             * Initialize from using current atomic coordinates.
             */
            potential.getCoordinates(x);
            /**
             * Initialize atomic velocities from a Maxwell-Boltzmann
             * distribution or set to 0.
             */
            if (initVelocities) {
                thermostat.maxwell(targetTemperature);
            } else {
                fill(v, 0.0);
            }
        }
    } else {
        /**
         * If MD has already been run (ie. Annealing or RepEx), then
         * initialize velocities if requested.
         */
        if (initVelocities) {
            thermostat.maxwell(targetTemperature);
        }
    }
    /**
     * Compute the current potential energy.
     */
    potential.setScaling(null);
    try {
        currentPotentialEnergy = potential.energyAndGradient(x, grad);
    } catch (EnergyException ex) {
        writeStoredSnapshots();
        throw ex;
    }
    /**
     * Initialize current and previous accelerations, unless they were just
     * loaded from a restart file.
     */
    if (!loadRestart || initialized) {
        for (int i = 0; i < numberOfVariables; i++) {
            a[i] = -Thermostat.convert * grad[i] / mass[i];
        }
        if (aPrevious != null) {
            arraycopy(a, 0, aPrevious, 0, numberOfVariables);
        }
    }
    /**
     * Compute the current kinetic energy.
     */
    thermostat.kineticEnergy();
    currentKineticEnergy = thermostat.getKineticEnergy();
    currentTemperature = thermostat.getCurrentTemperature();
    currentTotalEnergy = currentKineticEnergy + currentPotentialEnergy;
    initialized = true;
    logger.info(format("\n  %8s %12s %12s %12s %8s %8s", "Time", "Kinetic", "Potential", "Total", "Temp", "CPU"));
    logger.info(format("  %8s %12s %12s %12s %8s %8s", "psec", "kcal/mol", "kcal/mol", "kcal/mol", "K", "sec"));
    logger.info(format("  %8s %12.4f %12.4f %12.4f %8.2f", "", currentKineticEnergy, currentPotentialEnergy, currentTotalEnergy, currentTemperature));
    /**
     * Store the initialized state.
     */
    storeState();
    /**
     * Integrate Newton's equations of motion for the requested number of
     * steps, unless early termination is requested.
     */
    time = System.nanoTime();
    for (int step = 1; step <= nSteps; step++) {
        /* Notify MonteCarlo handlers such as PhMD or rotamer drivers. */
        if (monteCarloListener != null && mcNotification == MonteCarloNotification.EACH_STEP) {
            monteCarloListener.mcUpdate(thermostat.getCurrentTemperature());
        }
        /**
         * Do the half-step thermostat operation.
         */
        thermostat.halfStep(dt);
        /**
         * Do the half-step integration operation.
         */
        integrator.preForce(potential);
        double priorPE = currentPotentialEnergy;
        /**
         * Compute the potential energy and gradients.
         */
        try {
            currentPotentialEnergy = potential.energyAndGradient(x, grad);
        } catch (EnergyException ex) {
            writeStoredSnapshots();
            throw ex;
        }
        detectAtypicalEnergy(priorPE, defaultDeltaPEThresh);
        /**
         * Add the potential energy of the slow degrees of freedom.
         */
        if (integrator instanceof Respa) {
            Respa r = (Respa) integrator;
            currentPotentialEnergy += r.getHalfStepEnergy();
        }
        /**
         * Do the full-step integration operation.
         */
        integrator.postForce(grad);
        /**
         * Compute the full-step kinetic energy.
         */
        thermostat.kineticEnergy();
        /**
         * Do the full-step thermostat operation.
         */
        thermostat.fullStep(dt);
        /**
         * Recompute the kinetic energy after the full-step thermostat
         * operation.
         */
        thermostat.kineticEnergy();
        /**
         * Remove center of mass motion ever ~100 steps.
         */
        if (thermostat.getRemoveCenterOfMassMotion() && step % removeCOMMotionFrequency == 0) {
            thermostat.centerOfMassMotion(true, false);
        }
        /**
         * Collect current kinetic energy, temperature, and total energy.
         */
        currentKineticEnergy = thermostat.getKineticEnergy();
        currentTemperature = thermostat.getCurrentTemperature();
        currentTotalEnergy = currentKineticEnergy + currentPotentialEnergy;
        /**
         * Update atomic velocity, acceleration and previous acceleration.
         */
        potential.setVelocity(v);
        potential.setAcceleration(a);
        potential.setPreviousAcceleration(aPrevious);
        /**
         * Update extended system variables if present.
         */
        if (esvSystem != null) {
            esvSystem.propagateESVs(currentTemperature, dt, step * dt);
        }
        /**
         * Log the current state every printFrequency steps.
         */
        totalSimTime += dt;
        if (step % printFrequency == 0) {
            /**
             * Original print statement
             */
            time = System.nanoTime() - time;
            mdTime = time;
            logger.info(format(" %7.3e %12.4f %12.4f %12.4f %8.2f %8.3f", totalSimTime, currentKineticEnergy, currentPotentialEnergy, currentTotalEnergy, currentTemperature, time * NS2SEC));
            time = System.nanoTime();
            // Shirts et al. logging info
            Crystal crystal = molecularAssembly.getCrystal();
            double volume = crystal.getUnitCell().volume;
        // Shirts et al. print statement
        // time = System.nanoTime() - time;
        // logger.info(format("Shirts %7.3e %12.4f %12.4f %12.4f %12.4f %8.2f %8.3f",
        // totalSimTime, currentKineticEnergy, currentPotentialEnergy,
        // currentTotalEnergy, volume, currentTemperature, time * NS2SEC));
        // time = System.nanoTime();
        }
        if (step % printEsvFrequency == 0 && esvSystem != null) {
            logger.info(format(" %7.3e %s", totalSimTime, esvSystem.getLambdaList()));
        }
        /**
         * Write out snapshots in selected format every
         * saveSnapshotFrequency steps.
         */
        if (saveSnapshotFrequency > 0 && step % saveSnapshotFrequency == 0) {
            for (AssemblyInfo ai : assemblies) {
                if (ai.archiveFile != null && !saveSnapshotAsPDB) {
                    if (ai.xyzFilter.writeFile(ai.archiveFile, true)) {
                        logger.info(String.format(" Appended snap shot to %s", ai.archiveFile.getName()));
                    } else {
                        logger.warning(String.format(" Appending snap shot to %s failed", ai.archiveFile.getName()));
                    }
                } else if (saveSnapshotAsPDB) {
                    if (ai.pdbFilter.writeFile(ai.pdbFile, false)) {
                        logger.info(String.format(" Wrote PDB file to %s", ai.pdbFile.getName()));
                    } else {
                        logger.warning(String.format(" Writing PDB file to %s failed.", ai.pdbFile.getName()));
                    }
                }
            }
        }
        /**
         * Write out restart files every saveRestartFileFrequency steps.
         */
        if (saveRestartFileFrequency > 0 && step % saveRestartFileFrequency == 0) {
            if (dynFilter.writeDYN(restartFile, molecularAssembly.getCrystal(), x, v, a, aPrevious)) {
                logger.info(String.format(" Wrote dynamics restart file to " + restartFile.getName()));
            } else {
                logger.info(String.format(" Writing dynamics restart file to " + restartFile.getName() + " failed"));
            }
        }
        /**
         * Notify the algorithmListener.
         */
        if (algorithmListener != null && step % printFrequency == 0) {
            // algorithmListener.algorithmUpdate(molecularAssembly);
            for (AssemblyInfo assembly : assemblies) {
                // Probably unwise to parallelize this, so that it doesn't
                // hit the GUI with parallel updates.
                algorithmListener.algorithmUpdate(assembly.getAssembly());
            }
        }
        /**
         * Check for a termination request.
         */
        if (terminate) {
            logger.info(String.format("\n Terminating after %8d time steps\n", step));
            break;
        }
    }
    /**
     * Log normal completion.
     */
    if (!terminate) {
        logger.info(String.format(" Completed %8d time steps\n", nSteps));
    }
    /**
     * Reset the done and terminate flags.
     */
    done = true;
    terminate = false;
    if (monteCarloListener != null && mcNotification == MonteCarloNotification.AFTER_DYNAMICS) {
        monteCarloListener.mcUpdate(thermostat.getCurrentTemperature());
    }
}

Example 5 with Crystal

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

the class MolecularDynamicsOpenMM method init.

@Override
public void init(int numSteps, double timeStep, double printInterval, double saveInterval, String fileType, double restartFrequency, double temperature, boolean initVelocities, File dyn) {
    this.targetTemperature = temperature;
    this.dt = timeStep;
    this.printFrequency = (int) printInterval;
    this.restartFile = dyn;
    this.initVelocities = initVelocities;
    switch(thermostatType) {
        case BUSSI:
        case BERENDSEN:
            logger.info(String.format(" Replacing thermostat %s with OpenMM's Andersen thermostat", thermostatType));
            forceFieldEnergyOpenMM.addAndersenThermostat(temperature);
            break;
        default:
    }
    /**
     * Convert the print interval to a print frequency.
     */
    printFrequency = 100;
    // Time step is in fsec, so convert printInterval to fsec.
    printInterval *= 1000;
    if (printInterval >= dt) {
        printFrequency = (int) (printInterval / dt);
    }
    /**
     * Convert save interval to a save frequency.
     */
    saveSnapshotFrequency = 1000;
    if (saveInterval >= this.dt) {
        saveSnapshotFrequency = (int) (saveInterval / this.dt);
    }
    done = false;
    assemblyInfo();
    String firstFileName = FilenameUtils.removeExtension(molecularAssembly.getFile().getAbsolutePath());
    if (dyn == null) {
        restartFile = new File(firstFileName + ".dyn");
        loadRestart = false;
    } else {
        restartFile = dyn;
        loadRestart = true;
    }
    if (dynFilter == null) {
        dynFilter = new DYNFilter(molecularAssembly.getName());
    }
    dof = forceFieldEnergyOpenMM.calculateDegreesOfFreedom();
    if (!initialized) {
        if (loadRestart) {
            Crystal crystal = molecularAssembly.getCrystal();
            // possibly add check to see if OpenMM supports this space group.
            if (!dynFilter.readDYN(restartFile, crystal, x, v, a, aPrevious)) {
                String message = " Could not load the restart file - dynamics terminated.";
                logger.log(Level.WARNING, message);
                done = true;
            } else {
                forceFieldEnergyOpenMM.setCrystal(crystal);
                // Load positions into the main FFX data structure, move into primary unit cell, then load to OpenMM.
                Atom[] atoms = molecularAssembly.getAtomArray();
                double[] xyz = new double[3];
                double[] vel = new double[3];
                double[] acc = new double[3];
                double[] accPrev = new double[3];
                for (int i = 0; i < atoms.length; i++) {
                    Atom atom = atoms[i];
                    int i3 = i * 3;
                    for (int j = 0; j < 3; j++) {
                        xyz[j] = x[i3 + j];
                        vel[j] = v[i3 + j];
                        acc[j] = a[i3 + j];
                        accPrev[j] = aPrevious[i3 + j];
                    }
                    atom.setXYZ(xyz);
                    atom.setVelocity(vel);
                    atom.setAcceleration(acc);
                    atom.setPreviousAcceleration(accPrev);
                }
                molecularAssembly.moveAllIntoUnitCell();
            }
        } else {
            if (initVelocities) {
                getThermostat().setQuiet(true);
                getThermostat().maxwell(targetTemperature);
                forceFieldEnergyOpenMM.setOpenMMVelocities(v, numberOfVariables);
                Atom[] atoms = molecularAssembly.getAtomArray();
                double[] vel = new double[3];
                for (int i = 0; i < atoms.length; i++) {
                    Atom atom = atoms[i];
                    int i3 = i * 3;
                    for (int j = 0; j < 3; j++) {
                        vel[j] = v[i3 + j];
                    }
                    atom.setVelocity(vel);
                }
            }
        }
        setOpenMMState();
    }
    saveSnapshotAsPDB = true;
    if (fileType.equalsIgnoreCase("XYZ")) {
        saveSnapshotAsPDB = false;
        logger.info("Snapshots will be saved as XYZ");
    } else if (!fileType.equalsIgnoreCase("PDB")) {
        logger.warning("Snapshot file type unrecognized; saving snaphshots as PDB.\n");
    }
    Atom[] atoms = molecularAssembly.getAtomArray();
    natoms = atoms.length;
    int i = 0;
    running = false;
    // logger.info(" Calling OpenMM Update from MD Init.");
    updateFromOpenMM(i, running);
}