gaussiansbasis.java

化学图形处理软件

      {
        I[j] = new double[(nj+ni+1)-j];
        
        for(int i=0; i<=(nj+ni-j); i++)
          I[j][i] = I[nj-1][ni+1]+(xi-xj)*I[nj-1][ni];
      }   
    }

    return I[nj][ni];
  }

  public double calcS(int i, int j)
  {
    //logger.debug("S: i="+i+" j="+j+" r[i]="+r[i]+" r[j]="+r[j]);
    return calcD(norm[i], norm[j], alpha[i],alpha[j],r[i],r[j]) * 
           calcI(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0]) *
           calcI(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1]) * 
           calcI(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2]);
  }

  /**
   * Transfer equation the the calculation of the impulse
   */
  public double calcJ(int ni, int nj, double alphai, double alphaj, double xi, double xj)
  {
    if (ni>1)
      return -4d*alphai*alphai*   calcI(ni+2,nj,alphai,alphaj,xi,xj)
             +2d*alphai*(2*ni+1)* calcI(ni  ,nj,alphai,alphaj,xi,xj)
             -ni*(ni-1)*          calcI(ni-2,nj,alphai,alphaj,xi,xj);
    else if (ni==1)
      return -4d*alphai*alphai*   calcI(3   ,nj,alphai,alphaj,xi,xj)
             +6d*alphai*          calcI(1   ,nj,alphai,alphaj,xi,xj);
    else if (ni==0)
      return -4d*alphai*alphai*   calcI(2   ,nj,alphai,alphaj,xi,xj)
             +2d*alphai*          calcI(0   ,nj,alphai,alphaj,xi,xj);
             
    System.err.println("Error [Basis.calcJ]: ni="+ni);
    return Double.NaN; // Falls fehlerhafte Parameter
  }

  public double calcJ(int i, int j)
  {
    return calcD(norm[i], norm[j], alpha[i],alpha[j],r[i],r[j])*
           (calcJ(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0])*
            calcI(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1])*
            calcI(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2])+

            calcI(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0])*
            calcJ(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1])*
            calcI(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2])+

            calcI(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0])*
            calcI(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1])*
            calcJ(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2]));
  }

  /**
   * Transfer equation for the calculation of core potentials
   */
  public double calcG(int n, double t, double alphai, double alphaj, double xi, double xj, double xN)
  {
    if (n>1)
      return ((n-1)/(2*(alphai+alphaj)))*                       calcG(n-2, t, alphai, alphaj, xi, xj, xN)
             -(n-1)*t*t*                                        calcG(n-2, t, alphai, alphaj, xi, xj, xN)
             +(((alphai*xi+alphaj*xj)/(alphai+alphaj))-xi)*     calcG(n-1, t, alphai, alphaj, xi, xj, xN)
             +(((alphai*xi+alphaj*xj)/(alphai+alphaj))-xN)*t*t* calcG(n-1, t, alphai, alphaj, xi, xj, xN);
             
    else if (n==1)
      return (((alphai*xi+alphaj*xj)/(alphai+alphaj))-xi)*      calcG(0, t, alphai, alphaj, xi, xj, xN)
             +(((alphai*xi+alphaj*xj)/(alphai+alphaj))-xN)*t*t* calcG(0, t, alphai, alphaj, xi, xj, xN);
             
    else if (n==0)
      return Math.sqrt(Math.PI)/Math.sqrt(alphai+alphaj);
      
    System.err.println("Error [Basis.calcG]: n="+n);
    return Double.NaN; // Falls fehlerhafte Parameter
  }

  /**
   * Transfer equation for the calculation of core potentials.
   */
  private double calcI(int ni, int nj, double t, double alphai, double alphaj, 
                        double xi, double xj, double xN)
  {
    if (nj>0)
      return calcI(ni+1, nj-1, t, alphai, alphaj, xi, xj, xN)+
             (xi-xj)*calcI(ni, nj-1, t, alphai, alphaj, xi, xj, xN);
             
    else if (nj==0)
      return calcG(ni, t, alphai, alphaj, xi, xj, xN);
      
    System.err.println("Error [Basis.calcI()]: nj="+nj);
    return Double.NaN; // Falls fehlerhafte Parameter
  }

  /**
   * Calculates the core potential.
   * It use a 10 point Simpson formula.
   *
   * @param i Index of the first base
   * @param j Index of the second base
   * @param rN Position the core potential
   * @param Z Atomic number of the nucleous
   */
  public double calcV(int i, int j, Vector rN, double Z)
  {
    double f,t;

    double sum1,sum2;
    double f1,f2;

    int steps = 10;
    double h = 1d/steps;

    double alphaij = alpha[i]+alpha[j];

    double rxij = (alpha[i]*r[i].vector[0]+alpha[j]*r[j].vector[0])/alphaij;
    double ryij = (alpha[i]*r[i].vector[1]+alpha[j]*r[j].vector[1])/alphaij;
    double rzij = (alpha[i]*r[i].vector[2]+alpha[j]*r[j].vector[2])/alphaij;

    double X = alphaij*((rxij-rN.vector[0])*(rxij-rN.vector[0]) +
                        (ryij-rN.vector[1])*(ryij-rN.vector[1]) +
                        (rzij-rN.vector[2])*(rzij-rN.vector[2]));

    double C = 2*calcD(norm[i], norm[j], 
              alpha[i], alpha[j], r[i], r[j])*Math.sqrt(alphaij)/Math.sqrt(Math.PI);

    sum1 = 0;
    for(f = 1; f<steps; f=f+2)
    {
      t = f*h;
      sum1 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j], 
                                     r[i].vector[0], r[j].vector[0], rN.vector[0]) *
                               calcI(ny[i], ny[j], t, alpha[i], alpha[j], 
                                     r[i].vector[1], r[j].vector[1], rN.vector[1]) *
                               calcI(nz[i], nz[j], t, alpha[i], alpha[j], 
                                     r[i].vector[2], r[j].vector[2], rN.vector[2]);
    }

    sum2 = 0;
    for(f = 2; f<steps; f=f+2)
    {
      t = f*h;
      sum2 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j],
                                     r[i].vector[0], r[j].vector[0], rN.vector[0]) *
                               calcI(ny[i], ny[j], t, alpha[i], alpha[j], 
                                     r[i].vector[1], r[j].vector[1], rN.vector[1]) *
                               calcI(nz[i], nz[j], t, alpha[i], alpha[j], 
                                     r[i].vector[2], r[j].vector[2], rN.vector[2]);
    }

    t = 0d;
    f1 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j],
                                     r[i].vector[0], r[j].vector[0], rN.vector[0]) *
                               calcI(ny[i], ny[j], t, alpha[i], alpha[j], 
                                     r[i].vector[1], r[j].vector[1], rN.vector[1]) *
                               calcI(nz[i], nz[j], t, alpha[i], alpha[j], 
                                     r[i].vector[2], r[j].vector[2], rN.vector[2]);

    t = 1d;
    f2 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j],
                                     r[i].vector[0], r[j].vector[0], rN.vector[0]) *
                               calcI(ny[i], ny[j], t, alpha[i], alpha[j], 
                                     r[i].vector[1], r[j].vector[1], rN.vector[1]) *
                               calcI(nz[i], nz[j], t, alpha[i], alpha[j], 
                                     r[i].vector[2], r[j].vector[2], rN.vector[2]);

    return (h/3)*(f1 + 4*sum1 + 2*sum2 + f2)*Z*C;
  }

  /**
   * Calculates the core potential.
   * It use a 10 point Simpson formula.
   *
   * @param i Index of the first base
   * @param j Index of the second base
   */
  public double calcV(int i, int j)
  {
    double result = 0d;
    for(int k=0; k<count_atoms; k++)
    { 
      //logger.debug("k="+k+" r="+r[k]);
      result += calcV(i,j, rN[k], oz[k]);
    }
    return -result; // Vorsicht negatives Vorzeichen
  }

  /**
   * Transfer equation for a four center integral.
   */
  public double calcG(int n, int m, double u, 
     double alphai, double alphaj, double alphak, double alphal, double xi, double xj, double xk, double xl)
  {
    if ((n<0) || (m<0))
    {
//      logger.debug("Error(CalcG):Bad parameter n="+n+" m="+m);
      return Double.NaN;
    }

    double alphaij = alphai+alphaj;
    double alphakl = alphak+alphal;

    double xij = (alphai*xi+alphaj*xj)/alphaij;
    double xkl = (alphak*xk+alphal*xl)/alphakl;

    double C00 = (xij-xi)-((u*u*alphakl*(xij-xkl))/(u*u*(alphaij+alphakl)+alphaij*alphakl));
    double Cs00 = (xkl-xk)+((u*u*alphaij*(xij-xkl))/(u*u*(alphaij+alphakl)+alphaij*alphakl));
    double B00 = u*u/(2*(u*u*(alphaij+alphakl)+alphaij*alphakl));
    double B10 = (u*u+alphakl)/(2*(u*u*(alphaij+alphakl)+alphaij*alphakl));
    double Bs01 = (u*u+alphaij)/(2*(u*u*(alphaij+alphakl)+alphaij*alphakl));

    double[][] G = new double[n+1][m+1];

    int i,j;

    G[0][0] = 1d;
    
    // Nach 1)
    if (n>0)
      G[1][0] = C00;

    // Nach 1)
    for(i=2; i<=n; i++)
      G[i][0] = (i-1)*B10   *G[i-2][0]+
                C00         *G[i-1][0];

    // Nach 2)
    if (m>0)
      G[0][1] = Cs00;

    // Nach 2)
    for(i=2; i<=m; i++)
      G[0][i] = (i-1)*Bs01 *G[0][i-2]+
                Cs00       *G[0][i-1];

    // Nach 1)
    if (n>0)
      for(i=1; i<=m; i++)
        G[1][i] = i*B00       *G[0][i-1]+
                  C00         *G[0][i];

    // Nach 1)
    for(i=2; i<=n; i++)
      for(j=1; j<=m; j++)
        G[i][j] = (i-1)*B10   *G[i-2][j]+
                  j*B00       *G[i-1][j-1]+
                  C00         *G[i-1][j];

    return G[n][m];
  }  

  /**
   * Transfer equation for a four center integral.
   */
  public double calcI(int ni, int nj, int nk, int nl, double u, 
                      double alphai, double alphaj, double alphak, double alphal,
                      double xi, double xj, double xk, double xl)
  {
    if (nj>0)
      return          calcI(ni+1,nj-1,nk  ,nl  ,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl)+
             (xj-xi)*calcI(ni  ,nj-1,nk  ,nl  ,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl);

    if (nl>0)
      return          calcI(ni  ,nj  ,nk+1,nl-1,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl)+
             (xl-xk)*calcI(ni  ,nj  ,nk  ,nl-1,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl);

    if ((ni==0) && (nj==0) && (nk==0) && (nl==0))
      return 1d;

    if ((nj==0) && (nl==0))
      return calcG(ni,nk,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl);

//    logger.debug("Error(CalcI):Bad parameter ni="+ni+" nj="+nj+" nk="+nk+" nl="+nl);
    return Double.NaN;
  }

  public double calcI(int i, int j, int k, int l)
  {
    double f,t;
    
    double sum1,sum2;
    double f1,f2;
    
    // Berechnen der Integration nach Simson
    int steps = 10;
    double h = 1d/steps;
    //double h2 = 2d*h;

    double alphaij = alpha[i]+alpha[j];
    double alphakl = alpha[k]+alpha[l];

    double rxij = (alpha[i]*r[i].vector[0]+alpha[j]*r[j].vector[0])/alphaij;
    double ryij = (alpha[i]*r[i].vector[1]+alpha[j]*r[j].vector[1])/alphaij;
    double rzij = (alpha[i]*r[i].vector[2]+alpha[j]*r[j].vector[2])/alphaij;

    double rxkl = (alpha[k]*r[k].vector[0]+alpha[l]*r[l].vector[0])/alphakl;
    double rykl = (alpha[k]*r[k].vector[1]+alpha[l]*r[l].vector[1])/alphakl;
    double rzkl = (alpha[k]*r[k].vector[2]+alpha[l]*r[l].vector[2])/alphakl;

    double alpha0 = alphaij*alphakl/(alphaij+alphakl);

    double X = alpha0*((rxij-rxkl)*(rxij-rxkl) +
                       (ryij-rykl)*(ryij-rzkl) +
                       (rzij-rzkl)*(rzij-rzkl));
    
    double C = (Math.PI*Math.PI*Math.PI/Math.pow((alpha[i]+alpha[j])*(alpha[k]+alpha[l]),1.5))*
               Math.sqrt(alpha0)*
               calcD(norm[i], norm[j], alpha[i], alpha[j], r[i], r[j])*
               calcD(norm[k], norm[l], alpha[k], alpha[l], r[k], r[l])*
               (2d/Math.sqrt(Math.PI));

    sum1 = 0;
    for(f = 1; f<steps; f=f+2)
    {
      t = f*h;
      sum1 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
                               calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
                               calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);

    }

    sum2 = 0;
    for(f = 2; f<steps; f=f+2)
    {
      t = f*h;
      sum2 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
                               calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
                               calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);
    }                                
    
    t = 0d;
    f1 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
                               calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
                               calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);
                                     
    t = 1d;                          
    f2 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
                               calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
                               calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
                 r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);

    return C * (h/3)*(f1 + 4*sum1 + 2*sum2 + f2);
  }
}