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_gpetite_h
00029 #define _chemistry_qc_basis_gpetite_h
00030
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034
00035 #include <stdexcept>
00036
00037 #include <scconfig.h>
00038 #include <util/misc/scint.h>
00039 #include <chemistry/qc/basis/basis.h>
00040 #include <chemistry/qc/basis/petite.h>
00041
00042 namespace sc {
00043
00047 class canonical_aaaa {
00048 public:
00049 canonical_aaaa();
00050 canonical_aaaa(const Ref<GaussianBasisSet> bi,
00051 const Ref<GaussianBasisSet> bj,
00052 const Ref<GaussianBasisSet> bk,
00053 const Ref<GaussianBasisSet> bl
00054 );
00055 sc_int_least64_t offset(int i, int j, int k, int l) {
00056 long ij = (i>j?(((i*long(i+1))>>1)+j):(((j*long(j+1))>>1)+i));
00057 long kl = (k>l?(((k*long(k+1))>>1)+l):(((l*long(l+1))>>1)+k));
00058 sc_int_least64_t
00059 off = (ij>kl?(((ij*sc_int_least64_t(ij+1))>>1)+kl)
00060 :(((kl*sc_int_least64_t(kl+1))>>1)+ij));
00061 return off;
00062 }
00063 };
00064
00069 class canonical_aabc {
00070 long nk_, nl_;
00071 public:
00072 canonical_aabc(const Ref<GaussianBasisSet> bi,
00073 const Ref<GaussianBasisSet> bj,
00074 const Ref<GaussianBasisSet> bk,
00075 const Ref<GaussianBasisSet> bl
00076 );
00077 sc_int_least64_t offset(int i, int j, int k, int l) {
00078 long ij = (i>j?(((i*long(i+1))>>1)+j):(((j*long(j+1))>>1)+i));
00079 return k + nk_*sc_int_least64_t(l + nl_*ij);
00080 }
00081 };
00082
00087 class canonical_aabb {
00088 long nij_;
00089 public:
00090 canonical_aabb(const Ref<GaussianBasisSet> bi,
00091 const Ref<GaussianBasisSet> bj,
00092 const Ref<GaussianBasisSet> bk,
00093 const Ref<GaussianBasisSet> bl
00094 );
00095 sc_int_least64_t offset(int i, int j, int k, int l) {
00096 long ij = (i>j?(((i*long(i+1))>>1)+j):(((j*long(j+1))>>1)+i));
00097 long kl = (k>l?(((k*long(k+1))>>1)+l):(((l*long(l+1))>>1)+k));
00098 return ij + nij_*sc_int_least64_t(kl);
00099 }
00100 };
00101
00105 class canonical_abcd {
00106 int ni_, nj_, nk_;
00107 public:
00108 canonical_abcd(const Ref<GaussianBasisSet> bi,
00109 const Ref<GaussianBasisSet> bj,
00110 const Ref<GaussianBasisSet> bk,
00111 const Ref<GaussianBasisSet> bl
00112 );
00113 sc_int_least64_t offset(int i, int j, int k, int l) {
00114 return (i + ni_*sc_int_least64_t(j + nj_*long(k + nk_*l)));
00115 }
00116 };
00117
00124 template <class C4>
00125 class GPetite4: public RefCount {
00126 bool c1_;
00127 int ng_;
00128 C4 c_;
00129 int **shell_map_i_;
00130 int **shell_map_j_;
00131 int **shell_map_k_;
00132 int **shell_map_l_;
00133 Ref<GaussianBasisSet> b1_, b2_, b3_, b4_;
00134 public:
00135 GPetite4(const Ref<GaussianBasisSet> &b1,
00136 const Ref<GaussianBasisSet> &b2,
00137 const Ref<GaussianBasisSet> &b3,
00138 const Ref<GaussianBasisSet> &b4,
00139 const C4& c): c_(c) {
00140 int **atom_map;
00141 b1_ = b1;
00142 b2_ = b2;
00143 b3_ = b3;
00144 b4_ = b4;
00145
00146 ng_ = b1->molecule()->point_group()->char_table().order();
00147 if (b2->molecule()->point_group()->char_table().order() != ng_
00148 || b3->molecule()->point_group()->char_table().order() != ng_
00149 || b4->molecule()->point_group()->char_table().order() != ng_) {
00150 throw
00151 std::runtime_error("GPetite4: not all point groups are the same");
00152 }
00153 c1_ = (ng_ == 1);
00154
00155 atom_map = compute_atom_map(b1);
00156 shell_map_i_ = compute_shell_map(atom_map,b1);
00157 delete_atom_map(atom_map,b1);
00158
00159 atom_map = compute_atom_map(b2);
00160 shell_map_j_ = compute_shell_map(atom_map,b2);
00161 delete_atom_map(atom_map,b2);
00162
00163 atom_map = compute_atom_map(b3);
00164 shell_map_k_ = compute_shell_map(atom_map,b3);
00165 delete_atom_map(atom_map,b3);
00166
00167 atom_map = compute_atom_map(b4);
00168 shell_map_l_ = compute_shell_map(atom_map,b4);
00169 delete_atom_map(atom_map,b4);
00170 }
00171 ~GPetite4() {
00172 delete_shell_map(shell_map_i_,b1_);
00173 delete_shell_map(shell_map_j_,b2_);
00174 delete_shell_map(shell_map_k_,b3_);
00175 delete_shell_map(shell_map_l_,b4_);
00176 }
00177 int in_p4(int i, int j, int k, int l);
00178 };
00179
00180 template <class C4>
00181 inline int
00182 GPetite4<C4>::in_p4(int i, int j, int k, int l)
00183 {
00184 if (c1_) return 1;
00185
00186 sc_int_least64_t ijkl = c_.offset(i,j,k,l);
00187 int nijkl = 1;
00188
00189 for (int g=1; g < ng_; g++) {
00190 int gi = shell_map_i_[i][g];
00191 int gj = shell_map_j_[j][g];
00192 int gk = shell_map_k_[k][g];
00193 int gl = shell_map_l_[l][g];
00194 sc_int_least64_t gijkl = c_.offset(gi,gj,gk,gl);
00195
00196 if (gijkl > ijkl) return 0;
00197 else if (gijkl == ijkl) nijkl++;
00198 }
00199
00200 return ng_/nijkl;
00201 }
00202
00203 }
00204
00205
00206
00207 #endif
00208
00209
00210
00211
00212