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