1 //
2 // clhf.cc --- implementation of the closed shell Hartree-Fock SCF class
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Edward Seidl <seidl@janed.com>
7 // Maintainer: LPS
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 __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <math.h>
33
34 #include <util/misc/timer.h>
35 #include <util/misc/formio.h>
36 #include <util/state/stateio.h>
37
38 #include <chemistry/qc/basis/petite.h>
39
40 #include <chemistry/qc/scf/clhf.h>
41 #include <chemistry/qc/scf/lgbuild.h>
42 #include <chemistry/qc/scf/clhftmpl.h>
43
44 using namespace std;
45
46 ///////////////////////////////////////////////////////////////////////////
47 // CLHF
48
49 static ClassDesc CLHF_cd(
50 typeid(CLHF),"CLHF",1,"public CLSCF",
51 0, create<CLHF>, create<CLHF>);
52
53 CLHF::CLHF(StateIn& s) :
54 SavableState(s),
55 CLSCF(s)
56 {
57 }
58
59 CLHF::CLHF(const Ref<KeyVal>& keyval) :
60 CLSCF(keyval)
61 {
62 }
63
64 CLHF::~CLHF()
65 {
66 }
67
68 void
69 CLHF::save_data_state(StateOut& s)
70 {
71 CLSCF::save_data_state(s);
72 }
73
74 int
75 CLHF::value_implemented() const
76 {
77 return 1;
78 }
79
80 int
81 CLHF::gradient_implemented() const
82 {
83 return 1;
84 }
85
86 void
87 CLHF::print(ostream&o) const
88 {
89 CLSCF::print(o);
90 }
91
92 //////////////////////////////////////////////////////////////////////////////
93
94 void
95 CLHF::ao_fock(double accuracy)
96 1 0 0 0 0 0 0 0 0 0 {
97 int i;
98 int nthread = threadgrp_->nthread();
99
100 Ref<PetiteList> pl = integral()->petite_list(basis());
101
102 // calculate G. First transform cl_dens_diff_ to the AO basis, then
103 // scale the off-diagonal elements by 2.0
104 tim_enter("setup");
105 RefSymmSCMatrix dd = cl_dens_diff_;
106 cl_dens_diff_ = pl->to_AO_basis(dd);
107 cl_dens_diff_->scale(2.0);
108 cl_dens_diff_->scale_diagonal(0.5);
109 tim_exit("setup");
110
111 // now try to figure out the matrix specialization we're dealing with
112 // if we're using Local matrices, then there's just one subblock, or
113 // see if we can convert G and P to local matrices
114
115 if (debug_>1) {
116 cl_gmat_.print("cl_gmat before build");
117 cl_dens_diff_.print("cl_dens_diff before build");
118 }
119
120 if (local_ || local_dens_) {
121 // grab the data pointers from the G and P matrices
122 double *gmat, *pmat;
123 tim_enter("local data");
124 RefSymmSCMatrix gtmp = get_local_data(cl_gmat_, gmat, SCF::Accum);
125 RefSymmSCMatrix ptmp = get_local_data(cl_dens_diff_, pmat, SCF::Read);
126 tim_exit("local data");
127
128 tim_enter("init pmax");
129 signed char * pmax = init_pmax(pmat);
130 tim_exit("init pmax");
131
132 tim_enter("ao_gmat");
133 LocalGBuild<LocalCLHFContribution> **gblds =
134 new LocalGBuild<LocalCLHFContribution>*[nthread];
135 LocalCLHFContribution **conts = new LocalCLHFContribution*[nthread];
136
137 double **gmats = new double*[nthread];
138 gmats[0] = gmat;
139
140 Ref<GaussianBasisSet> bs = basis();
141 int ntri = i_offset(bs->nbasis());
142
143 double gmat_accuracy = accuracy;
144 if (min_orthog_res() < 1.0) { gmat_accuracy *= min_orthog_res(); }
145
146 for (i=0; i < nthread; i++) {
147 if (i) {
148 gmats[i] = new double[ntri];
149 memset(gmats[i], 0, sizeof(double)*ntri);
150 }
151 conts[i] = new LocalCLHFContribution(gmats[i], pmat);
152 gblds[i] = new LocalGBuild<LocalCLHFContribution>(*conts[i], tbis_[i],
153 pl, bs, scf_grp_, pmax, gmat_accuracy, nthread, i
154 );
155
156 threadgrp_->add_thread(i, gblds[i]);
157 }
158
159 tim_enter("start thread");
160 if (threadgrp_->start_threads() < 0) {
161 ExEnv::err() << node0 << indent
162 << "CLHF: error starting threads" << endl;
163 abort();
164 }
165 tim_exit("start thread");
166
167 tim_enter("stop thread");
168 if (threadgrp_->wait_threads() < 0) {
169 ExEnv::err() << node0 << indent
170 << "CLHF: error waiting for threads" << endl;
171 abort();
172 }
173 tim_exit("stop thread");
174
175 double tnint=0;
176 for (i=0; i < nthread; i++) {
177 tnint += gblds[i]->tnint;
178
179 if (i) {
180 for (int j=0; j < ntri; j++)
181 gmat[j] += gmats[i][j];
182 delete[] gmats[i];
183 }
184
185 delete gblds[i];
186 delete conts[i];
187 }
188
189 delete[] gmats;
190 delete[] gblds;
191 delete[] conts;
192 delete[] pmax;
193
194 scf_grp_->sum(&tnint, 1, 0, 0);
195 ExEnv::out() << node0 << indent << scprintf("%20.0f integrals\n", tnint);
196
197 tim_exit("ao_gmat");
198
199 // if we're running on multiple processors, then sum the G matrix
200 tim_enter("sum");
201 if (scf_grp_->n() > 1)
202 scf_grp_->sum(gmat, i_offset(basis()->nbasis()));
203 tim_exit("sum");
204
205 // if we're running on multiple processors, or we don't have local
206 // matrices, then accumulate gtmp back into G
207 tim_enter("accum");
208 if (!local_ || scf_grp_->n() > 1)
209 cl_gmat_->convert_accumulate(gtmp);
210 tim_exit("accum");
211 }
212
213 // for now quit
214 else {
215 ExEnv::err() << node0 << indent << "Cannot yet use anything but Local matrices\n";
216 abort();
217 }
218
219 tim_enter("symm");
220 // get rid of AO delta P
221 cl_dens_diff_ = dd;
222 dd = cl_dens_diff_.clone();
223
224 // now symmetrize the skeleton G matrix, placing the result in dd
225 RefSymmSCMatrix skel_gmat = cl_gmat_.copy();
226 skel_gmat.scale(1.0/(double)pl->order());
227 if (debug_>1) {
228 skel_gmat.print("skel_gmat before symmetrize");
229 }
230 pl->symmetrize(skel_gmat,dd);
231 if (debug_>1) {
232 dd.print("dd after symmetrize");
233 }
234 tim_exit("symm");
235
236 // F = H+G
237 cl_fock_.result_noupdate().assign(hcore_);
238 cl_fock_.result_noupdate().accumulate(dd);
239 accumddh_->accum(cl_fock_.result_noupdate());
240 cl_fock_.computed()=1;
241 0 0 1 0 0 0 0 0 0 0 }
242
243 /////////////////////////////////////////////////////////////////////////////
244
245 void
246 CLHF::two_body_energy(double &ec, double &ex)
247 {
248 tim_enter("clhf e2");
249 ec = 0.0;
250 ex = 0.0;
251
252 if (local_ || local_dens_) {
253 // grab the data pointers from the G and P matrices
254 double *pmat;
255 tim_enter("local data");
256 RefSymmSCMatrix dens = ao_density();
257 dens->scale(2.0);
258 dens->scale_diagonal(0.5);
259 RefSymmSCMatrix ptmp = get_local_data(dens, pmat, SCF::Read);
260 tim_exit("local data");
261
262 // initialize the two electron integral classes
263 Ref<TwoBodyInt> tbi = integral()->electron_repulsion();
264 tbi->set_integral_storage(0);
265
266 tim_enter("init pmax");
267 signed char * pmax = init_pmax(pmat);
268 tim_exit("init pmax");
269
270 LocalCLHFEnergyContribution lclc(pmat);
271 Ref<PetiteList> pl = integral()->petite_list();
272 LocalGBuild<LocalCLHFEnergyContribution>
273 gb(lclc, tbi, pl, basis(), scf_grp_, pmax,
274 1.e-20/*desired_value_accuracy()/100.0*/);
275 gb.run();
276
277 delete[] pmax;
278
279 ec = lclc.ec;
280 ex = lclc.ex;
281 }
282 else {
283 ExEnv::err() << node0 << indent << "Cannot yet use anything but Local matrices\n";
284 abort();
285 }
286 tim_exit("clhf e2");
287 }
288
289 /////////////////////////////////////////////////////////////////////////////
290
291 void
292 CLHF::two_body_deriv(double * tbgrad)
293 {
294 two_body_deriv_hf(tbgrad, 1.0);
295 }
296
297 /////////////////////////////////////////////////////////////////////////////
298
299 // Local Variables:
300 // mode: c++
301 // c-file-style: "ETS"
302 // End:
303