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