1 //
2 // tcscf.cc --- implementation of the two-configuration 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 <math/scmat/block.h>
39 #include <math/scmat/blocked.h>
40 #include <math/scmat/blkiter.h>
41 #include <math/scmat/local.h>
42
43 #include <math/optimize/scextrapmat.h>
44
45 #include <chemistry/qc/basis/petite.h>
46
47 #include <chemistry/qc/scf/scflocal.h>
48 #include <chemistry/qc/scf/scfops.h>
49 #include <chemistry/qc/scf/effh.h>
50
51 #include <chemistry/qc/scf/tcscf.h>
52
53 using namespace std;
54 using namespace sc;
55
56 ///////////////////////////////////////////////////////////////////////////
57 // TCSCF
58
59 static ClassDesc TCSCF_cd(
60 typeid(TCSCF),"TCSCF",1,"public SCF",
61 0, 0, 0);
62
TCSCF(StateIn & s)63 TCSCF::TCSCF(StateIn& s) :
64 SavableState(s),
65 SCF(s),
66 focka_(this),
67 fockb_(this),
68 ka_(this),
69 kb_(this)
70 {
71 focka_.result_noupdate() = basis_matrixkit()->symmmatrix(so_dimension());
72 focka_.restore_state(s);
73 focka_.result_noupdate().restore(s);
74
75 fockb_.result_noupdate() = basis_matrixkit()->symmmatrix(so_dimension());
76 fockb_.restore_state(s);
77 fockb_.result_noupdate().restore(s);
78
79 ka_.result_noupdate() = basis_matrixkit()->symmmatrix(so_dimension());
80 ka_.restore_state(s);
81 ka_.result_noupdate().restore(s);
82
83 kb_.result_noupdate() = basis_matrixkit()->symmmatrix(so_dimension());
84 kb_.restore_state(s);
85 kb_.result_noupdate().restore(s);
86
87 s.get(user_occupations_);
88 s.get(tndocc_);
89 s.get(nirrep_);
90 s.get(ndocc_);
91 s.get(osa_);
92 s.get(osb_);
93 s.get(occa_);
94 s.get(occb_);
95 s.get(ci1_);
96 s.get(ci2_);
97
98 // now take care of memory stuff
99 init_mem(8);
100 }
101
TCSCF(const Ref<KeyVal> & keyval)102 TCSCF::TCSCF(const Ref<KeyVal>& keyval) :
103 SCF(keyval),
104 focka_(this),
105 fockb_(this),
106 ka_(this),
107 kb_(this)
108 {
109 focka_.compute()=0;
110 focka_.computed()=0;
111
112 fockb_.compute()=0;
113 fockb_.computed()=0;
114
115 ka_.compute()=0;
116 ka_.computed()=0;
117
118 kb_.compute()=0;
119 kb_.computed()=0;
120
121 // calculate the total nuclear charge
122 double Znuc=molecule()->nuclear_charge();
123
124 // check to see if this is to be a charged molecule
125 double charge = keyval->doublevalue("total_charge");
126 int nelectrons = (int)(Znuc-charge+1.0e-4);
127
128 // figure out how many doubly occupied shells there are
129 if (keyval->exists("ndocc")) {
130 tndocc_ = keyval->intvalue("ndocc");
131 } else {
132 tndocc_ = (nelectrons-2)/2;
133 if ((nelectrons-2)%2) {
134 ExEnv::err0() << endl << indent
135 << "TCSCF::init: Warning, there's a leftover electron.\n"
136 << incindent
137 << indent << "total_charge = " << charge << endl
138 << indent << "total nuclear charge = " << Znuc << endl
139 << indent << "ndocc_ = " << tndocc_ << endl << decindent;
140 }
141 }
142
143 ExEnv::out0() << endl << indent << "TCSCF::init: total charge = "
144 << Znuc-2*tndocc_-2 << endl << endl;
145
146 nirrep_ = molecule()->point_group()->char_table().ncomp();
147
148 if (nirrep_==1) {
149 ExEnv::err0() << indent << "TCSCF::init: cannot do C1 symmetry\n";
150 abort();
151 }
152
153 occa_=occb_=1.0;
154 ci1_=ci2_ = 0.5*sqrt(2.0);
155
156 if (keyval->exists("ci1")) {
157 ci1_ = keyval->doublevalue("ci1");
158 ci2_ = sqrt(1.0 - ci1_*ci1_);
159 occa_ = 2.0*ci1_*ci1_;
160 occb_ = 2.0*ci2_*ci2_;
161 }
162
163 if (keyval->exists("occa")) {
164 occa_ = keyval->doublevalue("occa");
165 ci1_ = sqrt(occa_/2.0);
166 ci2_ = sqrt(1.0 - ci1_*ci1_);
167 occb_ = 2.0*ci2_*ci2_;
168 }
169
170 osa_=-1;
171 osb_=-1;
172
173 ndocc_ = read_occ(keyval, "docc", nirrep_);
174 int *nsocc = read_occ(keyval, "socc", nirrep_);
175 if (ndocc_ && nsocc) {
176 user_occupations_=1;
177 for (int i=0; i < nirrep_; i++) {
178 int nsi = nsocc[i];
179 if (nsi && osa_<0)
180 osa_=i;
181 else if (nsi && osb_<0)
182 osb_=i;
183 else if (nsi) {
184 ExEnv::err0() << indent << "TCSCF::init: too many open shells\n";
185 abort();
186 }
187 }
188 delete[] nsocc;
189 }
190 else if (ndocc_ && !nsocc || !ndocc_ && nsocc) {
191 ExEnv::outn() << "ERROR: TCSCF: only one of docc and socc specified: "
192 << "give both or none" << endl;
193 abort();
194 }
195 else {
196 ndocc_=0;
197 user_occupations_=0;
198 set_occupations(0);
199 }
200
201 int i;
202 ExEnv::out0() << indent << "docc = [";
203 for (i=0; i < nirrep_; i++)
204 ExEnv::out0() << " " << ndocc_[i];
205 ExEnv::out0() << " ]\n";
206
207 ExEnv::out0() << indent << "socc = [";
208 for (i=0; i < nirrep_; i++)
209 ExEnv::out0() << " " << (i==osa_ || i==osb_) ? 1 : 0;
210 ExEnv::out0() << " ]\n";
211
212 // check to see if this was done in SCF(keyval)
213 if (!keyval->exists("maxiter"))
214 maxiter_ = 200;
215
216 if (!keyval->exists("level_shift"))
217 level_shift_ = 0.25;
218
219 // now take care of memory stuff
220 init_mem(8);
221 }
222
~TCSCF()223 TCSCF::~TCSCF()
224 {
225 if (ndocc_) {
226 delete[] ndocc_;
227 ndocc_=0;
228 }
229 }
230
231 void
save_data_state(StateOut & s)232 TCSCF::save_data_state(StateOut& s)
233 {
234 SCF::save_data_state(s);
235
236 focka_.save_data_state(s);
237 focka_.result_noupdate().save(s);
238
239 fockb_.save_data_state(s);
240 fockb_.result_noupdate().save(s);
241
242 ka_.save_data_state(s);
243 ka_.result_noupdate().save(s);
244
245 kb_.save_data_state(s);
246 kb_.result_noupdate().save(s);
247
248 s.put(user_occupations_);
249 s.put(tndocc_);
250 s.put(nirrep_);
251 s.put(ndocc_,nirrep_);
252 s.put(osa_);
253 s.put(osb_);
254 s.put(occa_);
255 s.put(occb_);
256 s.put(ci1_);
257 s.put(ci2_);
258 }
259
260 double
occupation(int ir,int i)261 TCSCF::occupation(int ir, int i)
262 {
263 if (i < ndocc_[ir])
264 return 2.0;
265 else if (ir==osa_ && i==ndocc_[ir])
266 return occa_;
267 else if (ir==osb_ && i==ndocc_[ir])
268 return occb_;
269 else
270 return 0.0;
271 }
272
273 double
alpha_occupation(int ir,int i)274 TCSCF::alpha_occupation(int ir, int i)
275 {
276 if (i < ndocc_[ir])
277 return 1.0;
278 else if (ir==osa_ && i==ndocc_[ir])
279 return 0.5*occa_;
280 return 0.0;
281 }
282
283 double
beta_occupation(int ir,int i)284 TCSCF::beta_occupation(int ir, int i)
285 {
286 if (i < ndocc_[ir])
287 return 1.0;
288 else if (ir==osb_ && i==ndocc_[ir])
289 return 0.5*occb_;
290 return 0.0;
291 }
292
293 int
n_fock_matrices() const294 TCSCF::n_fock_matrices() const
295 {
296 return 4;
297 }
298
299 RefSymmSCMatrix
fock(int n)300 TCSCF::fock(int n)
301 {
302 if (n > 3) {
303 ExEnv::err0() << indent
304 << "TCSCF::fock: there are only four fock matrices, "
305 << scprintf("but fock(%d) was requested\n", n);
306 abort();
307 }
308
309 if (n==0)
310 return focka_.result();
311 else if (n==1)
312 return fockb_.result();
313 else if (n==2)
314 return ka_.result();
315 else
316 return kb_.result();
317 }
318
319 int
spin_polarized()320 TCSCF::spin_polarized()
321 {
322 return 1;
323 }
324
325 void
print(ostream & o) const326 TCSCF::print(ostream&o) const
327 {
328 int i;
329
330 SCF::print(o);
331
332 o << indent << "TCSCF Parameters:\n" << incindent
333 << indent << "ndocc = " << tndocc_ << endl
334 << indent << scprintf("occa = %f", occa_) << endl
335 << indent << scprintf("occb = %f", occb_) << endl
336 << indent << scprintf("ci1 = %9.6f", ci1_) << endl
337 << indent << scprintf("ci2 = %9.6f", ci2_) << endl
338 << indent << "docc = [";
339 for (i=0; i < nirrep_; i++)
340 o << " " << ndocc_[i];
341 o << " ]" << endl
342 << indent << "socc = [";
343 for (i=0; i < nirrep_; i++)
344 o << " " << (i==osa_ || i==osb_) ? 1 : 0;
345 o << " ]" << endl << decindent << endl;
346 }
347
348 //////////////////////////////////////////////////////////////////////////////
349
350 void
set_occupations(const RefDiagSCMatrix & ev)351 TCSCF::set_occupations(const RefDiagSCMatrix& ev)
352 {
353 if (user_occupations_)
354 return;
355
356 int i,j;
357
358 RefDiagSCMatrix evals;
359
360 if (ev.null()) {
361 initial_vector(0);
362 evals = eigenvalues_.result_noupdate();
363 }
364 else
365 evals = ev;
366
367 // first convert evals to something we can deal with easily
368 BlockedDiagSCMatrix *evalsb = require_dynamic_cast<BlockedDiagSCMatrix*>(evals,
369 "TCSCF::set_occupations");
370
371 double **vals = new double*[nirrep_];
372 for (i=0; i < nirrep_; i++) {
373 int nf=oso_dimension()->blocks()->size(i);
374 if (nf) {
375 vals[i] = new double[nf];
376 evalsb->block(i)->convert(vals[i]);
377 } else {
378 vals[i] = 0;
379 }
380 }
381
382 // now loop to find the tndocc_ lowest eigenvalues and populate those
383 // MO's
384 int *newdocc = new int[nirrep_];
385 memset(newdocc,0,sizeof(int)*nirrep_);
386
387 for (i=0; i < tndocc_; i++) {
388 // find lowest eigenvalue
389 int lir=0,ln=0;
390 double lowest=999999999;
391
392 for (int ir=0; ir < nirrep_; ir++) {
393 int nf=oso_dimension()->blocks()->size(ir);
394 if (!nf)
395 continue;
396 for (j=0; j < nf; j++) {
397 if (vals[ir][j] < lowest) {
398 lowest=vals[ir][j];
399 lir=ir;
400 ln=j;
401 }
402 }
403 }
404 vals[lir][ln]=999999999;
405 newdocc[lir]++;
406 }
407
408 int osa=-1, osb=-1;
409
410 for (i=0; i < 2; i++) {
411 // find lowest eigenvalue
412 int lir=0,ln=0;
413 double lowest=999999999;
414
415 for (int ir=0; ir < nirrep_; ir++) {
416 int nf=oso_dimension()->blocks()->size(ir);
417 if (!nf)
418 continue;
419 for (j=0; j < nf; j++) {
420 if (vals[ir][j] < lowest) {
421 lowest=vals[ir][j];
422 lir=ir;
423 ln=j;
424 }
425 }
426 }
427 vals[lir][ln]=999999999;
428
429 if (!i) {
430 osa=lir;
431 } else {
432 if (lir==osa) {
433 i--;
434 continue;
435 }
436 osb=lir;
437 }
438 }
439
440 if (osa > osb) {
441 int tmp=osa;
442 osa=osb;
443 osb=tmp;
444 }
445
446 // get rid of vals
447 for (i=0; i < nirrep_; i++)
448 if (vals[i])
449 delete[] vals[i];
450 delete[] vals;
451
452 if (!ndocc_) {
453 ndocc_=newdocc;
454 osa_=osa;
455 osb_=osb;
456 } else {
457 // test to see if newocc is different from ndocc_
458 for (i=0; i < nirrep_; i++) {
459 if (ndocc_[i] != newdocc[i]) {
460 ExEnv::err0() << indent << "TCSCF::set_occupations: WARNING!!!!\n"
461 << incindent << indent
462 << scprintf("occupations for irrep %d have changed\n", i+1)
463 << indent
464 << scprintf("ndocc was %d, changed to %d", ndocc_[i], newdocc[i])
465 << endl << decindent;
466 }
467 if (((osa != osa_ && osa != osb_) || (osb != osb_ && osb != osa_))) {
468 ExEnv::err0() << indent << "TCSCF::set_occupations: WARNING!!!!\n"
469 << incindent << indent << "open shell occupations have changed"
470 << endl << decindent;
471 osa_=osa;
472 osb_=osb;
473 reset_density();
474 }
475 }
476
477 memcpy(ndocc_,newdocc,sizeof(int)*nirrep_);
478
479 delete[] newdocc;
480 }
481 }
482
483 void
symmetry_changed()484 TCSCF::symmetry_changed()
485 {
486 SCF::symmetry_changed();
487 focka_.result_noupdate()=0;
488 fockb_.result_noupdate()=0;
489 ka_.result_noupdate()=0;
490 kb_.result_noupdate()=0;
491 nirrep_ = molecule()->point_group()->char_table().ncomp();
492 set_occupations(0);
493 }
494
495 //////////////////////////////////////////////////////////////////////////////
496 //
497 // scf things
498 //
499
500 void
init_vector()501 TCSCF::init_vector()
502 {
503 init_threads();
504
505 // allocate storage for other temp matrices
506 cl_dens_ = hcore_.clone();
507 cl_dens_.assign(0.0);
508
509 cl_dens_diff_ = hcore_.clone();
510 cl_dens_diff_.assign(0.0);
511
512 op_densa_ = hcore_.clone();
513 op_densa_.assign(0.0);
514
515 op_densa_diff_ = hcore_.clone();
516 op_densa_diff_.assign(0.0);
517
518 op_densb_ = hcore_.clone();
519 op_densb_.assign(0.0);
520
521 op_densb_diff_ = hcore_.clone();
522 op_densb_diff_.assign(0.0);
523
524 // gmat is in AO basis
525 ao_gmata_ = basis()->matrixkit()->symmmatrix(basis()->basisdim());
526 ao_gmata_.assign(0.0);
527
528 ao_gmatb_ = ao_gmata_.clone();
529 ao_gmatb_.assign(0.0);
530
531 ao_ka_ = ao_gmata_.clone();
532 ao_ka_.assign(0.0);
533
534 ao_kb_ = ao_gmata_.clone();
535 ao_kb_.assign(0.0);
536
537 // test to see if we need a guess vector
538 if (focka_.result_noupdate().null()) {
539 focka_ = hcore_.clone();
540 focka_.result_noupdate().assign(0.0);
541 fockb_ = hcore_.clone();
542 fockb_.result_noupdate().assign(0.0);
543 ka_ = hcore_.clone();
544 ka_.result_noupdate().assign(0.0);
545 kb_ = hcore_.clone();
546 kb_.result_noupdate().assign(0.0);
547 }
548
549 // set up trial vector
550 initial_vector(1);
551
552 oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
553 }
554
555 void
done_vector()556 TCSCF::done_vector()
557 {
558 done_threads();
559
560 cl_dens_ = 0;
561 cl_dens_diff_ = 0;
562 op_densa_ = 0;
563 op_densa_diff_ = 0;
564 op_densb_ = 0;
565 op_densb_diff_ = 0;
566
567 ao_gmata_ = 0;
568 ao_gmatb_ = 0;
569 ao_ka_ = 0;
570 ao_kb_ = 0;
571
572 oso_scf_vector_ = 0;
573 }
574
575 ////////////////////////////////////////////////////////////////////////////
576
577 RefSymmSCMatrix
density()578 TCSCF::density()
579 {
580 if (!density_.computed()) {
581 RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
582 RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());
583
584 so_density(dens, 2.0);
585 dens.scale(2.0);
586
587 so_density(dens1, occa_);
588 dens1.scale(occa_);
589 dens.accumulate(dens1);
590
591 so_density(dens1, occb_);
592 dens1.scale(occb_);
593 dens.accumulate(dens1);
594
595 dens1=0;
596
597 density_ = dens;
598 // only flag the density as computed if the calc is converged
599 if (!value_needed()) density_.computed() = 1;
600 }
601
602 return density_.result_noupdate();
603 }
604
605 RefSymmSCMatrix
alpha_density()606 TCSCF::alpha_density()
607 {
608 RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
609 RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());
610
611 so_density(dens, 2.0);
612 so_density(dens1, occa_);
613 dens.accumulate(dens1);
614
615 dens.scale(2.0);
616 return dens;
617 }
618
619 RefSymmSCMatrix
beta_density()620 TCSCF::beta_density()
621 {
622 RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
623 RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());
624
625 so_density(dens, 2.0);
626 so_density(dens1, occb_);
627 dens.accumulate(dens1);
628
629 dens.scale(2.0);
630 return dens;
631 }
632
633 void
reset_density()634 TCSCF::reset_density()
635 {
636 cl_dens_diff_.assign(cl_dens_);
637
638 ao_gmata_.assign(0.0);
639 op_densa_diff_.assign(op_densa_);
640
641 ao_gmatb_.assign(0.0);
642 op_densb_diff_.assign(op_densb_);
643
644 ao_ka_.assign(0.0);
645 ao_kb_.assign(0.0);
646 }
647
648 double
new_density()649 TCSCF::new_density()
650 {
651 // copy current density into density diff and scale by -1. later we'll
652 // add the new density to this to get the density difference.
653 cl_dens_diff_.assign(cl_dens_);
654 cl_dens_diff_.scale(-1.0);
655
656 op_densa_diff_.assign(op_densa_);
657 op_densa_diff_.scale(-1.0);
658
659 op_densb_diff_.assign(op_densb_);
660 op_densb_diff_.scale(-1.0);
661
662 so_density(cl_dens_, 2.0);
663 cl_dens_.scale(2.0);
664
665 so_density(op_densa_, occa_);
666 dynamic_cast<BlockedSymmSCMatrix*>(op_densa_.pointer())->block(osb_)->assign(0.0);
667 op_densa_.scale(2.0);
668
669 so_density(op_densb_, occb_);
670 dynamic_cast<BlockedSymmSCMatrix*>(op_densb_.pointer())->block(osa_)->assign(0.0);
671 op_densb_.scale(2.0);
672
673 cl_dens_diff_.accumulate(cl_dens_);
674 op_densa_diff_.accumulate(op_densa_);
675 op_densb_diff_.accumulate(op_densb_);
676
677 RefSymmSCMatrix del = cl_dens_diff_.copy();
678 del.accumulate(op_densa_diff_);
679 del.accumulate(op_densb_diff_);
680
681 Ref<SCElementScalarProduct> sp(new SCElementScalarProduct);
682 del.element_op(sp.pointer(), del);
683
684 double delta = sp->result();
685 delta = sqrt(delta/i_offset(cl_dens_diff_.n()));
686
687 return delta;
688 }
689
690 double
scf_energy()691 TCSCF::scf_energy()
692 {
693 // first calculate the elements of the CI matrix
694 SCFEnergy *eop = new SCFEnergy;
695 eop->reference();
696 Ref<SCElementOp2> op = eop;
697
698 RefSymmSCMatrix t = focka_.result_noupdate().copy();
699 t.accumulate(hcore_);
700
701 RefSymmSCMatrix d = cl_dens_.copy();
702 d.accumulate(op_densa_);
703
704 t.element_op(op, d);
705 double h11 = eop->result();
706
707 t.assign(fockb_.result_noupdate().copy());
708 t.accumulate(hcore_);
709
710 d.assign(cl_dens_);
711 d.accumulate(op_densb_);
712
713 eop->reset();
714 t.element_op(op, d);
715 double h22 = eop->result();
716
717 //t = ka_.result_noupdate();
718 //eop->reset();
719 //t.element_op(op, op_densb_);
720 //double h21 = eop->result();
721
722 t = kb_.result_noupdate();
723 eop->reset();
724 t.element_op(op, op_densa_);
725 double h12 = eop->result();
726
727 op=0;
728 eop->dereference();
729 delete eop;
730
731 // now diagonalize the CI matrix to get the coefficients
732 RefSCDimension l2 = new SCDimension(2);
733 Ref<SCMatrixKit> lkit = new LocalSCMatrixKit;
734 RefSymmSCMatrix h = lkit->symmmatrix(l2);
735 RefSCMatrix hv = lkit->matrix(l2,l2);
736 RefDiagSCMatrix hl = lkit->diagmatrix(l2);
737
738 h.set_element(0,0,h11);
739 h.set_element(1,1,h22);
740 h.set_element(1,0,h12);
741 h.diagonalize(hl,hv);
742
743 ci1_ = hv.get_element(0,0);
744 ci2_ = hv.get_element(1,0);
745 double c1c2 = ci1_*ci2_;
746
747 ExEnv::out0() << indent
748 << scprintf("c1 = %10.7f c2 = %10.7f", ci1_, ci2_)
749 << endl;
750
751 occa_ = 2*ci1_*ci1_;
752 occb_ = 2*ci2_*ci2_;
753
754 double eelec = 0.5*occa_*h11 + 0.5*occb_*h22 + 2.0*c1c2*h12;
755
756 return eelec;
757 }
758
759 Ref<SCExtrapData>
extrap_data()760 TCSCF::extrap_data()
761 {
762 RefSymmSCMatrix *m = new RefSymmSCMatrix[4];
763 m[0] = focka_.result_noupdate();
764 m[1] = fockb_.result_noupdate();
765 m[2] = ka_.result_noupdate();
766 m[3] = kb_.result_noupdate();
767
768 Ref<SCExtrapData> data = new SymmSCMatrixNSCExtrapData(4, m);
769 delete[] m;
770
771 return data;
772 }
773
774 RefSymmSCMatrix
effective_fock()775 TCSCF::effective_fock()
776 {
777 // use fock() instead of cl_fock_ just in case this is called from
778 // someplace outside SCF::compute_vector()
779 RefSymmSCMatrix mofocka(oso_dimension(), basis_matrixkit());
780 mofocka.assign(0.0);
781
782 RefSymmSCMatrix mofockb(oso_dimension(), basis_matrixkit());
783 mofockb.assign(0.0);
784
785 RefSymmSCMatrix moka = mofocka.clone();
786 moka.assign(0.0);
787
788 RefSymmSCMatrix mokb = mofocka.clone();
789 mokb.assign(0.0);
790
791 // use eigenvectors if oso_scf_vector_ is null
792 RefSCMatrix vec;
793 if (oso_scf_vector_.null()) {
794 vec = eigenvectors();
795 } else {
796 vec = so_to_orthog_so().t() * oso_scf_vector_;
797 }
798 mofocka.accumulate_transform(vec, fock(0),
799 SCMatrix::TransposeTransform);
800 mofockb.accumulate_transform(vec, fock(1),
801 SCMatrix::TransposeTransform);
802 moka.accumulate_transform(vec, fock(2),
803 SCMatrix::TransposeTransform);
804 mokb.accumulate_transform(vec, fock(3),
805 SCMatrix::TransposeTransform);
806
807 mofocka.scale(ci1_*ci1_);
808 mofockb.scale(ci2_*ci2_);
809 moka.scale(ci1_*ci2_);
810 mokb.scale(ci1_*ci2_);
811
812 RefSymmSCMatrix mofock = mofocka.copy();
813 mofock.accumulate(mofockb);
814
815 BlockedSymmSCMatrix *F = dynamic_cast<BlockedSymmSCMatrix*>(mofock.pointer());
816 BlockedSymmSCMatrix *Fa = dynamic_cast<BlockedSymmSCMatrix*>(mofocka.pointer());
817 BlockedSymmSCMatrix *Fb = dynamic_cast<BlockedSymmSCMatrix*>(mofockb.pointer());
818 BlockedSymmSCMatrix *Ka = dynamic_cast<BlockedSymmSCMatrix*>(moka.pointer());
819 BlockedSymmSCMatrix *Kb = dynamic_cast<BlockedSymmSCMatrix*>(mokb.pointer());
820
821 double scalea = (fabs(ci1_) < fabs(ci2_)) ? 1.0/(ci1_*ci1_ + 0.05) : 1.0;
822 double scaleb = (fabs(ci2_) < fabs(ci1_)) ? 1.0/(ci2_*ci2_ + 0.05) : 1.0;
823
824 for (int b=0; b < Fa->nblocks(); b++) {
825 if (b==osa_) {
826 RefSymmSCMatrix f = F->block(b);
827 RefSymmSCMatrix fa = Fa->block(b);
828 RefSymmSCMatrix fb = Fb->block(b);
829 RefSymmSCMatrix kb = Kb->block(b);
830
831 int i,j;
832
833 i=ndocc_[b];
834 for (j=0; j < ndocc_[b]; j++)
835 f->set_element(i,j,
836 scaleb*(fb->get_element(i,j)-kb->get_element(i,j)));
837
838 j=ndocc_[b];
839 for (i=ndocc_[b]+1; i < f->n(); i++)
840 f->set_element(i,j,
841 scalea*(fa->get_element(i,j)+kb->get_element(i,j)));
842
843 } else if (b==osb_) {
844 RefSymmSCMatrix f = F->block(b);
845 RefSymmSCMatrix fa = Fa->block(b);
846 RefSymmSCMatrix fb = Fb->block(b);
847 RefSymmSCMatrix ka = Ka->block(b);
848
849 int i,j;
850
851 i=ndocc_[b];
852 for (j=0; j < ndocc_[b]; j++)
853 f->set_element(i,j,
854 scalea*(fa->get_element(i,j)-ka->get_element(i,j)));
855
856 j=ndocc_[b];
857 for (i=ndocc_[b]+1; i < f->n(); i++)
858 f->set_element(i,j,
859 scaleb*(fb->get_element(i,j)+ka->get_element(i,j)));
860 }
861 }
862
863 return mofock;
864 }
865
866 /////////////////////////////////////////////////////////////////////////////
867
868 void
init_gradient()869 TCSCF::init_gradient()
870 {
871 // presumably the eigenvectors have already been computed by the time
872 // we get here
873 oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
874 }
875
876 void
done_gradient()877 TCSCF::done_gradient()
878 {
879 cl_dens_=0;
880 op_densa_=0;
881 op_densb_=0;
882 oso_scf_vector_ = 0;
883 }
884
885 /////////////////////////////////////////////////////////////////////////////
886
887 // MO lagrangian
888 // c o v
889 // c |2*FC|2*FC|0|
890 // -------------
891 // o |2*FC| FO |0|
892 // -------------
893 // v | 0 | 0 |0|
894 //
895 RefSymmSCMatrix
lagrangian()896 TCSCF::lagrangian()
897 {
898 RefSCMatrix vec = so_to_orthog_so().t() * oso_scf_vector_;
899
900 RefSymmSCMatrix mofocka = focka_.result_noupdate().clone();
901 mofocka.assign(0.0);
902 mofocka.accumulate_transform(vec, focka_.result_noupdate(),
903 SCMatrix::TransposeTransform);
904 mofocka.scale(ci1_*ci1_);
905
906 RefSymmSCMatrix mofockb = mofocka.clone();
907 mofockb.assign(0.0);
908 mofockb.accumulate_transform(vec, fockb_.result_noupdate(),
909 SCMatrix::TransposeTransform);
910 mofockb.scale(ci2_*ci2_);
911
912 // FOa = c1^2*Fa + c1c2*Kb
913 RefSymmSCMatrix moka = mofocka.clone();
914 moka.assign(0.0);
915 moka.accumulate_transform(vec, kb_.result_noupdate(),
916 SCMatrix::TransposeTransform);
917 moka.scale(ci1_*ci2_);
918 moka.accumulate(mofocka);
919
920 // FOb = c1^2*Fb + c1c2*Ka
921 RefSymmSCMatrix mokb = mofocka.clone();
922 mokb.assign(0.0);
923 mokb.accumulate_transform(vec, ka_.result_noupdate(),
924 SCMatrix::TransposeTransform);
925 mokb.scale(ci1_*ci2_);
926 mokb.accumulate(mofockb);
927
928 dynamic_cast<BlockedSymmSCMatrix*>(moka.pointer())->block(osb_)->assign(0.0);
929 dynamic_cast<BlockedSymmSCMatrix*>(mokb.pointer())->block(osa_)->assign(0.0);
930
931 moka.accumulate(mokb);
932 mokb=0;
933
934 // FC = c1^2*Fa + c2^2*Fb
935 mofocka.accumulate(mofockb);
936 mofockb=0;
937
938 Ref<SCElementOp2> op = new MOLagrangian(this);
939 mofocka.element_op(op, moka);
940 moka=0;
941 mofocka.scale(2.0);
942
943 // transform MO lagrangian to SO basis
944 RefSymmSCMatrix so_lag(so_dimension(), basis_matrixkit());
945 so_lag.assign(0.0);
946 so_lag.accumulate_transform(vec, mofocka);
947
948 // and then from SO to AO
949 Ref<PetiteList> pl = integral()->petite_list();
950 RefSymmSCMatrix ao_lag = pl->to_AO_basis(so_lag);
951
952 ao_lag.scale(-1.0);
953
954 return ao_lag;
955 }
956
957 RefSymmSCMatrix
gradient_density()958 TCSCF::gradient_density()
959 {
960 cl_dens_ = basis_matrixkit()->symmmatrix(so_dimension());
961 op_densa_ = cl_dens_.clone();
962 op_densb_ = cl_dens_.clone();
963
964 so_density(cl_dens_, 2.0);
965 cl_dens_.scale(2.0);
966
967 so_density(op_densa_, occa_);
968 op_densa_.scale(occa_);
969
970 so_density(op_densb_, occb_);
971 op_densb_.scale(occb_);
972
973 dynamic_cast<BlockedSymmSCMatrix*>(op_densa_.pointer())->block(osb_)->assign(0.0);
974 dynamic_cast<BlockedSymmSCMatrix*>(op_densb_.pointer())->block(osa_)->assign(0.0);
975
976 Ref<PetiteList> pl = integral()->petite_list(basis());
977
978 cl_dens_ = pl->to_AO_basis(cl_dens_);
979 op_densa_ = pl->to_AO_basis(op_densa_);
980 op_densb_ = pl->to_AO_basis(op_densb_);
981
982 RefSymmSCMatrix tdens = cl_dens_.copy();
983 tdens.accumulate(op_densa_);
984 tdens.accumulate(op_densb_);
985
986 op_densa_.scale(2.0/occa_);
987 op_densb_.scale(2.0/occb_);
988
989 return tdens;
990 }
991
992 /////////////////////////////////////////////////////////////////////////////
993
994 void
init_hessian()995 TCSCF::init_hessian()
996 {
997 }
998
999 void
done_hessian()1000 TCSCF::done_hessian()
1001 {
1002 }
1003
1004 /////////////////////////////////////////////////////////////////////////////
1005
1006 // Local Variables:
1007 // mode: c++
1008 // c-file-style: "ETS"
1009 // End:
1010