← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/mpqc/oneiric

« back to all changes in this revision

Viewing changes to src/lib/chemistry/qc/mbptr12/dualbasis_mp2.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2005-11-27 11:41:49 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051127114149-zgz9r3gk50w8ww2q
Tags: 2.3.0-1
http://bugs.debian.org/339232
* New upstream release.
* debian/rules (SONAME): Activate awk snippet for automatic so-name
  detection again, resulting in a bump to `7' and making a `c2a' for
  the C++ allocator change unnecessary; closes: #339232.
* debian/patches/00list (08_gcc-4.0_fixes): Removed, no longer needed.
* debian/rules (test): Remove workarounds, do not abort build if tests
  fail.
* debian/ref: Removed.
* debian/control.in (libsc): Added Conflict against libsc6c2.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/chemistry/qc/mbptr12/dualbasis_mp2.cc
 
1
//
 
2
// dualbasis_mp2.cc
 
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
#include <stdexcept>
 
29
#include <stdlib.h>
 
30
#include <math.h>
 
31
#include <limits.h>
 
32
 
 
33
#include <scconfig.h>
 
34
#include <util/misc/formio.h>
 
35
#include <util/misc/timer.h>
 
36
#include <util/class/class.h>
 
37
#include <util/state/state.h>
 
38
#include <util/state/state_text.h>
 
39
#include <util/state/state_bin.h>
 
40
#include <math/scmat/matrix.h>
 
41
#include <chemistry/molecule/molecule.h>
 
42
#include <chemistry/qc/basis/integral.h>
 
43
#include <chemistry/qc/mbpt/bzerofast.h>
 
44
#include <chemistry/qc/mbptr12/r12ia.h>
 
45
#include <chemistry/qc/mbptr12/vxb_eval_info.h>
 
46
#include <chemistry/qc/mbptr12/pairiter.h>
 
47
#include <chemistry/qc/mbptr12/r12int_eval.h>
 
48
 
 
49
using namespace std;
 
50
using namespace sc;
 
51
 
 
52
void
 
53
R12IntEval::compute_dualEmp2_()
 
54
{
 
55
  if (evaluated_)
 
56
    return;
 
57
  Ref<MessageGrp> msg = r12info()->msg();
 
58
  Ref<MemoryGrp> mem = r12info()->mem();
 
59
  Ref<ThreadGrp> thr = r12info()->thr();
 
60
  const int num_te_types = 1;
 
61
  enum te_types {eri=0};
 
62
 
 
63
  tim_enter("dual-basis MP2 energy");
 
64
  ExEnv::out0() << endl << indent
 
65
               << "Entered dual-basis MP2 energy evaluator" << endl;
 
66
  ExEnv::out0() << incindent;
 
67
 
 
68
  int me = msg->me();
 
69
  int nproc = msg->n();
 
70
  
 
71
  // Do the AO->MO transform
 
72
  form_canonvir_space_();
 
73
  Ref<MOIntsTransformFactory> tfactory = r12info_->tfactory();
 
74
  tfactory->set_spaces(r12info_->act_occ_space(),canonvir_space_,
 
75
                       r12info_->act_occ_space(),canonvir_space_);
 
76
  Ref<TwoBodyMOIntsTransform> ipjq_tform = tfactory->twobody_transform_13("(ix|jy)");
 
77
  ipjq_tform->set_num_te_types(num_te_types);
 
78
  ipjq_tform->compute();
 
79
  Ref<R12IntsAcc> ijpq_acc = ipjq_tform->ints_acc();
 
80
  if (num_te_types != ijpq_acc->num_te_types())
 
81
    throw std::runtime_error("R12IntEval::compute_dualEmp2_() -- number of MO integral types is wrong");
 
82
 
 
83
  int nocc_act = r12info()->nocc_act();
 
84
  int ncanonvir = canonvir_space_->rank();
 
85
 
 
86
  ExEnv::out0() << indent << "Begin computation of energies" << endl;
 
87
  SpatialMOPairIter_eq kl_iter(r12info_->act_occ_space());
 
88
  int naa = kl_iter.nij_aa();          // Number of alpha-alpha pairs (i > j)
 
89
  int nab = kl_iter.nij_ab();          // Number of alpha-beta pairs
 
90
  if (debug_) {
 
91
    ExEnv::out0() << indent << "naa = " << naa << endl;
 
92
    ExEnv::out0() << indent << "nab = " << nab << endl;
 
93
  }
 
94
 
 
95
  // Compute the number of tasks that have full access to the integrals
 
96
  // and split the work among them
 
97
  vector<int> proc_with_ints;
 
98
  int nproc_with_ints = tasks_with_ints_(ijpq_acc,proc_with_ints);
 
99
 
 
100
  RefDiagSCMatrix act_occ_evals = r12info_->act_occ_space()->evals();
 
101
  RefDiagSCMatrix canonvir_evals = canonvir_space_->evals();
 
102
 
 
103
  if (ijpq_acc->has_access(me)) {
 
104
 
 
105
    for(kl_iter.start();int(kl_iter);kl_iter.next()) {
 
106
 
 
107
      const int kl = kl_iter.ij();
 
108
      // Figure out if this task will handle this kl
 
109
      int kl_proc = kl%nproc_with_ints;
 
110
      if (kl_proc != proc_with_ints[me])
 
111
        continue;
 
112
      const int k = kl_iter.i();
 
113
      const int l = kl_iter.j();
 
114
      const int kl_aa = kl_iter.ij_aa();
 
115
      const int kl_ab = kl_iter.ij_ab();
 
116
      const int lk_ab = kl_iter.ij_ba();
 
117
 
 
118
      if (debug_)
 
119
        ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
 
120
 
 
121
      // Get (|1/r12|) integrals
 
122
      tim_enter("MO ints retrieve");
 
123
      double *klxy_buf_eri = ijpq_acc->retrieve_pair_block(k,l,R12IntsAcc::eri);
 
124
      tim_exit("MO ints retrieve");
 
125
 
 
126
      if (debug_)
 
127
        ExEnv::outn() << indent << "task " << me << ": obtained kl blocks" << endl;
 
128
 
 
129
      // Compute MP2 energies
 
130
      double emp2_aa = 0.0;
 
131
      double emp2_ab = 0.0;
 
132
      for(int a=0; a<ncanonvir; a++) {
 
133
        for(int b=0; b<ncanonvir; b++) {
 
134
          const int ab_offset = a*ncanonvir+b;
 
135
          const int ba_offset = b*ncanonvir+a;
 
136
          const double oo_delta_ijab = 1.0/(act_occ_evals(k)+act_occ_evals(l)-canonvir_evals(a)-canonvir_evals(b));
 
137
          const double eri_kalb = klxy_buf_eri[ab_offset];
 
138
          const double eri_kbla = klxy_buf_eri[ba_offset];
 
139
          emp2_ab += 0.5*(eri_kalb * eri_kalb + eri_kbla * eri_kbla) * oo_delta_ijab;
 
140
          if (kl_aa != -1) {
 
141
            emp2_aa += (eri_kalb - eri_kbla) * (eri_kalb - eri_kbla) * oo_delta_ijab;
 
142
          }
 
143
        }
 
144
      }
 
145
      emp2pair_ab_.set_element(kl_ab,emp2_ab);
 
146
      if (kl_ab != lk_ab)
 
147
        emp2pair_ab_.set_element(lk_ab,emp2_ab);
 
148
      if (kl_aa != -1)
 
149
        emp2pair_aa_.set_element(kl_aa,emp2_aa);
 
150
 
 
151
      ijpq_acc->release_pair_block(k,l,R12IntsAcc::eri);
 
152
    }
 
153
  }
 
154
 
 
155
  // Tasks that don't do any work here still need to create these timers
 
156
  tim_enter("MO ints retrieve");
 
157
  tim_exit("MO ints retrieve");
 
158
 
 
159
  ExEnv::out0() << indent << "End of computation of energies" << endl;
 
160
  ijpq_acc->deactivate();
 
161
  
 
162
  globally_sum_intermeds_();
 
163
 
 
164
  ExEnv::out0() << decindent;
 
165
  ExEnv::out0() << indent << "Exited dual-basis MP2 energy evaluator" << endl;
 
166
 
 
167
  tim_exit("dual-basis MP2 energy");
 
168
  checkpoint_();
 
169
  
 
170
  return;
 
171
}
 
172
 
 
173
void
 
174
R12IntEval::compute_dualEmp1_()
 
175
{
 
176
  if (evaluated_)
 
177
    return;
 
178
  Ref<MessageGrp> msg = r12info()->msg();
 
179
  Ref<MemoryGrp> mem = r12info()->mem();
 
180
  Ref<ThreadGrp> thr = r12info()->thr();
 
181
  const int num_te_types = 1;
 
182
  enum te_types {eri=0};
 
183
 
 
184
  tim_enter("dual-basis MP1 energy");
 
185
  ExEnv::out0() << endl << indent
 
186
               << "Entered dual-basis MP1 energy evaluator" << endl;
 
187
  ExEnv::out0() << incindent;
 
188
 
 
189
  int me = msg->me();
 
190
  int nproc = msg->n();
 
191
  
 
192
  // Compute act.occ./aux.virt. Fock matrix
 
193
  form_canonvir_space_();
 
194
  Ref<MOIndexSpace> occ_space = r12info_->occ_space();
 
195
  RefSCMatrix F_aocc_canonvir = fock_(occ_space,occ_space,canonvir_space_);
 
196
 
 
197
  int nocc = r12info()->nocc();
 
198
  int ncanonvir = canonvir_space_->rank();
 
199
  RefDiagSCMatrix occ_evals = r12info_->occ_space()->evals();
 
200
  RefDiagSCMatrix canonvir_evals = canonvir_space_->evals();
 
201
 
 
202
  double emp1 = 0.0;
 
203
  for(int i=0; i<nocc; i++) {
 
204
    for(int a=0; a<ncanonvir; a++) {
 
205
      const double Fia = F_aocc_canonvir.get_element(i,a);
 
206
      emp1 += Fia*Fia/(-occ_evals(i)+canonvir_evals(a));
 
207
    }
 
208
  }
 
209
  ExEnv::out0() << indent << "MP1 energy correction to HF energy [au] :   "
 
210
                << 2.0*emp1 << endl;
 
211
  ExEnv::out0() << indent << "HF energy estimated in new basis [au]   :   "
 
212
                << r12info_->ref()->energy() - 2.0*emp1 << endl;
 
213
 
 
214
  ExEnv::out0() << decindent;
 
215
  ExEnv::out0() << endl << "Exited dual-basis MP1 energy evaluator" << endl;
 
216
 
 
217
  tim_exit("dual-basis MP1 energy");
 
218
}
 
219
 
 
220
////////////////////////////////////////////////////////////////////////////
 
221
 
 
222
// Local Variables:
 
223
// mode: c++
 
224
// c-file-style: "CLJ-CONDENSED"
 
225
// End: