gaussiansbasis.java

化学图形处理软件

/* $RCSfile$
 * $Author: egonw $
 * $Date: 2007-01-04 18:46:10 +0100 (Thu, 04 Jan 2007) $
 * $Revision: 7636 $
 * 
 * Copyright (C) 2001-2007  The Chemistry Development Kit (CDK) project
 * 
 * Contact: [e-mail]
 * 
 * This program 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.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program 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 Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 *  
 */
package org.openscience.cdk.math.qm;
 
import org.openscience.cdk.math.Matrix;
import org.openscience.cdk.math.Vector;

/**
 * This class contains the information to use gauss function as a base
 * for calculation of quantum mechanics. The function is defined as:
 * <pre>
 * f(x,y,z) = (x-rx)^nx * (y-ry)^ny * (z-rz)^nz * exp(-alpha*(r-ri)^2)
 * </pre>
 *
 * <p>
 * S = &lt;phi_i|phi_j><br>
 * J = &lt;d/dr phi_i | d/dr phi_j><br>
 * V = &lt;phi_i | 1/r | phi_j><br>
 * 
 * @author  Stephan Michels <[e-mail]>
 * @cdk.created 2001-06-14
 *
 * @cdk.keyword Gaussian basis set
 */ 
public class GaussiansBasis implements IBasis
{
  private int count; // number of the basis functions
  private int[] nx; // [base]
  private int[] ny; // [base]
  private int[] nz; // [base]
  private double[] alpha; // [base]
  private double[] norm; // Normalize the bases
  private Vector[] r; // [basis] Positions of base functions

  private int count_atoms; // number of the atoms
  private Vector[] rN; // [atom] Positions of the atoms
  private int[] oz; // [atom] Atomic numbers of the atoms
  
  // For the volume
  private double minx = 0; private double maxx = 0;
  private double miny = 0; private double maxy = 0;
  private double minz = 0; private double maxz = 0;

  public GaussiansBasis()
  {
  }

  /**
   * Set up basis with gauss funktions
   * f(x,y,z) = (x-rx)^nx * (y-ry)^ny * (z-rz)^nz * exp(-alpha*(r-ri)^2).
   *
   * @param atoms The atom will need to calculate the core potential
   */
  public GaussiansBasis(int[] nx, int[] ny, int[] nz, double[] alpha, Vector[] r, org.openscience.cdk.interfaces.IAtom[] atoms)
  {
    setBasis(nx, ny, nz, alpha, r, atoms);
  }

  /**
   * Set up basis with gauss funktions
   * f(x,y,z) = (x-rx)^nx * (y-ry)^ny * (z-rz)^nz * exp(-alpha*(r-ri)^2).
   *
   * @param atoms The atom will need to calculate the core potential
   */
  protected void setBasis(int[] nx, int[] ny, int[] nz, double[] alpha, Vector[] r, org.openscience.cdk.interfaces.IAtom[] atoms)
  {
    int i;

    //count_atoms = molecule.getSize();
    count_atoms = atoms.length;
    
//    logger.debug("Count of atoms: "+count_atoms);
    
    this.rN = new Vector[count_atoms];
    this.oz = new int[count_atoms];
    for(i=0; i<count_atoms; i++)
    { 
      this.rN[i] = (new Vector(atoms[i].getPoint3d())).mul(1.8897);
      this.oz[i] = atoms[i].getAtomicNumber();
//      logger.debug((i+1)+".Atom Z="+this.oz[i]+" r="+(new Vector(atoms[i].getPoint3d()))+"[angstrom]");
    }
//    logger.debug();

    count = Math.min(nx.length,
            Math.min(ny.length,
            Math.min(nz.length,alpha.length)));

//    logger.debug("Count of bases: "+count);

    this.nx = new int[count];
    this.ny = new int[count];
    this.nz = new int[count];
    this.alpha = new double[count];
    //this.atoms = new int[count];
    this.norm = new double[count];
    this.r = new Vector[count];
    this.oz = new int[count];

    for(i=0; i<count; i++)
    {
      this.nx[i] = nx[i];
      this.ny[i] = ny[i];
      this.nz[i] = nz[i];
      this.alpha[i] = alpha[i];
      //this.atoms[i] = atoms[i];
      //this.r[i] = new Vector(atoms[i].getPoint3d()).mul(1.8897);
      this.r[i] = r[i].mul(1.8897);

      norm[i] = Math.sqrt(calcS(i,i));
      if (norm[i]==0d) 
        norm[i] = 1d;
      else
        norm[i] = 1/norm[i];

//      logger.debug((i+1)+".Base nx="+nx[i]+" ny="+ny[i]+" nz="+nz[i]+" alpha="+
//              alpha[i]+" r="+r[i]+" norm="+norm[i]);

      if (i>0)
      {
        minx = Math.min(minx,this.r[i].vector[0]); maxx = Math.max(maxx,this.r[i].vector[0]);
        miny = Math.min(miny,this.r[i].vector[1]); maxy = Math.max(maxy,this.r[i].vector[1]);
        minz = Math.min(minz,this.r[i].vector[2]); maxz = Math.max(maxz,this.r[i].vector[2]);
      }
      else
      {
        minx = r[0].vector[0]; maxx = r[0].vector[0];
        miny = r[0].vector[1]; maxy = r[0].vector[1];
        minz = r[0].vector[2]; maxz = r[0].vector[2];
      }
    }

    minx -= 2; maxx += 2;
    miny -= 2; maxy += 2;
    minz -= 2; maxz += 2;

//    logger.debug();
  }

  /**
   * Gets the number of base vectors.
   */
  public int getSize()
  {
    return count;
  }

  /**
   * Gets the dimension of the volume, which describes the base.
   */
  public double getMinX()
  {
    return minx;
  }

  /**
   * Gets the dimension of the volume, which describes the base.
   */
  public double getMaxX()
  {
    return maxx;
  }

  /**
   * Gets the dimension of the volume, which describes the base.
   */
  public double getMinY()
  {
    return miny;
  }
  
  /**
   * Gets the dimension of the volume, which describes the base.
   */
  public double getMaxY()
  {
    return maxy;
  }

  /**
   * Gets the dimension of the volume, which describes the base.
   */
  public double getMinZ()
  {
    return minz;
  }
  
  /**
   * Gets the dimension of the volume, which describes the base.
   */
  public double getMaxZ()
  {
    return maxz;
  }

  /**
   * Calculates the function value an (x,y,z).
   * @param index The number of the base 
   */
  public double getValue(int index, double x, double y, double z)
  {
    double dx = (x*1.8897)-r[index].vector[0];
    double dy = (y*1.8897)-r[index].vector[1];
    double dz = (z*1.8897)-r[index].vector[2];
    double result = 1d;
    int i;
    for(i=0; i<nx[index]; i++)
      result *= dx;
    for(i=0; i<ny[index]; i++)
      result *= dy;
    for(i=0; i<nz[index]; i++)
      result *= dz;
    return result*Math.exp(-alpha[index]*(dx*dx+dy*dy+dz*dz));
  }

  /**
   * Calculates the function values.
   * @param index The number of the base 
   */
  public Vector getValues(int index, Matrix m)
  {
    if (m.rows!=3)
      return null;

    double x,y,z,dx,dy,dz,mx,my,mz;

    x = m.matrix[0][0]; y = m.matrix[1][0]; z = m.matrix[2][0];
            
    dx = (x*1.8897)-r[index].vector[0];
    dy = (y*1.8897)-r[index].vector[1];
    dz = (z*1.8897)-r[index].vector[2];
    
    Vector result = new Vector(m.columns);
    int i,j; 

    mx = 1d;
    for(i=0; i<nx[index]; i++)
      mx *= dx;

    my = 1d;
    for(i=0; i<ny[index]; i++)
      my *= dy;

    mz = 1d;
    for(i=0; i<nz[index]; i++)
      mz *= dz;

    dx *= dx; dy *= dy; dz *= dz;

    result.vector[0] = mx*my*mz*Math.exp(-alpha[index]*(dx+dy+dz));

    for(j=1; j<m.columns; j++)
    {
      if (x!=m.matrix[0][j])
      {
        x = m.matrix[0][j];
        dx = (x*1.8897)-r[index].vector[0];
        mx = 1d;
        for(i=0; i<nx[index]; i++)
          mx *= dx;
        dx *= dx;
      }

      if (y!=m.matrix[1][j])
      {
        y = m.matrix[1][j];
        dy = (y*1.8897)-r[index].vector[1];
        my = 1d; 
        for(i=0; i<ny[index]; i++)
          my *= dy;
        dy *= dy;
      }

      if (z!=m.matrix[2][j])
      {
        z = m.matrix[2][j];
        dz = (z*1.8897)-r[index].vector[2];
        mz = 1d; 
        for(i=0; i<nz[index]; i++)
          mz *= dz;
        dz *= dz;
      }

      result.vector[j] = mx*my*mz*Math.exp(-alpha[index]*(dx+dy+dz));
    }

    return result;
  }

  /**
   * Gets the position of a base.
   */
  public Vector getPosition(int index)
  {
    return r[index].duplicate().mul(0.52918);
  }

  public double calcD(double normi, double normj, double alphai, double alphaj, Vector ri, Vector rj)
  {
    double dx = ri.vector[0]-rj.vector[0];
    double dy = ri.vector[1]-rj.vector[1];
    double dz = ri.vector[2]-rj.vector[2];
    return Math.exp(-((alphai*alphaj)/(alphai+alphaj))*(dx*dx+dy*dy+dz*dz))*normi*normj;
  }

  /**
   * Transfer equation for the calculation of the overlap integral..
   */
  private double calcI(int ni, int nj, double alphai, double alphaj, double xi, double xj)
  {
    if ((ni<0) || (nj<0))
    {
      System.err.println("Error [Basis.calcI()]: nj="+nj);
      return Double.NaN; // Falls fehlerhafte Parameter
    } 
    double[][] I = new double[nj+1][];
    
    double alphaij = alphai+alphaj;
    double xij = (alphai*xi+alphaj*xj)/alphaij;
    
    I[0] = new double[(nj+ni+1)];
    I[0][0] = Math.sqrt(Math.PI)/Math.sqrt(alphaij); // I(0,0)=G(0)
    
    if ((nj+ni+1)>1)
    {
      I[0][1] = -(xi-xij)*I[0][0]; // I(0,1)=G(1)
      
      // I(0,n)=G(n)
      for(int i=2; i<=(nj+ni); i++)
        I[0][i] = ((i-1)/(2*alphaij))*I[0][i-2]-(xi-xij)*I[0][i-1];
        
      for(int j=1; j<=nj; j++)