00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef __GNUG__
00029 #pragma interface
00030 #endif
00031
00032 #ifndef _chemistry_qc_intv3_int2e_h
00033 #define _chemistry_qc_intv3_int2e_h
00034
00035 #include <limits.h>
00036
00037 #include <util/ref/ref.h>
00038 #include <chemistry/qc/basis/basis.h>
00039 #include <chemistry/qc/oint3/build.h>
00040 #include <chemistry/qc/intv3/fjt.h>
00041 #include <chemistry/qc/intv3/types.h>
00042 #include <chemistry/qc/intv3/storage.h>
00043 #include <chemistry/qc/intv3/array.h>
00044 #include <chemistry/qc/intv3/macros.h>
00045
00046 namespace sc {
00047
00048 class Integral;
00049
00050 #define CHECK_INTEGRAL_ALGORITHM 0
00051
00055 class Int2eV3: public RefCount {
00056 protected:
00057 Integral *integral_;
00058
00059 BuildIntV3 build;
00060 Ref<IntegralStorer> storer;
00061
00062 Ref<GaussianBasisSet> bs1_;
00063 Ref<GaussianBasisSet> bs2_;
00064 Ref<GaussianBasisSet> bs3_;
00065 Ref<GaussianBasisSet> bs4_;
00066
00067 Ref<MessageGrp> grp_;
00068
00069 int bs1_shell_offset_;
00070 int bs2_shell_offset_;
00071 int bs3_shell_offset_;
00072 int bs4_shell_offset_;
00073 int bs1_func_offset_;
00074 int bs2_func_offset_;
00075 int bs3_func_offset_;
00076 int bs4_func_offset_;
00077 int bs1_prim_offset_;
00078 int bs2_prim_offset_;
00079 int bs3_prim_offset_;
00080 int bs4_prim_offset_;
00081
00082
00083 public:
00084 enum { STORAGE_CHUNK = 4096 };
00085 protected:
00086 struct store_list {
00087 void* data[STORAGE_CHUNK];
00088 struct store_list* p;
00089 };
00090 typedef struct store_list store_list_t;
00091 int n_store_last;
00092 store_list_t* store;
00093 typedef int (BuildIntV3::*intfunc)();
00094 intfunc build_routine[4][4][4][4][2];
00095
00096 int osh1, osh2, osh3, osh4;
00097
00098 int opr1, opr2, opr3, opr4;
00099
00100 int saved_am12,saved_am34,saved_ncon;
00101
00102 IntV3Arrayint3 contract_length;
00103
00104
00105 protected:
00106
00107 int g1,g2,g3,g4;
00108
00109 double AmB[3];
00110
00111 double CmD[3];
00112 int eAB;
00113 double *buf34;
00114 double *buf12;
00115 double *bufshared;
00116
00117 int redundant_;
00118 int permute_;
00119
00120 protected:
00121 Ref<FJT> fjt_;
00122
00123 int *int_shell_to_prim;
00124 IntV3Arraydouble2 int_shell_r;
00125 IntV3Arraydouble2 int_prim_zeta;
00126 IntV3Arraydouble2 int_prim_k;
00127 IntV3Arraydouble2 int_prim_oo2zeta;
00128 IntV3Arraydouble3 int_prim_p;
00129
00130 double *int_buffer;
00131 double *int_derint_buffer;
00132
00133 Ref<GaussianBasisSet> int_cs1;
00134 Ref<GaussianBasisSet> int_cs2;
00135 Ref<GaussianBasisSet> int_cs3;
00136 Ref<GaussianBasisSet> int_cs4;
00137
00138 GaussianShell *int_shell1;
00139 GaussianShell *int_shell2;
00140 GaussianShell *int_shell3;
00141 GaussianShell *int_shell4;
00142
00143 IntV3Arraydoublep2 ****e0f0_con_ints_array;
00144
00145 int int_expweight1;
00146 int int_expweight2;
00147 int int_expweight3;
00148 int int_expweight4;
00149
00150
00151
00152
00153
00154
00155
00156 int int_unit2;
00157 int int_unit4;
00158 GaussianShell* int_unit_shell;
00159
00160 int int_integral_storage;
00161 int int_store1;
00162 int int_store2;
00163 int int_derivative_bounds;
00164
00165
00166 protected:
00167 void add_store(void *p);
00168 void free_store();
00169 void _free_store(store_list_t* s, int n);
00170 void build_not_using_gcs(int nc1, int nc2, int nc3, int nc4,
00171 int minam1, int minam3, int maxam12, int maxam34,
00172 int dam1, int dam2, int dam3, int dam4, int eAB);
00173 void build_using_gcs(int nc1, int nc2, int nc3, int nc4,
00174 int minam1, int minam3, int maxam12, int maxam34,
00175 int dam1, int dam2, int dam3, int dam4, int eAB);
00176 void gen_prim_intermediates(int pr1, int pr2, int pr3, int pr4, int am);
00177 void gen_prim_intermediates_with_norm(int pr1, int pr2, int pr3, int pr4,
00178 int am, double norm);
00179 void gen_shell_intermediates(int sh1, int sh2, int sh3, int sh4);
00180 void blockbuildprim(int minam1, int maxam12, int minam3, int maxam34);
00181 void blockbuildprim_1(int am12min, int am12max, int am34, int m);
00182 void blockbuildprim_3(int am34min, int am34max, int m);
00183
00184
00185 protected:
00186 void int_init_buildgc(int order,
00187 int am1, int am2, int am3, int am4,
00188 int nc1, int nc2, int nc3, int nc4);
00189 void int_done_buildgc();
00190 void int_buildgcam(int minam1, int minam2, int minam3, int minam4,
00191 int maxam1, int maxam2, int maxam3, int maxam4,
00192 int dam1, int dam2, int dam3, int dam4,
00193 int sh1, int sh2, int sh3, int sh4,
00194 int eAB);
00195
00196
00197 protected:
00198 void int_offset_print(std::ostream &,
00199 double *buffer,
00200 Ref<GaussianBasisSet> c1, int s1,
00201 Ref<GaussianBasisSet> c2, int s2,
00202 Ref<GaussianBasisSet> c3, int s3,
00203 Ref<GaussianBasisSet> c4, int s4);
00204 void int_offset_print_n(std::ostream &, double *buffer,
00205 int n1, int n2, int n3, int n4,
00206 int o1, int o2, int o3, int o4,
00207 int e12, int e13e24, int e34);
00208 void int_print(std::ostream &, double *buffer,
00209 Ref<GaussianBasisSet> c1, int s1,
00210 Ref<GaussianBasisSet> c2, int s2,
00211 Ref<GaussianBasisSet> c3, int s3,
00212 Ref<GaussianBasisSet> c4, int s4);
00213 void int_print_n(std::ostream &, double *buffer,
00214 int n1, int n2, int n3, int n4,
00215 int e12, int e13e24, int e34);
00216 void int_print_intermediates(std::ostream &);
00217
00218
00219 protected:
00220 void shiftam_12(double *I0100, double *I1000, double *I0000,
00221 int am1, int am2, int am3, int am4);
00222 void shiftam_12eAB(double *I0100, double *I1000, double *I0000,
00223 int am1, int am2, int am3, int am4);
00224 void shiftam_34(double *I0001, double *I0010, double *I0000,
00225 int am1, int am2, int am3, int am4);
00226
00227
00228 protected:
00229 void int_init_shiftgc(int order, int am1, int am2, int am3, int am4);
00230 void int_done_shiftgc();
00231 double *int_shiftgcam(int gc1, int gc2, int gc3, int gc4,
00232 int tam1, int tam2, int tam3, int tam4, int peAB);
00233
00234
00235 protected:
00236 void alloc_inter(int nprim,int nshell);
00237 void compute_shell_1(Ref<GaussianBasisSet> cs, int, int);
00238 void compute_prim_1(Ref<GaussianBasisSet> cs1);
00239 void compute_shell_2(Ref<GaussianBasisSet> cs1,Ref<GaussianBasisSet> cs2);
00240 void compute_prim_2(Ref<GaussianBasisSet> cs1,Ref<GaussianBasisSet> cs2);
00241
00242
00243
00244 protected:
00245 double *int_initialize_erep(size_t storage, int order,
00246 const Ref<GaussianBasisSet> &cs1,
00247 const Ref<GaussianBasisSet> &cs2,
00248 const Ref<GaussianBasisSet> &cs3,
00249 const Ref<GaussianBasisSet> &cs4);
00250 void int_done_erep();
00251
00252
00253 protected:
00254 double *source;
00255 double *target;
00256 double *scratch;
00257 int nsourcemax;
00258
00259 void transform_init();
00260 void transform_done();
00261 void source_space(int nsource);
00262 void copy_to_source(double *integrals, int nsource);
00263 void do_gencon_sparse_transform_2e(Integral*integ,
00264 double *integrals, double *target,
00265 int index,
00266 GaussianShell *sh1, GaussianShell *sh2,
00267 GaussianShell *sh3, GaussianShell *sh4);
00268
00269
00270 void transform_2e_slow(Integral *,
00271 double *integrals, double *target,
00272 GaussianShell *sh1, GaussianShell *sh2,
00273 GaussianShell *sh3, GaussianShell *sh4);
00274 void transform_2e(Integral *,
00275 double *integrals, double *target,
00276 GaussianShell *sh1, GaussianShell *sh2,
00277 GaussianShell *sh3, GaussianShell *sh4);
00278
00279
00280 protected:
00281 void compute_erep(int flags, int *psh1, int *psh2, int *psh3, int *psh4,
00282 int dam1, int dam2, int dam3, int dam4);
00283 void compute_erep_1der(int flags, double *buffer,
00284 int *psh1, int *psh2, int *psh3, int *psh4,
00285 int dercenter);
00286 void nonredundant_erep(double *buffer, int e12, int e34, int e13e24,
00287 int n1, int n2, int n3, int n4,
00288 int *red_off, int *nonred_off);
00289 void compute_erep_bound1der(int flags, double *buffer,
00290 int *psh1, int *psh2, int *psh3, int *psh4);
00291
00292
00293 protected:
00294 void int_erep_bound1der(int flags, int bsh1, int bsh2, int *size);
00295
00296
00297
00298 protected:
00299 typedef signed char int_bound_t;
00300 enum { int_bound_min = SCHAR_MIN, int_bound_max = SCHAR_MAX };
00301 int_bound_t int_Q;
00302 int_bound_t int_R;
00303 int_bound_t *int_Qvec;
00304 int_bound_t *int_Rvec;
00305
00306
00307 protected:
00308 void int_init_bounds_nocomp();
00309 void int_init_bounds_1der_nocomp();
00310 void int_bounds_comp(int s1, int s2);
00311 void int_bounds_1der_comp(int s1, int s2);
00312 int int_erep_2bound(int s1, int s2);
00313 int int_erep_0bound_1der();
00314 int int_erep_2bound_1der(int s1, int s2);
00315
00316
00317 protected:
00318 void compute_bounds(int_bound_t *overall, int_bound_t *vec, int flag);
00319 void compute_bounds_shell(int_bound_t *overall, int_bound_t *vec,
00320 int flag, int sh1, int sh2);
00321
00322
00323 protected:
00324 int int_have_stored_integral(int sh1,int sh2,int sh3,int sh4,
00325 int p12,int p34,int p13p24);
00326 void int_store_integral(int sh1,int sh2,int sh3,int sh4,
00327 int p12,int p34,int p13p24,
00328 int size);
00329
00330
00331 protected:
00332 void int_initialize_offsets2();
00333 void int_done_offsets2();
00334
00335
00336 protected:
00337 void make_int_unit_shell();
00338 void delete_int_unit_shell();
00339
00340 protected:
00341
00342 int used_storage_;
00343 int used_storage_build_;
00344 int used_storage_shift_;
00345
00346 public:
00347 Int2eV3(Integral *,
00348 const Ref<GaussianBasisSet>&,
00349 const Ref<GaussianBasisSet>&,
00350 const Ref<GaussianBasisSet>&,
00351 const Ref<GaussianBasisSet>&,
00352 int order, size_t storage);
00353 ~Int2eV3();
00354
00355
00356 void init_storage(int size);
00357 void done_storage();
00358
00359
00360 int storage_used() { return used_storage_; }
00361
00362
00363 void init_bounds();
00364 void init_bounds_1der();
00365 void done_bounds();
00366 void done_bounds_1der();
00367
00368
00369
00370
00371
00372 static double logbound_to_bound(int);
00373 static int bound_to_logbound(double value);
00374
00375
00376
00377
00378 int redundant() { return redundant_; }
00379 void set_redundant(int i) { redundant_ = i; }
00380
00381
00382
00383 int permute() { return permute_; }
00384 void set_permute(int i) { permute_ = i; }
00385
00386 int used_storage() const { return used_storage_; }
00387
00388
00389 void erep(int &psh1, int &psh2, int &psh3, int &psh4);
00390 void erep(int *shells, int *sizes);
00391 void erep_all1der(int &psh1, int &psh2, int &psh3, int &psh4,
00392 der_centersv3_t *der_centers);
00393 void erep_all1der(int *shells, int *sizes,
00394 der_centersv3_t *dercenters);
00395
00396
00397 void erep_2center(int &psh1, int &psh2);
00398 void erep_2center(int *shells, int *sizes);
00399 void erep_3center(int &psh1, int &psh2, int &psh3);
00400 void erep_3center(int *shells, int *sizes);
00401
00402
00403 int erep_4bound(int s1, int s2, int s3, int s4);
00404 int erep_4bound_1der(int s1, int s2, int s3, int s4);
00405
00406 double *buffer() { return int_buffer; }
00407
00408 Ref<GaussianBasisSet> basis()
00409 {
00410 if (bs1_==bs2_ && bs1_ == bs3_ && bs1_ == bs4_) return bs1_;
00411 return 0;
00412 }
00413 Ref<GaussianBasisSet> basis1() { return bs1_; }
00414 Ref<GaussianBasisSet> basis2() { return bs2_; }
00415 Ref<GaussianBasisSet> basis3() { return bs3_; }
00416 Ref<GaussianBasisSet> basis4() { return bs4_; }
00417
00418 Ref<GaussianBasisSet> cs1() const { return int_cs1; }
00419 Ref<GaussianBasisSet> cs2() const { return int_cs2; }
00420 Ref<GaussianBasisSet> cs3() const { return int_cs3; }
00421 Ref<GaussianBasisSet> cs4() const { return int_cs4; }
00422
00423 GaussianBasisSet * pcs1() const { return int_cs1.pointer(); }
00424 GaussianBasisSet * pcs2() const { return int_cs2.pointer(); }
00425 GaussianBasisSet * pcs3() const { return int_cs3.pointer(); }
00426 GaussianBasisSet * pcs4() const { return int_cs4.pointer(); }
00427 };
00428
00429 }
00430
00431 #endif
00432
00433
00434
00435
00436