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_mbptr12_pairiter_h
00029 #define _chemistry_qc_mbptr12_pairiter_h
00030
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034
00035 #include <stdexcept>
00036 #include <chemistry/qc/mbptr12/moindexspace.h>
00037
00038 namespace sc {
00039
00041 class MOPairIter : public RefCount {
00042
00043 protected:
00044 bool i_eq_j_;
00045 int ni_;
00046 int nj_;
00047 int i_;
00048 int j_;
00049 int nij_;
00050 int ij_;
00051
00052 public:
00054 MOPairIter(const Ref<MOIndexSpace>& space_i, const Ref<MOIndexSpace>& space_j);
00055 virtual ~MOPairIter();
00056
00058 virtual void start(const int first_ij =0) =0;
00060 virtual void next() =0;
00062 virtual operator int() const =0;
00063
00065 int ni() const { return ni_; }
00067 int nj() const { return nj_; }
00069 int i() const { return i_; }
00071 int j() const { return j_; }
00073 int nij() const { return nij_; }
00075 int ij() const { return ij_; }
00076 };
00077
00078
00081 class SpatialMOPairIter : public MOPairIter {
00082
00083 public:
00085 SpatialMOPairIter(const Ref<MOIndexSpace>& space_i, const Ref<MOIndexSpace>& space_j) :
00086 MOPairIter(space_i,space_j) {};
00087 ~SpatialMOPairIter() {};
00088
00090 virtual int nij_aa() const =0;
00092 virtual int nij_ab() const =0;
00095 virtual int ij_aa() const =0;
00097 virtual int ij_ab() const =0;
00099 virtual int ij_ba() const =0;
00100 };
00101
00105 class SpatialMOPairIter_eq : public SpatialMOPairIter {
00106
00107 int nij_aa_;
00108 int nij_ab_;
00109 int ij_aa_;
00110 int ij_ab_;
00111 int ji_ab_;
00112
00113 void init_ij(const int ij) {
00114
00115 if (ij<0)
00116 throw std::runtime_error("SpatialMOPairIter_eq::start() -- argument ij out of range");
00117
00118 ij_ = 0;
00119 const int renorm_ij = ij%nij_;
00120
00121 i_ = (int)floor((sqrt(1.0+8.0*renorm_ij) - 1.0)/2.0);
00122 const int i_off = i_*(i_+1)/2;
00123 j_ = renorm_ij - i_off;
00124
00125 ij_ab_ = i_*nj_ + j_;
00126 ji_ab_ = j_*ni_ + i_;
00127
00128 if (i_ != 0) {
00129 const int i_off = i_*(i_-1)/2;
00130 ij_aa_ = i_off + j_;
00131 if (i_ == j_)
00132 ij_aa_--;
00133 }
00134 else {
00135 ij_aa_ = -1;
00136 }
00137 };
00138
00139 void inc_ij() {
00140 ij_++;
00141 if (ij_ab_ == nij_ab_-1) {
00142 i_ = 0;
00143 j_ = 0;
00144 ij_ab_ = 0;
00145 ji_ab_ = 0;
00146 ij_aa_ = -1;
00147 }
00148 else {
00149 if (i_ == j_) {
00150 i_++;
00151 j_ = 0;
00152 ji_ab_ = i_;
00153 ij_ab_ = i_*nj_;
00154 ij_aa_ += (i_ == j_) ? 0 : 1;
00155 }
00156 else {
00157 j_++;
00158 ji_ab_ += ni_;
00159 ij_ab_++;
00160 ij_aa_ += (i_ == j_) ? 0 : 1;
00161 }
00162 }
00163 };
00164
00165 public:
00167 SpatialMOPairIter_eq(const Ref<MOIndexSpace>& space1);
00168 ~SpatialMOPairIter_eq();
00169
00171 void start(const int ij_offset=0)
00172 {
00173 ij_ = 0;
00174 init_ij(ij_offset);
00175 };
00176
00178 void next() {
00179 inc_ij();
00180 };
00182 operator int() const { return (nij_ > ij_);};
00183
00185 int nij_aa() const { return nij_aa_; }
00187 int nij_ab() const { return nij_ab_; }
00190 int ij_aa() const { return (i_ == j_) ? -1 : ij_aa_; }
00192 int ij_ab() const { return ij_ab_; }
00194 int ij_ba() const { return ji_ab_; }
00195 };
00196
00197
00200 class SpatialMOPairIter_neq : public SpatialMOPairIter {
00201
00202 int IJ_;
00203
00204 void init_ij(const int ij) {
00205
00206 if (ij<0)
00207 throw std::runtime_error("SpatialMOPairIter_neq::start() -- argument ij out of range");
00208
00209 IJ_ = 0;
00210 const int renorm_ij = ij%nij_;
00211
00212 i_ = renorm_ij/nj_;
00213 j_ = renorm_ij - i_*nj_;
00214
00215 IJ_ = i_*nj_ + j_;
00216
00217 };
00218
00219 void inc_ij() {
00220 ij_++;
00221 IJ_++;
00222 if (IJ_ == nij_) {
00223 i_ = 0;
00224 j_ = 0;
00225 IJ_ = 0;
00226 }
00227 else {
00228 if (j_ == nj_-1) {
00229 i_++;
00230 j_ = 0;
00231 }
00232 else {
00233 j_++;
00234 }
00235 }
00236 };
00237
00238 public:
00240 SpatialMOPairIter_neq(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2);
00241 ~SpatialMOPairIter_neq();
00242
00244 void start(const int ij_offset=0)
00245 {
00246 ij_ = 0;
00247 init_ij(ij_offset);
00248 };
00249
00251 void next() {
00252 inc_ij();
00253 };
00255 operator int() const { return (nij_ > ij_);};
00256
00258 int nij_aa() const { return nij_; }
00260 int nij_ab() const { return nij_; }
00262 int ij_aa() const { return IJ_; }
00264 int ij_ab() const { return IJ_; }
00266 int ij_ba() const { return IJ_; }
00267 };
00268
00269
00272 class MOPairIterFactory {
00273
00274 public:
00275 MOPairIterFactory() {}
00276 ~MOPairIterFactory() {}
00277
00279 Ref<SpatialMOPairIter> mopairiter(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2);
00281 RefSCDimension scdim_aa(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2);
00283 RefSCDimension scdim_ab(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2);
00284 };
00285
00286 }
00287
00288 #endif
00289
00290
00291
00292
00293