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