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 #ifndef _chemistry_qc_basis_petite_h
00029 #define _chemistry_qc_basis_petite_h
00030
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034
00035 #include <iostream>
00036
00037 #include <util/misc/scint.h>
00038 #include <util/ref/ref.h>
00039 #include <util/container/array.h>
00040 #include <math/scmat/blocked.h>
00041 #include <math/scmat/offset.h>
00042 #include <chemistry/molecule/molecule.h>
00043 #include <chemistry/qc/basis/gaussbas.h>
00044 #include <chemistry/qc/basis/integral.h>
00045
00046
00047
00048 namespace sc {
00049
00050 inline sc_int_least64_t
00051 ij_offset64(sc_int_least64_t i, sc_int_least64_t j)
00052 {
00053 return (i>j) ? (((i*(i+1)) >> 1) + j) : (((j*(j+1)) >> 1) + i);
00054 }
00055
00056 inline sc_int_least64_t
00057 i_offset64(sc_int_least64_t i)
00058 {
00059 return ((i*(i+1)) >> 1);
00060 }
00061
00062
00063
00064 struct contribution {
00065 int bfn;
00066 double coef;
00067
00068 contribution();
00069 contribution(int b, double c);
00070 ~contribution();
00071 };
00072
00073 struct SO {
00074 int len;
00075 int length;
00076 contribution *cont;
00077
00078 SO();
00079 SO(int);
00080 ~SO();
00081
00082 SO& operator=(const SO&);
00083
00084 void set_length(int);
00085 void reset_length(int);
00086
00087
00088 int equiv(const SO& so);
00089 };
00090
00091 struct SO_block {
00092 int len;
00093 SO *so;
00094
00095 SO_block();
00096 SO_block(int);
00097 ~SO_block();
00098
00099 void set_length(int);
00100 void reset_length(int);
00101
00102 int add(SO& s, int i);
00103 void print(const char *title);
00104 };
00105
00106
00107
00108
00109 class PetiteList : public RefCount {
00110 private:
00111 int natom_;
00112 int nshell_;
00113 int ng_;
00114 int nirrep_;
00115 int nblocks_;
00116 int c1_;
00117
00118 Ref<GaussianBasisSet> gbs_;
00119 Ref<Integral> ints_;
00120
00121 char *p1_;
00122 int **atom_map_;
00123
00124 int **shell_map_;
00125
00126 char *lamij_;
00127
00128 int *nbf_in_ir_;
00129
00130 void init();
00131
00132 public:
00133 PetiteList(const Ref<GaussianBasisSet>&, const Ref<Integral>&);
00134 ~PetiteList();
00135
00136 int nirrep() const { return nirrep_; }
00137 int order() const { return ng_; }
00138 int atom_map(int n, int g) const { return (c1_) ? n : atom_map_[n][g]; }
00139 int shell_map(int n, int g) const { return (c1_) ? n : shell_map_[n][g]; }
00140 int lambda(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
00141 int lambda(int i, int j) const
00142 { return (c1_) ? 1 : (int) lamij_[ij_offset(i,j)]; }
00143
00144 int in_p1(int n) const { return (c1_) ? 1 : (int) p1_[n]; }
00145 int in_p2(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
00146 int in_p4(int ij, int kl, int i, int j, int k, int l) const;
00147
00148 int nfunction(int i) const
00149 { return (c1_) ? gbs_->nbasis() : nbf_in_ir_[i]; }
00150
00151 int nblocks() const { return nblocks_; }
00152
00153 void print(std::ostream& =ExEnv::out0(), int verbose=1);
00154
00155
00156 RefSCDimension AO_basisdim();
00157 RefSCDimension SO_basisdim();
00158
00159
00160 RefSCMatrix r(int g);
00161
00162
00163 SO_block * aotoso_info();
00164 RefSCMatrix aotoso();
00165 RefSCMatrix sotoao();
00166
00167
00168 void symmetrize(const RefSymmSCMatrix& skel, const RefSymmSCMatrix& sym);
00169
00170
00171
00172 RefSymmSCMatrix to_SO_basis(const RefSymmSCMatrix&);
00173
00174
00175 RefSymmSCMatrix to_AO_basis(const RefSymmSCMatrix&);
00176
00177
00178
00179 RefSCMatrix evecs_to_AO_basis(const RefSCMatrix&);
00180
00181 RefSCMatrix evecs_to_SO_basis(const RefSCMatrix&);
00182 };
00183
00184 inline int
00185 PetiteList::in_p4(int ij, int kl, int i, int j, int k, int l) const
00186 {
00187 if (c1_)
00188 return 1;
00189
00190 sc_int_least64_t ijkl = i_offset64(ij)+kl;
00191 int nijkl=0;
00192
00193 for (int g=0; g < ng_; g++) {
00194 int gij = ij_offset(shell_map_[i][g],shell_map_[j][g]);
00195 int gkl = ij_offset(shell_map_[k][g],shell_map_[l][g]);
00196 sc_int_least64_t gijkl = ij_offset64(gij,gkl);
00197
00198 if (gijkl > ijkl)
00199 return 0;
00200 else if (gijkl == ijkl)
00201 nijkl++;
00202 }
00203
00204 return ng_/nijkl;
00205 }
00206
00207 }
00208
00209 #endif
00210
00211
00212
00213
00214