← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/input/build_cart2pureang.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/input/build_cart2pureang.c
 
1
#define EXTERN
 
2
#include <stdio.h>
 
3
#include <stdlib.h>
 
4
#include <libciomr/libciomr.h>
 
5
#include <math.h>
 
6
#include <libqt/qt.h>
 
7
#include "input.h"
 
8
#include "global.h"
 
9
#include "defines.h"
 
10
 
 
11
 
 
12
double xyz2lm_Coeff(int,int,int,int,int);
 
13
static double **bc;
 
14
static double *fac;
 
15
 
 
16
/*-----------------------------------------------------------------------------------------------------------------
 
17
  This function builds cartesian to pure angular momentum transformation matrix
 
18
 -----------------------------------------------------------------------------------------------------------------*/
 
19
 
 
20
void build_cart2pureang()
 
21
{
 
22
  int i,j,k,m;
 
23
  int n;
 
24
  int class,lmax,irr,ao,atom,l;
 
25
  int shell, shell_first, shell_last;
 
26
  int ao_max;
 
27
  int cnt, so_cnt;
 
28
  int ao_offset;
 
29
  int n_max;
 
30
  double min,tmp;
 
31
 
 
32
 
 
33
    /*------------------------
 
34
      Initialize global arrays
 
35
     -------------------------*/
 
36
  
 
37
    num_so_per_irrep = init_int_array(nirreps);
 
38
    cart2pureang = (double ***) malloc(sizeof(double **)*(max_angmom+1));
 
39
    for(l=0;l<=max_angmom;l++)
 
40
      cart2pureang[l] = init_matrix(2*l+1,ioff[l+1]);
 
41
    fac = init_array(2*max_angmom+1);
 
42
    fac[0] = 1.0;
 
43
    for(i=1;i<=2*max_angmom;i++)
 
44
      fac[i] = i*fac[i-1];
 
45
    bc = init_matrix(max_angmom+1,max_angmom+1);
 
46
    for(i=0;i<=max_angmom;i++)
 
47
      for(j=0;j<=i;j++)
 
48
        bc[i][j] = combinations(i,j);
 
49
  
 
50
    /*------------------------
 
51
      Initialize local arrays
 
52
     -------------------------*/
 
53
  
 
54
    
 
55
    /*--------------------------------------------------------------------------------------------------
 
56
      Compute number of pure angular momentum SOs in each symmetry block for each angular momentum type
 
57
     --------------------------------------------------------------------------------------------------*/
 
58
 
 
59
    num_pureang_so = init_int_matrix(max_angmom+1,nirreps);
 
60
    num_redun_so = init_int_matrix(max_angmom+1,nirreps);
 
61
    for(l=0;l<=max_angmom;l++)
 
62
      for(irr=0;irr<nirreps;irr++)
 
63
        num_pureang_so[l][irr] = num_cart_so[l][irr];
 
64
    for(i=2;i<=max_angmom;i++)
 
65
      for(l=i%2;l<i;l+=2)
 
66
        for(irr=0;irr<nirreps;irr++) {
 
67
          num_redun_so[i][irr] += num_pureang_so[l][irr];
 
68
          num_pureang_so[i][irr] -= num_pureang_so[l][irr];
 
69
        }
 
70
 
 
71
 
 
72
    /*--------------------------------------------------------------------------------------------------
 
73
      Compute cartesian to pure angular momentum transformation matrices for each angular momentum type
 
74
     --------------------------------------------------------------------------------------------------*/
 
75
 
 
76
    for(l=0;l<=max_angmom;l++) {
 
77
      ao_max = ioff[l+1];
 
78
      for(m=-l;m<=l;m++)
 
79
        for(ao=0;ao<ao_max;ao++)
 
80
          cart2pureang[l][m+l][ao] = xyz2lm_Coeff(l,m,xexp_ao[l][ao],yexp_ao[l][ao],zexp_ao[l][ao]);
 
81
    }
 
82
 
 
83
 
 
84
    /*-----------------------------
 
85
      Remove after testing is over
 
86
     -----------------------------*/
 
87
    if (print_lvl >= DEBUGPRINT)
 
88
      for(l=0;l<=max_angmom;l++) {
 
89
        fprintf(outfile,"    -Cart2PureAng matrix for l=%d:\n",l);
 
90
        print_mat(cart2pureang[l],2*l+1,ioff[l+1],outfile);
 
91
        fprintf(outfile,"\n");
 
92
      }
 
93
 
 
94
    free(fac);
 
95
    free_matrix(bc,max_angmom+1);
 
96
    return;
 
97
}
 
98
 
 
99
 
 
100
 
 
101
/*---------------------------------------------------------------------------------------------
 
102
  Computes transformation coefficients from cartesian to real pure angular momentum functions.
 
103
  See IJQC 54, 83 (1995), eqn (15). If m is negative, imaginary part is computed, whereas
 
104
  a positive m indicates that the real part of spherical harmonic Ylm is requested.
 
105
 ---------------------------------------------------------------------------------------------*/
 
106
 
 
107
double xyz2lm_Coeff(int l, int m, int lx, int ly, int lz)
 
108
{
 
109
  int i,j,k,i_max;
 
110
  int k_min, k_max;
 
111
  int abs_m;
 
112
  int comp;
 
113
  double pfac, pfac1, sum, sum1;
 
114
 
 
115
  abs_m = abs(m);
 
116
  if ((lx + ly - abs(m))%2)
 
117
    return 0.0;
 
118
  else
 
119
    j = (lx + ly - abs(m))/2;
 
120
 
 
121
  if (j < 0)
 
122
    return 0.0;
 
123
 
 
124
  /*----------------------------------------------------------------------------------------
 
125
    Checking whether the cartesian polynomial contributes to the requested component of Ylm
 
126
   ----------------------------------------------------------------------------------------*/
 
127
  comp = (m >= 0) ? 1 : -1;
 
128
  i = abs_m-lx;
 
129
  if (comp != parity(i))
 
130
    return 0.0;
 
131
 
 
132
  pfac = sqrt(fac[2*lx]*fac[2*ly]*fac[2*lz]*fac[l-abs_m]/(fac[2*l]*fac[l]*fac[lx]*fac[ly]*fac[lz]*fac[l+abs_m]));
 
133
  pfac /= (1 << l);
 
134
 
 
135
  if (m < 0)
 
136
    pfac *= parity((i-1)/2);
 
137
  else
 
138
    pfac *= parity(i/2);
 
139
  
 
140
  i_max = (l-abs_m)/2;
 
141
  sum = 0.0;
 
142
  for(i=0;i<=i_max;i++) {
 
143
    pfac1 = bc[l][i]*bc[i][j];
 
144
    if (pfac1 == 0.0)
 
145
      continue;
 
146
    else
 
147
      pfac1 *= (parity(i)*fac[2*(l-i)]/fac[l-abs_m-2*i]);
 
148
    sum1 = 0.0;
 
149
    k_min = MAX((lx-abs_m)/2,0);
 
150
    k_max = MIN(j,lx/2);
 
151
    for(k=k_min;k<=k_max;k++)
 
152
      sum1 += bc[j][k]*bc[abs_m][lx-2*k]*parity(k);
 
153
    sum += pfac1*sum1;
 
154
  }
 
155
  if (m == 0)
 
156
    return pfac*sum;
 
157
  else
 
158
    return M_SQRT2*pfac*sum;
 
159
}
 
160