MPQC 2.3.1

- Main Page
- Related Pages
- Classes
- Files

- File List

intv3/int2e.h
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
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
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.
386 int redundant() { return redundant_; }
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.
391 int permute() { return permute_; }
392 void set_permute(int i) { permute_ = i; }
393
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
414 double *buffer() { return int_buffer; }
415
416 Ref<GaussianBasisSet> basis()
417 {
418 if (bs1_==bs2_ && bs1_ == bs3_ && bs1_ == bs4_) return bs1_;
419 return 0;
420 }
421 Ref<GaussianBasisSet> basis1() { return bs1_; }
422 Ref<GaussianBasisSet> basis2() { return bs2_; }
423 Ref<GaussianBasisSet> basis3() { return bs3_; }
424 Ref<GaussianBasisSet> basis4() { return bs4_; }
425
426 Ref<GaussianBasisSet> cs1() const { return int_cs1; }
427 Ref<GaussianBasisSet> cs2() const { return int_cs2; }
428 Ref<GaussianBasisSet> cs3() const { return int_cs3; }
429 Ref<GaussianBasisSet> cs4() const { return int_cs4; }
430
431 GaussianBasisSet * pcs1() const { return int_cs1.pointer(); }
432 GaussianBasisSet * pcs2() const { return int_cs2.pointer(); }
433 GaussianBasisSet * pcs3() const { return int_cs3.pointer(); }
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: