← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/oeprop/grid_dens_3d.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/oeprop/grid_dens_3d.c
1
 
#define EXTERN
2
 
#include "includes.h"
3
 
#include "globals.h"
4
 
#include "prototypes.h"
5
 
 
6
 
void compute_grid_dens_3d()
7
 
{
8
 
  int i,j,k,l,ig,jg,ibf,jbf,ib,jb,jlim,kk,ll;
9
 
  int iatom, jatom, iang, jang, i_1stbf, j_1stbf, nbfi, nbfj;
10
 
  int igmin, igmax;
11
 
  int ixm, iym, izm, iind;
12
 
  int ix, iy, iz;
13
 
  int lx, ly, lz;
14
 
  double ax, ay, az, xa, ya, za, ra2;
15
 
  double ai, norm_pf, ang_pf, exponent;
16
 
  double x,y,z;
17
 
  double r,r2,tmp;
18
 
  double *bf_values, *values;
19
 
  double *Density;
20
 
 
21
 
  if (spin_prop)    /* If spin_prop is set, compute the spin density... */
22
 
    Density = Pspin;
23
 
  else              /* .. otherwise compute the electron density */
24
 
    Density = Ptot;
25
 
 
26
 
  /* Initialize some intermediates */
27
 
  grid3d_pts = (double****) malloc(1*sizeof(double***));
28
 
  grid3d_pts[0] = init_box(nix+1,niy+1,niz+1);
29
 
  bf_values = init_array(nbfao);
30
 
 
31
 
  for(ix=0;ix<=nix;ix++) {
32
 
      for(iy=0;iy<=niy;iy++) {
33
 
        x = grid_origin[0] + grid_step_x[0]*ix + grid_step_y[0]*iy;
34
 
        y = grid_origin[1] + grid_step_x[1]*ix + grid_step_y[1]*iy;
35
 
        z = grid_origin[2] + grid_step_x[2]*ix + grid_step_y[2]*iy;
36
 
 
37
 
        for(iz=0;iz<=niz;iz++) {
38
 
          /* Shell loop */
39
 
          for(i=0;i<nshell;i++) {
40
 
              iatom = snuc[i] - 1;
41
 
              ax = geom[iatom][0];
42
 
              ay = geom[iatom][1];
43
 
              az = geom[iatom][2];
44
 
              iang = stype[i]-1;
45
 
              izm = 1;
46
 
              iym = iang+1;
47
 
              ixm = iym*iym;
48
 
              i_1stbf = sloc[i] - 1;
49
 
              nbfi = (iang+2)*(iang+1)/2;
50
 
              igmin = sprim[i] - 1;
51
 
              igmax = igmin + snumg[i] - 1;
52
 
    
53
 
              xa = x - ax;
54
 
              ya = y - ay;
55
 
              za = z - az;
56
 
              ra2 = xa*xa + ya*ya + za*za;
57
 
 
58
 
              /*--- Compute exponentional factor - loop over primitives ---*/
59
 
              exponent = 0.0;
60
 
              for(ig=igmin;ig<=igmax;ig++) {
61
 
                  ai = exps[ig];
62
 
                  exponent += contr[ig]*exp(-ai*ra2);
63
 
              }
64
 
 
65
 
              /*--- Loop over basis functions in the shell ---*/
66
 
              values = &(bf_values[i_1stbf]);
67
 
              for(ibf=0;ibf<nbfi;ibf++,values++) {
68
 
                  norm_pf = norm_bf[iang][ibf];
69
 
                  lx = xpow_bf[iang][ibf];
70
 
                  ly = ypow_bf[iang][ibf];
71
 
                  lz = zpow_bf[iang][ibf];
72
 
                  tmp = 1.0;
73
 
                  switch (lx) {
74
 
                  case 0:
75
 
                      break;
76
 
                  case 1:
77
 
                      tmp *= xa;
78
 
                      break;
79
 
                  case 2:
80
 
                      tmp *= xa*xa;
81
 
                      break;
82
 
                  case 3:
83
 
                      tmp *= xa*xa*xa;
84
 
                      break;
85
 
                  case 4:
86
 
                      tmp *= (xa*xa)*(xa*xa);
87
 
                      break;
88
 
                  default:
89
 
                      tmp *= pow(xa,lx);
90
 
                  }
91
 
                  ang_pf = tmp;
92
 
                  tmp = 1.0;
93
 
                  switch (ly) {
94
 
                  case 0:
95
 
                      break;
96
 
                  case 1:
97
 
                      tmp *= ya;
98
 
                      break;
99
 
                  case 2:
100
 
                      tmp *= ya*ya;
101
 
                      break;
102
 
                  case 3:
103
 
                      tmp *= ya*ya*ya;
104
 
                      break;
105
 
                  case 4:
106
 
                      tmp *= (ya*ya)*(ya*ya);
107
 
                      break;
108
 
                  default:
109
 
                      tmp *= pow(ya,ly);
110
 
                  }
111
 
                  ang_pf *= tmp;
112
 
                  tmp = 1.0;
113
 
                  switch (lz) {
114
 
                  case 0:
115
 
                      break;
116
 
                  case 1:
117
 
                      tmp *= za;
118
 
                      break;
119
 
                  case 2:
120
 
                      tmp *= za*za;
121
 
                      break;
122
 
                  case 3:
123
 
                      tmp *= za*za*za;
124
 
                      break;
125
 
                  case 4:
126
 
                      tmp *= (za*za)*(za*za);
127
 
                      break;
128
 
                  default:
129
 
                      tmp *= pow(za,lz);
130
 
                  }
131
 
                  ang_pf *= tmp;
132
 
                  
133
 
                  *values = norm_pf * ang_pf * exponent;
134
 
              }
135
 
          } /*--- end of shell loop ---*/
136
 
          tmp = 0;
137
 
          for(i=0;i<nbfao;i++)
138
 
            for(j=0;j<=i;j++)
139
 
              tmp += bf_values[i] * Density[ioff[i]+j] * bf_values[j] * (i!=j ? 2.0 : 1.0);
140
 
          grid3d_pts[0][ix][iy][iz] = tmp;
141
 
          
142
 
          x += grid_step_z[0];
143
 
          y += grid_step_z[1];
144
 
          z += grid_step_z[2];
145
 
        }
146
 
      }
147
 
  }
148
 
}
149