Package org.jmol.jvxl.readers

Source Code of org.jmol.jvxl.readers.IsoSolventReader

/* $RCSfile$
* $Author: hansonr $
* $Date: 2007-03-30 11:40:16 -0500 (Fri, 30 Mar 2007) $
* $Revision: 7273 $
*
* Copyright (C) 2007 Miguel, Bob, Jmol Development
*
* Contact: hansonr@stolaf.edu
*
*  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.jvxl.readers;

import java.util.BitSet;

import javax.vecmath.Point3f;
import javax.vecmath.Point3i;
import javax.vecmath.Point4f;

import org.jmol.util.Logger;

import org.jmol.api.AtomIndexIterator;
import org.jmol.jvxl.data.MeshData;

class IsoSolventReader extends AtomDataReader {

  IsoSolventReader(SurfaceGenerator sg) {
    super(sg);
  }

  ///// solvent-accessible, solvent-excluded surface //////

  /*
   * The surface fragment idea:
   *
   * isosurface solvent|sasurface both work on the SELECTED atoms, thus
   * allowing for a subset of the molecule to be involved. But in that
   * case we don't want to be creating a surface that goes right through
   * another atom. Rather, what we want (probably) is just the portion
   * of the OVERALL surface that involves these atoms.
   *
   * The addition of Mesh.voxelValue[] means that we can specify any
   * voxel we want to NOT be excluded (NaN). Here we first exclude any
   * voxel that would have been INSIDE a nearby atom. This will take care
   * of any portion of the vanderwaals surface that would be there. Then
   * we exclude any special-case voxel that is between two nearby atoms.
   * 
   *  Bob Hanson 13 Jul 2006
   * 
   */

  private float cavityRadius;
  private float envelopeRadius;
  private Point3f[] dots;

  private boolean doCalculateTroughs;
  private boolean isCavity, isPocket;
  private float solventRadius;
 
  protected void setup() {
    super.setup();
    cavityRadius = params.cavityRadius;
    envelopeRadius = params.envelopeRadius;
    solventRadius = params.solventRadius;
    point = params.point;

    isCavity = (params.isCavity && meshDataServer != null); // Jvxl cannot do this calculation on its own.
    isPocket = (params.pocket != null && meshDataServer != null);

    doCalculateTroughs = (atomDataServer != null && !isCavity // Jvxl needs an atom iterator to do this.
        && solventRadius > 0 && (dataType == Parameters.SURFACE_SOLVENT || dataType == Parameters.SURFACE_MOLECULAR));
    doUseIterator = doCalculateTroughs;
    getAtoms(Float.NaN, false, true);
    if (isCavity || isPocket)
      dots = meshDataServer.calculateGeodesicSurface(bsMySelected,
          envelopeRadius);

    setHeader("solvent/molecular surface", params.calculationType);
    setRangesAndAddAtoms(params.solvent_ptsPerAngstrom, params.solvent_gridMax,
        params.thePlane != null ? Integer.MAX_VALUE : Math.min(firstNearbyAtom,
            100));
  }

  //////////// meshData extensions ////////////

  public void selectPocket(boolean doExclude) {
    if (meshDataServer != null)
      meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null);
    //mark VERTICES for proximity to surface
    Point3f[] v = meshData.vertices;
    int nVertices = meshData.vertexCount;
    float[] vv = meshData.vertexValues;
    int nDots = dots.length;
    for (int i = 0; i < nVertices; i++) {
      for (int j = 0; j < nDots; j++) {
        if (dots[j].distance(v[i]) < envelopeRadius) {
          vv[i] = Float.NaN;
          continue;
        }
      }
    }
    meshData.getSurfaceSet();
    int nSets = meshData.nSets;
    BitSet pocketSet = new BitSet(nSets);
    BitSet ss;
    for (int i = 0; i < nSets; i++)
      if ((ss = meshData.surfaceSet[i]) != null)
        for (int j = ss.nextSetBit(0); j >= 0; j = ss.nextSetBit(j + 1))
          if (Float.isNaN(meshData.vertexValues[j])) {
            pocketSet.set(i);
            //System.out.println("pocket " + i + " " + j + " " + surfaceSet[i]);
            break;
          }
    //now clear all vertices that match the pocket toggle
    //"POCKET"   --> pocket TRUE means "show just the pockets"
    //"INTERIOR" --> pocket FALSE means "show everything that is not a pocket"
    for (int i = 0; i < nSets; i++)
      if (meshData.surfaceSet[i] != null && pocketSet.get(i) == doExclude)
        meshData.invalidateSurfaceSet(i);
    updateSurfaceData();
    if (!doExclude)
      meshData.surfaceSet = null;
    if (meshDataServer != null) {
      meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_SETS, null);
      meshData = new MeshData();
    }
  }
 
 
  /////////////// calculation methods //////////////
 
  protected void generateCube() {
    /*
     *
     * Jmol cavity rendering. Tim Driscoll suggested "filling a
     * protein with foam. Here you go...
     *
     * 1) Use a dot-surface extended x.xx Angstroms to define the
     *    outer envelope of the protein.
     * 2) Identify all voxel points outside the protein surface (v > 0)
     *    but inside the envelope (nearest distance to a dot > x.xx).
     * 3) First pass -- create the protein surface.
     * 4) Replace solvent atom set with "foam" ball of the right radius
     *    at the voxel vertex points.
     * 5) Run through a second time using these "atoms" to generate
     *    the surface around the foam spheres.
     *   
     *    Bob Hanson 3/19/07
     *
     */

    volumeData.voxelData = voxelData = new float[nPointsX][nPointsY][nPointsZ];
    if (isCavity && params.theProperty != null) {
      /*
       * couldn't get this to work -- we only have half of the points
       for (int x = 0, i = 0, ipt = 0; x < nPointsX; ++x)
       for (int y = 0; y < nPointsY; ++y)
       for (int z = 0; z < nPointsZ; ++z)
       voxelData[x][y][z] = (
       bs.get(i++) ?
       surface_data[ipt++]
       : Float.NaN);
       mappedDataMin = 0;
       mappedDataMax = 0;
       for (int i = 0; i < nAtoms; i++)
       if (surface_data[i] > mappedDataMax)
       mappedDataMax = surface_data[i];     
       */
      return;
    }
    generateSolventCube(true);
    if (isCavity && dataType != Parameters.SURFACE_NOMAP
        && dataType != Parameters.SURFACE_PROPERTY) {
      generateSolventCavity();
      generateSolventCube(false);
    }
    if (params.cappingObject instanceof Point4f) { // had a check for mapping here, turning this off
      //Logger.info("capping isosurface using " + params.cappingPlane);
      volumeData.capData((Point4f)params.cappingObject, params.cutoff);
      params.cappingObject = null;
    }
  }

private void generateSolventCavity() {
    //we have a ring of dots around the model.
    //1) identify which voxelData points are > 0 and within this volume
    //2) turn these voxel points into atoms with given radii
    //3) rerun the calculation to mark a solvent around these!
    BitSet bs = new BitSet(nPointsX * nPointsY * nPointsZ);
    int i = 0;
    int nDots = dots.length;
    int n = 0;
    //surface_data = new float[1000];

    // for (int j = 0; j < nDots; j++) {
    //   System.out.println("draw pt"+j+" {"+dots[j].x + " " +dots[j].y + " " +dots[j].z + "} " );
    // }
    /*
     Point3f pt = new Point3f(15.3375f, 1.0224999f, 5.1125f);
     for (int x = 0; x < nPointsX; ++x)
     for (int y = 0; y < nPointsY; ++y) {
     out: for (int z = 0; z < nPointsZ; ++z, ++i) {
     float d = voxelData[x][y][z];
     volumeData.voxelPtToXYZ(x, y, z, ptXyzTemp);
     if (d < Float.MAX_VALUE && d >= cavityRadius)
     if (ptXyzTemp.distance(pt) < 3f)
     System.out.println("draw pt"+(n++)+" {"+ptXyzTemp.x + " " +ptXyzTemp.y + " " +ptXyzTemp.z + "} # " +x + " " + y + " " + z + " " + voxelData[x][y][z]);
     if (false && d < Float.MAX_VALUE &&
     d >= cavityRadius) {
     volumeData.voxelPtToXYZ(x, y, z, ptXyzTemp);
     if (ptXyzTemp.x > 0 && ptXyzTemp.y > 0 && ptXyzTemp.z > 0)
     System.out.println("draw pt"+(n++)+" {"+ptXyzTemp.x + " " +ptXyzTemp.y + " " +ptXyzTemp.z + "} " );
     }
     }
     }
     n = 0;
     i = 0;
     */
    float d;
    float r2 = envelopeRadius;// - cavityRadius;

    for (int x = 0; x < nPointsX; ++x)
      for (int y = 0; y < nPointsY; ++y) {
        out: for (int z = 0; z < nPointsZ; ++z, ++i)
          if ((d = voxelData[x][y][z]) < Float.MAX_VALUE && d >= cavityRadius) {
            volumeData.voxelPtToXYZ(x, y, z, ptXyzTemp);
            for (int j = 0; j < nDots; j++) {
              if (dots[j].distance(ptXyzTemp) < r2)
                continue out;
            }
            bs.set(i);
            n++;
          }
      }
    Logger.info("cavities include " + n + " voxel points");
    atomRadius = new float[n];
    atomXyz = new Point3f[n];
    for (int x = 0, ipt = 0, apt = 0; x < nPointsX; ++x)
      for (int y = 0; y < nPointsY; ++y)
        for (int z = 0; z < nPointsZ; ++z)
          if (bs.get(ipt++)) {
            volumeData.voxelPtToXYZ(x, y, z, (atomXyz[apt] = new Point3f()));
            atomRadius[apt++] = voxelData[x][y][z];
          }
    myAtomCount = firstNearbyAtom = n;
  }

  final Point3f ptXyzTemp = new Point3f();

  void generateSolventCube(boolean isFirstPass) {
    float distance = params.distance;
    float rA, rB;
    Point3f ptA;
    Point3f ptY0 = new Point3f(), ptZ0 = new Point3f();
    Point3i pt0 = new Point3i(), pt1 = new Point3i();
    float value = Float.MAX_VALUE;
    if (Logger.debugging)
      Logger.startTimer();
    for (int x = 0; x < nPointsX; ++x)
      for (int y = 0; y < nPointsY; ++y)
        for (int z = 0; z < nPointsZ; ++z)
          voxelData[x][y][z] = value;
    if (dataType == Parameters.SURFACE_NOMAP)
      return;
    int atomCount = myAtomCount;
    float maxRadius = 0;
    float r0 = (isFirstPass && isCavity ? cavityRadius : 0);
    boolean isWithin = (isFirstPass && distance != Float.MAX_VALUE && point != null);
    AtomIndexIterator iter = (doCalculateTroughs ?
        atomDataServer.getSelectedAtomIterator(bsMySelected, true, false) : null);
    for (int iAtom = 0; iAtom < atomCount; iAtom++) {
      ptA = atomXyz[iAtom];
      rA = atomRadius[iAtom];
      if (rA > maxRadius)
        maxRadius = rA;
      if (isWithin && ptA.distance(point) > distance + rA + 0.5)
        continue;
      boolean isNearby = (iAtom >= firstNearbyAtom);
      float rA0 = rA + r0;
      setGridLimitsForAtom(ptA, rA0, pt0, pt1);
      volumeData.voxelPtToXYZ(pt0.x, pt0.y, pt0.z, ptXyzTemp);
      for (int i = pt0.x; i < pt1.x; i++) {
        ptY0.set(ptXyzTemp);
        for (int j = pt0.y; j < pt1.y; j++) {
          ptZ0.set(ptXyzTemp);
          for (int k = pt0.z; k < pt1.z; k++) {
            float v = ptXyzTemp.distance(ptA) - rA;
            if ((r0 == 0 || v <= rA0) && v < voxelData[i][j][k]) {
              voxelData[i][j][k] = (isNearby || isWithin
                  && ptXyzTemp.distance(point) > distance ? Float.NaN : v);
            }
            ptXyzTemp.add(volumetricVectors[2]);
          }
          ptXyzTemp.set(ptZ0);
          ptXyzTemp.add(volumetricVectors[1]);
        }
        ptXyzTemp.set(ptY0);
        ptXyzTemp.add(volumetricVectors[0]);
      }
    }
    if (isCavity && isFirstPass)
      return;
    if (doCalculateTroughs) {
      Point3i ptA0 = new Point3i();
      Point3i ptB0 = new Point3i();
      Point3i ptA1 = new Point3i();
      Point3i ptB1 = new Point3i();
      for (int iAtom = 0; iAtom < firstNearbyAtom - 1; iAtom++)
        if (atomNo[iAtom] > 0) {
          ptA = atomXyz[iAtom];
          rA = atomRadius[iAtom] + solventRadius;
          int iatomA = atomIndex[iAtom];
          if (isWithin && ptA.distance(point) > distance + rA + 0.5)
            continue;
          setGridLimitsForAtom(ptA, rA - solventRadius, ptA0, ptA1);
          atomDataServer.setIteratorForAtom(iter, iatomA, rA + solventRadius + maxRadius);
          //true ==> only atom index > this atom accepted
          while (iter.hasNext()) {
            int iatomB = iter.next();
            Point3f ptB = atomXyz[myIndex[iatomB]];
            rB = atomData.atomRadius[iatomB] + solventRadius;
            if (isWithin && ptB.distance(point) > distance + rB + 0.5)
              continue;
            if (params.thePlane != null
                && Math.abs(volumeData.distancePointToPlane(ptB)) > 2 * rB)
              continue;

            float dAB = ptA.distance(ptB);
            if (dAB >= rA + rB)
              continue;
            //defining pt0 and pt1 very crudely -- this could be refined
            setGridLimitsForAtom(ptB, rB - solventRadius, ptB0, ptB1);
            pt0.x = Math.min(ptA0.x, ptB0.x);
            pt0.y = Math.min(ptA0.y, ptB0.y);
            pt0.z = Math.min(ptA0.z, ptB0.z);
            pt1.x = Math.max(ptA1.x, ptB1.x);
            pt1.y = Math.max(ptA1.y, ptB1.y);
            pt1.z = Math.max(ptA1.z, ptB1.z);
            volumeData.voxelPtToXYZ(pt0.x, pt0.y, pt0.z, ptXyzTemp);
            for (int i = pt0.x; i < pt1.x; i++) {
              ptY0.set(ptXyzTemp);
              for (int j = pt0.y; j < pt1.y; j++) {
                ptZ0.set(ptXyzTemp);
                for (int k = pt0.z; k < pt1.z; k++) {
                  float dVS = checkSpecialVoxel(ptA, rA, ptB, rB, dAB,
                      ptXyzTemp);
                  if (!Float.isNaN(dVS)) {
                    float v = solventRadius - dVS;
                    if (v < voxelData[i][j][k]) {
                      voxelData[i][j][k] = (isWithin
                          && ptXyzTemp.distance(point) > distance ? Float.NaN
                          : v);
                    }
                  }
                  ptXyzTemp.add(volumetricVectors[2]);
                }
                ptXyzTemp.set(ptZ0);
                ptXyzTemp.add(volumetricVectors[1]);
              }
              ptXyzTemp.set(ptY0);
              ptXyzTemp.add(volumetricVectors[0]);
            }
          }
        }
      iter.release();
      iter = null;
    }
    if (params.thePlane == null) {
      for (int x = 0; x < nPointsX; ++x)
        for (int y = 0; y < nPointsY; ++y)
          for (int z = 0; z < nPointsZ; ++z)
            if (voxelData[x][y][z] == Float.MAX_VALUE)
              voxelData[x][y][z] = Float.NaN;
    } else { //solvent planes just focus on negative values
      value = 0.001f;
      for (int x = 0; x < nPointsX; ++x)
        for (int y = 0; y < nPointsY; ++y)
          for (int z = 0; z < nPointsZ; ++z)
            if (voxelData[x][y][z] < value) {
              // Float.NaN will also match ">=" this way
            } else {
              voxelData[x][y][z] = value;
            }
    }
    if (Logger.debugging)
      Logger.checkTimer("solvent surface time");
  }

  void setGridLimitsForAtom(Point3f ptA, float rA, Point3i pt0, Point3i pt1) {
    int n = 1;
    volumeData.xyzToVoxelPt(ptA.x - rA, ptA.y - rA, ptA.z - rA, pt0);
    pt0.x -= n;
    pt0.y -= n;
    pt0.z -= n;
    if (pt0.x < 0)
      pt0.x = 0;
    if (pt0.y < 0)
      pt0.y = 0;
    if (pt0.z < 0)
      pt0.z = 0;
    volumeData.xyzToVoxelPt(ptA.x + rA, ptA.y + rA, ptA.z + rA, pt1);
    pt1.x += n + 1;
    pt1.y += n + 1;
    pt1.z += n + 1;
    if (pt1.x >= nPointsX)
      pt1.x = nPointsX;
    if (pt1.y >= nPointsY)
      pt1.y = nPointsY;
    if (pt1.z >= nPointsZ)
      pt1.z = nPointsZ;
  }

  final Point3f ptS = new Point3f();

  float checkSpecialVoxel(Point3f ptA, float rAS, Point3f ptB, float rBS,
                          float dAB, Point3f ptV) {
    /*
     * Checking here for voxels that are in the situation:
     *
     * A------)(-----S-----)(------B  (not actually linear)
     * |-----rAS-----|-----rBS-----|
     * |-----------dAB-------------|
     *         ptV
     * |--dAV---|---------dBV------|
     *
     * A and B are the two atom centers; S is a hypothetical
     * PROJECTED solvent center based on the position of ptV
     * in relation to first A, then B.
     *
     * 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
     *
     *  ++ Here represents the van der Waals radius for each atom.
     *  ***** is the key "trough" location. The objective is to calculate dSV.
     * 
     * 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).
     *
     */
    float dAV = ptA.distance(ptV);
    float dBV = ptB.distance(ptV);
    float dVS;
    float f = rAS / dAV;
    if (f > 1) {
      // within solvent sphere of atom A
      // calculate point on solvent sphere A projected through ptV
      ptS.set(ptA.x + (ptV.x - ptA.x) * f, ptA.y + (ptV.y - ptA.y) * f, ptA.z
          + (ptV.z - ptA.z) * f);
      // If the distance of this point to B is less than the distance
      // of S to B, then we need to check this point
      if (ptB.distance(ptS) >= rBS)
        return Float.NaN;
      // we are somewhere in the triangle ASB, within the solvent sphere of A
      dVS = solventDistance(rAS, rBS, dAB, dAV, dBV);
      return (voxelIsInTrough(dVS, rAS * rAS, rBS, dAB, dAV) ? dVS : Float.NaN);
    }
    f = rBS / dBV;
    if (f <= 1)
      return Float.NaN; // not within solvent sphere of A or B
    // calculate point on solvent sphere B projected through ptV
    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)
      return Float.NaN;
    // we are somewhere in the triangle ASB, within the solvent sphere of B
    dVS = solventDistance(rBS, rAS, dAB, dBV, dAV);
    return (voxelIsInTrough(dVS, rAS * rAS, rBS, dAB, dAV) ? dVS : Float.NaN);
  }

  private boolean voxelIsInTrough(float dVS, float rAS2, float rBS, float dAB,
                          float dAV) {
    //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);
  }

  private float solventDistance(float rAS, float rBS, float dAB, float dAV, float dBV) {
    double dAV2 = dAV * dAV;
    double rAS2 = rAS * rAS;
    double dAB2 = dAB * dAB;
    double angleVAB = Math.acos((dAV2 + dAB2 - dBV * dBV)
        / (2 * dAV * dAB));
    double angleBAS = Math.acos((dAB2 + rAS2 - rBS * rBS)
        / (2 * dAB * rAS));
    float dVS = (float) Math.sqrt(rAS2 + dAV2 - 2 * rAS * dAV
        * Math.cos(angleBAS - angleVAB));
    return dVS;
  }
}
TOP

Related Classes of org.jmol.jvxl.readers.IsoSolventReader

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.