← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/cusp/phi.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cusp/phi.c
1
 
/*
2
 
** phi(): Compute the values of the atomic orbitals at a given point.
3
 
** -TDC, June 2001
4
 
*/
5
 
 
6
 
#include <stdio.h>
7
 
#include <string.h>
8
 
#include <math.h>
9
 
#include <libipv1/ip_lib.h>
10
 
#include <libciomr/libciomr.h>
11
 
#include <libchkpt/chkpt.h>
12
 
#include <libint/libint.h>  /* for the maximum angluar momentum, LIBINT_MAX_AM */
13
 
#include <psifiles.h>
14
 
 
15
 
#define MAXFACT 100
16
 
 
17
 
int nao, natom, nprim, nshell;
18
 
int *stype, *snuc, *snumg, *sprim;
19
 
double *a, *c, **geom, df[MAXFACT], **norm;
20
 
int **xexp, **yexp, **zexp, *l_length;
21
 
 
22
 
void setup_phi(void)
23
 
{
24
 
  static int done=0;
25
 
  int i,l,j,ao;
26
 
 
27
 
  if(done) return;
28
 
 
29
 
  chkpt_init(PSIO_OPEN_OLD);
30
 
  nao = chkpt_rd_nao();
31
 
  natom = chkpt_rd_natom();
32
 
  nprim = chkpt_rd_nprim();
33
 
  nshell = chkpt_rd_nshell();
34
 
  stype = chkpt_rd_stype();
35
 
  snuc = chkpt_rd_snuc();
36
 
  snumg = chkpt_rd_snumg();
37
 
  sprim = chkpt_rd_sprim();
38
 
  a = chkpt_rd_exps();
39
 
  c = chkpt_rd_contr();
40
 
  geom = chkpt_rd_geom();
41
 
  chkpt_close();
42
 
 
43
 
  /* compute double factorial --- df[n] = (n-1)!! */
44
 
  df[0] = 1.0;  df[1] = 1.0;  df[2] = 1.0;
45
 
  for(i=3; i < MAXFACT; i++) df[i] = (i-1) * df[i-2];
46
 
  
47
 
  /* Compute the AO ordering within a shell and the angle-independent
48
 
     normalization contants */
49
 
  xexp = (int **) malloc(LIBINT_MAX_AM * sizeof(int *));
50
 
  yexp = (int **) malloc(LIBINT_MAX_AM * sizeof(int *));
51
 
  zexp = (int **) malloc(LIBINT_MAX_AM * sizeof(int *));
52
 
  norm = (double **) malloc(LIBINT_MAX_AM * sizeof(double *));
53
 
  l_length = init_int_array(LIBINT_MAX_AM);
54
 
  l_length[0] = 1;
55
 
  for(l=0; l < (LIBINT_MAX_AM); l++) {
56
 
 
57
 
      if(l) l_length[l] = l_length[l-1] + l + 1;
58
 
      
59
 
      xexp[l] = init_int_array(l_length[l]);
60
 
      yexp[l] = init_int_array(l_length[l]);
61
 
      zexp[l] = init_int_array(l_length[l]);
62
 
      norm[l] = init_array(l_length[l]);
63
 
 
64
 
      for(i=0,ao=0; i <= l; i++) {
65
 
 
66
 
          for(j=0; j <= i; j++) {
67
 
 
68
 
              xexp[l][ao] = l - i;
69
 
              yexp[l][ao] = i - j;
70
 
              zexp[l][ao] = j;
71
 
 
72
 
              norm[l][ao] = sqrt(df[2*l]/(df[2*(l-i)]*df[2*(i-j)]*df[2*j]));
73
 
 
74
 
/*            printf("%d %d %20.10f\n", l, ao, norm[l][ao]); */
75
 
 
76
 
              ao++;
77
 
            }
78
 
        }
79
 
    }
80
 
 
81
 
  done = 1;
82
 
 
83
 
}
84
 
 
85
 
void compute_phi(double *phi, double x, double y, double z)
86
 
{
87
 
  int i, j, l, firstp, lastp;
88
 
  int am, atom, ao;
89
 
  double cexpr, dx, dy, dz, rr;
90
 
 
91
 
  setup_phi();
92
 
 
93
 
  /* Loop over the shells */
94
 
  for(i=0,ao=0; i < nshell; i++) {
95
 
 
96
 
      firstp = sprim[i]-1;
97
 
      lastp = firstp + snumg[i];
98
 
 
99
 
      atom = snuc[i] - 1;
100
 
      dx = x - geom[atom][0];
101
 
      dy = y - geom[atom][1];
102
 
      dz = z - geom[atom][2];
103
 
      
104
 
      /* Compute r^2 for this center */
105
 
      rr = dx*dx + dy*dy + dz*dz;
106
 
 
107
 
      /* Angular momentum level for this shell */
108
 
      am = stype[i] - 1;
109
 
 
110
 
      if(am > LIBINT_MAX_AM) {
111
 
          printf("Angular momentum max. exceeded.\n");
112
 
          exit(PSI_RETURN_FAILURE);
113
 
        }
114
 
 
115
 
      /* Loop over the primitive Gaussians in this shell */
116
 
      cexpr = 0;
117
 
      for(j=firstp; j < lastp; j++) {
118
 
          cexpr += c[j] * exp(-a[j] * rr);
119
 
        }
120
 
 
121
 
      for(l=0; l < l_length[am]; l++) {
122
 
          phi[ao+l] += pow(dx,(double) xexp[am][l]) *
123
 
                       pow(dy,(double) yexp[am][l]) *
124
 
                       pow(dz,(double) zexp[am][l]) *
125
 
                       cexpr * norm[am][l];
126
 
        }
127
 
 
128
 
      ao += l_length[am];
129
 
 
130
 
    }
131
 
}
132