1 // 2 // r12int_eval.h 3 // 4 // Copyright (C) 2004 Edward Valeev 5 // 6 // Author: Edward Valeev <edward.valeev@chemistry.gatech.edu> 7 // Maintainer: EV 8 // 9 // This file is part of the SC Toolkit. 10 // 11 // The SC Toolkit is free software; you can redistribute it and/or modify 12 // it under the terms of the GNU Library General Public License as published by 13 // the Free Software Foundation; either version 2, or (at your option) 14 // any later version. 15 // 16 // The SC Toolkit is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU Library General Public License for more details. 20 // 21 // You should have received a copy of the GNU Library General Public License 22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to 23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 24 // 25 // The U.S. Government is granted a limited license as per AL 91-7. 26 // 27 28 #ifdef __GNUG__ 29 #pragma interface 30 #endif 31 32 #ifndef _chemistry_qc_mbptr12_r12inteval_h 33 #define _chemistry_qc_mbptr12_r12inteval_h 34 35 #include <util/ref/ref.h> 36 #include <chemistry/qc/mbptr12/vxb_eval_info.h> 37 #include <chemistry/qc/mbptr12/linearr12.h> 38 #include <chemistry/qc/mbptr12/r12_amps.h> 39 40 namespace sc { 41 42 /** R12IntEval is the top-level class which computes intermediates occuring in linear R12 theories. 43 This class is used by all Wavefunction classes that implement linear R12 methods. 44 */ 45 46 class R12IntEval : virtual public SavableState { 47 48 bool evaluated_; 49 50 // Calculation information (number of basis functions, R12 approximation, etc.) 51 Ref<R12IntEvalInfo> r12info_; 52 53 // Note that intermediate B is symmetric but is stored as a full matrix to simplify the code 54 // that computes asymmetric form of B 55 RefSCMatrix Vaa_, Vab_, Xaa_, Xab_, Baa_, Bab_, Aaa_, Aab_, T2aa_, T2ab_; 56 RefSCMatrix Ac_aa_, Ac_ab_; 57 RefSCMatrix Raa_, Rab_; // Not sure if I'll compute and keep these explicitly later 58 Ref<R12Amplitudes> Amps_; // Amplitudes of various R12-contributed terms in pair functions 59 RefSCVector emp2pair_aa_, emp2pair_ab_; 60 RefSCDimension dim_ij_aa_, dim_ij_ab_, dim_ij_s_, dim_ij_t_; 61 RefSCDimension dim_ab_aa_, dim_ab_ab_; 62 63 bool gbc_; 64 bool ebc_; 65 LinearR12::ABSMethod abs_method_; 66 LinearR12::StandardApproximation stdapprox_; 67 bool spinadapted_; 68 bool include_mp1_; 69 /// should I follow Klopper-Samson approach in the intermediates formulation for the EBC-free method? 70 bool follow_ks_ebcfree_; 71 int debug_; 72 73 // Map to TwoBodyMOIntsTransform objects that have been computed previously 74 typedef std::map<std::string, Ref<TwoBodyMOIntsTransform> > TformMap; 75 TformMap tform_map_; 76 // Returns pointer to the appropriate transform. 77 // If the transform is not found then throw runtime_error 78 Ref<TwoBodyMOIntsTransform> get_tform_(const std::string&); 79 80 /// Fock-weighted occupied space |i_f> = f_i^R |R>, where R is a function in RI-BS 81 Ref<MOIndexSpace> focc_space_; 82 /// Form Fock-weighted occupied space 83 void form_focc_space_(); 84 /// Fock-weighted active occupied space |i_f> = f_i^R |R>, where R is a function in RI-BS 85 Ref<MOIndexSpace> factocc_space_; 86 /// Form Fock-weighted active occupied space 87 void form_factocc_space_(); 88 /// Space of canonical virtual MOs 89 Ref<MOIndexSpace> canonvir_space_; 90 /// Form space of auxiliary virtuals 91 void form_canonvir_space_(); 92 93 /// Initialize standard transforms 94 void init_tforms_(); 95 /// Set intermediates to zero + add the "diagonal" contributions 96 void init_intermeds_(); 97 /// Compute r^2 contribution to X 98 void r2_contrib_to_X_orig_(); 99 /// Compute r^2 contribution to X using compute_r2_() 100 void r2_contrib_to_X_new_(); 101 /// Compute <space1 space1|r_{12}^2|space1 space2> matrix 102 RefSCMatrix compute_r2_(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2); 103 /** Compute the Fock matrix between 2 spaces. scale_J and scale_K are used to scale Coulomb 104 and exchange contributions */ 105 RefSCMatrix fock_(const Ref<MOIndexSpace>& occ_space, const Ref<MOIndexSpace>& bra_space, 106 const Ref<MOIndexSpace>& ket_space, double scale_J = 1.0, double scale_K = 1.0); 107 /// Compute the coulomb matrix between 2 spaces 108 RefSCMatrix coulomb_(const Ref<MOIndexSpace>& occ_space, const Ref<MOIndexSpace>& bra_space, 109 const Ref<MOIndexSpace>& ket_space); 110 /// Compute the exchange matrix between 2 spaces 111 RefSCMatrix exchange_(const Ref<MOIndexSpace>& occ_space, const Ref<MOIndexSpace>& bra_space, 112 const Ref<MOIndexSpace>& ket_space); 113 114 /// Checkpoint the top-level molecular energy 115 void checkpoint_() const; 116 117 /** Compute the number tasks which have access to the integrals. 118 map_to_twi is a vector<int> each element of which will hold the 119 number of tasks with access to the integrals of lower rank than that task 120 (or -1 if the task doesn't have access to the integrals) */ 121 const int tasks_with_ints_(const Ref<R12IntsAcc> ints_acc, vector<int>& map_to_twi); 122 123 /** Compute contribution to V, X, and B of the following form: 124 0.5 * \bar{g}_{ij}^{pq} * \bar{r}_{pq}^{kl}, where p and q span mospace. 125 tform_name is the name of the transform to be used to get the integrals. 126 mospace is either space2() or space4() of that transform 127 */ 128 void contrib_to_VXB_a_symm_(const std::string& tform_name); 129 130 /** Compute contribution to V, X, and B of the following form: 131 \bar{g}_{ij}^{am} * \bar{r}_{am}^{kl}, where m and a span space1 and space2, respectively. 132 tform_name is the name of the transform to be used to get the integrals. 133 mospace1 and mospace2 are space2() and space4() of that transform, respectively 134 */ 135 void contrib_to_VXB_a_asymm_(const std::string& tform_name); 136 137 /// Compute OBS contribution to V, X, and B (these contributions are independent of the method) 138 void obs_contrib_to_VXB_gebc_vbseqobs_(); 139 140 /// Compute 1-ABS contribution to V, X, and B (these contributions are independent of the method) 141 void abs1_contrib_to_VXB_gebc_(); 142 143 /// Equiv to the sum of above, except for this doesn't assume that VBS is the same as OBS 144 void contrib_to_VXB_gebc_vbsneqobs_(); 145 146 /// Compute A using the "simple" formula obtained using direct substitution alpha'->a' 147 void compute_A_simple_(); 148 149 /// Compute A using the standard commutator approach 150 void compute_A_via_commutator_(); 151 152 /// Compute MP2 T2 153 void compute_T2_(); 154 155 /// Compute R "intermediate" (r12 integrals in occ-pair/vir-pair basis) 156 void compute_R_(); 157 158 /// Compute A*T2 contribution to V (needed if EBC is not assumed) 159 void AT2_contrib_to_V_(); 160 161 /// Compute -2*A*R contribution to B (needed if EBC is not assumed) 162 void AR_contrib_to_B_(); 163 164 /** Compute the first (r<sub>kl</sub>^<sup>AB</sup> f<sub>A</sub><sup>m</sup> r<sub>mB</sub>^<sup>ij</sup>) 165 contribution to B that vanishes under GBC */ 166 void compute_B_gbc_1_(); 167 168 /** Compute the second (r<sub>kl</sub>^<sup>AB</sup> r<sub>AB</sub>^<sup>Kj</sup> f<sub>K</sub><sup>i</sup>) 169 contribution to B that vanishes under GBC */ 170 void compute_B_gbc_2_(); 171 172 /// Compute dual-basis MP2 energy (contribution from doubles to MP2 energy) 173 void compute_dualEmp2_(); 174 175 /// Compute dual-basis MP1 energy (contribution from singles to HF energy) 176 void compute_dualEmp1_(); 177 178 /// This function computes T2 amplitudes 179 void compute_T2_vbsneqobs_(); 180 181 /** New general function to compute <ij|r<sub>12</sub>|pq> integrals. ipjq_tform 182 is the source of the integrals.*/ 183 void compute_R_vbsneqobs_(const Ref<TwoBodyMOIntsTransform>& ipjq_tform, 184 RefSCMatrix& Raa, RefSCMatrix& Rab); 185 186 /** Initialize amplitude objects */ 187 void compute_amps_(); 188 189 /** Sum contributions to SCMatrix A from all nodes and broadcast so 190 every node has the correct SCMatrix. If to_all_tasks is false, then 191 collect all contributions to task 0 and zero out the matrix on other tasks, 192 otherwise distribute the sum to every task. If to_average is true then 193 each result is scaled by the inverse of the number of tasks. */ 194 void globally_sum_scmatrix_(RefSCMatrix& A, bool to_all_tasks = false, bool to_average = false); 195 196 /** Sum contributions to SCVector A from all nodes and broadcast so 197 every node has the correct SCVector. If to_all_tasks is false, then 198 collect all contributions to task 0 and zero out the matrix on other tasks, 199 otherwise distribute the sum to every task. If to_average is true then 200 each result is scaled by the inverse of the number of tasks. */ 201 void globally_sum_scvector_(RefSCVector& A, bool to_all_tasks = false, bool to_average = false); 202 203 /** Sum contributions to the intermediates from all nodes and broadcast so 204 every node has the correct matrices. If to_all_tasks is false, then 205 collect all contributions to task 0 and zero out the matrix on other tasks, 206 otherwise distribute the sum to every task. */ 207 void globally_sum_intermeds_(bool to_all_tasks = false); 208 209 public: 210 R12IntEval(StateIn&); 211 /** Constructs R12IntEval. If follow_ks_ebcfree is true then follow formalism of Klopper and Samson 212 to compute EBC-free MP2-R12 energy. */ 213 R12IntEval(const Ref<R12IntEvalInfo>& info, bool gbc = true, bool ebc = true, 214 LinearR12::ABSMethod abs_method = LinearR12::ABS_CABSPlus, 215 LinearR12::StandardApproximation stdapprox = LinearR12::StdApprox_Ap, 216 bool follow_ks_ebcfree = false); 217 ~R12IntEval(); 218 219 void save_data_state(StateOut&); 220 virtual void obsolete(); 221 222 void include_mp1(bool include_mp1); 223 void set_debug(int debug); 224 void set_dynamic(bool dynamic); 225 void set_print_percent(double print_percent); 226 void set_memory(size_t nbytes); 227 gbc()228 const bool gbc() const { return gbc_; } ebc()229 const bool ebc() const { return ebc_; } stdapprox()230 const LinearR12::StandardApproximation stdapprox() const { return stdapprox_; } follow_ks_ebcfree()231 bool follow_ks_ebcfree() const { return follow_ks_ebcfree_; } 232 233 Ref<R12IntEvalInfo> r12info() const; 234 RefSCDimension dim_oo_aa() const; 235 RefSCDimension dim_oo_ab() const; 236 RefSCDimension dim_oo_s() const; 237 RefSCDimension dim_oo_t() const; 238 RefSCDimension dim_vv_aa() const; 239 RefSCDimension dim_vv_ab() const; 240 241 /// This function causes the intermediate matrices to be computed. 242 virtual void compute(); 243 244 /// Returns alpha-alpha block of the V intermediate matrix. 245 RefSCMatrix V_aa(); 246 /// Returns alpha-alpha block of the X intermediate matrix. 247 RefSCMatrix X_aa(); 248 /// Returns alpha-alpha block of the B intermediate matrix. 249 RefSymmSCMatrix B_aa(); 250 /// Returns alpha-alpha block of the A intermediate matrix. Returns 0 if EBC is assumed. 251 RefSCMatrix A_aa(); 252 /// Returns alpha-alpha block of the A intermediate matrix. Returns 0 if EBC is assumed. 253 RefSCMatrix Ac_aa(); 254 /// Returns alpha-alpha block of the MP2 T2 matrix. Returns 0 if EBC is assumed. 255 RefSCMatrix T2_aa(); 256 /// Returns alpha-beta block of the V intermediate matrix. 257 RefSCMatrix V_ab(); 258 /// Returns alpha-beta block of the X intermediate matrix. 259 RefSCMatrix X_ab(); 260 /// Returns alpha-beta block of the B intermediate matrix. 261 RefSymmSCMatrix B_ab(); 262 /// Returns alpha-beta block of the A intermediate matrix. Returns 0 if EBC is assumed 263 RefSCMatrix A_ab(); 264 /// Returns alpha-beta block of the A intermediate matrix. Returns 0 if EBC is assumed 265 RefSCMatrix Ac_ab(); 266 /// Returns alpha-beta block of the MP2 T2 matrix. Returns 0 if EBC is assumed 267 RefSCMatrix T2_ab(); 268 /// Returns alpha-alpha MP2 pair energies. 269 RefSCVector emp2_aa(); 270 /// Returns alpha-beta MP2 pair energies. 271 RefSCVector emp2_ab(); 272 /// Returns amplitudes of pair correlation functions 273 Ref<R12Amplitudes> amps(); 274 275 RefDiagSCMatrix evals() const; 276 }; 277 278 } 279 280 #endif 281 282 // Local Variables: 283 // mode: c++ 284 // c-file-style: "CLJ" 285 // End: 286 287 288