Revision: 18531
http://sourceforge.net/p/jmol/code/18531
Author: hansonr
Date: 2013-08-08 23:30:00 +0000 (Thu, 08 Aug 2013)
Log Message:
-----------
___JmolVersion="13.3.4_dev_2013.08.08"
code: incommensurate crystal work -- occupancy Crenel, displacive sawtooths
Modified Paths:
--------------
trunk/Jmol/src/org/jmol/adapter/readers/cif/CifReader.java
trunk/Jmol/src/org/jmol/adapter/readers/cif/ModulationReader.java
trunk/Jmol/src/org/jmol/adapter/readers/xtal/JanaReader.java
trunk/Jmol/src/org/jmol/symmetry/SpaceGroup.java
trunk/Jmol/src/org/jmol/symmetry/SymmetryOperation.java
trunk/Jmol/src/org/jmol/util/Modulation.java
trunk/Jmol/src/org/jmol/viewer/Jmol.properties
Modified: trunk/Jmol/src/org/jmol/adapter/readers/cif/CifReader.java
===================================================================
--- trunk/Jmol/src/org/jmol/adapter/readers/cif/CifReader.java 2013-08-08
14:18:07 UTC (rev 18530)
+++ trunk/Jmol/src/org/jmol/adapter/readers/cif/CifReader.java 2013-08-08
23:30:00 UTC (rev 18531)
@@ -257,7 +257,6 @@
@Override
protected void finalizeReader() throws Exception {
- finalizeIncommensurate();
if (atomSetCollection.getAtomCount() == nAtoms)
atomSetCollection.removeCurrentAtomSet();
else
@@ -1803,236 +1802,6 @@
}
////////////////////////////////////////////////////////////////
- // incommensurate modulation
- ////////////////////////////////////////////////////////////////
-
-// The occupational distortion of a given atom or rigid group is
-// usually parameterized by Fourier series. Each term of the series
-// commonly adopts two different representations: the sine-cosine
-// form,
-// Pc cos(2\p k r)+Ps sin(2\p k r),
-// and the modulus-argument form,
-// |P| cos(2\p k r+\d),
-// where k is the wave vector of the term and r is the atomic
-// average position. _atom_site_occ_Fourier_param_phase is the phase
-// (\d/2\p) in cycles corresponding to the Fourier term defined by
-// _atom_site_occ_Fourier_atom_site_label and
-// _atom_site_occ_Fourier_wave_vector_seq_id.
-
- private final static int WV_ID = 0;
- private final static int WV_X = 1;
- private final static int WV_Y = 2;
- private final static int WV_Z = 3;
- private final static int FWV_ID = 4;
- private final static int FWV_X = 5;
- private final static int FWV_Y = 6;
- private final static int FWV_Z = 7;
- private final static int FWV_Q1_COEF = 8;
- private final static int FWV_Q2_COEF = 9;
- private final static int FWV_Q3_COEF = 10;
-
- private final static int FWV_DISP_LABEL = 11;
- private final static int FWV_DISP_AXIS = 12;
- private final static int FWV_DISP_ID = 13;
- private final static int FWV_DISP_COS = 14;
- private final static int FWV_DISP_SIN = 15;
- private final static int FWV_DISP_MODULUS = 16;
- private final static int FWV_DISP_PHASE = 17;
-
- private final static int FWV_OCC_LABEL = 18;
- private final static int FWV_OCC_ID = 19;
- private final static int FWV_OCC_COS = 20;
- private final static int FWV_OCC_SIN = 21;
- private final static int FWV_OCC_MODULUS = 22;
- private final static int FWV_OCC_PHASE = 23;
-
- private final static int DISP_SPEC_LABEL = 24;
- private final static int DISP_SAW_AX = 25;
- private final static int DISP_SAW_AY = 26;
- private final static int DISP_SAW_AZ = 27;
- private final static int DISP_SAW_C = 28;
- private final static int DISP_SAW_W = 29;
-
- private final static int OCC_SPECIAL_LABEL = 30;
- private final static int OCC_CRENEL_C = 31;
- private final static int OCC_CRENEL_W = 32;
-
-
- final private static String[] modulationFields = {
- "_cell_wave_vector_seq_id",
- "_cell_wave_vector_x",
- "_cell_wave_vector_y",
- "_cell_wave_vector_z",
- "_atom_site_fourier_wave_vector_seq_id",
- "_atom_site_fourier_wave_vector_x",
- "_atom_site_fourier_wave_vector_y",
- "_atom_site_fourier_wave_vector_z",
- "_jana_atom_site_fourier_wave_vector_q1_coeff",
- "_jana_atom_site_fourier_wave_vector_q2_coeff",
- "_jana_atom_site_fourier_wave_vector_q3_coeff",
- "_atom_site_displace_fourier_atom_site_label",
- "_atom_site_displace_fourier_axis",
- "_atom_site_displace_fourier_wave_vector_seq_id",
- "_atom_site_displace_fourier_param_cos",
- "_atom_site_displace_fourier_param_sin",
- "_atom_site_displace_fourier_param_modulus",
- "_atom_site_displace_fourier_param_phase",
- "_atom_site_occ_fourier_atom_site_label",
- "_atom_site_occ_fourier_wave_vector_seq_id",
- "_atom_site_occ_fourier_param_cos",
- "_atom_site_occ_fourier_param_sin",
- "_atom_site_occ_fourier_param_modulus",
- "_atom_site_occ_fourier_param_phase",
- "_atom_site_displace_special_func_atom_site_label",
- "_atom_site_displace_special_func_sawtooth_ax",
- "_atom_site_displace_special_func_sawtooth_ay",
- "_atom_site_displace_special_func_sawtooth_az",
- "_atom_site_displace_special_func_sawtooth_c",
- "_atom_site_displace_special_func_sawtooth_w",
- "_atom_site_occ_special_func_atom_site_label",
- "_atom_site_occ_special_func_crenel_c",
- "_atom_site_occ_special_func_crenel_w"
- };
-
- /**
- * creates entries in htModulation with a key of
- * the form: type_id_axis;atomLabel where
- * type = W|F|D|O (wave vector, Fourier index, displacement, occupancy);
- * id = 1|2|3|S (Fourier index or special -- sawtooth or crenel);
- * axis (optional) = o|x|y|z (o indicates irrelevant -- occupancy);
- * and ;atomLabel is only for D and O.
- *
- * @throws Exception
- */
- private void processModulationLoopBlock() throws Exception {
- parseLoopParameters(modulationFields);
- int tok;
- while (tokenizer.getData()) {
- boolean ignore = false;
- String id = null;
- String atomLabel = null;
- String axis = null;
- P3 pt = P3.new3(Float.NaN, Float.NaN, Float.NaN);
- float c = Float.NaN;
- float w = Float.NaN;
- for (int i = 0; i < tokenizer.fieldCount; ++i) {
- switch (tok = fieldProperty(i)) {
- case WV_ID:
- case FWV_ID:
- case FWV_DISP_ID:
- case FWV_OCC_ID:
- switch (tok) {
- case WV_ID:
- id = "W_" + field;
- break;
- case FWV_ID:
- id = "F_" + field;
- break;
- case FWV_DISP_ID:
- id = "D_" + field;
- break;
- case FWV_OCC_ID:
- id = "O_" + field;
- break;
- }
- break;
- case DISP_SPEC_LABEL:
- id = "D_S";
- atomLabel = field;
- axis = "0";
- break;
- case OCC_SPECIAL_LABEL:
- id = "O_S";
- axis = "0";
- //$FALL-THROUGH$
- case FWV_DISP_LABEL:
- case FWV_OCC_LABEL:
- atomLabel = field;
- pt.z = 0;
- break;
- case FWV_DISP_AXIS:
- axis = field;
- if (modAxes != null
- && modAxes.indexOf(axis.toUpperCase()) < 0)
- ignore = true;
- break;
- case WV_X:
- case FWV_X:
- case FWV_DISP_COS:
- case FWV_OCC_COS:
- case OCC_CRENEL_C:
- case DISP_SAW_AX:
- pt.x = parseFloatStr(field);
- break;
- case FWV_Q1_COEF:
- id += "_q_";
- pt.x = parseFloatStr(field);
- switch (modDim) {
- case 1:
- pt.y = 0;
- //$FALL-THROUGH$
- case 2:
- pt.z = 0;
- }
- break;
- case FWV_DISP_MODULUS:
- case FWV_OCC_MODULUS:
- pt.x = parseFloatStr(field);
- pt.z = 1;
- break;
- case WV_Y:
- case FWV_Y:
- case FWV_Q2_COEF:
- case FWV_DISP_SIN:
- case FWV_DISP_PHASE:
- case FWV_OCC_SIN:
- case FWV_OCC_PHASE:
- case OCC_CRENEL_W:
- case DISP_SAW_AY:
- pt.y = parseFloatStr(field);
- break;
- case WV_Z:
- case FWV_Z:
- case FWV_Q3_COEF:
- case DISP_SAW_AZ:
- pt.z = parseFloatStr(field);
- break;
- case DISP_SAW_C:
- c = parseFloatStr(field);
- break;
- case DISP_SAW_W:
- w = parseFloatStr(field);
- break;
- }
- if (ignore || Float.isNaN(pt.x + pt.y + pt.z) || id == null)
- continue;
- switch (id.charAt(0)) {
- case 'W':
- case 'F':
- break;
- case 'D':
- case 'O':
- if (atomLabel == null || axis == null)
- continue;
- if (id.equals("D_S")) {
- // saw tooth displacement center/width/Axyz
- if (pt.x != 0)
- addModulation(null, "D_Sx;" + atomLabel, P3.new3(c, w, pt.x));
- if (pt.y != 0)
- addModulation(null, "D_Sy;" + atomLabel, P3.new3(c, w, pt.y));
- if (pt.z != 0)
- addModulation(null, "D_Sz;" + atomLabel, P3.new3(c, w, pt.z));
- continue;
- }
- id += axis + ";" + atomLabel;
- break;
- }
- addModulation(null, id, pt);
- }
- }
- }
-
- ////////////////////////////////////////////////////////////////
// symmetry operations
////////////////////////////////////////////////////////////////
@@ -2430,4 +2199,238 @@
setBs(atoms, i, bsBonds, bs);
}
}
+
+ ////////////////////////////////////////////////////////////////
+ // incommensurate modulation
+ ////////////////////////////////////////////////////////////////
+
+// The occupational distortion of a given atom or rigid group is
+// usually parameterized by Fourier series. Each term of the series
+// commonly adopts two different representations: the sine-cosine
+// form,
+// Pc cos(2\p k r)+Ps sin(2\p k r),
+// and the modulus-argument form,
+// |P| cos(2\p k r+\d),
+// where k is the wave vector of the term and r is the atomic
+// average position. _atom_site_occ_Fourier_param_phase is the phase
+// (\d/2\p) in cycles corresponding to the Fourier term defined by
+// _atom_site_occ_Fourier_atom_site_label and
+// _atom_site_occ_Fourier_wave_vector_seq_id.
+
+ private final static int WV_ID = 0;
+ private final static int WV_X = 1;
+ private final static int WV_Y = 2;
+ private final static int WV_Z = 3;
+ private final static int FWV_ID = 4;
+ private final static int FWV_X = 5;
+ private final static int FWV_Y = 6;
+ private final static int FWV_Z = 7;
+ private final static int FWV_Q1_COEF = 8;
+ private final static int FWV_Q2_COEF = 9;
+ private final static int FWV_Q3_COEF = 10;
+
+ private final static int FWV_DISP_LABEL = 11;
+ private final static int FWV_DISP_AXIS = 12;
+ private final static int FWV_DISP_ID = 13;
+ private final static int FWV_DISP_COS = 14;
+ private final static int FWV_DISP_SIN = 15;
+ private final static int FWV_DISP_MODULUS = 16;
+ private final static int FWV_DISP_PHASE = 17;
+
+ private final static int FWV_OCC_LABEL = 18;
+ private final static int FWV_OCC_ID = 19;
+ private final static int FWV_OCC_COS = 20;
+ private final static int FWV_OCC_SIN = 21;
+ private final static int FWV_OCC_MODULUS = 22;
+ private final static int FWV_OCC_PHASE = 23;
+
+ private final static int DISP_SPEC_LABEL = 24;
+ private final static int DISP_SAW_AX = 25;
+ private final static int DISP_SAW_AY = 26;
+ private final static int DISP_SAW_AZ = 27;
+ private final static int DISP_SAW_C = 28;
+ private final static int DISP_SAW_W = 29;
+
+ private final static int OCC_SPECIAL_LABEL = 30;
+ private final static int OCC_CRENEL_C = 31;
+ private final static int OCC_CRENEL_W = 32;
+
+
+ final private static String[] modulationFields = {
+ "_cell_wave_vector_seq_id",
+ "_cell_wave_vector_x",
+ "_cell_wave_vector_y",
+ "_cell_wave_vector_z",
+ "_atom_site_fourier_wave_vector_seq_id",
+ "_atom_site_fourier_wave_vector_x",
+ "_atom_site_fourier_wave_vector_y",
+ "_atom_site_fourier_wave_vector_z",
+ "_jana_atom_site_fourier_wave_vector_q1_coeff",
+ "_jana_atom_site_fourier_wave_vector_q2_coeff",
+ "_jana_atom_site_fourier_wave_vector_q3_coeff",
+ "_atom_site_displace_fourier_atom_site_label",
+ "_atom_site_displace_fourier_axis",
+ "_atom_site_displace_fourier_wave_vector_seq_id",
+ "_atom_site_displace_fourier_param_cos",
+ "_atom_site_displace_fourier_param_sin",
+ "_atom_site_displace_fourier_param_modulus",
+ "_atom_site_displace_fourier_param_phase",
+ "_atom_site_occ_fourier_atom_site_label",
+ "_atom_site_occ_fourier_wave_vector_seq_id",
+ "_atom_site_occ_fourier_param_cos",
+ "_atom_site_occ_fourier_param_sin",
+ "_atom_site_occ_fourier_param_modulus",
+ "_atom_site_occ_fourier_param_phase",
+ "_atom_site_displace_special_func_atom_site_label",
+ "_atom_site_displace_special_func_sawtooth_ax",
+ "_atom_site_displace_special_func_sawtooth_ay",
+ "_atom_site_displace_special_func_sawtooth_az",
+ "_atom_site_displace_special_func_sawtooth_c",
+ "_atom_site_displace_special_func_sawtooth_w",
+ "_atom_site_occ_special_func_atom_site_label",
+ "_atom_site_occ_special_func_crenel_c",
+ "_atom_site_occ_special_func_crenel_w"
+ };
+
+ /**
+ * creates entries in htModulation with a key of
+ * the form: type_id_axis;atomLabel where
+ * type = W|F|D|O (wave vector, Fourier index, displacement, occupancy);
+ * id = 1|2|3|0|S (Fourier index, Crenel(0), sawtooth);
+ * axis (optional) = 0|x|y|z (0 indicates irrelevant -- occupancy);
+ * and ;atomLabel is only for D and O.
+ *
+ * @throws Exception
+ */
+ private void processModulationLoopBlock() throws Exception {
+ parseLoopParameters(modulationFields);
+ int tok;
+ while (tokenizer.getData()) {
+ boolean ignore = false;
+ String id = null;
+ String atomLabel = null;
+ String axis = null;
+ P3 pt = P3.new3(Float.NaN, Float.NaN, Float.NaN);
+ float c = Float.NaN;
+ float w = Float.NaN;
+ for (int i = 0; i < tokenizer.fieldCount; ++i) {
+ switch (tok = fieldProperty(i)) {
+ case WV_ID:
+ case FWV_ID:
+ pt.x = pt.y = pt.z = 0;
+ //$FALL-THROUGH$
+ case FWV_DISP_ID:
+ case FWV_OCC_ID:
+ switch (tok) {
+ case WV_ID:
+ id = "W_" + field;
+ break;
+ case FWV_ID:
+ id = "F_" + field;
+ break;
+ case FWV_DISP_ID:
+ id = "D_" + field;
+ break;
+ case FWV_OCC_ID:
+ id = "O_" + field;
+ break;
+ }
+ break;
+ case DISP_SPEC_LABEL:
+ id = "D_S";
+ atomLabel = field;
+ axis = "0";
+ break;
+ case OCC_SPECIAL_LABEL:
+ id = "O_0";
+ axis = "0";
+ //$FALL-THROUGH$
+ case FWV_DISP_LABEL:
+ case FWV_OCC_LABEL:
+ atomLabel = field;
+ pt.z = 0;
+ break;
+ case FWV_DISP_AXIS:
+ axis = field;
+ if (modAxes != null
+ && modAxes.indexOf(axis.toUpperCase()) < 0)
+ ignore = true;
+ break;
+ case WV_X:
+ case FWV_X:
+ case FWV_DISP_COS:
+ case FWV_OCC_COS:
+ case OCC_CRENEL_C:
+ case DISP_SAW_AX:
+ pt.x = parseFloatStr(field);
+ break;
+ case FWV_Q1_COEF:
+ id += "_q_";
+ pt.x = parseFloatStr(field);
+ switch (modDim) {
+ case 1:
+ pt.y = 0;
+ //$FALL-THROUGH$
+ case 2:
+ pt.z = 0;
+ }
+ break;
+ case FWV_DISP_MODULUS:
+ case FWV_OCC_MODULUS:
+ pt.x = parseFloatStr(field);
+ pt.z = 1;
+ break;
+ case WV_Y:
+ case FWV_Y:
+ case FWV_Q2_COEF:
+ case FWV_DISP_SIN:
+ case FWV_DISP_PHASE:
+ case FWV_OCC_SIN:
+ case FWV_OCC_PHASE:
+ case OCC_CRENEL_W:
+ case DISP_SAW_AY:
+ pt.y = parseFloatStr(field);
+ break;
+ case WV_Z:
+ case FWV_Z:
+ case FWV_Q3_COEF:
+ case DISP_SAW_AZ:
+ pt.z = parseFloatStr(field);
+ break;
+ case DISP_SAW_C:
+ c = parseFloatStr(field);
+ break;
+ case DISP_SAW_W:
+ w = parseFloatStr(field);
+ break;
+ }
+ if (ignore || Float.isNaN(pt.x + pt.y + pt.z) || id == null)
+ continue;
+ switch (id.charAt(0)) {
+ case 'W':
+ case 'F':
+ break;
+ case 'D':
+ case 'O':
+ if (atomLabel == null || axis == null)
+ continue;
+ if (id.equals("D_S")) {
+ // saw tooth displacement center/width/Axyz
+ if (pt.x != 0)
+ addModulation(null, "D_Sx;" + atomLabel, P3.new3(c, w, pt.x));
+ if (pt.y != 0)
+ addModulation(null, "D_Sy;" + atomLabel, P3.new3(c, w, pt.y));
+ if (pt.z != 0)
+ addModulation(null, "D_Sz;" + atomLabel, P3.new3(c, w, pt.z));
+ continue;
+ }
+ id += axis + ";" + atomLabel;
+ break;
+ }
+ addModulation(null, id, pt);
+ }
+ }
+ }
+
+
}
Modified: trunk/Jmol/src/org/jmol/adapter/readers/cif/ModulationReader.java
===================================================================
--- trunk/Jmol/src/org/jmol/adapter/readers/cif/ModulationReader.java
2013-08-08 14:18:07 UTC (rev 18530)
+++ trunk/Jmol/src/org/jmol/adapter/readers/cif/ModulationReader.java
2013-08-08 23:30:00 UTC (rev 18531)
@@ -30,6 +30,7 @@
import java.util.Map;
import java.util.Map.Entry;
+import org.jmol.util.BSUtil;
import org.jmol.util.Logger;
import org.jmol.util.Matrix3f;
import org.jmol.util.Modulation;
@@ -61,11 +62,11 @@
protected boolean incommensurate;
protected Atom[] atoms;
- private JmolList<float[]> lattvecs;
- private V3 modT;
+ private V3 q1Norm;
private Matrix3f rot;
private Map<String, P3> htModulation;
private Map<String, JmolList<Modulation>> htAtomMods;
+ private int modT;
protected void initializeMod() throws Exception {
modAxes = getFilter("MODAXES=");
@@ -76,28 +77,6 @@
allowRotations = !checkFilterKey("NOSYM");
}
-
- protected void addLatticeVector(String data) {
- if (lattvecs == null)
- lattvecs = new JmolList<float[]>();
- float[] a = getTokensFloat(data, null, modDim + 3);
- boolean isOK = false;
- for (int i = a.length; --i >= 0 && !isOK;) {
- if (Float.isNaN(a[i]))
- return;
- isOK |= (a[i] != 0);
- }
- if (isOK)
- lattvecs.addLast(a);
- }
-
- protected void finalizeIncommensurate() {
- if (incommensurate)
-
atomSetCollection.setBaseSymmetryAtomCount(atomSetCollection.getAtomCount());
- if (lattvecs != null)
- atomSetCollection.getSymmetry().addLatticeVectors(lattvecs);
- }
-
protected void setModDim(String token) {
modDim = parseIntStr(token);
if (modAverage)
@@ -117,19 +96,12 @@
protected P3 getModulationVector(String id) {
return htModulation.get(id);
}
-
+
protected void addModulation(Map<String, P3> map, String id, P3 pt) {
if (map == null)
map = htModulation;
map.put(id, pt);
Logger.info("Adding " + id + " " + pt);
- if (id.charAt(0) == 'W' || id.charAt(0) == 'F') {
- appendLoadNote("Wave vector " + id +" = " + pt);
- if (id.equals("W_1")) {
- modT = V3.newV(pt);
- modT.normalize();
- }
- }
}
/**
@@ -139,13 +111,25 @@
protected void setModulation() {
if (!incommensurate || htModulation == null)
return;
+ int n = atomSetCollection.getAtomCount();
Map<String, P3> map = new Hashtable<String, P3>();
for (Entry<String, P3> e : htModulation.entrySet()) {
String key = e.getKey();
P3 pt = e.getValue();
switch (key.charAt(0)) {
+ case 'W':
+ appendLoadNote("Wave vector " + key + " = " + pt);
+
+ if (key.equals("W_1")) {
+ q1Norm = V3.newV(pt);
+ q1Norm.normalize();
+ }
+ break;
+ case 'O':
+ if (!modVib && atomSetCollection.bsAtoms == null)
+ atomSetCollection.bsAtoms = BSUtil.newBitSet2(0, n);
+ //$FALL-THROUGH$
case 'D':
- case 'O':
// fix modulus/phase option only for non-special modulations;
if (pt.z == 1 && key.charAt(2) != 'S') {
// modulus/phase M cos(2pi(q.r) + 2pi(p))
@@ -171,6 +155,8 @@
pf.scaleAdd2(pt.z, htModulation.get("W_3"), pf);
key = TextFormat.simpleReplace(key, "_q_", "");
addModulation(map, key, pf);
+ } else {
+ appendLoadNote("Wave vector " + key + " = " + pt);
}
break;
}
@@ -180,30 +166,31 @@
boolean haveAtomMods = false;
for (Entry<String, P3> e : htModulation.entrySet()) {
String key = e.getKey();
- switch (key.charAt(0)) {
+ P3 params = e.getValue();
+ String label = key.substring(key.indexOf(";") + 1);
+ System.out.println(key);
+ int type = key.charAt(0);
+ switch (type) {
case 'O':
- // TODO
- break;
case 'D':
+ char id = key.charAt(2);
char axis = key.charAt(3);
- if (key.charAt(2) == 'S') {
- // TODO -- sawtooth, so now what?
- } else {
- if (htAtomMods == null)
- htAtomMods = new Hashtable<String, JmolList<Modulation>>();
- P3 coefs = e.getValue();
- String label = key.substring(key.indexOf(";") + 1);
- int n = key.charAt(2) - '0';
- key = "F_" + n;
- //TODO -- THIS IS WRONG. n is just a label. It will break in 2D
- // but I don't know to determine "n" any other way in a CIF
file.
- P3 nq = htModulation.get(key);
- JmolList<Modulation> list = htAtomMods.get(label);
- if (list == null)
- htAtomMods.put(label, list = new JmolList<Modulation>());
- list.addLast(new Modulation(nq, n, axis, coefs));
- haveAtomMods = true;
- }
+ type = (id == 'S' ? Modulation.TYPE_DISP_SAWTOOTH
+ : id == '0' ? Modulation.TYPE_OCC_CRENEL
+ : type == 'O' ? Modulation.TYPE_OCC_FOURIER
+ : Modulation.TYPE_DISP_FOURIER);
+ if (htAtomMods == null)
+ htAtomMods = new Hashtable<String, JmolList<Modulation>>();
+ int fn = (id == 'S' ? -1 : id - '0');
+ key = (fn <= 0 ? "W_1" : "F_" + fn);
+ //TODO -- THIS IS WRONG. n is just a label. It will break in 2D
+ // but I don't know to determine "n" any other way in a CIF
file.
+ P3 nq = htModulation.get(key);
+ JmolList<Modulation> list = htAtomMods.get(label);
+ if (list == null)
+ htAtomMods.put(label, list = new JmolList<Modulation>());
+ list.addLast(new Modulation(nq, axis, type, fn, params));
+ haveAtomMods = true;
break;
}
}
@@ -213,16 +200,17 @@
symmetry = atomSetCollection.getSymmetry();
rot = new Matrix3f();
SB sb = new SB();
- int n = atomSetCollection.getAtomCount();
for (int i = 0; i < n; i++)
modulateAtom(atoms[i], sb);
atomSetCollection.setAtomSetAtomProperty("modt", sb.toString(), -1);
}
/**
- * Modulation generally involves x4 = q.r + t. Here we arbitrarily set t =
0,
- * but t could be a FILTER option MODT=n. There would need to be one t per
dimension.
- * The displacement will be set as the atom vibration vector.
+ * The displacement will be set as the atom vibration vector; the string
buffer
+ * will be appended with the t value for a given unit cell.
+ *
+ * Modulation generally involves x4 = q.r + t. Here we arbitrarily set t =
modT = 0,
+ * but modT could be a FILTER option MODT=n. There would need to be one modT
per dimension.
*
* @param a
* @param sb
@@ -236,21 +224,30 @@
iop = 0;
float epsilon = symmetry.getModParam(iop, 0);
float delta = symmetry.getModParam(iop, 1);
+ delta -= modT;
symmetry.getSpaceGroupOperation(iop).getRotationScale(rot);
+ System.out.println("=========MR i=" + a.index + " " + a.atomName + " " + a
+ " " + a.occupancy);
+ System.out.println("op=" + (iop + 1) + " " +
symmetry.getSpaceGroupXyz(iop, false) + " ep=" + epsilon + " de=" + delta);
a.vib = new V3();
- Modulation.modulateAtom(list, a, epsilon, delta, rot, a.vib);
- //System.out.println("=========MR i=" + i + " " + a.atomName + " " + a);
- //System.out.println("op=" + (iop + 1) + " " +
symmetry.getSpaceGroupXyz(iop, false) + " ep=" + epsilon + " de=" + delta);
- //System.out.println("a.vib(abc)=" + a.vib);
+ float val = Modulation.modulateAtom(list, a, epsilon, delta, rot, a.vib);
+ if (val < 0.5f) {
+ a.occupancy = 0;
+ System.out.println("occup 0 for " + a.index + " a=" + a + " " +
a.atomName);
+ if (!modVib)
+ atomSetCollection.bsAtoms.clear(a.index);
+ }
+ System.out.println("a.vib(abc)=" + a.vib);
+
+
// set property_modT to be Math.floor (q.r/|q|) -- really only for d=1
- float t = modT.dot(a);
+ float t = q1Norm.dot(a);
if (Math.abs(t - (int) t) > 0.001f)
t = (int) Math.floor(t);
sb.append(((int) t) + "\n");
- // displace the atom if not filter "MODVIB"
+ // displace the atom and reverse the vector only if not filter "MODVIB"
if (!modVib) {
a.add(a.vib);
a.vib.scale(-1);
Modified: trunk/Jmol/src/org/jmol/adapter/readers/xtal/JanaReader.java
===================================================================
--- trunk/Jmol/src/org/jmol/adapter/readers/xtal/JanaReader.java
2013-08-08 14:18:07 UTC (rev 18530)
+++ trunk/Jmol/src/org/jmol/adapter/readers/xtal/JanaReader.java
2013-08-08 23:30:00 UTC (rev 18531)
@@ -29,6 +29,7 @@
import org.jmol.adapter.readers.cif.ModulationReader;
import org.jmol.adapter.smarter.Atom;
import org.jmol.io.JmolBinary;
+import org.jmol.util.JmolList;
import org.jmol.util.Logger;
import org.jmol.util.P3;
import org.jmol.util.TextFormat;
@@ -41,7 +42,9 @@
public class JanaReader extends ModulationReader {
private String m40Data;
+ private JmolList<float[]> lattvecs;
private String[] tokens;
+
@Override
public void initializeReader() throws Exception {
setFractionalCoordinates(true);
@@ -105,7 +108,7 @@
ndim();
break;
case LATT:
- lattvec();
+ lattvec(line.substring(8));
break;
case SPG:
setSpaceGroupName(getTokens()[1]);
@@ -127,12 +130,68 @@
@Override
public void finalizeReader() throws Exception {
readM40Data();
- finalizeIncommensurate();
+ if (lattvecs != null)
+ atomSetCollection.getSymmetry().addLatticeVectors(lattvecs);
applySymmetryAndSetTrajectory();
setModulation();
finalizeReaderASCR();
}
+ private void cell() throws Exception {
+ for (int ipt = 0; ipt < 6; ipt++)
+ setUnitCellItem(ipt, parseFloat());
+ }
+
+ private void ndim() {
+ setModDim(line.substring(line.length() - 1));
+ }
+
+ private int qicount;
+
+ private void qi() {
+ P3 pt = P3.new3(parseFloat(), parseFloat(), parseFloat());
+ if (qicount == 0)
+ addModulation(null, "W_1", pt);
+ addModulation(null, "F_" + (++qicount), pt);
+ }
+ private void lattvec(String data) throws Exception {
+ float[] a;
+ char c = data.charAt(0);
+ switch(c) {
+ case 'P':
+ case 'X':
+ return;
+ case 'A':
+ case 'B':
+ case 'C':
+ case 'I':
+ a = new float[] {0.5f, 0.5f, 0.5f};
+ if (c != 'I')
+ a[c - 'A'] = 0;
+ break;
+ case 'F':
+ lattvec("A");
+ lattvec("B");
+ lattvec("C");
+ return;
+ case '0': // X explicit
+ if (data.indexOf(".") < 0)
+ return; // lattvec 0 0 0 unnecessary
+ a = getTokensFloat(data, null, modDim + 3);
+ break;
+ default:
+ appendLoadNote(line + " not supported");
+ return;
+ }
+ if (lattvecs == null)
+ lattvecs = new JmolList<float[]>();
+ lattvecs.addLast(a);
+ }
+
+ private void symmetry() throws Exception {
+ setSymmetryOperator(TextFormat.simpleReplace(line.substring(9).trim()," ",
","));
+ }
+
private final String labels = "xyz";
// 12 0 0 1
@@ -194,31 +253,5 @@
return line;
}
- private int qicount;
- private void qi() {
- P3 pt = P3.new3(parseFloat(), parseFloat(), parseFloat());
- if (qicount == 0)
- addModulation(null, "W_1", pt);
- addModulation(null, "F_" + (++qicount), pt);
- }
-
- private void ndim() {
- setModDim(line.substring(line.length() - 1));
- }
-
- private void lattvec() throws Exception {
- addLatticeVector(line.substring(8));
- }
-
- private void symmetry() throws Exception {
- setSymmetryOperator(TextFormat.simpleReplace(line.substring(9).trim()," ",
","));
- }
-
- private void cell() throws Exception {
- for (int ipt = 0; ipt < 6; ipt++)
- setUnitCellItem(ipt, parseFloat());
- }
-
-
}
Modified: trunk/Jmol/src/org/jmol/symmetry/SpaceGroup.java
===================================================================
--- trunk/Jmol/src/org/jmol/symmetry/SpaceGroup.java 2013-08-08 14:18:07 UTC
(rev 18530)
+++ trunk/Jmol/src/org/jmol/symmetry/SpaceGroup.java 2013-08-08 23:30:00 UTC
(rev 18531)
@@ -1385,13 +1385,14 @@
int nOps = operationCount;
for (int j = 0; j < lattvecs.size(); j++) {
float[] data = lattvecs.get(j);
- if (data.length != modulationDimension + 3)
+ if (data.length > modulationDimension + 3)
return;
for (int i = 0; i < nOps; i++) {
SymmetryOperation op = operations[i];
float[] rotTrans = op.rotTransMatrix;
SymmetryOperation newOp = new SymmetryOperation(null, null, 0, 0,
doNormalize);
+ newOp.modDim = modulationDimension;
newOp.rotTransMatrix = ArrayUtil.arrayCopyF(rotTrans, -1);
newOp.setFromMatrix(data, false);
newOp.xyzOriginal = newOp.xyz;
Modified: trunk/Jmol/src/org/jmol/symmetry/SymmetryOperation.java
===================================================================
--- trunk/Jmol/src/org/jmol/symmetry/SymmetryOperation.java 2013-08-08
14:18:07 UTC (rev 18530)
+++ trunk/Jmol/src/org/jmol/symmetry/SymmetryOperation.java 2013-08-08
23:30:00 UTC (rev 18531)
@@ -71,7 +71,7 @@
V3 originalTranslation;
private String[] myLabels;
- private int modDim;
+ int modDim;
float[] rotTransMatrix;
@@ -245,8 +245,6 @@
boolean setFromMatrix(float[] offset, boolean isReverse) {
float v = 0;
int pt = 0;
- if (offset != null)
- modDim = offset.length - 3;
myLabels = (modDim == 0 ? labelsXYZ : labelsX1_6);
for (int i = 0; i < rotTransMatrix.length; i++) {
if (Float.isNaN(rotTransMatrix[i]))
@@ -255,7 +253,7 @@
if (Math.abs(v) < 0.00001f)
v = 0;
if (i % 4 == 3) {
- if (offset != null)
+ if (offset != null && pt < offset.length)
v = v / 12 + offset[pt++];
setOrigTranslation(i, v);
v = normalizeTwelfths((v < 0 ? -1 : 1) * Math.round(Math.abs(v * 12))/
12f, doNormalize);
Modified: trunk/Jmol/src/org/jmol/util/Modulation.java
===================================================================
--- trunk/Jmol/src/org/jmol/util/Modulation.java 2013-08-08 14:18:07 UTC
(rev 18530)
+++ trunk/Jmol/src/org/jmol/util/Modulation.java 2013-08-08 23:30:00 UTC
(rev 18531)
@@ -1,110 +1,257 @@
package org.jmol.util;
/**
- * A class to allow for more complex vibrations and associated
- * phenomena, such as modulated crystals.
+ * A class to allow for more complex vibrations and associated phenomena, such
+ * as modulated crystals, including Fourier series, Crenel functions, and
+ * sawtooth functions
*
- * @author Bob Hanson hans...@stolaf.edu
+ * @author Bob Hanson hans...@stolaf.edu 8/8/2013
*
*/
public class Modulation extends V3 {
- private final static double TWOPI = 2 * Math.PI;
-
- private double ccos;
- private double csin;
+ private static final double TWOPI = 2 * Math.PI;
+
+ private double a1;
+ private double a2;
+ private double center;
+ private double left, right;
+
private V3 nq; // wave vector
- private int n; // power
+ private int fn; // power
private char axis;
+ private final int type;
-
+ public static final int TYPE_DISP_FOURIER = 0;
+ public static final int TYPE_DISP_SAWTOOTH = 1;
+ public static final int TYPE_OCC_FOURIER = 2;
+ public static final int TYPE_OCC_CRENEL = 3;
+
public V3 getWaveVector() {
return nq;
}
-
+
/**
* Each atomic modulation involves a fractional coordinate wave vector q, a
- * Fourier power n, a modulation axis (x, y, or, z), and specified
- * coefficients for cos and sin.
+ * Fourier power n, a modulation axis (x, y, or, z), and specified parameters
+ * that depend upon the type of function.
*
*
* @param nq
- * @param n
* @param axis
- * @param coefs
+ * @param type
+ * @param fn
+ * @param params
*/
- public Modulation(P3 nq, int n, char axis, P3 coefs) {
+ public Modulation(P3 nq, char axis, int type, int fn, P3 params) {
this.nq = V3.newV(nq);
- //this.q.scale(1f/n); leave this nq.
- this.n = n;
this.axis = axis;
- this.ccos = coefs.x;
- this.csin = coefs.y;
+ this.type = type;
+ switch (type) {
+ case TYPE_DISP_FOURIER:
+ case TYPE_OCC_FOURIER:
+ this.fn = fn;
+ a1 = params.x;
+ a2 = params.y;
+ break;
+ case TYPE_DISP_SAWTOOTH:
+ case TYPE_OCC_CRENEL:
+ this.fn = 1;
+ center = params.x;
+ float width = params.y;
+ if (width > 1)
+ width = 1; //
http://b-incstrdb.ehu.es/incstrdb/CIFFile.php?RefCode=Bi-Sr-Ca-Cu-O_rNdCbetq
+ left = center - width / 2;
+ right = center + width / 2;
+ if (left < 0)
+ left += 1;
+ if (right > 1)
+ right -= 1;
+ if (left >= right && left - right < 0.01f)
+ left = right + 0.01f;
+ a1 = 2 * params.z / width;
+ break;
+ }
}
-
+
/**
- * Starting with fractional coordinates, determine the overall modulation
displacement.
+ * Starting with fractional coordinates, determine the overall modulation.
*
- * @param mods a given atom's modulations
- * @param r fractional xyz
- * @param epsilon as in x4´ = epsilon x4 + delta
- * @param delta as in x4´ = epsilon x4 + delta
- * @param rot symmetry rotation to be applied
- * @param d displacement return
+ * @param mods
+ * a given atom's modulations
+ * @param r
+ * fractional xyz
+ * @param epsilon
+ * as in x4´ = epsilon x4 + delta
+ * @param delta
+ * as in x4´ = epsilon x4 + delta
+ * @param rot
+ * symmetry rotation to be applied
+ * @param d
+ * displacement return
+ * @return crenel value
*/
- public static void modulateAtom(JmolList<Modulation> mods, Tuple3f r, float
epsilon, float delta, Matrix3f rot, V3 d) {
+ public static float modulateAtom(JmolList<Modulation> mods, Tuple3f r,
+ float epsilon, float delta, Matrix3f rot,
+ V3 d) {
d.set(0, 0, 0);
- for (int i = mods.size(); --i >= 0;)
- mods.get(i).apply(r, epsilon, delta, d);
- rot.transform(d);
+ float v = Float.NaN;
+ for (int i = mods.size(); --i >= 0;) {
+ float f = mods.get(i).apply(r, epsilon, delta, d);
+ if (!Float.isNaN(f)) {
+ if (Float.isNaN(v))
+ v = 0;
+ v += f;
+ }
+ }
+ rot.transform(d);
+ return v;
}
/**
*
* In general, we have:
- *
- * u_axis(x) = sum[A1 cos(theta) + B1 sin(theta)]
*
+ * u_axis(x) = sum[A1 cos(theta) + B1 sin(theta)]
+ *
* where axis is x, y, or z, and theta = 2n pi x
*
- * However, for symmetry-related atoms, we need to do a 4D transformation,
- * not just a 3D one. We need to operate on x first BEFORE applying u(x):
+ * However, for symmetry-related atoms, we need to do a 4D transformation,
not
+ * just a 3D one. We need to operate on x first BEFORE applying u(x):
*
- * u(x4') = Ru(x4)
+ * u(x4') = Ru(x4)
*
- * where we need to express x4 in terms of a "transformed" x4', because
- * x4' is for our rotated point.
+ * where we need to express x4 in terms of a "transformed" x4', because x4'
is
+ * for our rotated point.
*
- * x4' = epsilon x4 + delta
+ * x4' = epsilon x4 + delta
*
* where epsilon = +/-1, so
*
- * x4 = epsilon (x4' - delta)
+ * x4 = epsilon (x4' - delta)
*
* More generally, we might have something like:
*
- * x4' = x5 + 1/2; x5' = x4 - 1/2
- *
+ * x4' = x5 + 1/2; x5' = x4 - 1/2
+ *
* Will have to work on that later!
*
*
- * @param r
- * @param epsilon
- * @param delta
- * @param vecMod
+ * @param r
+ * @param epsilon
+ * @param delta
+ * @param vecMod
+ * @return crenel value if that type
*
*/
- private void apply(Tuple3f r, float epsilon, float delta, V3 vecMod) {
-
+ private float apply(Tuple3f r, float epsilon, float delta, V3 vecMod) {
+
// TODO: must be adapted for d > 1 modulation
-
- double theta = TWOPI * (epsilon * (nq.dot(r) + n * delta));
+
double v = 0;
- if (ccos != 0)
- v += ccos * Math.cos(theta);
- if (csin != 0)
- v += csin * Math.sin(theta);
+ //if (type == TYPE_OCC_CRENEL)
+ //delta = 0;
+
+ double x4 = (epsilon * (nq.dot(r) - fn * delta));
+
+ switch (type) {
+ case TYPE_DISP_FOURIER:
+ case TYPE_OCC_FOURIER:
+ double theta = TWOPI * x4;
+ if (a1 != 0)
+ v += a1 * Math.cos(theta);
+ if (a2 != 0)
+ v += a2 * Math.sin(theta);
+ break;
+ case TYPE_OCC_CRENEL:
+
+ // An occupational crenel function along the internal space is
+ // defined as follows:
+ //
+ // p(x4)=1 if x4 belongs to the interval [c-w/2,c+w/2]
+ // p(x4)=0 if x4 is outside the interval [c-w/2,c+w/2],
+
+ x4 -= Math.floor(x4);
+ return (range(x4) ? 1 : 0);
+
+ case TYPE_DISP_SAWTOOTH:
+
+ // _atom_site_displace_special_func_sawtooth_ items are the
+ // adjustable parameters of a sawtooth function. A displacive sawtooth
+ // function along the internal space is defined as follows:
+ //
+ // u_x = 2a_x[(x4 − c)/w]
+ //
+ // for x4 belonging to the interval [c − (w/2), c + (w/2)], where ax,
+ // ay and az are the amplitudes (maximum displacements) along each
+ // crystallographic axis, w is its width, x4 is the internal coordinate
+ // and c is the centre of the function in internal space. ux, uy and
+ // uz must be expressed in relative units.
+
+ // here we have set a1 = 2a_xyz/w
+
+ x4 -= Math.floor(x4);
+ if (!range(x4))
+ return Float.NaN;
+
+ // x < L < c
+ //
+ // /|
+ // / |
+ // / x------------->
+ // | | L /|
+ // --------+--|---|----c-+------
+ // 0 R | / 1
+ // | /
+ // | /
+ // |/
+
+ // becomes
+
+ // /|
+ // / |
+ // / x|
+ // L / |
+ // --------+------|----c-+--|----
+ // 0 | / 1 R
+ // | /
+ // | /
+ // |/
+
+ // x > R > c
+ //
+ // /|
+ // / |
+ // / |
+ // / | L
+ // --------+c----|---------|---+-------
+ // 0 R | 1
+ // <-----------------x /
+ // | /
+ // |/
+
+ // becomes
+
+ // /|
+ // / |
+ // / |
+ // L / |
+ // ----|---+c----|-------------+--------
+ // | 0 R 1
+ // |x /
+ // | /
+ // |/
+
+ if (left > right) {
+ if (x < left && left < center)
+ x += 1;
+ else if (x > right && right > center)
+ x -= 1;
+ }
+ v = a1 * (x4 - center);
+ break;
+ }
switch (axis) {
case 'x':
vecMod.x += v;
@@ -115,8 +262,23 @@
case 'z':
vecMod.z += v;
break;
+ default:
+ return (float) v;
}
+ return Float.NaN;
//System.out.println("MOD q=" + nq + " r=" + r + " axis=" + axis + "
theta=" + theta + " ccos=" + ccos + " csin=" + csin + " delta=" + delta + " v="
+ v);
}
+ /**
+ * Check that left < x4 < right, but allow for folding
+ *
+ * @param x4
+ * @return true only if x4 is in the (possibly folded) range of left and
right
+ *
+ */
+ private boolean range(double x4) {
+ return (left < right ? left <= x4 && x4 <= right : left <= x4
+ || x4 <= right);
+ }
+
}
Modified: trunk/Jmol/src/org/jmol/viewer/Jmol.properties
===================================================================
--- trunk/Jmol/src/org/jmol/viewer/Jmol.properties 2013-08-08 14:18:07 UTC
(rev 18530)
+++ trunk/Jmol/src/org/jmol/viewer/Jmol.properties 2013-08-08 23:30:00 UTC
(rev 18531)
@@ -13,7 +13,7 @@
___JmolVersion="13.3.4_dev_2013.08.08"
-code: incommensurate crystal work
+code: incommensurate crystal work -- occupancy Crenel, displacive sawtooths
bug fix: shapeInfo not reporting visibility of isosurface
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
------------------------------------------------------------------------------
Get 100% visibility into Java/.NET code with AppDynamics Lite!
It's a free troubleshooting tool designed for production.
Get down to code-level detail for bottlenecks, with <2% overhead.
Download for free and get started troubleshooting in minutes.
http://pubads.g.doubleclick.net/gampad/clk?id=48897031&iu=/4140/ostg.clktrk
_______________________________________________
Jmol-commits mailing list
Jmol-commits@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/jmol-commits
