Revision: 5206
Author: hansonr
Date: 2006-06-06 15:57:31 -0700 (Tue, 06 Jun 2006)
ViewCVS: http://svn.sourceforge.net/jmol/?rev=5206&view=rev
Log Message:
-----------
bob200603 real, calculated orbitals!!!! :)
Modified Paths:
--------------
branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java
branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java
Modified: branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java
===================================================================
--- branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java 2006-06-06
19:36:48 UTC (rev 5205)
+++ branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java 2006-06-06
22:57:31 UTC (rev 5206)
@@ -3980,6 +3980,8 @@
int fileIndexPt = 0;
int signPt = 0;
int contourPt = 0;
+ int orbitalPt = 0;
+ float[] nlmZR = new float[5];
viewer.loadShape(JmolConstants.SHAPE_ISOSURFACE);
viewer.setShapeProperty(JmolConstants.SHAPE_ISOSURFACE, "init", null);
boolean colorSeen = false;
@@ -4019,6 +4021,11 @@
colorRangeStage = 1;
break;
} else if (str.equalsIgnoreCase("phase")) {
+ if (orbitalPt != 0) {
+ propertyName = "phase";
+ propertyValue = "_orb";
+ break;
+ }
colorByPhase = true;
break;
} else if (str.equalsIgnoreCase("fileIndex")) {
@@ -4039,6 +4046,9 @@
} else if (str.equalsIgnoreCase("gridPoints")) {
propertyName = "gridPoints";
break;
+ } else if (str.equalsIgnoreCase("orbital")) {
+ orbitalPt = i + nlmZR.length;
+ break;
} else if (str.equalsIgnoreCase("contour")) {
propertyName = "contour";
contourPt = i + 1;
@@ -4073,13 +4083,26 @@
}
break;
case Token.decimal:
- if (colorRangeStage == 0) {
+ if (colorRangeStage == 0 && orbitalPt == 0) {
propertyName = "cutoff";
propertyValue = token.value;
break;
}
// fall into
case Token.integer:
+ if (orbitalPt > 0 && i <= orbitalPt) {
+ int ipt = nlmZR.length - (orbitalPt - i) - 1;
+ if (ipt < 3 && token.tok != Token.integer)
+ integerExpected();
+ // n l m Znuc ptsPerAngstrom
+ nlmZR[ipt] = floatParameter(i);
+ if (i != orbitalPt)
+ break;
+ propertyName = "orbital";
+ propertyValue = nlmZR;
+ orbitalPt = Integer.MIN_VALUE;
+ break;
+ }
if (i == contourPt) {
if (token.tok != Token.integer)
integerExpected();
Modified: branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java
===================================================================
--- branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java 2006-06-06
19:36:48 UTC (rev 5205)
+++ branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java 2006-06-06
22:57:31 UTC (rev 5206)
@@ -101,6 +101,7 @@
final static float defaultMappedDataMin = 0f;
final static float defaultMappedDataMax = 1.0f;
final static float defaultCutoff = 0.02f;
+ final static float defaultOrbitalCutoff = 0.02f;
final static int defaultContourCount = 10;
final static String colorScheme = "roygb";
final static int nContourMax = 100;
@@ -113,6 +114,7 @@
int colorPtr;
boolean colorByPhase;
int colorPhase;
+ boolean isOrbital;
int edgeFractionBase;
int edgeFractionRange;
@@ -154,7 +156,6 @@
int nContours;
int thisContour;
boolean dotsOnly;
- boolean loadAll;
BufferedReader br;
@@ -181,7 +182,7 @@
colorNeg = defaultColorNegative;
colorPos = defaultColorPositive;
colorByPhase = false;
- loadAll = false;
+ isOrbital = false;
cutoff = Float.MAX_VALUE;
super.setProperty("meshID", null, null);
return;
@@ -267,13 +268,32 @@
return;
}
applyColorScale();
- currentMesh.visible = true;
+ if (currentMesh != null)
+ currentMesh.visible = true;
return;
}
+ if ("orbital" == propertyName) {
+ isOrbital = true;
+ float[] nlmZR = (float[]) value;
+ psi_n = (int) nlmZR[0];
+ psi_l = (int) nlmZR[1];
+ psi_m = (int) nlmZR[2];
+ psi_Znuc = nlmZR[3];
+ psi_ptsPerAngstrom = nlmZR[4];
+ if (psi_ptsPerAngstrom <= 0 || psi_Znuc <= 0 || Math.abs(psi_m) > psi_l
+ || psi_l >= psi_n) {
+ System.out.println("invalid n, l, m, Znuc, or ptsPerAngstrom : "
+ + psi_n + " " + psi_l + " " + psi_m + " " + psi_Znuc + " "
+ + psi_ptsPerAngstrom);
+ return;
+ }
+ propertyName = "bufferedReader";
+ }
if ("bufferedReader" == propertyName) {
+ if (!isOrbital)
+ br = (BufferedReader) value;
System.out.println("loading voxel data...");
- br = (BufferedReader) value;
checkFlags();
createIsosurface();
initializeMesh(force2SidedTriangles);
@@ -405,7 +425,7 @@
currentMesh.hasGridPoints = true;
}
if (cutoff == Float.MAX_VALUE)
- cutoff = defaultCutoff;
+ cutoff = (isOrbital ? defaultOrbitalCutoff : defaultCutoff);
currentMesh.jvxlSurfaceData = "";
currentMesh.jvxlEdgeData = "";
currentMesh.jvxlColorData = "";
@@ -473,11 +493,15 @@
int readVolumetricHeader() {
try {
- readTitleLines();
- readAtomCountAndOrigin();
- for (int i = 0; i < 3; ++i) {
- readVoxelVector(i);
- System.out.println("cube axes vector:" + volumetricVectors[i]);
+ if (isOrbital) {
+ setupOrbital();
+ } else {
+ readTitleLines();
+ readAtomCountAndOrigin();
+ for (int i = 0; i < 3; ++i) {
+ readVoxelVector(i);
+ System.out.println("cube axes vector:" + volumetricVectors[i]);
+ }
}
setupMatrix(volumetricMatrix, volumetricVectors);
readAtoms();
@@ -662,7 +686,9 @@
*
*/
- if (isMapData || thePlane == null)
+ if (isOrbital)
+ voxelValue = getPsi2(x, y, z);
+ else if (isMapData || thePlane == null)
voxelValue = getNextVoxelValue();
else
voxelValue = calcVoxelPlaneDistance(x, y, z);
@@ -1073,7 +1099,7 @@
}
final static String[] colorPhases = { "x", "y", "z", "xy", "yz", "xz",
- "x2-y2", "z2" };
+ "x2-y2", "z2", "_orb" };
float getPhase(Point3f pt) {
switch (colorPhase) {
@@ -1092,7 +1118,9 @@
case 6:
return (pt.x * pt.x - pt.y * pt.y > 0 ? 1 : -1);
case 7:
- return (pt.z * pt.z - pt.x * pt.x - pt.y * pt.y > 0 ? 1 : -1);
+ return (pt.z * pt.z * 2f - pt.x * pt.x - pt.y * pt.y > 0 ? 1 : -1);
+ case 8:
+ return (hydrogenAtomPsiAt(pt, psi_n, psi_l, psi_m) > 0 ? 1 : -1);
}
return 1;
}
@@ -1428,7 +1456,7 @@
compressedData += jvxlCompressString(mesh.jvxlEdgeData
+ mesh.jvxlColorData);
int nColor = (mesh.jvxlColorData.length() - 1);
- if (isJvxlColorPrecision && nColor != -1)
+ if (isJvxlColorPrecision && nColor != -1)
nColor = -nColor;
data += mesh.jvxlSurfaceData.length() + " "
+ (mesh.jvxlEdgeData.length() - 1) + " " + nColor;
@@ -1444,8 +1472,8 @@
data += " CONTOUR PLANE " + mesh.jvxlPlane;
if (!isJvxl)
data += " approximate compressionRatio="
- + (int) (((float)mesh.nBytes + mesh.jvxlFileHeader.length()) /
- (data.length() + compressedData.length())) + ":1";
+ + (int) (((float) mesh.nBytes + mesh.jvxlFileHeader.length()) / (data
+ .length() + compressedData.length())) + ":1";
data += "\n" + compressedData;
if (msg != null)
data += "#-------end of jvxl file data-------\n";
@@ -1602,7 +1630,8 @@
continue;
}
++surfaceCount;
- processOneVoxel(insideMask, cutoff, voxelPointIndexes, x, y, z,
contouring);
+ processOneVoxel(insideMask, cutoff, voxelPointIndexes, x, y, z,
+ contouring);
}
}
}
@@ -2758,4 +2787,141 @@
vector.y /= planarVectorLengths[1];
}
+ /////// orbitals ///////
+
+ //just for small factorials
+ static float factorial(int n) {
+ if (n == 0)
+ return 1;
+ return n * factorial(n - 1);
+ }
+
+ static float[] fact = new float[20];
+ static double[] rfactor = new double[10];
+ static double[] pfactor = new double[10];
+ int lastFactorial = -1;
+
+ void calcFactors(int n, int el, int m) {
+ int abm = Math.abs(m);
+ if (lastFactorial < n + el) {
+ for (int i = lastFactorial + 1; i <= n + el; i++)
+ fact[i] = factorial(i);
+ lastFactorial = n + el;
+ }
+ double Nnl = Math.pow(2 * psi_Znuc / n / A0, 1.5)
+ * Math.sqrt(fact[n - el - 1] / 2 / n / Math.pow(fact[n + el], 3));
+ double Lnl = fact[n + el] * fact[n + el];
+ double Plm = Math.pow(2, -el) * fact[el] * fact[el + abm]
+ * Math.sqrt((2 * el + 1) * fact[el - abm] / 2 / fact[el + abm]);
+
+ for (int p = 0; p <= n - el - 1; p++)
+ rfactor[p] = Nnl * Lnl / fact[p] / fact[n - el - p - 1]
+ / fact[2 * el + p + 1];
+ for (int p = abm; p <= el; p++)
+ pfactor[p] = Math.pow(-1, el - p) * Plm / fact[p] / fact[el + abm - p]
+ / fact[el - p] / fact[p - abm];
+ }
+
+ Point3f psiOrigin = new Point3f(0, 0, 0);
+ final static double A0 = 0.52918f; //x10^-10 meters
+ final static double ROOT2 = 1.414214;
+
+ float getPsi2(int x, int y, int z) {
+ Point3f ptPsi = voxelPtToXYZ(x, y, z);
+ float value = (float) hydrogenAtomPsiAt(ptPsi, psi_n, psi_l, psi_m);
+ return value * value;
+ }
+
+ double hydrogenAtomPsiAt(Point3f pt, int n, int el, int m) {
+ // ref: http://www.stolaf.edu/people/hansonr/imt/concept/schroed.pdf
+ int abm = Math.abs(m);
+ float x = pt.x - psiOrigin.x;
+ float y = pt.y - psiOrigin.y;
+ float z = pt.z - psiOrigin.z;
+ double x2y2 = x * x + y * y;
+ double r2 = x2y2 + z * z;
+ double r = Math.sqrt(r2);
+ double rho = 2d * psi_Znuc * r / n / A0;
+ double ph, th, cth, sth;
+ double theta_lm = 0;
+ double phi_m = 0;
+ double sum = 0;
+ for (int p = 0; p <= n - el - 1; p++)
+ sum += Math.pow(-rho, p) * rfactor[p];
+ double rnl = Math.exp(-rho / 2) * Math.pow(rho, el) * sum;
+ ph = Math.atan2(y, x);
+ th = Math.atan2(Math.sqrt(x2y2), z);
+ cth = Math.cos(th);
+ sth = Math.sin(th);
+ sum = 0;
+ for (int p = abm; p <= el; p++)
+ sum += Math.pow(1 + cth, p - abm) * Math.pow(1 - cth, el - p)
+ * pfactor[p];
+ theta_lm = Math.abs(Math.pow(sth, abm)) * sum;
+ if (m == 0)
+ phi_m = 1;
+ else if (m > 0)
+ phi_m = Math.cos(m * ph) * ROOT2;
+ else
+ phi_m = Math.sin(-m * ph) * ROOT2;
+ if (Math.abs(phi_m) < 0.0000000001)
+ phi_m = 0;
+ return rnl * theta_lm * phi_m;
+ }
+
+ int psi_n = 2;
+ int psi_l = 1;
+ int psi_m = 1;
+ float psi_ptsPerAngstrom = 5f;
+ float psi_radiusAngstroms;
+ float psi_Znuc = 1; // hydrogen
+
+ float autoScaleOrbital() {
+ float w = (psi_n * (psi_n + 3) - 5f) / psi_Znuc;
+ if (w < 1)
+ w = 1;
+ if (psi_n < 3)
+ w += 1;
+ return w;
+ }
+
+ int psi_gridMax = 61;
+
+ void setupOrbital() {
+ psi_radiusAngstroms = autoScaleOrbital();
+ int nGrid = (int) (psi_radiusAngstroms * psi_ptsPerAngstrom) * 2 + 1;
+ if (nGrid > psi_gridMax) {
+ nGrid = psi_gridMax;
+ psi_ptsPerAngstrom = ((nGrid - 1) / 2) / psi_radiusAngstroms;
+ }
+ jvxlFileHeader = "orbital psi squared\nn=" + psi_n + ", l=" + psi_l
+ + ", m=" + psi_m + " Znuc=" + psi_Znuc + " res=" + psi_ptsPerAngstrom
+ + " rad=" + psi_radiusAngstroms + "\n";
+ voxelCounts[0] = voxelCounts[1] = voxelCounts[2] = nGrid;
+ float dx = 1f / psi_ptsPerAngstrom;
+ volumetricVectors[0].set(dx, 0, 0);
+ volumetricVectors[1].set(0, dx, 0);
+ volumetricVectors[2].set(0, 0, dx);
+
+ volumetricVectorLengths[0] = volumetricVectorLengths[1] =
volumetricVectorLengths[2] = dx;
+ unitVolumetricVectors[0].set(1, 0, 0);
+ unitVolumetricVectors[1].set(0, 1, 0);
+ unitVolumetricVectors[2].set(0, 0, 1);
+
+ float origin = dx * (1 - nGrid) / 2;
+
+ volumetricOrigin.set(origin, origin, origin);
+
+ jvxlFileHeader += "-1 " + (origin / ANGSTROMS_PER_BOHR) + " "
+ + (origin / ANGSTROMS_PER_BOHR) + " " + (origin / ANGSTROMS_PER_BOHR)
+ + "\n";
+ for (int i = 0; i < 3; i++)
+ jvxlFileHeader += nGrid + " "
+ + (volumetricVectors[i].x / ANGSTROMS_PER_BOHR) + " "
+ + (volumetricVectors[i].y / ANGSTROMS_PER_BOHR) + " "
+ + (volumetricVectors[i].z / ANGSTROMS_PER_BOHR) + "\n";
+ atomCount = 0;
+ negativeAtomCount = false;
+ calcFactors(psi_n, psi_l, psi_m);
+ }
}
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
_______________________________________________
Jmol-commits mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/jmol-commits
