3
\brief Enter brief description of file here
11
#include<libciomr/libciomr.h>
12
#include<libchkpt/chkpt.h>
14
#include<libiwl/iwl.h>
15
#include<libint/libint.h>
23
#include"quartet_data.h"
24
#include"norm_quartet.h"
25
#include "data_structs.h"
27
#include"taylor_fm_eval.h"
32
#include"quartet_permutations.h"
33
#include"mkpt2_ints.h"
40
/*-------------------------------
41
Explicit function declarations
42
-------------------------------*/
43
extern void *mkpt2_ints_thread(void *);
45
/*--------------------------------
46
Varixbles common to all threads
47
--------------------------------*/
48
pthread_mutex_t mkpt2_energy_mutex;
49
pthread_mutex_t *mkpt2_sindex_mutex;
50
pthread_cond_t mkpt2_energy_cond;
51
MkPT2_Status_t MkPT2_Status;
56
unsigned long int n_ij;
57
unsigned long int n_xy;
58
unsigned long int n_ab;
68
pthread_attr_t thread_attr;
71
/*--- Various data structures ---*/
73
long int libint_memory;
75
int max_num_prim_comb;
80
int num_ibatch, num_i_per_ibatch, ibatch, ibatch_length;
82
int mo_i, mo_j, mo_a, mo_b, mo_ij;
83
int rs_offset, rsi_offset, rsp_offset;
87
double *raw_data; /* pointer to the unnormalized taregt quartet of integrals */
88
double *data; /* pointer to the transformed normalized target quartet of integrals */
90
REALTYPE *target_ints; /* Pointer to the location of the target quartet on the stack of
91
integrals quartets if libint.a is using other than regular doubles */
97
double *rsiq_buf; /* buffer for (rs|iq) integrals, where r,s run over shell sets,
98
i runs over I-batch, q runs over all AOs */
100
double **ix_buf; /* buffer for one |ix) ket */
104
double ixjy, iyjx, pfac;
106
double temp1,temp2,*iq_row,*ip_row;
113
init_Taylor_Fm_Eval(BasisSet.max_am*4-4,UserOptions.cutoff);
115
init_fjt(BasisSet.max_am*4);
119
timer_on("Schwartz");
121
timer_off("Schwartz");
122
MkPT2_Status.num_arrived = 0;
125
/*-------------------------
126
Allocate data structures
127
-------------------------*/
128
n_ij = MOInfo.ndocc*(MOInfo.ndocc+1)/2;
129
n_ab = MOInfo.nuocc*(MOInfo.nuocc+1)/2;
130
n_xy = MOInfo.num_mo*(MOInfo.num_mo+1)/2;
131
max_bf_per_shell = ioff[BasisSet.max_am];
132
/*--- Use this dirty trick to get how much memory integrals library needs ---*/
133
max_num_prim_comb = (BasisSet.max_num_prims*
134
BasisSet.max_num_prims)*
135
(BasisSet.max_num_prims*
136
BasisSet.max_num_prims);
137
libint_memory = libint_storage_required(BasisSet.max_am-1,max_num_prim_comb);
138
UserOptions.memory -= libint_memory*UserOptions.num_threads;
140
Minimum number of I-batches -
141
take sizes of rsiq_buf, rsix_buf, jsix_buf,
142
jyix_buf ,abij_buf, asij_buf, rsij_buf, xy_buf into account
144
fprintf(outfile,"\n Computing MkPT2 integrals\n");
145
num_i_per_ibatch = UserOptions.memory / (UserOptions.num_threads*
147
BasisSet.num_ao*max_bf_per_shell*max_bf_per_shell + /*rsiq*/
148
MOInfo.num_mo*max_bf_per_shell*max_bf_per_shell + /*rsix*/
149
MOInfo.ndocc*max_bf_per_shell*max_bf_per_shell /*rsij*/
151
MOInfo.num_mo*MOInfo.ndocc*BasisSet.num_ao + /*jsix*/
152
MOInfo.nuocc*MOInfo.ndocc*BasisSet.num_ao + /*asij*/
153
n_ab*MOInfo.ndocc + /*abij*/
154
MOInfo.num_mo*MOInfo.ndocc*MOInfo.num_mo /*jyix*/
156
+ MOInfo.num_mo*MOInfo.num_mo /*xy_buf*/
159
if (num_i_per_ibatch > MOInfo.ndocc)
160
num_i_per_ibatch = MOInfo.ndocc;
161
if (num_i_per_ibatch < 1)
162
throw std::domain_error("Not enough memory for direct MkPT2 integrals");
163
num_ibatch = (MOInfo.ndocc + num_i_per_ibatch - 1) / num_i_per_ibatch;
164
/*--- Recompute number of MOs per I-batch ---*/
165
num_i_per_ibatch = (MOInfo.ndocc + num_ibatch - 1) / num_ibatch;
166
MkPT2_Status.num_ibatch = num_ibatch;
167
MkPT2_Status.num_i_per_ibatch = num_i_per_ibatch;
168
jsix_buf = init_array(MOInfo.ndocc*BasisSet.num_ao* num_i_per_ibatch*MOInfo.num_mo);
169
jyix_buf = init_array(MOInfo.ndocc*MOInfo.num_mo* num_i_per_ibatch*MOInfo.num_mo);
170
asij_buf = init_array(MOInfo.nuocc*BasisSet.num_ao*num_i_per_ibatch*MOInfo.ndocc);
171
abij_buf = init_array(n_ab* num_i_per_ibatch*MOInfo.ndocc);
173
iwl_buf_init(&ERIOUT,PSIF_MO_TEI,UserOptions.cutoff,0,0);
175
xy_buf = init_array(MOInfo.num_mo*MOInfo.num_mo);
177
fprintf(outfile," Using %d %s\n\n",num_ibatch, (num_ibatch == 1) ? "pass" : "passes");
179
/*--------------------------
180
Start compute threads now
181
--------------------------*/
182
thread_id = (pthread_t *) malloc(UserOptions.num_threads*sizeof(pthread_t));
183
pthread_attr_init(&thread_attr);
184
pthread_attr_setscope(&thread_attr,
185
PTHREAD_SCOPE_SYSTEM);
186
pthread_mutex_init(&mkpt2_energy_mutex,NULL);
187
pthread_cond_init(&mkpt2_energy_cond,NULL);
189
mkpt2_sindex_mutex = (pthread_mutex_t *) malloc(ioff[BasisSet.num_shells]*sizeof(pthread_mutex_t));
190
for(i=0;i<ioff[BasisSet.num_shells];i++)
192
mkpt2_sindex_mutex = (pthread_mutex_t *) malloc(BasisSet.num_ao*sizeof(pthread_mutex_t));
193
for(i=0;i<BasisSet.num_ao;i++)
195
pthread_mutex_init(&(mkpt2_sindex_mutex[i]),NULL);
196
for(long int i=0;i<UserOptions.num_threads-1;i++)
197
pthread_create(&(thread_id[i]),&thread_attr,
198
mkpt2_ints_thread,(void *)i);
199
mkpt2_ints_thread( (void *) (UserOptions.num_threads - 1) );
200
for(i=0;i<UserOptions.num_threads-1;i++)
201
pthread_join(thread_id[i], NULL);
203
pthread_mutex_destroy(&mkpt2_energy_mutex);
205
for(i=0;i<ioff[BasisSet.num_shells];i++)
207
for(i=0;i<BasisSet.num_ao;i++)
209
pthread_mutex_destroy(&(mkpt2_sindex_mutex[i]));
212
double *temp_arr = init_array(MOInfo.num_mo*MOInfo.num_mo*MOInfo.num_mo*MOInfo.num_mo);
213
iwl_buf_init(&ERIOUT,PSIF_MO_TEI,UserOptions.cutoff,1,1);
214
iwl_buf_rd_all(&ERIOUT,temp_arr,ioff,ioff,1,ioff,1,outfile);
215
iwl_buf_close(&ERIOUT, 1);
222
free(mkpt2_sindex_mutex);
224
iwl_buf_flush(&ERIOUT, 1);
225
iwl_buf_close(&ERIOUT, 1);
234
free_Taylor_Fm_Eval();