use of ffx.crystal.Crystal in project ffx by mjschnie.
the class MolecularAssembly method moveIntoUnitCell.
/**
* Move the center of each listed chemical entity into the primary unit cell.
* @param groups
*/
private void moveIntoUnitCell(List<MSNode> groups) {
Crystal cryst = getCrystal();
if (cryst.aperiodic()) {
return;
}
for (MSNode group : groups) {
double[] com = new double[3];
double[] xyz = new double[3];
double[] translate = new double[3];
List<Atom> atoms = group.getAtomList();
double totMass = 0;
for (Atom atom : atoms) {
double mass = atom.getMass();
totMass += mass;
xyz = atom.getXYZ(xyz);
for (int i = 0; i < 3; i++) {
double diff = xyz[i] - com[i];
diff /= totMass;
com[i] += diff;
}
}
cryst.toPrimaryCell(com, translate);
for (int i = 0; i < 3; i++) {
translate[i] -= com[i];
}
for (Atom atom : atoms) {
atom.move(translate);
}
}
}
use of ffx.crystal.Crystal in project ffx by mjschnie.
the class RotamerOptimization method boxOptimization.
/**
* Breaks down a structure into a number of overlapping boxes for
* optimization.
*
* @return Potential energy of final structure.
*/
private double boxOptimization(ArrayList<Residue> residueList) {
this.usingBoxOptimization = true;
long beginTime = -System.nanoTime();
Residue[] residues = residueList.toArray(new Residue[residueList.size()]);
boolean firstCellSaved = false;
/*
* A new dummy Crystal will be constructed for an aperiodic system. The
* purpose is to avoid using the overly large dummy Crystal used for
* Ewald purposes. Atoms are not and should not be moved into the dummy
* Crystal boundaries; to check if an Atom is inside a cell, use an
* array of coordinates adjusted to be 0 < coordinate < 1.0.
*/
Crystal crystal = generateSuperbox(residueList);
// Cells indexed by x*(YZ divisions) + y*(Z divisions) + z.
// Also initializes cell count if using -bB
int totalCells = getTotalCellCount(crystal);
if (boxStart > totalCells - 1) {
logger.severe(format(" FFX shutting down: Box optimization start is out of range of total boxes: %d > %d", (boxStart + 1), totalCells));
}
if (boxEnd > totalCells - 1) {
boxEnd = totalCells - 1;
logIfMaster(" Final box out of range: reset to last possible box.");
} else if (boxEnd < 0) {
boxEnd = totalCells - 1;
}
BoxOptCell[] cells = loadCells(crystal, residues);
int numCells = cells.length;
logIfMaster(format(" Optimizing boxes %d to %d", (boxStart + 1), (boxEnd + 1)));
for (int i = 0; i < numCells; i++) {
BoxOptCell celli = cells[i];
ArrayList<Residue> residuesList = celli.getResiduesAsList();
int[] cellIndices = celli.getXYZIndex();
logIfMaster(format("\n Iteration %d of the box optimization.", (i + 1)));
logIfMaster(format(" Cell index (linear): %d", (celli.getLinearIndex() + 1)));
logIfMaster(format(" Cell xyz indices: x = %d, y = %d, z = %d", cellIndices[0] + 1, cellIndices[1] + 1, cellIndices[2] + 1));
int nResidues = residuesList.size();
if (nResidues > 0) {
readyForSingles = false;
finishedSelf = false;
readyFor2Body = false;
finished2Body = false;
readyFor3Body = false;
finished3Body = false;
energiesToWrite = Collections.synchronizedList(new ArrayList<String>());
receiveThread = new ReceiveThread(residuesList.toArray(new Residue[1]));
receiveThread.start();
if (master && writeEnergyRestart && printFiles) {
if (energyWriterThread != null) {
int waiting = 0;
while (energyWriterThread.getState() != java.lang.Thread.State.TERMINATED) {
try {
if (waiting++ > 20) {
logger.log(Level.SEVERE, " ReceiveThread/EnergyWriterThread from previous box locked up.");
}
logIfMaster(" Waiting for previous iteration's communication threads to shut down... ");
Thread.sleep(10000);
} catch (InterruptedException ex) {
}
}
}
energyWriterThread = new EnergyWriterThread(receiveThread, i + 1, cellIndices);
energyWriterThread.start();
}
if (loadEnergyRestart) {
boxLoadIndex = i + 1;
boxLoadCellIndices = new int[3];
boxLoadCellIndices[0] = cellIndices[0];
boxLoadCellIndices[1] = cellIndices[1];
boxLoadCellIndices[2] = cellIndices[2];
}
long boxTime = -System.nanoTime();
Residue firstResidue = residuesList.get(0);
Residue lastResidue = residuesList.get(nResidues - 1);
if (firstResidue != lastResidue) {
logIfMaster(format(" Residues %s ... %s", firstResidue.toString(), lastResidue.toString()));
} else {
logIfMaster(format(" Residue %s", firstResidue.toString()));
}
if (revert) {
ResidueState[] coordinates = ResidueState.storeAllCoordinates(residuesList);
// If x has not yet been constructed, construct it.
if (x == null) {
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
x = new double[nAtoms * 3];
}
double startingEnergy = 0;
double finalEnergy = 0;
try {
startingEnergy = currentEnergy(residuesList);
} catch (ArithmeticException ex) {
logger.severe(String.format(" Exception %s in calculating starting energy of a box; FFX shutting down", ex.toString()));
}
globalOptimization(residuesList);
try {
finalEnergy = currentEnergy(residuesList);
} catch (ArithmeticException ex) {
logger.severe(String.format(" Exception %s in calculating starting energy of a box; FFX shutting down", ex.toString()));
}
if (startingEnergy <= finalEnergy) {
logger.warning("Optimization did not yield a better energy. Reverting to orginal coordinates.");
ResidueState.revertAllCoordinates(residuesList, coordinates);
}
long currentTime = System.nanoTime();
boxTime += currentTime;
logIfMaster(format(" Time elapsed for this iteration: %11.3f sec", boxTime * 1.0E-9));
logIfMaster(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
} else {
globalOptimization(residuesList);
long currentTime = System.nanoTime();
boxTime += currentTime;
logIfMaster(format(" Time elapsed for this iteration: %11.3f sec", boxTime * 1.0E-9));
logIfMaster(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
}
if (master && printFiles) {
String filename = FilenameUtils.removeExtension(molecularAssembly.getFile().getAbsolutePath()) + ".partial";
File file = new File(filename);
if (firstCellSaved) {
file.delete();
}
// Don't write a file if it's the final iteration.
if (i == (numCells - 1)) {
continue;
}
PDBFilter windowFilter = new PDBFilter(file, molecularAssembly, null, null);
try {
windowFilter.writeFile(file, false);
if (firstResidue != lastResidue) {
logIfMaster(format(" File with residues %s ... %s in window written.", firstResidue.toString(), lastResidue.toString()));
} else {
logIfMaster(format(" File with residue %s in window written.", firstResidue.toString()));
}
firstCellSaved = true;
} catch (Exception e) {
logger.warning(format("Exception writing to file: %s", file.getName()));
}
}
/*for (Residue residue : residueList) {
if (residue instanceof MultiResidue) {
((MultiResidue) residue).setDefaultResidue();
residue.reInitOriginalAtomList();
}
}*/
} else {
logIfMaster(format(" Empty box: no residues found."));
}
}
return 0.0;
}
use of ffx.crystal.Crystal in project ffx by mjschnie.
the class RotamerOptimization method distanceMatrix.
/**
* Calculates a residue-residue distance matrix.
* <p>
* Residue-residue distance is defined as the shortest atom-atom distance in
* any possible rotamer-rotamer pair if the residues are neighbors (central
* atom-central atom distances are within a cutoff). Otherewise, distances
* are set to a default of Double.MAX_VALUE.
* </p>
* <p>
* The intent of using a neighbor list is to avoid tediously searching
* rotamer- rotamer pairs when two residues are so far apart we will never
* need the exact distance. We use the distance matrix for adding residues
* to the sliding window and determining whether to set 3-body energy to
* 0.0.
* </p>
* <p>
* If the central atoms are too distant from each other, we can safely
* assume no atoms will ever be close enough for addition to sliding window
* or to cause explicit calculation of 3-body energy.
* </p>
*/
private void distanceMatrix() {
distanceMatrix = new double[numResidues - 1][][][];
long numDistances = 0L;
for (int i = 0; i < (numResidues - 1); i++) {
Residue residuei = allResiduesArray[i];
int lengthRi;
try {
if (checkIfForced(residuei)) {
lengthRi = 1;
} else {
lengthRi = residuei.getRotamers(library).length;
}
} catch (IndexOutOfBoundsException ex) {
if (useForcedResidues) {
logger.warning(ex.toString());
} else {
logIfMaster(format(" Non-forced Residue i %s has null rotamers.", residuei.toFormattedString(false, true)), Level.WARNING);
}
continue;
}
distanceMatrix[i] = new double[lengthRi][][];
for (int ri = 0; ri < lengthRi; ri++) {
distanceMatrix[i][ri] = new double[numResidues][];
for (int j = (i + 1); j < numResidues; j++) {
Residue residuej = allResiduesArray[j];
int lengthRj;
try {
if (checkIfForced(residuej)) {
lengthRj = 1;
} else {
lengthRj = residuej.getRotamers(library).length;
}
} catch (IndexOutOfBoundsException ex) {
if (useForcedResidues) {
logger.warning(ex.toString());
} else {
logIfMaster(format(" Residue j %s has null rotamers.", residuej.toFormattedString(false, true)));
}
continue;
}
distanceMatrix[i][ri][j] = new double[lengthRj];
numDistances += lengthRj;
if (!lazyMatrix) {
fill(distanceMatrix[i][ri][j], Double.MAX_VALUE);
} else {
fill(distanceMatrix[i][ri][j], -1.0);
}
}
}
}
logger.info(format(" Number of pairwise distances: %d", numDistances));
if (!lazyMatrix) {
ResidueState[] orig = ResidueState.storeAllCoordinates(allResiduesList);
int nMultiRes = 0;
/**
* Build a list that contains one atom from each Residues: CA from
* amino acids, C1 from nucleic acids, or the first atom otherwise.
*/
Atom[] atoms = new Atom[numResidues];
for (int i = 0; i < numResidues; i++) {
Residue residuei = allResiduesArray[i];
atoms[i] = residuei.getReferenceAtom();
if (residuei instanceof MultiResidue) {
++nMultiRes;
}
}
/**
* Use of the pre-existing ParallelTeam causes a conflict when
* MultiResidues must re-init the force field. Temporary solution
* for sequence optimization: if > 1 residue optimized, run on only
* one thread.
*/
int nThreads = 1;
if (molecularAssembly.getPotentialEnergy().getParallelTeam() != null) {
nThreads = (nMultiRes > 1) ? 1 : molecularAssembly.getPotentialEnergy().getParallelTeam().getThreadCount();
} else {
// Suggested: nThreads = (nMultiRes > 1) ? 1 : ParallelTeam.getDefaultThreadCount();
nThreads = 16;
}
ParallelTeam parallelTeam = new ParallelTeam(nThreads);
Crystal crystal = molecularAssembly.getCrystal();
int nSymm = crystal.spaceGroup.getNumberOfSymOps();
logger.info("\n Computing Residue Distance Matrix");
double nlistCutoff = Math.max(Math.max(distance, twoBodyCutoffDist), threeBodyCutoffDist);
/**
* I think this originated from the fact that side-chain
* (and later nucleic acid) atoms could be fairly distant
* from the reference atom.
*/
double magicNumberBufferOfUnknownOrigin = 25.0;
nlistCutoff += magicNumberBufferOfUnknownOrigin;
NeighborList neighborList = new NeighborList(null, crystal, atoms, nlistCutoff, 0.0, parallelTeam);
// Expand coordinates
double[][] xyz = new double[nSymm][3 * numResidues];
double[] in = new double[3];
double[] out = new double[3];
for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
for (int i = 0; i < numResidues; i++) {
int i3 = i * 3;
int iX = i3 + 0;
int iY = i3 + 1;
int iZ = i3 + 2;
Atom atom = atoms[i];
in[0] = atom.getX();
in[1] = atom.getY();
in[2] = atom.getZ();
crystal.applySymOp(in, out, symOp);
xyz[iSymOp][iX] = out[0];
xyz[iSymOp][iY] = out[1];
xyz[iSymOp][iZ] = out[2];
}
}
// Build the residue neighbor-list.
int[][][] lists = new int[nSymm][numResidues][];
boolean[] use = new boolean[numResidues];
fill(use, true);
boolean forceRebuild = true;
boolean printLists = false;
long neighborTime = -System.nanoTime();
neighborList.buildList(xyz, lists, use, forceRebuild, printLists);
neighborTime += System.nanoTime();
logger.info(format(" Built residue neighbor list: %8.3f sec", neighborTime * 1.0e-9));
DistanceRegion distanceRegion = new DistanceRegion(parallelTeam.getThreadCount(), numResidues, crystal, lists, neighborList.getPairwiseSchedule());
long parallelTime = -System.nanoTime();
try {
parallelTeam.execute(distanceRegion);
} catch (Exception e) {
String message = " Exception compting residue distance matrix.";
logger.log(Level.SEVERE, message, e);
}
parallelTime += System.nanoTime();
logger.info(format(" Pairwise distance matrix: %8.3f sec\n", parallelTime * 1.0e-9));
ResidueState.revertAllCoordinates(allResiduesList, orig);
try {
parallelTeam.shutdown();
} catch (Exception ex) {
logger.warning(format(" Exception shutting down parallel team for the distance matrix: %s", ex.toString()));
}
}
}
use of ffx.crystal.Crystal in project ffx by mjschnie.
the class MainPanel method saveAsP1.
/**
* Save the currently selected FFXSystem to disk.
*
* @param file File to save the system to.
* @since 1.0
*/
public void saveAsP1(File file) {
FFXSystem system = hierarchy.getActive();
if (system != null && !system.isClosing()) {
File saveFile = file;
if (saveFile == null) {
resetFileChooser();
fileChooser.setCurrentDirectory(pwd);
fileChooser.setFileFilter(xyzFileFilter);
fileChooser.setAcceptAllFileFilterUsed(false);
int result = fileChooser.showSaveDialog(this);
if (result == JFileChooser.APPROVE_OPTION) {
saveFile = fileChooser.getSelectedFile();
pwd = saveFile.getParentFile();
}
}
if (saveFile != null) {
XYZFilter filter = new XYZFilter(saveFile, system, null, null);
ForceField forceField = system.getForceField();
final double a = forceField.getDouble(ForceFieldDouble.A_AXIS, 10.0);
final double b = forceField.getDouble(ForceFieldDouble.B_AXIS, a);
final double c = forceField.getDouble(ForceFieldDouble.C_AXIS, a);
final double alpha = forceField.getDouble(ForceFieldDouble.ALPHA, 90.0);
final double beta = forceField.getDouble(ForceFieldDouble.BETA, 90.0);
final double gamma = forceField.getDouble(ForceFieldDouble.GAMMA, 90.0);
final String spacegroup = forceField.getString(ForceFieldString.SPACEGROUP, "P1");
Crystal crystal = new Crystal(a, b, c, alpha, beta, gamma, spacegroup);
if (filter.writeFileAsP1(saveFile, false, crystal)) {
// Refresh Panels with the new System name
hierarchy.setActive(system);
}
activeFilter = filter;
}
}
}
use of ffx.crystal.Crystal in project ffx by mjschnie.
the class CIFFilter method getReflectionList.
/**
* {@inheritDoc}
*/
@Override
public ReflectionList getReflectionList(File cifFile, CompositeConfiguration properties) {
try {
BufferedReader br = new BufferedReader(new FileReader(cifFile));
String string;
while ((string = br.readLine()) != null) {
// Reached reflections, break.
if (string.startsWith("_refln.")) {
break;
}
String[] strArray = string.split("\\s+");
if (strArray[0].startsWith("_reflns") || strArray[0].startsWith("_cell") || strArray[0].startsWith("_symmetry")) {
String[] cifArray = strArray[0].split("\\.+");
switch(Header.toHeader(cifArray[1])) {
case d_resolution_high:
resHigh = Double.parseDouble(strArray[1]);
break;
case number_all:
nAll = Integer.parseInt(strArray[1]);
break;
case number_obs:
nObs = Integer.parseInt(strArray[1]);
break;
case length_a:
cell[0] = Double.parseDouble(strArray[1]);
break;
case length_b:
cell[1] = Double.parseDouble(strArray[1]);
break;
case length_c:
cell[2] = Double.parseDouble(strArray[1]);
break;
case angle_alpha:
cell[3] = Double.parseDouble(strArray[1]);
break;
case angle_beta:
cell[4] = Double.parseDouble(strArray[1]);
break;
case angle_gamma:
cell[5] = Double.parseDouble(strArray[1]);
break;
case Int_Tables_number:
spacegroupNum = Integer.parseInt(strArray[1]);
break;
case space_group_name_H_M:
String[] spacegroupNameArray = string.split("'+");
if (spacegroupNameArray.length > 1) {
spacegroupName = spacegroupNameArray[1];
} else if (strArray.length > 1) {
spacegroupName = strArray[1];
}
break;
}
}
}
br.close();
} catch (IOException e) {
String message = " CIF IO Exception.";
logger.log(Level.WARNING, message, e);
return null;
}
if (spacegroupNum < 0 && spacegroupName != null) {
spacegroupNum = SpaceGroup.spaceGroupNumber(SpaceGroup.pdb2ShortName(spacegroupName));
}
if (spacegroupNum < 0 || resHigh < 0 || cell[0] < 0) {
logger.info(" The CIF header contains insufficient information to generate the reflection list.");
return null;
}
if (logger.isLoggable(Level.INFO)) {
StringBuilder sb = new StringBuilder();
sb.append(String.format("\nOpening %s\n", cifFile.getName()));
sb.append(String.format("setting up Reflection List based on CIF:\n"));
sb.append(String.format(" spacegroup #: %d (name: %s)\n", spacegroupNum, SpaceGroup.spaceGroupNames[spacegroupNum - 1]));
sb.append(String.format(" resolution: %8.3f\n", 0.999999 * resHigh));
sb.append(String.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]));
sb.append(String.format("\n CIF # HKL (observed): %d\n", nObs));
sb.append(String.format(" CIF # HKL (all): %d\n", nAll));
logger.info(sb.toString());
}
Crystal crystal = new Crystal(cell[0], cell[1], cell[2], cell[3], cell[4], cell[5], SpaceGroup.spaceGroupNames[spacegroupNum - 1]);
double sampling = 1.0 / 1.5;
if (properties != null) {
sampling = properties.getDouble("sampling", 1.0 / 1.5);
}
Resolution resolution = new Resolution(0.999999 * resHigh, sampling);
ReflectionList reflectionlist = new ReflectionList(crystal, resolution, properties);
return reflectionlist;
}
Aggregations