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