← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/MkPT2/mkpt2_ints.cc

  • 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/cints/MkPT2/mkpt2_ints.cc
 
1
/*! \file
 
2
    \ingroup CINTS
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include<cmath>
 
6
#include<cstdio>
 
7
#include<cstring>
 
8
#include<cstdlib>
 
9
#include<memory.h>
 
10
#include<pthread.h>
 
11
#include<libciomr/libciomr.h>
 
12
#include<libchkpt/chkpt.h>
 
13
#include<libqt/qt.h>
 
14
#include<libiwl/iwl.h>
 
15
#include<libint/libint.h>
 
16
 
 
17
#include"defines.h"
 
18
#include"psifiles.h"
 
19
#define EXTERN
 
20
#include"global.h"
 
21
#include <stdexcept>
 
22
#include"schwartz.h"
 
23
#include"quartet_data.h"
 
24
#include"norm_quartet.h"
 
25
#include "data_structs.h"
 
26
#ifdef USE_TAYLOR_FM
 
27
  #include"taylor_fm_eval.h"
 
28
#else
 
29
  #include"int_fjt.h"
 
30
  #include"fjt.h"
 
31
#endif
 
32
#include"quartet_permutations.h"
 
33
#include"mkpt2_ints.h"
 
34
 
 
35
 
 
36
namespace psi { 
 
37
  namespace CINTS {
 
38
    namespace mkpt2 {
 
39
    
 
40
    /*-------------------------------
 
41
      Explicit function declarations
 
42
      -------------------------------*/
 
43
    extern void *mkpt2_ints_thread(void *);
 
44
    
 
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;
 
52
    double *jsix_buf;
 
53
    double *jyix_buf;
 
54
    double *asij_buf;
 
55
    double *abij_buf;
 
56
    unsigned long int n_ij;
 
57
    unsigned long int n_xy;
 
58
    unsigned long int n_ab;
 
59
#if MkPT2_USE_IWL
 
60
    iwlbuf ERIOUT;
 
61
#else
 
62
    double * xy_buf;
 
63
#endif
 
64
    
 
65
    
 
66
    void mkpt2_ints()
 
67
    {
 
68
      pthread_attr_t thread_attr;
 
69
      pthread_t *thread_id;
 
70
  
 
71
      /*--- Various data structures ---*/
 
72
      Libint_t Libint;
 
73
      long int libint_memory;
 
74
      int max_bf_per_shell;
 
75
      int max_num_prim_comb;
 
76
      
 
77
      int i;
 
78
      
 
79
      
 
80
      int num_ibatch, num_i_per_ibatch, ibatch, ibatch_length;
 
81
      int imin, imax, jmin;
 
82
      int mo_i, mo_j, mo_a, mo_b, mo_ij;
 
83
      int rs_offset, rsi_offset, rsp_offset;
 
84
      
 
85
      double AB2, CD2;
 
86
      
 
87
      double *raw_data;             /* pointer to the unnormalized taregt quartet of integrals */
 
88
      double *data;                 /* pointer to the transformed normalized target quartet of integrals */
 
89
#ifdef NONDOUBLE_INTS
 
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 */
 
92
#endif
 
93
      
 
94
      double *rspq_ptr;
 
95
      double temp;
 
96
      double *mo_vec;
 
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 */
 
99
      double *rsi_row;
 
100
      double **ix_buf;              /* buffer for one |ix) ket */
 
101
      double *i_row;
 
102
      double *jsi_row;
 
103
      double *jbi_row;
 
104
      double ixjy, iyjx, pfac;
 
105
      
 
106
      double temp1,temp2,*iq_row,*ip_row;
 
107
      int rs,qrs;
 
108
      
 
109
      /*---------------
 
110
        Initixlization
 
111
        ---------------*/
 
112
#ifdef USE_TAYLOR_FM
 
113
      init_Taylor_Fm_Eval(BasisSet.max_am*4-4,UserOptions.cutoff);
 
114
#else
 
115
      init_fjt(BasisSet.max_am*4);
 
116
#endif
 
117
      init_libint_base();
 
118
      timer_init();
 
119
      timer_on("Schwartz");
 
120
      schwartz_eri();
 
121
      timer_off("Schwartz");
 
122
      MkPT2_Status.num_arrived = 0;
 
123
      
 
124
      
 
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;
 
139
      /*---
 
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
 
143
      ---*/
 
144
      fprintf(outfile,"\n  Computing MkPT2 integrals\n");
 
145
      num_i_per_ibatch = UserOptions.memory / (UserOptions.num_threads*
 
146
                              (
 
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*/
 
150
                              ) +
 
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*/
 
155
#if !MkPT2_USE_IWL
 
156
                             + MOInfo.num_mo*MOInfo.num_mo /*xy_buf*/
 
157
#endif
 
158
                             );
 
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);
 
172
#if MkPT2_USE_IWL
 
173
      iwl_buf_init(&ERIOUT,PSIF_MO_TEI,UserOptions.cutoff,0,0);
 
174
#else
 
175
      xy_buf = init_array(MOInfo.num_mo*MOInfo.num_mo);
 
176
#endif
 
177
      fprintf(outfile,"  Using %d %s\n\n",num_ibatch, (num_ibatch == 1) ? "pass" : "passes");
 
178
      
 
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);
 
188
#if LOCK_RS_SHELL
 
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++)
 
191
#else
 
192
        mkpt2_sindex_mutex = (pthread_mutex_t *) malloc(BasisSet.num_ao*sizeof(pthread_mutex_t));
 
193
      for(i=0;i<BasisSet.num_ao;i++)
 
194
#endif
 
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);
 
202
    free(thread_id);
 
203
    pthread_mutex_destroy(&mkpt2_energy_mutex);
 
204
#if LOCK_RS_SHELL
 
205
      for(i=0;i<ioff[BasisSet.num_shells];i++)
 
206
#else
 
207
      for(i=0;i<BasisSet.num_ao;i++)
 
208
#endif
 
209
        pthread_mutex_destroy(&(mkpt2_sindex_mutex[i]));
 
210
      
 
211
#if MkPT2_TEST
 
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);  
 
216
     free(temp_arr);
 
217
#endif
 
218
 
 
219
      /*---------
 
220
        Clean-up
 
221
        ---------*/
 
222
      free(mkpt2_sindex_mutex);
 
223
#if MkPT2_USE_IWL
 
224
  iwl_buf_flush(&ERIOUT, 1);
 
225
  iwl_buf_close(&ERIOUT, 1);  
 
226
#else
 
227
      free(xy_buf);
 
228
#endif
 
229
      free(jyix_buf);
 
230
      free(abij_buf);
 
231
      free(jsix_buf);
 
232
      free(asij_buf);
 
233
#ifdef USE_TAYLOR_FM
 
234
      free_Taylor_Fm_Eval();
 
235
#else
 
236
      free_fjt();
 
237
#endif
 
238
      
 
239
      timer_done();
 
240
      
 
241
      return;
 
242
    }
 
243
  }
 
244
}
 
245
}
 
246