1 //
2 // int2e.h
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifdef __GNUG__
29 #pragma interface
30 #endif
31 
32 #ifndef _chemistry_qc_intv3_int2e_h
33 #define _chemistry_qc_intv3_int2e_h
34 
35 #include <limits.h>
36 
37 #include <util/ref/ref.h>
38 #include <chemistry/qc/basis/basis.h>
39 #include <chemistry/qc/oint3/build.h>
40 #include <chemistry/qc/intv3/fjt.h>
41 #include <chemistry/qc/intv3/types.h>
42 #include <chemistry/qc/intv3/storage.h>
43 #include <chemistry/qc/intv3/array.h>
44 #include <chemistry/qc/intv3/macros.h>
45 
46 namespace sc {
47 
48 class Integral;
49 
50 #define CHECK_INTEGRAL_ALGORITHM 0
51 
52 /** Int2eV3 is a class wrapper for the two body part of the C language
53     IntV3 library.  It is used by TwoBodyIntV3 and TwoBodyDerivIntV3 to
54     implement IntegralV3. */
55 class Int2eV3: public RefCount {
56   protected:
57     Integral *integral_;
58 
59     BuildIntV3 build;
60     Ref<IntegralStorer> storer;
61 
62     Ref<GaussianBasisSet> bs1_;
63     Ref<GaussianBasisSet> bs2_;
64     Ref<GaussianBasisSet> bs3_;
65     Ref<GaussianBasisSet> bs4_;
66 
67     // the permuted bases
68     Ref<GaussianBasisSet> pbs1_;
69     Ref<GaussianBasisSet> pbs2_;
70     Ref<GaussianBasisSet> pbs3_;
71     Ref<GaussianBasisSet> pbs4_;
72 
73     Ref<MessageGrp> grp_;
74 
75     int bs1_shell_offset_;
76     int bs2_shell_offset_;
77     int bs3_shell_offset_;
78     int bs4_shell_offset_;
79     int bs1_func_offset_;
80     int bs2_func_offset_;
81     int bs3_func_offset_;
82     int bs4_func_offset_;
83     int bs1_prim_offset_;
84     int bs2_prim_offset_;
85     int bs3_prim_offset_;
86     int bs4_prim_offset_;
87 
88     // statics from vrr.cc
89   public:
90     enum { STORAGE_CHUNK = 4096 };
91   protected:
92     struct store_list {
93         void* data[STORAGE_CHUNK];
94         struct store_list* p;
95     };
96     typedef struct store_list store_list_t;
97     int n_store_last;
98     store_list_t* store;
99     typedef int (BuildIntV3::*intfunc)();
100     intfunc build_routine[4][4][4][4][2];
101     /* Offset shell numbers. */
102     int osh1, osh2, osh3, osh4;
103     /* Offset primitive numbers. */
104     int opr1, opr2, opr3, opr4;
105     /* Saved initialization parameters used to free data. */
106     int saved_am12,saved_am34,saved_ncon;
107     /* Stores the length of the inner loop for integral contraction. */
108     IntV3Arrayint3 contract_length;
109 
110     // statics from hrr.cc
111   protected:
112     /* The general contraction numbers. */
113     int g1,g2,g3,g4;
114     /* A[] - B[] */
115     double AmB[3];
116     /* C[] - D[] */
117     double CmD[3];
118     int eAB;
119     double *buf34;
120     double *buf12;
121     double *bufshared;
122 
123     int redundant_;
124     int permute_;
125 
126   protected:
127     Ref<FJT> fjt_;
128 
129     int *int_shell_to_prim;
130     IntV3Arraydouble2 int_shell_r;
131     IntV3Arraydouble2 int_prim_zeta;
132     IntV3Arraydouble2 int_prim_k;
133     IntV3Arraydouble2 int_prim_oo2zeta;
134     IntV3Arraydouble3 int_prim_p;
135 
136     double *int_buffer;
137     double *int_derint_buffer;
138 
139     Ref<GaussianBasisSet> int_cs1;
140     Ref<GaussianBasisSet> int_cs2;
141     Ref<GaussianBasisSet> int_cs3;
142     Ref<GaussianBasisSet> int_cs4;
143 
144     GaussianShell *int_shell1;
145     GaussianShell *int_shell2;
146     GaussianShell *int_shell3;
147     GaussianShell *int_shell4;
148 
149     IntV3Arraydoublep2 ****e0f0_con_ints_array;  /* The contr. int. inter. */
150 
151     int int_expweight1; // For exponent weighted contractions.
152     int int_expweight2; // For exponent weighted contractions.
153     int int_expweight3; // For exponent weighted contractions.
154     int int_expweight4; // For exponent weighted contractions.
155 
156     // These are used to compute two and three center electron repulsion
157     // integrals.  int_unit2 is 1 if shell 2 is to have value one everywhere
158     // and int_unit4 is 1 if shell4 is to be a unit function.  Otherwise,
159     // they should be zero.
160     //
161 
162     int int_unit2;
163     int int_unit4;
164     GaussianShell* int_unit_shell;
165 
166     int int_integral_storage;
167     int int_store1;
168     int int_store2;
169     int int_derivative_bounds;
170 
171     // locals from vrr.cc
172   protected:
173     void add_store(void *p);
174     void free_store();
175     void _free_store(store_list_t* s, int n);
176     void build_not_using_gcs(int nc1, int nc2, int nc3, int nc4,
177                              int minam1, int minam3, int maxam12, int maxam34,
178                              int dam1, int dam2, int dam3, int dam4, int eAB);
179     void build_using_gcs(int nc1, int nc2, int nc3, int nc4,
180                          int minam1, int minam3, int maxam12, int maxam34,
181                          int dam1, int dam2, int dam3, int dam4, int eAB);
182     void gen_prim_intermediates(int pr1, int pr2, int pr3, int pr4, int am);
183     void gen_prim_intermediates_with_norm(int pr1, int pr2, int pr3, int pr4,
184                                  int am, double norm);
185     void gen_shell_intermediates(int sh1, int sh2, int sh3, int sh4);
186     void blockbuildprim(int minam1, int maxam12, int minam3, int maxam34);
187     void blockbuildprim_1(int am12min, int am12max, int am34, int m);
188     void blockbuildprim_3(int am34min, int am34max, int m);
189 
190     // globals from vrr.cc
191   protected:
192     void int_init_buildgc(int order,
193                           int am1, int am2, int am3, int am4,
194                           int nc1, int nc2, int nc3, int nc4);
195     void int_done_buildgc();
196     void int_buildgcam(int minam1, int minam2, int minam3, int minam4,
197                        int maxam1, int maxam2, int maxam3, int maxam4,
198                        int dam1, int dam2, int dam3, int dam4,
199                        int sh1, int sh2, int sh3, int sh4,
200                        int eAB);
201 
202     // globals from print2e.cc
203   protected:
204     void int_offset_print(std::ostream &,
205                           double *buffer,
206                           Ref<GaussianBasisSet> c1, int s1,
207                           Ref<GaussianBasisSet> c2, int s2,
208                           Ref<GaussianBasisSet> c3, int s3,
209                           Ref<GaussianBasisSet> c4, int s4);
210     void int_offset_print_n(std::ostream &, double *buffer,
211                             int n1, int n2, int n3, int n4,
212                             int o1, int o2, int o3, int o4,
213                             int e12, int e13e24, int e34);
214     void int_print(std::ostream &, double *buffer,
215                    Ref<GaussianBasisSet> c1, int s1,
216                    Ref<GaussianBasisSet> c2, int s2,
217                    Ref<GaussianBasisSet> c3, int s3,
218                    Ref<GaussianBasisSet> c4, int s4);
219     void int_print_n(std::ostream &, double *buffer,
220                      int n1, int n2, int n3, int n4,
221                      int e12, int e13e24, int e34);
222     void int_print_intermediates(std::ostream &);
223 
224     // locals from hrr.cc
225   protected:
226     void shiftam_12(double *I0100, double *I1000, double *I0000,
227                     int am1, int am2, int am3, int am4);
228     void shiftam_12eAB(double *I0100, double *I1000, double *I0000,
229                        int am1, int am2, int am3, int am4);
230     void shiftam_34(double *I0001, double *I0010, double *I0000,
231                     int am1, int am2, int am3, int am4);
232 
233     // globals from hrr.cc
234   protected:
235     void int_init_shiftgc(int order, int am1, int am2, int am3, int am4);
236     void int_done_shiftgc();
237     double *int_shiftgcam(int gc1, int gc2, int gc3, int gc4,
238                           int tam1, int tam2, int tam3, int tam4, int peAB);
239 
240     // locals from init2e.cc
241   protected:
242     void alloc_inter(int nprim,int nshell);
243     void compute_shell_1(Ref<GaussianBasisSet> cs, int, int);
244     void compute_prim_2(Ref<GaussianBasisSet> cs1,int,int,
245                         Ref<GaussianBasisSet> cs2,int,int);
246 
247 
248     // globals from init2e.cc
249   protected:
250     double *int_initialize_erep(size_t storage, int order,
251                                 const Ref<GaussianBasisSet> &cs1,
252                                 const Ref<GaussianBasisSet> &cs2,
253                                 const Ref<GaussianBasisSet> &cs3,
254                                 const Ref<GaussianBasisSet> &cs4);
255     void int_done_erep();
256 
257     // from tformv3.cc
258   protected:
259     double *source;
260     double *target;
261     double *scratch;
262     int nsourcemax;
263     // transform implementation functions:
264     void transform_init();
265     void transform_done();
266     void source_space(int nsource);
267     void copy_to_source(double *integrals, int nsource);
268     void do_gencon_sparse_transform_2e(Integral*integ,
269                                        double *integrals, double *target,
270                                        int index,
271                                        GaussianShell *sh1, GaussianShell *sh2,
272                                        GaussianShell *sh3, GaussianShell *sh4);
273     // functions for general use outside of tformv3.cc:
274     // integrals and target may overlap
275     void transform_2e_slow(Integral *,
276                       double *integrals, double *target,
277                       GaussianShell *sh1, GaussianShell *sh2,
278                       GaussianShell *sh3, GaussianShell *sh4);
279     void transform_2e(Integral *,
280                       double *integrals, double *target,
281                       GaussianShell *sh1, GaussianShell *sh2,
282                       GaussianShell *sh3, GaussianShell *sh4);
283 
284     // locals from comp2e.cc
285   protected:
286     void compute_erep(int flags, int *psh1, int *psh2, int *psh3, int *psh4,
287                       int dam1, int dam2, int dam3, int dam4);
288     void compute_erep_1der(int flags, double *buffer,
289                            int *psh1, int *psh2, int *psh3, int *psh4,
290                            int dercenter);
291     void nonredundant_erep(double *buffer, int e12, int e34, int e13e24,
292                            int n1, int n2, int n3, int n4,
293                            int *red_off, int *nonred_off);
294     void compute_erep_bound1der(int flags, double *buffer,
295                                 int *psh1, int *psh2, int *psh3, int *psh4);
296 
297     // globals from comp2e.cc
298   protected:
299     void int_erep_bound1der(int flags, int bsh1, int bsh2, int *size);
300 
301 
302     // global vars from bounds.h
303   protected:
304     typedef signed char int_bound_t;
305     enum { int_bound_min = SCHAR_MIN, int_bound_max = SCHAR_MAX };
306     int_bound_t int_Q;
307     int_bound_t int_R;
308     int_bound_t *int_Qvec;
309     int_bound_t *int_Rvec;
310 
311     // global routines from bounds.cc
312   protected:
313     void int_init_bounds_nocomp();
314     void int_init_bounds_1der_nocomp();
315     void int_bounds_comp(int s1, int s2);
316     void int_bounds_1der_comp(int s1, int s2);
317     int int_erep_2bound(int s1, int s2);
318     int int_erep_0bound_1der();
319     int int_erep_2bound_1der(int s1, int s2);
320 
321     // local routines from bounds.cc
322   protected:
323     void compute_bounds(int_bound_t *overall, int_bound_t *vec, int flag);
324     void compute_bounds_shell(int_bound_t *overall, int_bound_t *vec,
325                               int flag, int sh1, int sh2);
326 
327     // global routines from storage.cc
328   protected:
329     int int_have_stored_integral(int sh1,int sh2,int sh3,int sh4,
330                                  int p12,int p34,int p13p24);
331     void int_store_integral(int sh1,int sh2,int sh3,int sh4,
332                             int p12,int p34,int p13p24,
333                             int size);
334 
335     // from offsets.cc
336   protected:
337     void int_initialize_offsets2();
338     void int_done_offsets2();
339 
340     // from comp2e3c.cc
341   protected:
342     void make_int_unit_shell();
343     void delete_int_unit_shell();
344 
345   protected:
346     // for intermediate storage:
347     int used_storage_;
348     int used_storage_build_;
349     int used_storage_shift_;
350 
351   public:
352     // bs4 must be null for 3 center integrals
353     // bs2 must be null for 2 center integrals
354     // bs1 and bs3 must be nonnull.
355     Int2eV3(Integral *,
356             const Ref<GaussianBasisSet>&bs1,
357             const Ref<GaussianBasisSet>&bs2,
358             const Ref<GaussianBasisSet>&bs3,
359             const Ref<GaussianBasisSet>&bs4,
360             int order, size_t storage);
361     ~Int2eV3();
362 
363     // storage.cc: for the storage of integrals
364     void init_storage(int size);
365     void done_storage();
366 
367     // for intermediate storage
storage_used()368     int storage_used() { return used_storage_; }
369 
370     // bounds.cc
371     void init_bounds();
372     void init_bounds_1der();
373     void done_bounds();
374     void done_bounds_1der();
375     // Covert a bound to/from the log of the bound (returns 2^bound)
376     // replace:
377     //double int_bound_to_double(int bound);
378     //double int_bound_double(int value);
379     //int int_bound_log(double value);
380     static double logbound_to_bound(int);
381     static int bound_to_logbound(double value);
382 
383     // If redundant is false the redundant integrals that arise
384     // when a shell index is repeated are stored.
385     // The default is true.
redundant()386     int redundant() { return redundant_; }
set_redundant(int i)387     void set_redundant(int i) { redundant_ = i; }
388 
389     // If permute is true the routines are allowed to permute indices.
390     // The default is false.
permute()391     int permute() { return permute_; }
set_permute(int i)392     void set_permute(int i) { permute_ = i; }
393 
used_storage()394     int used_storage() const { return used_storage_; }
395 
396     // from comp2e.cc
397     void erep(int &psh1, int &psh2, int &psh3, int &psh4);
398     void erep(int *shells, int  *sizes);
399     void erep_all1der(int &psh1, int &psh2, int &psh3, int &psh4,
400                       der_centersv3_t *der_centers);
401     void erep_all1der(int *shells, int  *sizes,
402                       der_centersv3_t *dercenters);
403 
404     // from comp2e3c.cc
405     void erep_2center(int &psh1, int &psh2);
406     void erep_2center(int *shells, int  *sizes);
407     void erep_3center(int &psh1, int &psh2, int &psh3);
408     void erep_3center(int *shells, int  *sizes);
409 
410     // from bounds.cc
411     int erep_4bound(int s1, int s2, int s3, int s4);
412     int erep_4bound_1der(int s1, int s2, int s3, int s4);
413 
buffer()414     double *buffer() { return int_buffer; }
415 
basis()416     Ref<GaussianBasisSet> basis()
417     {
418       if (bs1_==bs2_ && bs1_ == bs3_ && bs1_ == bs4_) return bs1_;
419       return 0;
420     }
basis1()421     Ref<GaussianBasisSet> basis1() { return bs1_; }
basis2()422     Ref<GaussianBasisSet> basis2() { return bs2_; }
basis3()423     Ref<GaussianBasisSet> basis3() { return bs3_; }
basis4()424     Ref<GaussianBasisSet> basis4() { return bs4_; }
425 
cs1()426     Ref<GaussianBasisSet> cs1() const { return int_cs1; }
cs2()427     Ref<GaussianBasisSet> cs2() const { return int_cs2; }
cs3()428     Ref<GaussianBasisSet> cs3() const { return int_cs3; }
cs4()429     Ref<GaussianBasisSet> cs4() const { return int_cs4; }
430 
pcs1()431     GaussianBasisSet * pcs1() const { return int_cs1.pointer(); }
pcs2()432     GaussianBasisSet * pcs2() const { return int_cs2.pointer(); }
pcs3()433     GaussianBasisSet * pcs3() const { return int_cs3.pointer(); }
pcs4()434     GaussianBasisSet * pcs4() const { return int_cs4.pointer(); }
435 };
436 
437 }
438 
439 #endif
440 
441 // Local Variables:
442 // mode: c++
443 // c-file-style: "CLJ"
444 // End:
445