org.jmol.adapter.smarter.XtalSymmetry Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of jmol Show documentation
Show all versions of jmol Show documentation
Jmol: an open-source Java viewer for chemical structures in 3D
/* $RCSfile$
* $Author: hansonr $
* $Date: 2014-01-10 09:19:33 -0600 (Fri, 10 Jan 2014) $
* $Revision: 19162 $
*
* Copyright (C) 2003-2005 Miguel, Jmol Development, www.jmol.org
*
* Contact: [email protected]
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
*/
package org.jmol.adapter.smarter;
import java.util.Hashtable;
import java.util.Map;
import java.util.Map.Entry;
import javajs.util.Lst;
import javajs.util.M3;
import javajs.util.M4;
import javajs.util.P3;
import javajs.util.P3i;
import javajs.util.PT;
import javajs.util.SB;
import javajs.util.T3;
import javajs.util.V3;
import org.jmol.api.JmolModulationSet;
import org.jmol.api.SymmetryInterface;
import org.jmol.java.BS;
import org.jmol.symmetry.Symmetry;
import org.jmol.symmetry.SymmetryOperation;
import org.jmol.util.BSUtil;
import org.jmol.util.Logger;
import org.jmol.util.SimpleUnitCell;
import org.jmol.util.Tensor;
import org.jmol.util.Vibration;
/**
*
* A class used by AtomSetCollection for building the symmetry of a model and
* generating new atoms based on that symmetry.
*
*/
public class XtalSymmetry {
private AtomSetCollection asc;
private AtomSetCollectionReader acr;
public XtalSymmetry() {
// for reflection
}
public XtalSymmetry set(AtomSetCollectionReader reader) {
this.acr = reader;
this.asc = reader.asc;
getSymmetry();
return this;
}
public SymmetryInterface symmetry;
SymmetryInterface getSymmetry() {
return (symmetry == null ?(symmetry = (Symmetry) acr.getInterface("org.jmol.symmetry.Symmetry")) : symmetry);
}
SymmetryInterface setSymmetry(SymmetryInterface symmetry) {
return (this.symmetry = symmetry);
}
private float[] unitCellParams = new float[6];
private float[] baseUnitCell;
// expands to 26 for cartesianToFractional matrix as array (PDB) and supercell
private float symmetryRange;
private boolean doCentroidUnitCell;
private boolean centroidPacked;
private float packingError;
private String filterSymop;
private void setSymmetryRange(float factor) {
symmetryRange = factor;
asc.setInfo("symmetryRange", Float.valueOf(factor));
}
private boolean applySymmetryToBonds = false;
private int[] latticeCells;
private void setLatticeCells() {
// int[] latticeCells, boolean applySymmetryToBonds,
// }
// boolean doPackUnitCell, boolean doCentroidUnitCell,
// boolean centroidPacked, String strSupercell,
// P3 ptSupercell) {
//set when unit cell is determined
// x <= 555 and y >= 555 indicate a range of cells to load
// AROUND the central cell 555 and that
// we should normalize (z = 1) or pack unit cells (z = -1) or not (z = 0)
// in addition (Jmol 11.7.36) z = -2 does a full 3x3x3 around the designated cells
// but then only delivers the atoms that are within the designated cells.
// Normalization is the moving of the center of mass into the unit cell.
// Starting with Jmol 12.0.RC23 we do not normalize a CIF file that
// is being loaded without {i j k} indicated.
latticeCells = acr.latticeCells;
boolean isLatticeRange = (latticeCells[0] <= 555 && latticeCells[1] >= 555 && (latticeCells[2] == 0
|| latticeCells[2] == 1 || latticeCells[2] == -1));
doNormalize = latticeCells[0] != 0
&& (!isLatticeRange || latticeCells[2] == 1);
applySymmetryToBonds = acr.applySymmetryToBonds;
doPackUnitCell = acr.doPackUnitCell;
doCentroidUnitCell = acr.doCentroidUnitCell;
centroidPacked = acr.centroidPacked;
filterSymop = acr.filterSymop;
if (acr.strSupercell == null)
setSupercellFromPoint(acr.ptSupercell);
}
private P3 ptSupercell;
//private float[] fmatSupercell;
private M4 matSupercell;
public void setSupercellFromPoint(P3 pt) {
ptSupercell = pt;
if (pt == null) {
matSupercell = null;
return;
}
matSupercell = new M4();
matSupercell.m00 = pt.x;
matSupercell.m11 = pt.y;
matSupercell.m22 = pt.z;
matSupercell.m33 = 1;
Logger.info("Using supercell \n" + matSupercell);
}
private Lst trajectoryUnitCells;
private void setUnitCell(float[] info, M3 matUnitCellOrientation,
P3 unitCellOffset) {
unitCellParams = new float[info.length];
//this.unitCellOffset = unitCellOffset;
for (int i = 0; i < info.length; i++)
unitCellParams[i] = info[i];
asc.haveUnitCell = true;
asc.setCurrentModelInfo("unitCellParams", unitCellParams);
if (asc.isTrajectory) {
if (trajectoryUnitCells == null) {
trajectoryUnitCells = new Lst();
asc.setInfo("unitCells", trajectoryUnitCells);
}
trajectoryUnitCells.addLast(unitCellParams);
}
asc.setGlobalBoolean(AtomSetCollection.GLOBAL_UNITCELLS);
getSymmetry().setUnitCell(unitCellParams, false);
// we need to set the auxiliary info as well, because
// ModelLoader creates a new symmetry object.
if (unitCellOffset != null) {
symmetry.setOffsetPt(unitCellOffset);
asc.setCurrentModelInfo("unitCellOffset", unitCellOffset);
}
if (matUnitCellOrientation != null) {
symmetry.initializeOrientation(matUnitCellOrientation);
asc.setCurrentModelInfo("matUnitCellOrientation",
matUnitCellOrientation);
}
}
int addSpaceGroupOperation(String xyz, boolean andSetLattice) {
if (andSetLattice)
setLatticeCells();
symmetry.setSpaceGroup(doNormalize);
return symmetry.addSpaceGroupOperation(xyz, 0);
}
public void setLatticeParameter(int latt) {
symmetry.setSpaceGroup(doNormalize);
symmetry.setLattice(latt);
}
private boolean doNormalize = true;
private boolean doPackUnitCell = false;
private SymmetryInterface baseSymmetry;
private SymmetryInterface sym2;
SymmetryInterface applySymmetryFromReader(SymmetryInterface readerSymmetry)
throws Exception {
asc.setCoordinatesAreFractional(acr.iHaveFractionalCoordinates);
setUnitCell(acr.unitCellParams, acr.matUnitCellOrientation,
acr.unitCellOffset);
setAtomSetSpaceGroupName(acr.sgName);
setSymmetryRange(acr.symmetryRange);
if (acr.doConvertToFractional || acr.fileCoordinatesAreFractional) {
setLatticeCells();
boolean doApplySymmetry = true;
if (acr.ignoreFileSpaceGroupName || !acr.iHaveSymmetryOperators) {
if (!acr.merging || readerSymmetry == null)
readerSymmetry = acr.getNewSymmetry();
doApplySymmetry = readerSymmetry.createSpaceGroup(
acr.desiredSpaceGroupIndex, (acr.sgName.indexOf("!") >= 0 ? "P1"
: acr.sgName), acr.unitCellParams);
} else {
acr.doPreSymmetry();
readerSymmetry = null;
}
packingError = acr.packingError;
if (doApplySymmetry) {
if (readerSymmetry != null)
setSpaceGroupFrom(readerSymmetry);
//parameters are counts of unit cells as [a b c]
applySymmetryLattice();
if (readerSymmetry != null && filterSymop == null)
setAtomSetSpaceGroupName(readerSymmetry.getSpaceGroupName());
}
}
if (acr.iHaveFractionalCoordinates && acr.merging && readerSymmetry != null) {
// when merging (with appendNew false), we must return cartesians
Atom[] atoms = asc.atoms;
for (int i = asc.getLastAtomSetAtomIndex(), n = asc.ac; i < n; i++)
readerSymmetry.toCartesian(atoms[i], true);
asc.setCoordinatesAreFractional(false);
// We no longer allow merging of multiple-model files
// when the file to be appended has fractional coordinates and vibrations
acr.addVibrations = false;
}
return symmetry;
}
public void setSpaceGroupFrom(SymmetryInterface readerSymmetry) {
getSymmetry().setSpaceGroupFrom(readerSymmetry);
}
private void setAtomSetSpaceGroupName(String spaceGroupName) {
asc.setCurrentModelInfo("spaceGroup", spaceGroupName + "");
}
private void applySymmetryLattice() throws Exception {
if (!asc.coordinatesAreFractional || symmetry.getSpaceGroup() == null)
return;
sym2 = null;
int maxX = latticeCells[0];
int maxY = latticeCells[1];
int maxZ = Math.abs(latticeCells[2]);
firstSymmetryAtom = asc.getLastAtomSetAtomIndex();
BS bsAtoms = asc.bsAtoms;
if (bsAtoms != null)
firstSymmetryAtom = bsAtoms.nextSetBit(firstSymmetryAtom);
rminx = rminy = rminz = Float.MAX_VALUE;
rmaxx = rmaxy = rmaxz = -Float.MAX_VALUE;
P3 pt0 = null;
if (acr.fillRange != null) {
if (bsAtoms == null)
asc.bsAtoms = bsAtoms = new BS();
bsAtoms.setBits(firstSymmetryAtom, asc.ac);
doPackUnitCell = false;
minXYZ = new P3i();
maxXYZ = P3i.new3(1, 1, 1);
P3[] oabc = new P3[4];
for (int i = 0; i < 4; i++)
oabc[i] = P3.newP(acr.fillRange[i]);
adjustRangeMinMax(oabc);
//Logger.info("setting min/max for original lattice to " + minXYZ + " and "
// + maxXYZ);
if (sym2 == null) {
sym2 = new Symmetry();
sym2.getUnitCell(acr.fillRange, false, null);
}
applyAllSymmetry(acr.ms, bsAtoms);
pt0 = new P3();
Atom[] atoms = asc.atoms;
for (int i = asc.ac; --i >= firstSymmetryAtom; ) {
pt0.setT(atoms[i]);
symmetry.toCartesian(pt0, false);
sym2.toFractional(pt0, false);
if (acr.fixJavaFloat)
PT.fixPtFloats(pt0, PT.FRACTIONAL_PRECISION);
if (!isWithinCell(dtype, pt0, 0, 1, 0, 1, 0, 1, packingError))
bsAtoms.clear(i);
}
return;
}
P3 offset = null;
nVib = 0;
T3 va = null, vb = null, vc = null;
baseSymmetry = symmetry;
String supercell = acr.strSupercell;
T3[] oabc = null;
boolean isSuper = (supercell != null && supercell.indexOf(",") >= 0);
if (isSuper) {
// expand range to accommodate this alternative cell
// oabc will be cartesian
oabc = symmetry.getV0abc(supercell);
if (oabc != null) {
// set the bounds for atoms in the new unit cell
// in terms of the old unit cell
minXYZ = new P3i();
maxXYZ = P3i.new3(maxX, maxY, maxZ);
symmetry.setMinMaxLatticeParameters(minXYZ, maxXYZ);
// base origin for new unit cell
pt0 = P3.newP(oabc[0]);
// base vectors for new unit cell
va = P3.newP(oabc[1]);
vb = P3.newP(oabc[2]);
vc = P3.newP(oabc[3]);
adjustRangeMinMax(oabc);
}
}
int iAtomFirst = asc.getLastAtomSetAtomIndex();
if (bsAtoms != null)
iAtomFirst = bsAtoms.nextSetBit(iAtomFirst);
if (rminx == Float.MAX_VALUE) {
matSupercell = null;
supercell = null;
oabc = null;
} else {
boolean doPack0 = doPackUnitCell;
doPackUnitCell = doPack0;//(doPack0 || oabc != null && acr.forcePacked);
if (asc.bsAtoms == null)
asc.bsAtoms = BSUtil.newBitSet2(0, asc.ac);
bsAtoms = asc.bsAtoms;
applyAllSymmetry(acr.ms, null);
doPackUnitCell = doPack0;
// 2) set all atom coordinates to Cartesians
Atom[] atoms = asc.atoms;
int atomCount = asc.ac;
for (int i = iAtomFirst; i < atomCount; i++) {
symmetry.toCartesian(atoms[i], true);
bsAtoms.set(i);
}
// 3) create the supercell unit cell
symmetry = null;
symmetry = getSymmetry();
setUnitCell(new float[] { 0, 0, 0, 0, 0, 0, va.x, va.y, va.z, vb.x, vb.y,
vb.z, vc.x, vc.y, vc.z }, null, offset);
setAtomSetSpaceGroupName(oabc == null || supercell == null ? "P1"
: "cell=" + supercell);
symmetry.setSpaceGroup(doNormalize);
symmetry.addSpaceGroupOperation("x,y,z", 0);
// 4) reset atoms to fractional values in this new system
if (pt0 != null)
symmetry.toFractional(pt0, true);
for (int i = iAtomFirst; i < atomCount; i++) {
symmetry.toFractional(atoms[i], true);
if (pt0 != null)
atoms[i].sub(pt0);
}
// 5) apply the full lattice symmetry now
asc.haveAnisou = false;
// ?? TODO
asc.setCurrentModelInfo("matUnitCellOrientation", null);
}
minXYZ = new P3i();
maxXYZ = P3i.new3(maxX, maxY, maxZ);
symmetry.setMinMaxLatticeParameters(minXYZ, maxXYZ);
if (oabc == null) {
applyAllSymmetry(acr.ms, bsAtoms);
return;
}
if (acr.forcePacked || doPackUnitCell) {
// trim atom set based on original unit cell
BS bs = asc.bsAtoms;
Atom[] atoms = asc.atoms;
if (bs == null)
bs = asc.bsAtoms = BSUtil.newBitSet2(0, asc.ac);
for (int i = bs.nextSetBit(iAtomFirst); i >= 0; i = bs.nextSetBit(i + 1)) {
if (!isWithinCell(dtype, atoms[i], minXYZ.x, maxXYZ.x, minXYZ.y,
maxXYZ.y, minXYZ.z, maxXYZ.z, packingError))
bs.clear(i);
}
}
// but we leave matSupercell, because we might need it for vibrations in CASTEP
}
private void adjustRangeMinMax(T3[] oabc) {
// be sure to add in packing error adjustments
// so that we include all needed atoms
// load "" packed x.x supercell...
P3 pa = new P3();
P3 pb = new P3();
P3 pc = new P3();
if (acr.forcePacked) {
pa.setT(oabc[1]);
pb.setT(oabc[2]);
pc.setT(oabc[3]);
pa.scale(packingError);
pb.scale(packingError);
pc.scale(packingError);
}
// account for lattice specification
// load "" {x y z} supercell...
oabc[0].scaleAdd2(minXYZ.x, oabc[1], oabc[0]);
oabc[0].scaleAdd2(minXYZ.y, oabc[2], oabc[0]);
oabc[0].scaleAdd2(minXYZ.z, oabc[3], oabc[0]);
// add in packing adjustment
oabc[0].sub(pa);
oabc[0].sub(pb);
oabc[0].sub(pc);
// fractionalize and adjust min/max
P3 pt = P3.newP(oabc[0]);
symmetry.toFractional(pt, true);
setSymmetryMinMax(pt);
// account for lattice specification
oabc[1].scale(maxXYZ.x - minXYZ.x);
oabc[2].scale(maxXYZ.y - minXYZ.y);
oabc[3].scale(maxXYZ.z - minXYZ.z);
// add in packing adjustment
oabc[1].scaleAdd2(2, pa, oabc[1]);
oabc[2].scaleAdd2(2, pb, oabc[2]);
oabc[3].scaleAdd2(2, pc, oabc[3]);
// run through six of the corners -- a, b, c, ab, ac, bc
for (int i = 0; i < 3; i++) {
for (int j = i + 1; j < 4; j++) {
pt.add2(oabc[i], oabc[j]);
if (i != 0)
pt.add(oabc[0]);
symmetry.toFractional(pt, false);
setSymmetryMinMax(pt);
}
}
// bc in the end, so we need abc
symmetry.toCartesian(pt, false);
pt.add(oabc[1]);
symmetry.toFractional(pt, false);
setSymmetryMinMax(pt);
// allow for some imprecision
minXYZ = P3i.new3((int) Math.min(0,Math.floor(rminx + 0.001f)),
(int) Math.min(0, Math.floor(rminy + 0.001f)),
(int) Math.min(0, Math.floor(rminz + 0.001f)));
maxXYZ = P3i.new3((int) Math.max(1,Math.ceil(rmaxx - 0.001f)),
(int) Math.max(1,Math.ceil(rmaxy - 0.001f)),
(int) Math.max(1,Math.ceil(rmaxz - 0.001f)));
}
/**
* range minima and maxima -- also usedf for cartesians comparisons
*
*/
private float rminx, rminy, rminz, rmaxx, rmaxy, rmaxz;
private void setSymmetryMinMax(P3 c) {
if (rminx > c.x)
rminx = c.x;
if (rminy > c.y)
rminy = c.y;
if (rminz > c.z)
rminz = c.z;
if (rmaxx < c.x)
rmaxx = c.x;
if (rmaxy < c.y)
rmaxy = c.y;
if (rmaxz < c.z)
rmaxz = c.z;
}
private final P3 ptOffset = new P3();
//private P3 unitCellOffset;
private P3i minXYZ, maxXYZ;
private P3 minXYZ0, maxXYZ0;
public boolean isWithinCell(int dtype, P3 pt, float minX, float maxX,
float minY, float maxY, float minZ, float maxZ,
float slop) {
return (pt.x > minX - slop && pt.x < maxX + slop
&& (dtype < 2 || pt.y > minY - slop && pt.y < maxY + slop)
&& (dtype < 3 || pt.z > minZ - slop && pt.z < maxZ + slop));
}
// /**
// * A problem arises when converting to JavaScript, because JavaScript numbers are all
// * doubles, while here we have floats. So what we do is to multiply by a number that
// * is beyond the precision of our data but within the range of floats -- namely, 100000,
// * and then integerize. This ensures that both doubles and floats compare the same number.
// * Unfortunately, it will break Java reading of older files, so we check for legacy versions.
// *
// * @param dtype
// * @param pt
// * @param minX
// * @param maxX
// * @param minY
// * @param maxY
// * @param minZ
// * @param maxZ
// * @param slop
// * @return true if within range
// */
// private boolean xxxisWithinCellInt(int dtype, P3 pt, float minX, float maxX,
// float minY, float maxY, float minZ, float maxZ,
// int slop) {
// switch (dtype) {
// case 3:
// if (Math.round((minZ - pt.z) * 100000) >= slop || Math.round((pt.z - maxZ) * 100000) >= slop)
// return false;
// //$FALL-THROUGH$
// case 2:
// if (Math.round((minY - pt.y) * 100000) >= slop || Math.round((pt.y - maxY) * 100000) >= slop)
// return false;
// //$FALL-THROUGH$
// case 1:
// if (Math.round((minX - pt.x) * 100000) >= slop || Math.round((pt.x - maxX) * 100000) >= slop)
// return false;
// break;
// }
// return true;
// }
/**
* @param ms
* modulated structure interface
* @param bsAtoms
* relating to supercells
* @throws Exception
*/
private void applyAllSymmetry(MSInterface ms, BS bsAtoms) throws Exception {
if (asc.ac == 0 || bsAtoms != null && bsAtoms.isEmpty())
return;
int n = noSymmetryCount = asc.baseSymmetryAtomCount > 0 ? n = asc.baseSymmetryAtomCount
: bsAtoms == null ? asc.getLastAtomSetAtomCount() : asc.ac
- bsAtoms.nextSetBit(asc.getLastAtomSetAtomIndex());
asc.setTensors();
bondCount0 = asc.bondCount;
finalizeSymmetry(symmetry);
int operationCount = symmetry.getSpaceGroupOperationCount();
BS excludedOps = (acr.thisBiomolecule == null ? null : new BS());
if (excludedOps != null)
asc.checkSpecial = true;
dtype = (int) symmetry.getUnitCellInfoType(SimpleUnitCell.INFO_DIMENSIONS);
symmetry.setMinMaxLatticeParameters(minXYZ, maxXYZ);
if (doCentroidUnitCell)
asc.setInfo("centroidMinMax", new int[] { minXYZ.x, minXYZ.y, minXYZ.z,
maxXYZ.x, maxXYZ.y, maxXYZ.z, (centroidPacked ? 1 : 0) });
if (ptSupercell != null) {
asc.setCurrentModelInfo("supercell", ptSupercell);
switch (dtype) {
case 3:
// standard
minXYZ.z *= (int) Math.abs(ptSupercell.z);
maxXYZ.z *= (int) Math.abs(ptSupercell.z);
//$FALL-THROUGH$;
case 2:
// slab or standard
minXYZ.y *= (int) Math.abs(ptSupercell.y);
maxXYZ.y *= (int) Math.abs(ptSupercell.y);
//$FALL-THROUGH$;
case 1:
// slab, polymer, or standard
minXYZ.x *= (int) Math.abs(ptSupercell.x);
maxXYZ.x *= (int) Math.abs(ptSupercell.x);
}
}
if (doCentroidUnitCell || doPackUnitCell || symmetryRange != 0
&& maxXYZ.x - minXYZ.x == 1 && maxXYZ.y - minXYZ.y == 1
&& maxXYZ.z - minXYZ.z == 1) {
// weird Mac bug does not allow Point3i.new3(minXYZ) !!
minXYZ0 = P3.new3(minXYZ.x, minXYZ.y, minXYZ.z);
maxXYZ0 = P3.new3(maxXYZ.x, maxXYZ.y, maxXYZ.z);
if (ms != null) {
ms.setMinMax0(minXYZ0, maxXYZ0);
minXYZ.set((int) minXYZ0.x, (int) minXYZ0.y, (int) minXYZ0.z);
maxXYZ.set((int) maxXYZ0.x, (int) maxXYZ0.y, (int) maxXYZ0.z);
}
switch (dtype) {
case 3:
// standard
minXYZ.z--;
maxXYZ.z++;
//$FALL-THROUGH$;
case 2:
// slab or standard
minXYZ.y--;
maxXYZ.y++;
//$FALL-THROUGH$;
case 1:
// slab, polymer, or standard
minXYZ.x--;
maxXYZ.x++;
}
}
int nCells = (maxXYZ.x - minXYZ.x) * (maxXYZ.y - minXYZ.y)
* (maxXYZ.z - minXYZ.z);
int cartesianCount = (asc.checkSpecial || acr.thisBiomolecule != null ? n
* operationCount * nCells : symmetryRange > 0 ? n * operationCount // checking against {1 1 1}
// : symmetryRange < 0 ? 1 // checking against symop=1555 set; just a box
: 1 // not checking
);
P3[] cartesians = new P3[cartesianCount];
Atom[] atoms = asc.atoms;
for (int i = 0; i < n; i++)
atoms[i + firstSymmetryAtom].bsSymmetry = BS.newN(operationCount
* (nCells + 1));
int pt = 0;
int[] unitCells = new int[nCells];
unitCellTranslations = new V3[nCells];
int iCell = 0;
int cell555Count = 0;
float absRange = Math.abs(symmetryRange);
boolean checkCartesianRange = (symmetryRange != 0);
boolean checkRangeNoSymmetry = (symmetryRange < 0);
boolean checkRange111 = (symmetryRange > 0);
if (checkCartesianRange) {
rminx = rminy = rminz = Float.MAX_VALUE;
rmaxx = rmaxy = rmaxz = -Float.MAX_VALUE;
}
// always do the 555 cell first
// incommensurate symmetry can have lattice centering, resulting in
// duplication of operators. There's a bug later on that requires we
// only do this with the first atom set for now, at least.
SymmetryInterface symmetry = this.symmetry;
SymmetryInterface lastSymmetry = symmetry;
latticeOp = symmetry.getLatticeOp();
checkAll = (asc.atomSetCount == 1 && asc.checkSpecial && latticeOp >= 0);
latticeOnly = (asc.checkLatticeOnly && latticeOp >= 0); // CrystalReader
P3 pttemp = null;
M4 op = symmetry.getSpaceGroupOperation(0);
if (doPackUnitCell) {
pttemp = new P3();
ptOffset.set(0, 0, 0);
}
int[] atomMap = (bondCount0 > asc.bondIndex0 && applySymmetryToBonds ? new int[n]
: null);
Lst lstNCS = acr.lstNCS;
if (lstNCS != null && lstNCS.get(0).m33 == 0) {
int nOp = symmetry.getSpaceGroupOperationCount();
int nn = lstNCS.size();
for (int i = nn; --i >= 0;) {
M4 m = lstNCS.get(i);
m.m33 = 1;
symmetry.toFractionalM(m);
}
for (int i = 1; i < nOp; i++) {
M4 m1 = symmetry.getSpaceGroupOperation(i);
for (int j = 0; j < nn; j++) {
M4 m = M4.newM4(lstNCS.get(j));
m.mul2(m1, m);
if (doNormalize)
SymmetryOperation.setOffset(m, atoms, firstSymmetryAtom,
noSymmetryCount);
lstNCS.addLast(m);
}
}
}
for (int tx = minXYZ.x; tx < maxXYZ.x; tx++)
for (int ty = minXYZ.y; ty < maxXYZ.y; ty++)
for (int tz = minXYZ.z; tz < maxXYZ.z; tz++) {
unitCellTranslations[iCell] = V3.new3(tx, ty, tz);
unitCells[iCell++] = 555 + tx * 100 + ty * 10 + tz;
if (tx != 0 || ty != 0 || tz != 0 || cartesians.length == 0)
continue;
// base cell only
for (pt = 0; pt < n; pt++) {
Atom atom = atoms[firstSymmetryAtom + pt];
if (ms != null) {
symmetry = ms.getAtomSymmetry(atom, this.symmetry);
if (symmetry != lastSymmetry) {
if (symmetry.getSpaceGroupOperationCount() == 0)
finalizeSymmetry(lastSymmetry = symmetry);
op = symmetry.getSpaceGroupOperation(0);
}
}
P3 c = P3.newP(atom);
op.rotTrans(c);
symmetry.toCartesian(c, false);
if (doPackUnitCell) {
symmetry.toUnitCell(c, ptOffset);
pttemp.setT(c);
symmetry.toFractional(pttemp, false);
if (acr.fixJavaFloat)
PT.fixPtFloats(pttemp, PT.FRACTIONAL_PRECISION);
// when bsAtoms != null, we are
// setting it to be correct for a
// second unit cell -- the supercell
if (bsAtoms == null)
atom.setT(pttemp);
else if (atom.distance(pttemp) < 0.0001f)
bsAtoms.set(atom.index);
else {// not in THIS unit cell
bsAtoms.clear(atom.index);
continue;
}
}
if (bsAtoms != null)
atom.bsSymmetry.clearAll();
atom.bsSymmetry.set(iCell * operationCount);
atom.bsSymmetry.set(0);
if (checkCartesianRange)
setSymmetryMinMax(c);
if (pt < cartesianCount)
cartesians[pt] = c;
}
if (checkRangeNoSymmetry) {
rminx -= absRange;
rminy -= absRange;
rminz -= absRange;
rmaxx += absRange;
rmaxy += absRange;
rmaxz += absRange;
}
cell555Count = pt = symmetryAddAtoms(0, 0, 0, 0, pt, iCell
* operationCount, cartesians, ms, excludedOps, atomMap);
}
if (checkRange111) {
rminx -= absRange;
rminy -= absRange;
rminz -= absRange;
rmaxx += absRange;
rmaxy += absRange;
rmaxz += absRange;
}
// now apply all the translations
iCell = 0;
for (int tx = minXYZ.x; tx < maxXYZ.x; tx++)
for (int ty = minXYZ.y; ty < maxXYZ.y; ty++)
for (int tz = minXYZ.z; tz < maxXYZ.z; tz++) {
iCell++;
if (tx != 0 || ty != 0 || tz != 0)
pt = symmetryAddAtoms(tx, ty, tz, cell555Count, pt, iCell
* operationCount, cartesians, ms, excludedOps, atomMap);
}
if (iCell * n == asc.ac - firstSymmetryAtom)
duplicateAtomProperties(iCell);
setSymmetryOps();
asc.setCurrentModelInfo("presymmetryAtomIndex",
Integer.valueOf(firstSymmetryAtom));
asc.setCurrentModelInfo("presymmetryAtomCount", Integer.valueOf(n));
asc.setCurrentModelInfo("latticeDesignation",
symmetry.getLatticeDesignation());
asc.setCurrentModelInfo("unitCellRange", unitCells);
asc.setCurrentModelInfo("unitCellTranslations", unitCellTranslations);
baseUnitCell = unitCellParams;
unitCellParams = new float[6];
reset();
}
private boolean checkAll;
private int bondCount0;
private int symmetryAddAtoms(int transX, int transY, int transZ,
int baseCount, int pt, int iCellOpPt,
P3[] cartesians, MSInterface ms, BS excludedOps,
int[] atomMap) throws Exception {
boolean isBaseCell = (baseCount == 0);
boolean addBonds = (atomMap != null);
if (doPackUnitCell)
ptOffset.set(transX, transY, transZ);
//symmetryRange < 0 : just check symop=1 set
//symmetryRange > 0 : check against {1 1 1}
// if we are not checking special atoms, then this is a PDB file
// and we return all atoms within a cubical volume around the
// target set. The user can later use select within() to narrow that down
// This saves immensely on time.
float range2 = symmetryRange * symmetryRange;
boolean checkRangeNoSymmetry = (symmetryRange < 0);
boolean checkRange111 = (symmetryRange > 0);
boolean checkSymmetryMinMax = (isBaseCell && checkRange111);
checkRange111 &= !isBaseCell;
int nOp = symmetry.getSpaceGroupOperationCount();
Lst lstNCS = acr.lstNCS;
int nNCS = (lstNCS == null ? 0 : lstNCS.size());
int nOperations = nOp + nNCS;
nNCS = nNCS / nOp;
boolean checkSpecial = (nOperations == 1 && !doPackUnitCell ? false
: asc.checkSpecial);
boolean checkSymmetryRange = (checkRangeNoSymmetry || checkRange111);
boolean checkDistances = (checkSpecial || checkSymmetryRange);
boolean checkOps = (excludedOps != null);
boolean addCartesian = (checkSpecial || checkSymmetryMinMax);
BS bsAtoms = (acr.isMolecular ? null : asc.bsAtoms);
SymmetryInterface symmetry = this.symmetry;
if (checkRangeNoSymmetry)
baseCount = noSymmetryCount;
int atomMax = firstSymmetryAtom + noSymmetryCount;
P3 ptAtom = new P3();
String code = null;
float d0 = (checkOps ? 0.01f : 0.0001f);
char subSystemId = '\0';
int j00 = (bsAtoms == null ? firstSymmetryAtom : bsAtoms.nextSetBit(0));
out: for (int iSym = 0; iSym < nOperations; iSym++) {
if (isBaseCell && iSym == 0 || latticeOnly && iSym > 0
&& iSym != latticeOp || excludedOps != null && excludedOps.get(iSym))
continue;
/* pt0 sets the range of points cross-checked.
* If we are checking special positions, then we have to check
* all previous atoms.
* If we are doing a symmetry range check relative to {1 1 1}, then
* we have to check only the base set. (checkRange111 true)
* If we are doing a symmetry range check on symop=1555 (checkRangeNoSymmetry true),
* then we don't check any atoms and just use the box.
*
*/
int pt0 = (checkSpecial || excludedOps != null ? pt
: checkRange111 ? baseCount : 0);
float spinOp = (iSym >= nOp ? 0 : asc.vibScale == 0 ? symmetry
.getSpinOp(iSym) : asc.vibScale);
int i0 = Math.max(firstSymmetryAtom,
(bsAtoms == null ? 0 : bsAtoms.nextSetBit(0)));
boolean checkDistance = checkDistances;
int spt = (iSym >= nOp ? (iSym - nOp) / nNCS : iSym);
int cpt = spt + iCellOpPt;
for (int i = i0; i < atomMax; i++) {
Atom a = asc.atoms[i];
if (a.ignoreSymmetry || bsAtoms != null && !bsAtoms.get(i))
continue;
if (ms == null) {
symmetry.newSpaceGroupPoint(iSym, a, ptAtom, transX, transY, transZ,
(iSym >= nOp ? lstNCS.get(iSym - nOp) : null));
} else {
symmetry = ms.getAtomSymmetry(a, this.symmetry);
symmetry.newSpaceGroupPoint(iSym, a, ptAtom, transX, transY, transZ,
null);
// COmmensurate structures may use a symmetry operator
// to changes space groups.
code = symmetry.getSpaceGroupOperationCode(iSym);
if (code != null) {
subSystemId = code.charAt(0);
symmetry = ms.getSymmetryFromCode(code);
if (symmetry.getSpaceGroupOperationCount() == 0)
finalizeSymmetry(symmetry);
}
}
if (acr.fixJavaFloat)
PT.fixPtFloats(ptAtom, PT.FRACTIONAL_PRECISION);
P3 c = P3.newP(ptAtom); // cartesian position
symmetry.toCartesian(c, false);
if (doPackUnitCell) {
// note that COmmensurate structures may need
// modulation at this point.
symmetry.toUnitCell(c, ptOffset);
ptAtom.setT(c);
symmetry.toFractional(ptAtom, false);
if (acr.fixJavaFloat)
PT.fixPtFloats(ptAtom, PT.FRACTIONAL_PRECISION);
if (!isWithinCell(dtype, ptAtom, minXYZ0.x, maxXYZ0.x, minXYZ0.y,
maxXYZ0.y, minXYZ0.z, maxXYZ0.z, packingError))
continue;
}
if (checkSymmetryMinMax)
setSymmetryMinMax(c);
Atom special = null;
if (checkDistance) {
// for range checking, we first make sure we are not out of range
// for the cartesian
if (checkSymmetryRange
&& (c.x < rminx || c.y < rminy || c.z < rminz || c.x > rmaxx
|| c.y > rmaxy || c.z > rmaxz))
continue;
float minDist2 = Float.MAX_VALUE;
int j0 = (checkAll ? asc.ac : pt0);
String name = a.atomName;
char id = (code == null ? a.altLoc : subSystemId);
for (int j = j00; j < j0; j++) {
if (bsAtoms != null && !bsAtoms.get(j))
continue;
P3 pc = cartesians[j];
if (pc == null)
continue;
float d2 = c.distanceSquared(pc);
//System.out.println(iSym + " " + i + " " + j + " " + pc + " " + c + " " + d2);
if (checkSpecial && d2 < d0) {
/* checkSpecial indicates that we are looking for atoms with (nearly) the
* same cartesian position.
*/
if (checkOps) {
// if a matching atom is found for a model built
// from a mix of crystallographic and noncrystallographic
// operators, we throw out the entire operation, not just this atom
excludedOps.set(iSym);
continue out;
}
special = asc.atoms[firstSymmetryAtom + j];
if ((special.atomName == null || special.atomName.equals(name))
&& special.altLoc == id)
break;
special = null;
}
if (checkRange111 && j < baseCount && d2 < minDist2)
minDist2 = d2;
}
if (checkRange111 && minDist2 > range2)
continue;
}
if (checkOps) {
// if we did not find a common atom for the first atom when checking operators,
// we do not have to check again for any other atom.
checkDistance = false;
}
int atomSite = a.atomSite;
if (special != null) {
if (addBonds)
atomMap[atomSite] = special.index;
special.bsSymmetry.set(cpt);
special.bsSymmetry.set(spt);
} else {
if (addBonds)
atomMap[atomSite] = asc.ac;
Atom atom1 = asc.newCloneAtom(a);
if (asc.bsAtoms != null)
asc.bsAtoms.set(atom1.index);
atom1.setT(ptAtom);
if (spinOp != 0 && atom1.vib != null) {
// spinOp is making the correction for spin being a pseudoVector, not a standard vector
symmetry.getSpaceGroupOperation(iSym).rotate(atom1.vib);
atom1.vib.scale(spinOp);
}
atom1.atomSite = atomSite;
if (code != null)
atom1.altLoc = subSystemId;
atom1.bsSymmetry = BSUtil.newAndSetBit(cpt);
atom1.bsSymmetry.set(spt);
if (addCartesian)
cartesians[pt++] = c;
Lst