1 //
2 // scfvector.cc --- implementation of SCF::compute_vector
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 #include <util/misc/timer.h>
29 #include <util/misc/formio.h>
30
31 #include <math/scmat/offset.h>
32 #include <math/scmat/blocked.h>
33 #include <math/scmat/blkiter.h>
34
35 #include <math/optimize/diis.h>
36 #include <math/optimize/scextrapmat.h>
37
38 #include <chemistry/qc/basis/symmint.h>
39
40 #include <chemistry/qc/scf/scf.h>
41 #include <chemistry/qc/scf/scfops.h>
42 #include <chemistry/qc/scf/scflocal.h>
43
44 using namespace std;
45
46 #undef GENERALIZED_EIGENSOLVER
47
48 ///////////////////////////////////////////////////////////////////////////
49
50 extern "C" {
51 void
52 dsygv_(int *ITYPE, const char *JOBZ, const char *UPLO,
53 int *N, double *A, int *LDA, double *B, int *LDB,
54 double *W, double *WORK, int *LWORK, int *INFO);
55 }
56
57 double
58 SCF::compute_vector(double& eelec)
59 {
60 tim_enter("vector");
61 int i;
62
63 // reinitialize the extrapolation object
64 extrap_->reinitialize();
65
66 // create level shifter
67 LevelShift *level_shift = new LevelShift(this);
68 level_shift->reference();
69
70 // calculate the core Hamiltonian
71 hcore_ = core_hamiltonian();
72
73 // add density independant contributions to Hcore
74 accumdih_->accum(hcore_);
75
76 // set up subclass for vector calculation
77 init_vector();
78
79 // calculate the nuclear repulsion energy
80 double nucrep = molecule()->nuclear_repulsion_energy();
81 ExEnv::out() << node0 << indent
82 << scprintf("nuclear repulsion energy = %15.10f", nucrep)
83 << endl << endl;
84
85 RefDiagSCMatrix evals(oso_dimension(), basis_matrixkit());
86
87 double delta = 1.0;
88 int iter, iter_since_reset = 0;
89 double accuracy = 1.0;
90 for (iter=0; iter < maxiter_; iter++, iter_since_reset++) {
91 // form the density from the current vector
92 tim_enter("density");
93 delta = new_density();
94 tim_exit("density");
95
96 // check convergence
97 if (delta < desired_value_accuracy()
98 && accuracy < desired_value_accuracy()) break;
99
100 // reset the density from time to time
101 if (iter_since_reset && !(iter_since_reset%dens_reset_freq_)) {
102 reset_density();
103 iter_since_reset = 0;
104 }
105
106 // form the AO basis fock matrix & add density dependant H
107 tim_enter("fock");
108 double base_accuracy = delta;
109 if (base_accuracy < desired_value_accuracy())
110 base_accuracy = desired_value_accuracy();
111 double new_accuracy = 0.01 * base_accuracy;
112 if (new_accuracy > 0.001) new_accuracy = 0.001;
113 if (iter == 0) accuracy = new_accuracy;
114 else if (new_accuracy < accuracy) {
115 accuracy = new_accuracy/10.0;
116 if (iter_since_reset > 0) {
117 reset_density();
118 iter_since_reset = 0;
119 }
120 }
121 ao_fock(accuracy);
122 tim_exit("fock");
123
124 // calculate the electronic energy
125 eelec = scf_energy();
126 double eother = 0.0;
127 if (accumddh_.nonnull()) eother = accumddh_->e();
128 ExEnv::out() << node0 << indent
129 << scprintf("iter %5d energy = %15.10f delta = %10.5e",
130 iter+1, eelec+eother+nucrep, delta)
131 << endl;
132
133 // now extrapolate the fock matrix
134 tim_enter("extrap");
135 Ref<SCExtrapData> data = extrap_data();
136 Ref<SCExtrapError> error = extrap_error();
137 extrap_->extrapolate(data,error);
138 data=0;
139 error=0;
140 tim_exit("extrap");
141
142 #ifdef GENERALIZED_EIGENSOLVER
143 // Get the fock matrix and overlap in the SO basis. The fock matrix
144 // used here works for CLOSED SHELL ONLY.
145 RefSymmSCMatrix bfmatref = fock(0);
146 RefSymmSCMatrix bsmatref = overlap();
147 BlockedSymmSCMatrix *bfmat
148 = dynamic_cast<BlockedSymmSCMatrix*>(bfmatref.pointer());
149 BlockedSymmSCMatrix *bsmat
150 = dynamic_cast<BlockedSymmSCMatrix*>(bsmatref.pointer());
151 BlockedDiagSCMatrix *bevals
152 = dynamic_cast<BlockedDiagSCMatrix*>(evals.pointer());
153 BlockedSCMatrix *bvec
154 = dynamic_cast<BlockedSCMatrix*>(oso_scf_vector_.pointer());
155
156 ExEnv::out() << node0 << indent
157 << "solving generalized eigenvalue problem" << endl;
158
159 for (int iblock=0; iblock<bfmat->nblocks(); iblock++) {
160 RefSymmSCMatrix fmat = bfmat->block(iblock);
161 RefSymmSCMatrix smat = bsmat->block(iblock);
162 RefDiagSCMatrix evalblock = bevals->block(iblock);
163 RefSCMatrix oso_scf_vector_block = bvec->block(iblock);
164 int nbasis = fmat.dim().n();
165 int nbasis2 = nbasis*nbasis;
166
167 if (!nbasis) continue;
168
169 // Convert to the lapack storage format.
170 double *fso = new double[nbasis2];
171 double *sso = new double[nbasis2];
172 int ij=0;
173 for (int i=0; i<nbasis; i++) {
174 for (int j=0; j<nbasis; j++,ij++) {
175 fso[ij] = fmat(i,j);
176 sso[ij] = smat(i,j);
177 }
178 }
179
180 // solve generalized eigenvalue problem with DSYGV
181 int itype = 1;
182 double *epsilon = new double[nbasis];
183 int lwork = -1;
184 int info;
185 double optlwork;
186 dsygv_(&itype,"V","U",&nbasis,fso,&nbasis,sso,&nbasis,
187 epsilon,&optlwork,&lwork,&info);
188 if (info) {
189 ExEnv::out() << "dsygv could not determine work size: info = "
190 << info << endl;
191 abort();
192 }
193 lwork = (int)optlwork;
194 double *work = new double[lwork];
195 dsygv_(&itype,"V","U",&nbasis,fso,&nbasis,sso,&nbasis,
196 epsilon,work,&lwork,&info);
197 if (info) {
198 ExEnv::out() << "dsygv could not diagonalize matrix: info = "
199 << info << endl;
200 abort();
201 }
202 double *z = fso; // the vector is placed in fso
203
204 // make sure everyone agrees on the new arrays
205 scf_grp_->bcast(z, nbasis2);
206 scf_grp_->bcast(epsilon, nbasis);
207
208 evalblock->assign(epsilon);
209 oso_scf_vector_block->assign(z);
210
211 // cleanup
212 delete[] fso;
213 delete[] work;
214 delete[] sso;
215 delete[] epsilon;
216 }
217
218 oso_scf_vector_ = (oso_scf_vector_ * so_to_orthog_so_inverse()).t();
219 #else
220 // diagonalize effective MO fock to get MO vector
221 tim_enter("evals");
222 RefSCMatrix nvector(oso_dimension(),oso_dimension(),basis_matrixkit());
223
224 0 0 0 0 1 0 0 0 0 0 RefSymmSCMatrix eff = effective_fock();
225
226 // level shift effective fock
227 level_shift->set_shift(level_shift_);
228 eff.element_op(level_shift);
229
230 if (debug_>1) {
231 eff.print("effective 1 body hamiltonian in current mo basis");
232 }
233
234 // transform eff to the oso basis to diagonalize it
235 RefSymmSCMatrix oso_eff(oso_dimension(), basis_matrixkit());
236 oso_eff.assign(0.0);
237 oso_eff.accumulate_transform(oso_scf_vector_,eff);
238 eff = 0;
239 oso_eff.diagonalize(evals, oso_scf_vector_);
240
241 tim_exit("evals");
242
243 1 0 0 0 0 0 0 0 0 0 if (debug_>0 && level_shift_ != 0.0) {
244 evals.print("level shifted scf eigenvalues");
245 }
246
247 // now un-level shift eigenvalues
248 level_shift->set_shift(-level_shift_);
249 evals.element_op(level_shift);
250 #endif
251
252 if (debug_>0) {
253 evals.print("scf eigenvalues");
254 }
255
256 if (reset_occ_)
257 set_occupations(evals);
258
259 if (debug_>1) {
260 oso_scf_vector_.print("OSO basis scf vector");
261 (oso_scf_vector_.t()*oso_scf_vector_).print(
262 "vOSO.t()*vOSO",ExEnv::out(),14);
263 }
264 }
265
266 eigenvalues_ = evals;
267 eigenvalues_.computed() = 1;
268 eigenvalues_.set_actual_accuracy(accuracy<delta?delta:accuracy);
269
270 // search for HOMO and LUMO
271 // first convert evals to something we can deal with easily
272 BlockedDiagSCMatrix *evalsb = require_dynamic_cast<BlockedDiagSCMatrix*>(evals,
273 "SCF::compute_vector");
274
275 CharacterTable ct = molecule()->point_group()->char_table();
276
277 int homo_ir=0, lumo_ir=0;
278 int homo_mo = -1, lumo_mo = -1;
279 double homo=-1e99, lumo=1e99;
280 for (i=0; i < oso_dimension()->blocks()->nblock(); i++) {
281 int nf=oso_dimension()->blocks()->size(i);
282 if (nf) {
283 double *vals = new double[nf];
284 evalsb->block(i)->convert(vals);
285
286 for (int mo=0; mo < nf; mo++) {
287 if (occupation(i, mo) > 0.0) {
288 if (vals[mo] > homo) {
289 homo = vals[mo];
290 homo_ir = i;
291 homo_mo = mo;
292 }
293 } else {
294 if (vals[mo] < lumo) {
295 lumo = vals[mo];
296 lumo_ir = i;
297 lumo_mo = mo;
298 }
299 }
300 }
301
302 if (print_all_evals_ || print_occ_evals_) {
303 ExEnv::out() << node0 << endl
304 << indent << ct.gamma(i).symbol() << endl << incindent;
305 for (int m=0; m < nf; m++) {
306 if (occupation(i,m) < 1e-8 && !print_all_evals_)
307 break;
308 ExEnv::out() << node0 << indent
309 << scprintf("%5d %10.5f %10.5f", m+1, vals[m], occupation(i,m))
310 << endl;
311 }
312 ExEnv::out() << node0 << decindent;
313 }
314
315 delete[] vals;
316 }
317 }
318
319 if (homo_mo >= 0) {
320 ExEnv::out() << node0 << endl << indent
321 << scprintf("HOMO is %5d %3s = %10.6f",
322 homo_mo+1,
323 ct.gamma(homo_ir).symbol(),
324 homo)
325 << endl;
326 }
327 if (lumo_mo >= 0) {
328 ExEnv::out() << node0 << indent
329 << scprintf("LUMO is %5d %3s = %10.6f",
330 lumo_mo+1,
331 ct.gamma(lumo_ir).symbol(),
332 lumo)
333 << endl;
334 }
335
336 // free up evals
337 evals = 0;
338
339 oso_eigenvectors_ = oso_scf_vector_;
340 oso_eigenvectors_.computed() = 1;
341 oso_eigenvectors_.set_actual_accuracy(delta);
342
343 // now clean up
344 done_vector();
345 hcore_ = 0;
346
347 level_shift->dereference();
348 delete level_shift;
349
350 tim_exit("vector");
351 //tim_print(0);
352
353 return delta;
354 }
355
356 ////////////////////////////////////////////////////////////////////////////
357
358 class ExtrapErrorOp : public BlockedSCElementOp {
359 private:
360 SCF *scf_;
361
362 public:
363 ExtrapErrorOp(SCF *s) : scf_(s) {}
364 ~ExtrapErrorOp() {}
365
366 int has_side_effects() { return 1; }
367
368 void process(SCMatrixBlockIter& bi) {
369 int ir=current_block();
370
371 for (bi.reset(); bi; bi++) {
372 7 0 2 0 4 0 0 0 0 0 int i=bi.i();
373 3 0 1 0 0 0 0 0 0 0 int j=bi.j();
374 10 0 9 2 7 2 0 0 0 0 if (scf_->occupation(ir,i) == scf_->occupation(ir,j))
375 2 0 3 0 3 0 0 0 0 0 bi.set(0.0);
376 }
377 }
378 };
379
380 Ref<SCExtrapError>
381 SCF::extrap_error()
382 {
383 RefSymmSCMatrix mofock = effective_fock();
384
385 Ref<SCElementOp> op = new ExtrapErrorOp(this);
386 mofock.element_op(op);
387
388 RefSymmSCMatrix aoerror(so_dimension(), basis_matrixkit());
389 aoerror.assign(0.0);
390 aoerror.accumulate_transform(so_to_orthog_so().t()*oso_scf_vector_, mofock);
391 0 0 0 0 0 0 1 0 0 0 mofock=0;
392
393 Ref<SCExtrapError> error = new SymmSCMatrixSCExtrapError(aoerror);
394 aoerror=0;
395
396 return error;
397 }
398
399 /////////////////////////////////////////////////////////////////////////////
400
401 // Local Variables:
402 // mode: c++
403 // c-file-style: "ETS"
404 // End:
405