3
\brief Enter brief description of file here
9
#include<libciomr/libciomr.h>
10
#include<libint/libint.h>
17
namespace psi { namespace CINTS {
19
static const int use_cca_integrals_standard = (PSI_INTEGRALS_STANDARD == 1);
21
/*-------------------------------
22
Explicit function declarations
23
-------------------------------*/
24
static double **init_bf_norm(int);
25
static double ***init_cart2pureang(int);
26
static double ****init_cc2pp(int);
27
static void init_sparse_cc2pp(int);
28
static double xyz2lm_Coeff(int, int, int, int, int);
32
GTOs.bf_norm = init_bf_norm(BasisSet.max_am);
33
if (BasisSet.puream && UserOptions.symm_ints) {
34
GTOs.cart2pureang = init_cart2pureang(BasisSet.max_am);
35
#if USE_MM && !SPARSE_C2P
36
GTOs.cc2pp = init_cc2pp(BasisSet.max_am);
38
#if USE_MM && SPARSE_C2P
39
init_sparse_cc2pp(BasisSet.max_am);
51
free_matrix(GTOs.bf_norm,BasisSet.max_am);
53
if (BasisSet.puream && UserOptions.symm_ints) {
54
for(l=0;l<BasisSet.max_am;l++)
55
free_block(GTOs.cart2pureang[l]);
56
free(GTOs.cart2pureang);
60
/*!---------------------------------------------------------
61
Computes normalization constants for cartesian Gaussians
62
---------------------------------------------------------*/
63
double **init_bf_norm(int max_am)
66
int am,bf,i,j,l1,m1,n1;
68
bf_norm = (double **) malloc(sizeof(double *)*max_am);
69
for(am=0; am<max_am; am++) {
71
bf_norm[am] = init_array(ioff[am+1]);
72
for(i=0; i<=am; i++) {
77
if (use_cca_integrals_standard)
78
bf_norm[am][bf++] = 1.0;
80
bf_norm[am][bf++] = sqrt(df[2*am]/(df[2*l1]*df[2*m1]*df[2*n1]));
89
double ***init_cart2pureang(int max_am)
92
double ***cart2pureang;
94
cart2pureang = (double ***) malloc(sizeof(double **)*(max_am));
96
cart2pureang[l] = block_matrix(2*l+1,ioff[l+1]);
97
for(l=0;l<max_am;l++) {
102
cart2pureang[l][m+l][bf++] = xyz2lm_Coeff(l,m,l-i,i-j,j);
104
if (UserOptions.print_lvl >= PRINT_DEBUG) {
105
fprintf(outfile," -cart->sph.harm matrix (l = %d):\n",l);
106
print_mat(cart2pureang[l],2*l+1,ioff[l+1],outfile);
114
double ****init_cc2pp(int max_am)
116
int i1,j1,l1,m1,bf1,i2,j2,l2,m2,bf2;
120
cc2pp = (double ****) malloc(sizeof(double ***)*max_am);
121
for(l1=0;l1<max_am;l1++)
122
cc2pp[l1] = (double ***) malloc(sizeof(double **)*max_am);
123
for(l1=0;l1<max_am;l1++)
124
for(l2=0;l2<max_am;l2++) {
127
cc2pp[l1][l2] = block_matrix((2*l1+1)*np2,ioff[l1+1]*nc2);
128
for(m1=-l1;m1<=l1;m1++)
129
for(m2=-l2;m2<=l2;m2++) {
131
for(i1=0;i1<=l1;i1++)
132
for(j1=0;j1<=i1;j1++,bf1++) {
134
for(i2=0;i2<=l2;i2++)
135
for(j2=0;j2<=i2;j2++,bf2++)
136
cc2pp[l1][l2][(m1+l1)*np2+m2+l2][bf1*nc2+bf2] = xyz2lm_Coeff(l1,m1,l1-i1,i1-j1,j1)*
137
xyz2lm_Coeff(l2,m2,l2-i2,i2-j2,j2);
146
void init_sparse_cc2pp(int max_am)
148
/*--- Used for convenience ---*/
149
mat_elem ****cc2pp_sparse;
150
int ***cc2pp_rowlength;
151
mat_elem ****pp2cc_sparse;
152
int ***pp2cc_rowlength;
155
int i1,j1,l1,m1,bf1,i2,j2,l2,m2,bf2;
157
int row_index, col_index;
158
int col_count, total_count, length;
162
/*---------------------------------------------------------------------
163
cc2pp matrices are computed for all l1 and l2, not only for l1 >= l2
165
---------------------------------------------------------------------*/
166
cc2pp_sparse = (mat_elem ****) malloc(sizeof(mat_elem ***)*max_am);
167
cc2pp_rowlength = (int ***) malloc(sizeof(int **)*max_am);
168
for(l1=0;l1<max_am;l1++) {
169
cc2pp_sparse[l1] = (mat_elem ***) malloc(sizeof(mat_elem **)*max_am);
170
cc2pp_rowlength[l1] = (int **) malloc(sizeof(int *)*max_am);
172
for(l1=0;l1<max_am;l1++)
173
for(l2=0;l2<max_am;l2++) {
176
cc2pp_sparse[l1][l2] = (mat_elem **) malloc(sizeof(mat_elem *)*ioff[l1+1]*nc2);
177
cc2pp_rowlength[l1][l2] = (int *) malloc(sizeof(int)*ioff[l1+1]*nc2);
178
/* allocate some space for now, realloc later */
179
increment = (2*l1+1)*np2*ioff[l1+1]*nc2;
180
scratch = (mat_elem *) malloc(sizeof(mat_elem)*increment);
184
for(i1=0;i1<=l1;i1++)
185
for(j1=0;j1<=i1;j1++,bf1++) {
187
for(i2=0;i2<=l2;i2++)
188
for(j2=0;j2<=i2;j2++,bf2++) {
189
row_index = bf1*nc2+bf2;
190
cc2pp_sparse[l1][l2][row_index] = &(scratch[total_count]);
192
for(m1=-l1;m1<=l1;m1++)
193
for(m2=-l2;m2<=l2;m2++) {
194
cc2pp_sparse[l1][l2][row_index][col_count].column = (m1+l1)*np2+m2+l2;
195
cc2pp_sparse[l1][l2][row_index][col_count].row = row_index;
196
value = xyz2lm_Coeff(l1,m1,l1-i1,i1-j1,j1) * xyz2lm_Coeff(l2,m2,l2-i2,i2-j2,j2);
197
if (fabs(value) > ZERO)
198
cc2pp_sparse[l1][l2][row_index][col_count].value = value;
203
if (total_count == length) {
204
scratch = (mat_elem *) realloc(scratch, sizeof(mat_elem)*increment);
208
cc2pp_rowlength[l1][l2][row_index] = col_count;
213
pp2cc_sparse = (mat_elem ****) malloc(sizeof(mat_elem ***)*max_am);
214
pp2cc_rowlength = (int ***) malloc(sizeof(int **)*max_am);
215
for(l1=0;l1<max_am;l1++) {
216
pp2cc_sparse[l1] = (mat_elem ***) malloc(sizeof(mat_elem **)*max_am);
217
pp2cc_rowlength[l1] = (int **) malloc(sizeof(int *)*max_am);
219
for(l1=0;l1<max_am;l1++)
220
for(l2=0;l2<max_am;l2++) {
223
pp2cc_sparse[l1][l2] = (mat_elem **) malloc(sizeof(mat_elem *)*(2*l1+1)*np2);
224
pp2cc_rowlength[l1][l2] = (int *) malloc(sizeof(int)*(2*l1+1)*np2);
225
/* allocate some space for now, realloc later */
226
increment = (2*l1+1)*np2*ioff[l1+1]*nc2;
227
scratch = (mat_elem *) malloc(sizeof(mat_elem)*increment);
230
for(m1=-l1;m1<=l1;m1++)
231
for(m2=-l2;m2<=l2;m2++) {
232
row_index = (m1+l1)*np2+m2+l2;
233
pp2cc_sparse[l1][l2][row_index] = &(scratch[total_count]);
236
for(i1=0;i1<=l1;i1++)
237
for(j1=0;j1<=i1;j1++,bf1++) {
239
for(i2=0;i2<=l2;i2++)
240
for(j2=0;j2<=i2;j2++,bf2++) {
241
pp2cc_sparse[l1][l2][row_index][col_count].row = row_index;
242
pp2cc_sparse[l1][l2][row_index][col_count].column = bf1*nc2+bf2;
243
value = xyz2lm_Coeff(l1,m1,l1-i1,i1-j1,j1) * xyz2lm_Coeff(l2,m2,l2-i2,i2-j2,j2);
244
if (fabs(value) > ZERO)
245
pp2cc_sparse[l1][l2][row_index][col_count].value = value;
250
if (total_count == length) {
251
scratch = (mat_elem *) realloc(scratch, sizeof(mat_elem)*increment);
255
pp2cc_rowlength[l1][l2][row_index] = col_count;
260
GTOs.cc2pp_sparse = cc2pp_sparse;
261
GTOs.cc2pp_rowlength = cc2pp_rowlength;
262
GTOs.pp2cc_sparse = pp2cc_sparse;
263
GTOs.pp2cc_rowlength = pp2cc_rowlength;
270
#define parity(m) ((m)%2 ? -1 : 1) /*Returns (-1)^m */
271
/*!---------------------------------------------------------------------------------------------
272
Computes transformation coefficients from cartesian to real pure angular momentum functions.
273
See IJQC 54, 83 (1995), eqn (15). If m is negative, imaginary part is computed, whereas
274
a positive m indicates that the real part of spherical harmonic Ylm is requested.
275
---------------------------------------------------------------------------------------------*/
276
double xyz2lm_Coeff(int l, int m, int lx, int ly, int lz)
283
double pfac, pfac1, sum, sum1;
286
if ((lx + ly - abs(m))%2)
289
j = (lx + ly - abs(m))/2;
294
/*----------------------------------------------------------------------------------------
295
Checking whether the cartesian polynomial contributes to the requested component of Ylm
296
----------------------------------------------------------------------------------------*/
297
comp = (m >= 0) ? 1 : -1;
298
/* if (comp != ((abs_m-lx)%2 ? -1 : 1))*/
300
if (comp != parity(abs(i)))
303
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]));
304
/* pfac = sqrt(fac[l-abs_m]/(fac[l]*fac[l]*fac[l+abs_m]));*/
307
pfac *= parity((i-1)/2);
313
for(i=0;i<=i_max;i++) {
314
pfac1 = bc[l][i]*bc[i][j];
318
pfac1 *= (parity(i)*fac[2*(l-i)]/fac[l-abs_m-2*i]);
320
k_min = MAX((lx-abs_m)/2,0);
322
for(k=k_min;k<=k_max;k++)
323
sum1 += bc[j][k]*bc[abs_m][lx-2*k]*parity(k);
327
if (use_cca_integrals_standard)
328
sum *= sqrt(df[2*l]/(df[2*lx]*df[2*ly]*df[2*lz]));
333
return M_SQRT2*pfac*sum;
added
removed
src/bin/cints/Tools/gto.cc
