Revision: 5208 Author: hansonr Date: 2006-06-08 21:09:22 -0700 (Thu, 08 Jun 2006) ViewCVS: http://svn.sourceforge.net/jmol/?rev=5208&view=rev
Log Message: ----------- bob200603 10.x.11 isosurface, phases IV (orbitals) and V (sasurfaces) adds orbitals and solvent-accessible surfaces as isosurfaces adds dots sasurface for testing and comparison, not for publication, I think maybe fixes normix problem with static object issue see http://www.stolaf.edu/people/hansonr/jmol/test/proto/isosurface.htm Modified Paths: -------------- branches/bob200603/Jmol/src/org/jmol/g3d/Normix3D.java branches/bob200603/Jmol/src/org/jmol/viewer/Dots.java branches/bob200603/Jmol/src/org/jmol/viewer/DotsRenderer.java branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java branches/bob200603/Jmol/src/org/jmol/viewer/JmolConstants.java Modified: branches/bob200603/Jmol/src/org/jmol/g3d/Normix3D.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/g3d/Normix3D.java 2006-06-07 19:38:55 UTC (rev 5207) +++ branches/bob200603/Jmol/src/org/jmol/g3d/Normix3D.java 2006-06-09 04:09:22 UTC (rev 5208) @@ -50,8 +50,8 @@ final int normixCount; final static int[] vertexCounts = {12, 42, 162, 642, 2562}; - final static short[][] faceNormixesArrays = - new short[NORMIX_GEODESIC_LEVEL + 1][]; + static short[][] faceNormixesArrays; + // not "final" = new short[NORMIX_GEODESIC_LEVEL + 1][]; private final static boolean TIMINGS = false; private final static boolean DEBUG_WITH_SEQUENTIAL_SEARCH = false; @@ -288,6 +288,8 @@ } short[] getFaceNormixes(int level) { + if (faceNormixesArrays == null) + faceNormixesArrays = new short[NORMIX_GEODESIC_LEVEL + 1][]; short[] faceNormixes = faceNormixesArrays[level]; if (faceNormixes != null) return faceNormixes; Modified: branches/bob200603/Jmol/src/org/jmol/viewer/Dots.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/viewer/Dots.java 2006-06-07 19:38:55 UTC (rev 5207) +++ branches/bob200603/Jmol/src/org/jmol/viewer/Dots.java 2006-06-09 04:09:22 UTC (rev 5208) @@ -143,13 +143,20 @@ } + boolean showSurface; void setSize(int size, BitSet bsSelected) { // miguel 9 Sep 2005 // this is a hack ... // if mad == 0 then turn it off // if mad > 0 then use vdw radius // if mad < 0 then use ionic radius + // 9999 turns crude solvent surface on short mad = (short) size; + if (mad == 9999) { + showSurface = true; + return; + } + showSurface = false; this.mad = mad; boolean newUseVanderwaalsRadius = (mad > 0); timeBeginExecution = System.currentTimeMillis(); Modified: branches/bob200603/Jmol/src/org/jmol/viewer/DotsRenderer.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/viewer/DotsRenderer.java 2006-06-07 19:38:55 UTC (rev 5207) +++ branches/bob200603/Jmol/src/org/jmol/viewer/DotsRenderer.java 2006-06-09 04:09:22 UTC (rev 5208) @@ -87,7 +87,7 @@ short[] colixesConvex = dots.colixesConvex; int myVisibilityFlag = dots.myVisibilityFlag; boolean isInMotion = viewer.getInMotion(); - boolean iShowSolid = viewer.getTestFlag3() && dots.useBobsAlgorithm; + boolean iShowSolid = (viewer.getTestFlag3()||dots.showSurface) && dots.useBobsAlgorithm; for (int i = dots.dotsConvexMax; --i >= 0;) { Atom atom = atoms[i]; if (atom.isShapeVisible(myVisibilityFlag)) { Modified: branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java 2006-06-07 19:38:55 UTC (rev 5207) +++ branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java 2006-06-09 04:09:22 UTC (rev 5208) @@ -2661,6 +2661,9 @@ // I don't know what to do with this thing yet mad = (short) dotsParam; break; + case Token.sasurface: + mad = 9999; + break; default: booleanOrNumberExpected(); } @@ -4012,6 +4015,11 @@ propertyValue = new Integer(getArgbParam(i)); colorRangeStage++; break; + case Token.solvent: + propertyName = "solvent"; + propertyValue = new Float(viewer.getSolventProbeRadius()); + System.out.println (propertyName+" "+propertyValue); + break; case Token.identifier: String str = (String) token.value; if (str.equalsIgnoreCase("sign")) { @@ -4083,7 +4091,7 @@ } break; case Token.decimal: - if (colorRangeStage == 0 && orbitalPt == 0) { + if (colorRangeStage == 0 && orbitalPt <= 0) { propertyName = "cutoff"; propertyValue = token.value; break; @@ -4100,7 +4108,7 @@ break; propertyName = "orbital"; propertyValue = nlmZR; - orbitalPt = Integer.MIN_VALUE; + orbitalPt = -999; break; } if (i == contourPt) { @@ -4108,6 +4116,7 @@ integerExpected(); propertyName = "contour"; propertyValue = new Integer(token.intValue); + contourPt = 0; break; } if (i == fileIndexPt) { Modified: branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java 2006-06-07 19:38:55 UTC (rev 5207) +++ branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java 2006-06-09 04:09:22 UTC (rev 5208) @@ -107,6 +107,7 @@ final static int nContourMax = 100; final static int defaultColorNegative = Graphics3D.getArgbFromString("red"); final static int defaultColorPositive = Graphics3D.getArgbFromString("blue"); + final static float defaultSolventRadius = 1.2f; boolean colorBySign; int colorNeg; @@ -114,7 +115,10 @@ int colorPtr; boolean colorByPhase; int colorPhase; + boolean isCalculation; boolean isOrbital; + boolean isSolvent; + float solventRadius; int edgeFractionBase; int edgeFractionRange; @@ -150,7 +154,7 @@ boolean associateNormals; boolean force2SidedTriangles; boolean iDoContourPlane; - boolean datasetIs2dContour; + boolean dataIsJvxl2dContour; Point4f thePlane; int nContours; @@ -177,12 +181,14 @@ mappedDataMin = Float.MAX_VALUE; thisContour = -1; iDoContourPlane = true; - datasetIs2dContour = false; + dataIsJvxl2dContour = false; colorBySign = false; colorNeg = defaultColorNegative; colorPos = defaultColorPositive; colorByPhase = false; isOrbital = false; + isSolvent = false; + isCalculation = false; cutoff = Float.MAX_VALUE; super.setProperty("meshID", null, null); return; @@ -273,8 +279,19 @@ return; } + if ("solvent" == propertyName) { + isSolvent = isCalculation = true; + solventRadius = ((Float) value).floatValue(); + if (solventRadius < 0) + solventRadius = defaultSolventRadius; + System.out.println("creating solvent-accessible surface with radius " + + solventRadius); + cutoff = 0.0f; + propertyName = "bufferedReader"; + } + if ("orbital" == propertyName) { - isOrbital = true; + isOrbital = isCalculation = true; float[] nlmZR = (float[]) value; psi_n = (int) nlmZR[0]; psi_l = (int) nlmZR[1]; @@ -291,16 +308,16 @@ propertyName = "bufferedReader"; } if ("bufferedReader" == propertyName) { - if (!isOrbital) + if (!isCalculation) br = (BufferedReader) value; System.out.println("loading voxel data..."); checkFlags(); createIsosurface(); initializeMesh(force2SidedTriangles); jvxlFileMessage = (isJvxlColorMapped ? "mapped" : ""); - if (datasetIs2dContour) + if (dataIsJvxl2dContour) colorIsosurface(); - discardTempData(datasetIs2dContour); + discardTempData(dataIsJvxl2dContour); currentMesh.nBytes = nBytes; if (thePlane == null) System.out.println("\n" + jvxlGetFile(currentMesh, jvxlFileMessage)); @@ -308,6 +325,8 @@ return; } if ("colorReader" == propertyName) { + isSolvent = false; + isOrbital = false; System.out.println("mapping data..."); br = (BufferedReader) value; checkFlags(); @@ -344,12 +363,11 @@ } void colorIsosurface() { + setMapRanges(); if (thePlane == null) { - setMapRanges(); applyColorScale(); } else { - setMapRanges(); - generateContourData(datasetIs2dContour); + generateContourData(dataIsJvxl2dContour); initializeMesh(true); if (!colorByContourOnly) applyColorScale(); @@ -360,6 +378,9 @@ } void setMapRanges() { + // ["colorReader" | "bufferedReader/jvxl-2dContour] --> colorIsosurface + // ["phase" | colorIsoSurface] --> applyColorScale + // colorIsosurface --> generateContourData --> createContours if (colorByPhase) { mappedDataMin = -1; mappedDataMax = 1; @@ -495,6 +516,8 @@ try { if (isOrbital) { setupOrbital(); + } else if (isSolvent) { + setupSolvent(); } else { readTitleLines(); readAtomCountAndOrigin(); @@ -624,6 +647,7 @@ int nBytes; int nDataPoints; String surfaceData; + int nPointsX, nPointsY, nPointsZ; void readVoxelData(boolean isMapData) throws Exception { /* @@ -641,9 +665,9 @@ boolean inside = false; int dataCount = 0; - int nPointsX = voxelCounts[0]; - int nPointsY = voxelCounts[1]; - int nPointsZ = voxelCounts[2]; + nPointsX = voxelCounts[0]; + nPointsY = voxelCounts[1]; + nPointsZ = voxelCounts[2]; int nPoints = nPointsX * nPointsY * nPointsZ; ichNextParse = 0; @@ -688,6 +712,8 @@ if (isOrbital) voxelValue = getPsi2(x, y, z); + else if (isSolvent) + voxelValue = getSolventValue(x, y, z); else if (isMapData || thePlane == null) voxelValue = getNextVoxelValue(); else @@ -772,7 +798,7 @@ isJvxlColorPrecision = true; } nEdgeData = -1; - datasetIs2dContour = (isJvxlColorMapped && iDoContourPlane); + dataIsJvxl2dContour = (isJvxlColorMapped && iDoContourPlane); try { thePlane = new Point4f(parseFloat(), parseFloat(), parseFloat(), parseFloat()); @@ -1074,7 +1100,7 @@ if (colorByPhase) datum = value = getPhase(mesh.vertices[vertexIndex]); - else if (datasetIs2dContour) + else if (dataIsJvxl2dContour) datum = value = getInterpolatedPixelValue(mesh.vertices[vertexIndex]); else datum = value = lookupInterpolatedVoxelValue(mesh.vertices[vertexIndex]); @@ -1146,13 +1172,14 @@ for (int i = 0; i < vertexCount; i += incr) if (mesh.firstViewableVertex == 0 || i < mesh.firstViewableVertex) { float challenger; - if (datasetIs2dContour) + if (dataIsJvxl2dContour) challenger = getInterpolatedPixelValue(vertexes[i]); else challenger = lookupInterpolatedVoxelValue(vertexes[i]); if (challenger < min) min = challenger; } + System.out.println("minimum mapped value: " + min); return min; } @@ -1165,7 +1192,6 @@ if (challenger > max) max = challenger; } - System.out.println("maximum mapped value: " + max); return max; } @@ -1177,13 +1203,14 @@ for (int i = 0; i < vertexCount; i += incr) if (mesh.firstViewableVertex == 0 || i < mesh.firstViewableVertex) { float challenger; - if (datasetIs2dContour) + if (dataIsJvxl2dContour) challenger = getInterpolatedPixelValue(vertexes[i]); else challenger = lookupInterpolatedVoxelValue(vertexes[i]); if (challenger > max) max = challenger; } + System.out.println("maximum mapped value: " + max); return max; } @@ -1212,7 +1239,7 @@ int yDown = indexDown(y, iMax = voxelCounts[1] - 1); int yUp = yDown + (y < 0 || yDown == iMax ? 0 : 1); int zDown = indexDown(z, iMax = voxelCounts[2] - 1); - int zUp = zDown + (z < 0 || zDown == iMax || datasetIs2dContour ? 0 : 1); + int zUp = zDown + (z < 0 || zDown == iMax || dataIsJvxl2dContour ? 0 : 1); float v1 = getFractional2DValue(x - xDown, y - yDown, voxelData[xDown][yDown][zDown], voxelData[xUp][yDown][zDown], voxelData[xDown][yUp][zDown], voxelData[xUp][yUp][zDown]); @@ -1825,7 +1852,7 @@ new Vector3f(0, 1, 0), new Vector3f(1, 1, 0), new Vector3f(1, 1, 1), new Vector3f(0, 1, 1) }; - /*static*/ Vector3f[] voxelVertexVectors = new Vector3f[8]; + Vector3f[] voxelVertexVectors = new Vector3f[8]; void calcVoxelVertexVectors() { for (int i = 8; --i >= 0;) @@ -2109,7 +2136,7 @@ final Vector3f[] planarVectors = new Vector3f[3]; final Vector3f[] unitPlanarVectors = new Vector3f[3]; final float[] planarVectorLengths = new float[2]; - final /*static*/ Matrix3f matXyzToPlane = new Matrix3f(); + final Matrix3f matXyzToPlane = new Matrix3f(); { planarVectors[0] = new Vector3f(); planarVectors[1] = new Vector3f(); @@ -2574,7 +2601,7 @@ return fraction; } - /*static*/ Vector3f[] pixelVertexVectors = new Vector3f[4]; + Vector3f[] pixelVertexVectors = new Vector3f[4]; void calcPixelVertexVectors() { for (int i = 4; --i >= 0;) @@ -2769,13 +2796,13 @@ final Point3i ptiTemp = new Point3i(); - Point3i locatePixel(Point3f ptXYZ) { - pointVector.set(ptXYZ); + Point3i locatePixel(Point3f ptXyz) { + pointVector.set(ptXyz); toPlanarCoord(pointVector); ptiTemp.x = (int) (pointVector.x + 0.5f); //NOTE: fails if negative -- (int) (-0.9 + 0.5) = (int) (-0.4) = 0 ptiTemp.y = (int) (pointVector.y + 0.5f); - //System.out.println("locate " + ptXYZ + " " + pointVector + " " + ptiTemp); + //System.out.println("locate " + ptXyz + " " + pointVector + " " + ptiTemp); return ptiTemp; } @@ -2787,8 +2814,72 @@ vector.y /= planarVectorLengths[1]; } - /////// orbitals /////// + /////// for any sort of functional mapping //////// + class Voxel extends Point3i { + Point3f ptXyz; + float value; + + void setValue(int x, int y, int z, float value) { + this.x = x; + this.y = y; + this.z = z; + this.value = value; + this.ptXyz = voxelPtToXYZ(x, y, z); + } + + void setValue(float value) { + if (this.value < value) + return; + if (logMessages) + System.out.println("voxel.setValue " + x + " " + y + " " + z + ptXyz + + ": " + value + " was " + this.value); + this.value = value; + } + } + + void setVoxelRange(int index, float min, float max, float ptsPerAngstrom, + int gridMax) { + float range = max - min; + int nGrid = (int) (range * ptsPerAngstrom); + if (nGrid > solvent_gridMax) { + nGrid = solvent_gridMax; + ptsPerAngstrom = nGrid / range; + } + + float d = volumetricVectorLengths[index] = 1f / ptsPerAngstrom; + voxelCounts[index] = nGrid + 3; + + switch (index) { + case 0: + unitVolumetricVectors[0].set(1, 0, 0); + volumetricVectors[0].set(d, 0, 0); + volumetricOrigin.x = min - d; + break; + case 1: + unitVolumetricVectors[1].set(0, 1, 0); + volumetricVectors[1].set(0, d, 0); + volumetricOrigin.y = min - d; + break; + case 2: + unitVolumetricVectors[2].set(0, 0, 1); + volumetricVectors[2].set(0, 0, d); + volumetricOrigin.z = min - d; + } + } + + String jvxlGetVolumeHeader() { + String str = "-1 " + (volumetricOrigin.x / ANGSTROMS_PER_BOHR) + " " + + (volumetricOrigin.y / ANGSTROMS_PER_BOHR) + " " + + (volumetricOrigin.z / ANGSTROMS_PER_BOHR) + "\n"; + for (int i = 0; i < 3; i++) + str += voxelCounts[i] + " " + + (volumetricVectors[i].x / ANGSTROMS_PER_BOHR) + " " + + (volumetricVectors[i].y / ANGSTROMS_PER_BOHR) + " " + + (volumetricVectors[i].z / ANGSTROMS_PER_BOHR) + "\n"; + return str; + } + //just for small factorials static float factorial(int n) { if (n == 0) @@ -2796,11 +2887,51 @@ return n * factorial(n - 1); } - /*static*/ float[] fact = new float[20]; - /*static*/ double[] rfactor = new double[10]; - /*static*/ double[] pfactor = new double[10]; + float[] fact = new float[20]; + + /////// orbitals /////// + + double[] rfactor = new double[10]; + double[] pfactor = new double[10]; int lastFactorial = -1; + int psi_gridMax = 60; + Point3f psiOrigin = new Point3f(0, 0, 0); + final static double A0 = 0.52918f; //x10^-10 meters + final static double ROOT2 = 1.414214; + int psi_n = 2; + int psi_l = 1; + int psi_m = 1; + float psi_ptsPerAngstrom = 5f; + float psi_radiusAngstroms; + float psi_Znuc = 1; // hydrogen + + void setupOrbital() { + psi_radiusAngstroms = autoScaleOrbital(); + + for (int i = 0; i < 3; i++) + setVoxelRange(i, -psi_radiusAngstroms, psi_radiusAngstroms, + psi_ptsPerAngstrom, psi_gridMax); + + jvxlFileHeader = "orbital psi squared\nn=" + psi_n + ", l=" + psi_l + + ", m=" + psi_m + " Znuc=" + psi_Znuc + " res=" + psi_ptsPerAngstrom + + " rad=" + psi_radiusAngstroms + "\n"; + + jvxlFileHeader += jvxlGetVolumeHeader(); + atomCount = 0; + negativeAtomCount = false; + calcFactors(psi_n, psi_l, psi_m); + } + + 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; + } + void calcFactors(int n, int el, int m) { int abm = Math.abs(m); if (lastFactorial < n + el) { @@ -2822,10 +2953,6 @@ / 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); @@ -2869,59 +2996,207 @@ 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 + ///// solvent-accessible surface ////// - 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; - } + float solvent_ptsPerAngstrom = 8f; + int solvent_gridMax = 40; + BitSet solvent_bsAtoms; + int solvent_modelIndex; + float[] solvent_atomRadius; + Point3f[] solvent_ptAtom; + int solvent_nAtoms; - int psi_gridMax = 61; + Voxel solvent_voxel = new Voxel(); - 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; + void setupSolvent() { + solvent_bsAtoms = viewer.getSelectionSet(); + Atom[] atoms = frame.atoms; + Point3f xyzMin = new Point3f(Integer.MAX_VALUE, Integer.MAX_VALUE, + Integer.MAX_VALUE); + Point3f xyzMax = new Point3f(-Integer.MAX_VALUE, -Integer.MAX_VALUE, + -Integer.MAX_VALUE); + solvent_modelIndex = -1; + int iAtom = 0; + for (int i = solvent_bsAtoms.size(); --i >= 0;) + if (solvent_bsAtoms.get(i)) { + if (solvent_modelIndex < 0) + solvent_modelIndex = frame.atoms[i].modelIndex; + if (solvent_modelIndex != frame.atoms[i].modelIndex) + continue; + ++iAtom; + } + solvent_nAtoms = iAtom; + System.out.println(iAtom + + " atoms will be used in the solvent-accessible surface calculation"); + if (iAtom > 0) { + solvent_atomRadius = new float[iAtom]; + solvent_ptAtom = new Point3f[iAtom]; } - 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); + iAtom = 0; + for (int i = solvent_bsAtoms.size(); --i >= 0;) { + if (!solvent_bsAtoms.get(i) + || solvent_modelIndex != frame.atoms[i].modelIndex) + continue; + Point3f pt = solvent_ptAtom[iAtom] = atoms[i].point3f; + float rA = solvent_atomRadius[iAtom++] = frame.atoms[i] + .getVanderwaalsRadiusFloat(); + if (pt.x - rA < xyzMin.x) + xyzMin.x = pt.x - rA; + if (pt.x + rA > xyzMax.x) + xyzMax.x = pt.x + rA; + if (pt.y - rA < xyzMin.y) + xyzMin.y = pt.y - rA; + if (pt.y + rA > xyzMax.y) + xyzMax.y = pt.y + rA; + if (pt.z - rA < xyzMin.z) + xyzMin.z = pt.z - rA; + if (pt.z + rA > xyzMax.z) + xyzMax.z = pt.z + rA; + } + System.out.println("surface range " + xyzMin + " to " + xyzMax); + jvxlFileHeader = "solvent accesible surface\nrange " + xyzMin + " to " + + xyzMax + "\n"; - 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); + int maxGrid = solvent_gridMax / (iAtom > 50 ? 2 : 1); + setVoxelRange(0, xyzMin.x, xyzMax.x, solvent_ptsPerAngstrom, + maxGrid); + setVoxelRange(1, xyzMin.y, xyzMax.y, solvent_ptsPerAngstrom, + maxGrid); + setVoxelRange(2, xyzMin.z, xyzMax.z, solvent_ptsPerAngstrom, + maxGrid); - 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"; + jvxlFileHeader += jvxlGetVolumeHeader(); atomCount = 0; negativeAtomCount = false; - calcFactors(psi_n, psi_l, psi_m); } + + float getSolventValue(int x, int y, int z) { + solvent_voxel.setValue(x, y, z, Float.MAX_VALUE); + float rA, rB; + Point3f ptA, ptB; + for (int i = 0; i < solvent_nAtoms && solvent_voxel.value >= -0.5; i++) { + ptA = solvent_ptAtom[i]; + rA = solvent_atomRadius[i]; + float v = solvent_voxel.ptXyz.distance(ptA) - rA; + if (v < solvent_voxel.value) + solvent_voxel.setValue(v); + } + if (solventRadius == 0) + return solvent_voxel.value; + for (int i = 0; i < solvent_nAtoms - 1 && solvent_voxel.value >= -0.5; i++) { + ptA = solvent_ptAtom[i]; + rA = solvent_atomRadius[i] + solventRadius; + for (int j = i + 1; j < solvent_nAtoms && solvent_voxel.value >= -0.5; j++) { + ptB = solvent_ptAtom[j]; + rB = solvent_atomRadius[j] + solventRadius; + float dAB = ptA.distance(ptB); + if (dAB >= rA + rB) + continue; + checkSpecialVoxel(ptA, rA, ptB, rB, dAB); + } + } + return solvent_voxel.value; + } + + final Point3f ptS = new Point3f(); + + void checkSpecialVoxel(Point3f ptA, float rAS, Point3f ptB, float rBS, + float dAB) { + /* + * Checking here for voxels that are in the situation: + * + * A------)-- V ---((--))-- S --(------B + * |----d--------| + * or + * + * B------)-- V ---((--))-- S --(------A + * |----d--------| + * + * A and B are the two atom centers; V is the voxel; S is a hypothetical + * PROJECTED solvent center based on the position of V in relation to + * first A, then B; ( and ) are atom radii and (( )) are the overlapping + * atom+solvent radii. + * + * That is, where the projected solvent location for one voxel is + * within the solvent radius sphere of another, this voxel should + * be checked in relation to solvent distance, not atom distance. + * + * + * S + *++ / \ ++ + * ++ / | \ ++ + * + V + x want V such that angle ASV < angle ASB + * / ****** \ + * A ----------B + * b + * + * A, B are atoms; S is solvent center; V is voxel point + * objective is to calculate dSV. ++ Here represents the van der Waals + * radius for each atom. ***** is the key "trough" location. + * + * Getting dVS: + * + * Known: rAB, rAS, rBS, giving angle BAS (theta) + * Known: rAB, rAV, rBV, giving angle VAB (alpha) + * Determined: angle VAS (theta - alpha), and from that, dSV, using + * the cosine law: + * + * a^2 + b^2 - 2ab Cos(theta) = c^2. + * + * The trough issue: + * + * Since the voxel might be at point x (above), outside the + * triangle, we have to test for that. What we will be looking + * for in the "trough" will be that angle ASV < angle ASB + * that is, cosASB < cosASV + * + * If we find the voxel in the "trough", then we set its value to + * (solvent radius - dVS). + * + */ + Point3f ptV = solvent_voxel.ptXyz; + float dAV = ptA.distance(ptV); + float dBV = ptB.distance(ptV); + float dVS; + float f = rAS / dAV; + if (f > 1) { + ptS.set(ptA.x + (ptV.x - ptA.x) * f, ptA.y + (ptV.y - ptA.y) * f, ptA.z + + (ptV.z - ptA.z) * f); + if (ptB.distance(ptS) < rBS) { + dVS = solventDistance(ptV, ptA, ptB, rAS, rBS, dAB, dAV, dBV); + if (voxelIsInTrough(dVS, rAS * rAS, rBS, dAB, dAV, dBV)) + solvent_voxel.setValue(solventRadius - dVS); + } + return; + } + f = rBS / dBV; + if (f <= 1) + return; + ptS.set(ptB.x + (ptV.x - ptB.x) * f, ptB.y + (ptV.y - ptB.y) * f, ptB.z + + (ptV.z - ptB.z) * f); + if (ptA.distance(ptS) < rAS) { + dVS = solventDistance(ptV, ptB, ptA, rBS, rAS, dAB, dBV, dAV); + if (voxelIsInTrough(dVS, rAS * rAS, rBS, dAB, dAV, dBV)) + solvent_voxel.setValue(solventRadius - dVS); + } + } + + boolean voxelIsInTrough(float dVS, float rAS2, float rBS, float dAB, + float dAV, float dBV) { + //only calculate what we need -- a factor proportional to cos + float cosASBf = (rAS2 + rBS * rBS - dAB * dAB) / rBS; // /2 /rAS); + float cosASVf = (rAS2 + dVS * dVS - dAV * dAV) / dVS; // /2 /rAS); + return (cosASBf < cosASVf); + } + + float solventDistance(Point3f ptV, Point3f ptA, Point3f ptB, float rAS, + float rBS, float dAB, float dAV, float dBV) { + double angleVAB = Math.acos((dAV * dAV + dAB * dAB - dBV * dBV) + / (2 * dAV * dAB)); + double angleBAS = Math.acos((dAB * dAB + rAS * rAS - rBS * rBS) + / (2 * dAB * rAS)); + float dVS = (float) Math.sqrt(rAS * rAS + dAV * dAV - 2 * rAS * dAV + * Math.cos(angleBAS - angleVAB)); + return dVS; + } } Modified: branches/bob200603/Jmol/src/org/jmol/viewer/JmolConstants.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/viewer/JmolConstants.java 2006-06-07 19:38:55 UTC (rev 5207) +++ branches/bob200603/Jmol/src/org/jmol/viewer/JmolConstants.java 2006-06-09 04:09:22 UTC (rev 5208) @@ -39,7 +39,7 @@ // for now, just update this by hand // perhaps use ant filter later ... but mth doesn't like it :-( public final static String copyright = "(C) 2006 Jmol Development"; - public final static String version = "10.x.10(branch bob200603)"; + public final static String version = "10.x.11(branch bob200603)"; public final static String cvsDate = "$Date$"; public final static String date = cvsDate.substring(7, 23); 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
