lebedev.c

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

/* These are the comments from the F77 translation by Dr. Christoph van Wuellen
chvd
chvd   This subroutine is part of a set of subroutines that generate
chvd   Lebedev grids [1-6] for integration on a sphere. The original 
chvd   C-code [1] was kindly provided by Dr. Dmitri N. Laikov and 
chvd   translated into fortran by Dr. Christoph van Wuellen.
chvd   This subroutine was translated from C to fortran77 by hand.
chvd
chvd   Users of this code are asked to include reference [1] in their
chvd   publications, and in the user- and programmers-manuals 
chvd   describing their codes.
chvd
chvd   This code was distributed through CCL (http://www.ccl.net/).
chvd
chvd   [1] V.I. Lebedev, and D.N. Laikov
chvd       "A quadrature formula for the sphere of the 131st
chvd        algebraic order of accuracy"
chvd       Doklady Mathematics, Vol. 59, No. 3, 1999, pp. 477-481.
chvd
chvd   [2] V.I. Lebedev
chvd       "A quadrature formula for the sphere of 59th algebraic
chvd        order of accuracy"
chvd       Russian Acad. Sci. Dokl. Math., Vol. 50, 1995, pp. 283-286. 
chvd
chvd   [3] V.I. Lebedev, and A.L. Skorokhodov
chvd       "Quadrature formulas of orders 41, 47, and 53 for the sphere"
chvd       Russian Acad. Sci. Dokl. Math., Vol. 45, 1992, pp. 587-592. 
chvd
chvd   [4] V.I. Lebedev
chvd       "Spherical quadrature formulas exact to orders 25-29"
chvd       Siberian Mathematical Journal, Vol. 18, 1977, pp. 99-107. 
chvd
chvd   [5] V.I. Lebedev
chvd       "Quadratures on a sphere"
chvd       Computational Mathematics and Mathematical Physics, Vol. 16,
chvd       1976, pp. 10-24. 
chvd
chvd   [6] V.I. Lebedev
chvd       "Values of the nodes and weights of ninth to seventeenth 
chvd        order Gauss-Markov quadrature formulae invariant under the
chvd        octahedron group with inversion"
chvd       Computational Mathematics and Mathematical Physics, Vol. 15,
chvd       1975, pp. 44-51.
chvd
*/  

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define INTHISFILE

#ifdef DEFINE_MAIN
int
main ()
{
#define NMAX 65
#define MMAX ((NMAX*2+3)*(NMAX*2+3)/3)
  int i, j, k, m, M, N;
  static int NW[NMAX+1] = {
       0,   6,  14,  26,  38,  50,  74,  86, 110, 146, 170, 194, 230, 266, 302,
     350, 386, 434, 482, 530, 590, 650, 698, 770, 830, 890, 974,1046,1118,1202,
    1274,1358,1454,1538,1622,1730,1814,1910,2030,2126,2222,2354,2450,2558,2702,
    2810,2930,3074,3182,3314,3470,3590,3722,3890,4010,4154,4334,4466,4610,4802,
    4934,5090,5294,5438,5606,5810
  };
  static double x[MMAX], y[MMAX], z[MMAX], w[MMAX];
  static double xn[MMAX*(NMAX+1)], yn[MMAX*(NMAX+1)], zn[MMAX*(NMAX+1)];
  static double s[NMAX+2];
  double S, S0, r, rmax;
  for (N = 1; N <= NMAX; N++) {
    M = NW[N];
    if (Lebedev_Laikov_sphere (M, x, y, z, w) == 0) continue;
    s[0] = 1.0;
    for (k = 1; k <= N+1; k++) {
      s[k] = (2.0*k-1.0)*s[k-1];
    }
    for (m = 0; m < M; m++) {
      xn[m*(N+1)] = 1.0;
      yn[m*(N+1)] = 1.0;
      zn[m*(N+1)] = 1.0;
      for (k = 1; k <= N; k++) {
        xn[k+m*(N+1)] = xn[k-1+m*(N+1)]*x[m]*x[m];
        yn[k+m*(N+1)] = yn[k-1+m*(N+1)]*y[m]*y[m];
        zn[k+m*(N+1)] = zn[k-1+m*(N+1)]*z[m]*z[m];
      }
    }
    rmax = 0.0;
    for (i = 0; i <= N; i++) {
      for (j = 0; j <= N-i; j++) {
        k = N-i-j;
        for (S = 0.0, m = 0; m < M; m++) {
          S += w[m]*xn[i+m*(N+1)]*yn[j+m*(N+1)]*zn[k+m*(N+1)];
        }
        S0 = s[i]*s[j]*s[k]/s[1+i+j+k];
        r = fabs ((S-S0)/S0);
        if (r > rmax) rmax = r;
#if 0
        printf (" %2d %2d %2d  %20.15lf  %20.15lf\n", i, j, k, S0, r);
#endif
      }
    }
    printf (" M = %4d  rmax = %8.1e\n", M, rmax);
  }
  return 0;
}
#endif

#ifdef INTHISFILE
/*
** Lebedev_Laikov_npoint
**
** lvalue : grid complete through this value of angular momentum quantum number l.
**
** return value : number of points in sought Lebedev-Laikov grid.
**
*/
int Lebedev_Laikov_npoint(int lvalue)
{

  int fraction, tmp_lvalue;
  tmp_lvalue = lvalue;
  
  if (lvalue <3) {
      return 6;
    }
  else if (lvalue <= 31) {
      /* grids complete through l = 2m+1 */
      tmp_lvalue -= 1;
      fraction = tmp_lvalue%2;
      tmp_lvalue /= 2;
      if (fraction) tmp_lvalue++; /* round up */
      switch (tmp_lvalue) {
      case 0:  return 6;
      case 1:  return 6;
      case 2:  return 14;
      case 3:  return 26;
      case 4:  return 38;
      case 5:  return 50;
      case 6:  return 74;
      case 7:  return 86;
      case 8:  return 110;
      case 9:  return 146;
      case 10: return 170;
      case 11: return 194;
      case 12: return 230;
      case 13: return 266;
      case 14: return 302;
      case 15: return 350;
        } 
    }
  else if (lvalue <= 131) {
      /* grids complete through l = 6m+5 */
      tmp_lvalue -= 5;
      fraction = tmp_lvalue%6;
      tmp_lvalue /= 6;
      if (fraction) tmp_lvalue++; /* round up */
      switch (tmp_lvalue) {
      case 5:  return 434;
      case 6:  return 590;
      case 7:  return 770;
      case 8:  return 974;
      case 9:  return 1202;
      case 10: return 1454;
      case 11: return 1730;
      case 12: return 2030;
      case 13: return 2354;
      case 14: return 2702;
      case 15: return 3074;
      case 16: return 3470;
      case 17: return 3890;
      case 18: return 4334;
      case 19: return 4802;
      case 20: return 5294;
      case 21: return 5810;
        }
    }
  else {
      printf(" Lebedev_Laikov_npoint: lvalue > 131.  No grids of this type available.\n");
      abort();
    }
  return 0;
}

int Lebedev_Laikov_Oh (int n, double a, double b, double v,
                       double *x, double *y, double *z, double *w);

int
Lebedev_Laikov_sphere (int N, double *X, double *Y, double *Z, double *W)
{
#define A n+=Lebedev_Laikov_Oh(
#define B ,X+n,Y+n,Z+n,W+n);
  int n;
  n = 0;
  switch (N) {
    case 6:
    A 1, 0.0                  , 0.0                  , 0.1666666666666667e+0 B
    break;
    case 14:
    A 1, 0.0                  , 0.0                  , 0.6666666666666667e-1 B