1 //
2 // hsosscf.cc --- implementation of the high-spin open shell 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 #include <chemistry/qc/scf/hsosscf.h>
51
52 #include <chemistry/qc/scf/ltbgrad.h>
53 #include <chemistry/qc/scf/hsoshftmpl.h>
54
55 using namespace std;
56 using namespace sc;
57
58 ///////////////////////////////////////////////////////////////////////////
59 // HSOSSCF
60
61 static ClassDesc HSOSSCF_cd(
62 typeid(HSOSSCF),"HSOSSCF",2,"public SCF",
63 0, 0, 0);
64
HSOSSCF(StateIn & s)65 HSOSSCF::HSOSSCF(StateIn& s) :
66 SavableState(s),
67 SCF(s),
68 cl_fock_(this),
69 op_fock_(this)
70 {
71 cl_fock_.result_noupdate() =
72 basis_matrixkit()->symmmatrix(so_dimension());
73 cl_fock_.restore_state(s);
74 cl_fock_.result_noupdate().restore(s);
75
76 op_fock_.result_noupdate() =
77 basis_matrixkit()->symmmatrix(so_dimension());
78 op_fock_.restore_state(s);
79 op_fock_.result_noupdate().restore(s);
80
81 s.get(user_occupations_);
82 s.get(tndocc_);
83 s.get(tnsocc_);
84 s.get(nirrep_);
85 s.get(ndocc_);
86 s.get(nsocc_);
87 if (s.version(::class_desc<HSOSSCF>()) >= 2) {
88 s.get(initial_ndocc_);
89 s.get(initial_nsocc_);
90 most_recent_pg_ << SavableState::restore_state(s);
91 } else {
92 initial_ndocc_ = new int[nirrep_];
93 memcpy(initial_ndocc_, ndocc_, sizeof(int)*nirrep_);
94 initial_nsocc_ = new int[nirrep_];
95 memcpy(initial_nsocc_, nsocc_, sizeof(int)*nirrep_);
96 }
97
98 // now take care of memory stuff
99 init_mem(4);
100 }
101
HSOSSCF(const Ref<KeyVal> & keyval)102 HSOSSCF::HSOSSCF(const Ref<KeyVal>& keyval) :
103 SCF(keyval),
104 cl_fock_(this),
105 op_fock_(this)
106 {
107 int i;
108
109 cl_fock_.compute()=0;
110 cl_fock_.computed()=0;
111
112 op_fock_.compute()=0;
113 op_fock_.computed()=0;
114
115 // calculate the total nuclear charge
116 double Znuc=molecule()->nuclear_charge();
117
118 // check to see if this is to be a charged molecule
119 double charge = keyval->doublevalue("total_charge");
120 int nelectrons = (int)(Znuc-charge+1.0e-4);
121
122 // first let's try to figure out how many open shells there are
123 if (keyval->exists("nsocc")) {
124 tnsocc_ = keyval->intvalue("nsocc");
125 } else if (keyval->exists("multiplicity")) {
126 tnsocc_ = keyval->intvalue("multiplicity")-1;
127 } else {
128 // if there's an odd number of electrons, then do a doublet, otherwise
129 // do a triplet
130 if (nelectrons%2)
131 tnsocc_=1;
132 else
133 tnsocc_=2;
134 }
135
136 // now do the same for the number of doubly occupied shells
137 if (keyval->exists("ndocc")) {
138 tndocc_ = keyval->intvalue("ndocc");
139 } else {
140 tndocc_ = (nelectrons-tnsocc_)/2;
141 if ((nelectrons-tnsocc_)%2) {
142 ExEnv::err0() << endl << indent
143 << "HSOSSCF::init: Warning, there's a leftover electron.\n"
144 << incindent << indent << "total_charge = " << charge << endl
145 << indent << "total nuclear charge = " << Znuc << endl
146 << indent << "ndocc_ = " << tndocc_ << endl
147 << indent << "nsocc_ = " << tnsocc_ << endl << decindent;
148 }
149 }
150
151 ExEnv::out0() << endl << indent
152 << "HSOSSCF::init: total charge = " << Znuc-2*tndocc_-tnsocc_
153 << endl << endl;
154
155 nirrep_ = molecule()->point_group()->char_table().ncomp();
156
157 ndocc_ = read_occ(keyval, "docc", nirrep_);
158 nsocc_ = read_occ(keyval, "socc", nirrep_);
159 if (ndocc_ && nsocc_) {
160 user_occupations_=1;
161 initial_ndocc_ = new int[nirrep_];
162 memcpy(initial_ndocc_, ndocc_, sizeof(int)*nirrep_);
163 initial_nsocc_ = new int[nirrep_];
164 memcpy(initial_nsocc_, nsocc_, sizeof(int)*nirrep_);
165 }
166 else if (ndocc_ && !nsocc_ || !ndocc_ && nsocc_) {
167 ExEnv::outn() << "ERROR: HSOSSCF: only one of docc and socc specified: "
168 << "give both or none" << endl;
169 abort();
170 }
171 else {
172 ndocc_=0;
173 nsocc_=0;
174 initial_ndocc_=0;
175 initial_nsocc_=0;
176 user_occupations_=0;
177 set_occupations(0);
178 }
179
180 ExEnv::out0() << indent << "docc = [";
181 for (i=0; i < nirrep_; i++)
182 ExEnv::out0() << " " << ndocc_[i];
183 ExEnv::out0() << " ]\n";
184
185 ExEnv::out0() << indent << "socc = [";
186 for (i=0; i < nirrep_; i++)
187 ExEnv::out0() << " " << nsocc_[i];
188 ExEnv::out0() << " ]\n";
189
190 // check to see if this was done in SCF(keyval)
191 if (!keyval->exists("maxiter"))
192 maxiter_ = 100;
193
194 if (!keyval->exists("level_shift"))
195 level_shift_ = 0.25;
196
197 // now take care of memory stuff
198 init_mem(4);
199 }
200
~HSOSSCF()201 HSOSSCF::~HSOSSCF()
202 {
203 if (ndocc_) {
204 delete[] ndocc_;
205 ndocc_=0;
206 }
207 if (nsocc_) {
208 delete[] nsocc_;
209 nsocc_=0;
210 }
211 delete[] initial_ndocc_;
212 delete[] initial_nsocc_;
213 }
214
215 void
save_data_state(StateOut & s)216 HSOSSCF::save_data_state(StateOut& s)
217 {
218 SCF::save_data_state(s);
219
220 cl_fock_.save_data_state(s);
221 cl_fock_.result_noupdate().save(s);
222
223 op_fock_.save_data_state(s);
224 op_fock_.result_noupdate().save(s);
225
226 s.put(user_occupations_);
227 s.put(tndocc_);
228 s.put(tnsocc_);
229 s.put(nirrep_);
230 s.put(ndocc_,nirrep_);
231 s.put(nsocc_,nirrep_);
232 s.put(initial_ndocc_,nirrep_);
233 s.put(initial_nsocc_,nirrep_);
234 SavableState::save_state(most_recent_pg_.pointer(),s);
235 }
236
237 double
occupation(int ir,int i)238 HSOSSCF::occupation(int ir, int i)
239 {
240 if (i < ndocc_[ir]) return 2.0;
241 else if (i < ndocc_[ir] + nsocc_[ir]) return 1.0;
242 return 0.0;
243 }
244
245 double
alpha_occupation(int ir,int i)246 HSOSSCF::alpha_occupation(int ir, int i)
247 {
248 if (i < ndocc_[ir] + nsocc_[ir]) return 1.0;
249 return 0.0;
250 }
251
252 double
beta_occupation(int ir,int i)253 HSOSSCF::beta_occupation(int ir, int i)
254 {
255 if (i < ndocc_[ir]) return 1.0;
256 return 0.0;
257 }
258
259 int
n_fock_matrices() const260 HSOSSCF::n_fock_matrices() const
261 {
262 return 2;
263 }
264
265 RefSymmSCMatrix
fock(int n)266 HSOSSCF::fock(int n)
267 {
268 if (n > 1) {
269 ExEnv::err0() << indent
270 << "HSOSSCF::fock: there are only two fock matrices, "
271 << scprintf("but fock(%d) was requested\n",n);
272 abort();
273 }
274
275 if (n==0)
276 return cl_fock_.result();
277 else
278 return op_fock_.result();
279 }
280
281 int
spin_polarized()282 HSOSSCF::spin_polarized()
283 {
284 return 1;
285 }
286
287 void
print(ostream & o) const288 HSOSSCF::print(ostream&o) const
289 {
290 int i;
291
292 SCF::print(o);
293 o << indent << "HSOSSCF Parameters:\n" << incindent
294 << indent << "charge = " << molecule()->nuclear_charge()
295 - 2*tndocc_ - tnsocc_ << endl
296 << indent << "ndocc = " << tndocc_ << endl
297 << indent << "nsocc = " << tnsocc_ << endl
298 << indent << "docc = [";
299 for (i=0; i < nirrep_; i++)
300 o << " " << ndocc_[i];
301 o << " ]" << endl;
302
303 o << indent << "socc = [";
304 for (i=0; i < nirrep_; i++)
305 o << " " << nsocc_[i];
306 o << " ]" << endl << decindent << endl;
307 }
308
309 //////////////////////////////////////////////////////////////////////////////
310
311 void
set_occupations(const RefDiagSCMatrix & ev)312 HSOSSCF::set_occupations(const RefDiagSCMatrix& ev)
313 {
314 if (user_occupations_ || (initial_ndocc_ && initial_nsocc_ && ev.null())) {
315 if (form_occupations(ndocc_, initial_ndocc_)
316 &&form_occupations(nsocc_, initial_nsocc_)) {
317 most_recent_pg_ = new PointGroup(molecule()->point_group());
318 return;
319 }
320 delete[] ndocc_; ndocc_ = 0;
321 delete[] nsocc_; nsocc_ = 0;
322 ExEnv::out0() << indent
323 << "HSOSSCF: WARNING: reforming occupation vectors from scratch"
324 << endl;
325 }
326
327 if (nirrep_==1) {
328 delete[] ndocc_;
329 ndocc_=new int[1];
330 ndocc_[0]=tndocc_;
331 if (!initial_ndocc_) {
332 initial_ndocc_=new int[1];
333 initial_ndocc_[0]=tndocc_;
334 }
335
336 delete[] nsocc_;
337 nsocc_=new int[1];
338 nsocc_[0]=tnsocc_;
339 if (!initial_nsocc_) {
340 initial_nsocc_=new int[1];
341 initial_nsocc_[0]=tnsocc_;
342 }
343
344 return;
345 }
346
347 int i,j;
348
349 RefDiagSCMatrix evals;
350
351 if (ev.null()) {
352 initial_vector(0);
353 evals = eigenvalues_.result_noupdate();
354 }
355 else
356 evals = ev;
357
358 // first convert evals to something we can deal with easily
359 BlockedDiagSCMatrix *evalsb = require_dynamic_cast<BlockedDiagSCMatrix*>(evals,
360 "HSOSSCF::set_occupations");
361
362 double **vals = new double*[nirrep_];
363 for (i=0; i < nirrep_; i++) {
364 int nf=oso_dimension()->blocks()->size(i);
365 if (nf) {
366 vals[i] = new double[nf];
367 evalsb->block(i)->convert(vals[i]);
368 } else {
369 vals[i] = 0;
370 }
371 }
372
373 // now loop to find the tndocc_ lowest eigenvalues and populate those
374 // MO's
375 int *newdocc = new int[nirrep_];
376 memset(newdocc,0,sizeof(int)*nirrep_);
377
378 for (i=0; i < tndocc_; i++) {
379 // find lowest eigenvalue
380 int lir=0,ln=0;
381 double lowest=999999999;
382
383 for (int ir=0; ir < nirrep_; ir++) {
384 int nf=oso_dimension()->blocks()->size(ir);
385 if (!nf)
386 continue;
387 for (j=0; j < nf; j++) {
388 if (vals[ir][j] < lowest) {
389 lowest=vals[ir][j];
390 lir=ir;
391 ln=j;
392 }
393 }
394 }
395 vals[lir][ln]=999999999;
396 newdocc[lir]++;
397 }
398
399 int *newsocc = new int[nirrep_];
400 memset(newsocc,0,sizeof(int)*nirrep_);
401
402 for (i=0; i < tnsocc_; i++) {
403 // find lowest eigenvalue
404 int lir=0,ln=0;
405 double lowest=999999999;
406
407 for (int ir=0; ir < nirrep_; ir++) {
408 int nf=oso_dimension()->blocks()->size(ir);
409 if (!nf)
410 continue;
411 for (j=0; j < nf; j++) {
412 if (vals[ir][j] < lowest) {
413 lowest=vals[ir][j];
414 lir=ir;
415 ln=j;
416 }
417 }
418 }
419 vals[lir][ln]=999999999;
420 newsocc[lir]++;
421 }
422
423 // get rid of vals
424 for (i=0; i < nirrep_; i++)
425 if (vals[i])
426 delete[] vals[i];
427 delete[] vals;
428
429 if (!ndocc_) {
430 ndocc_=newdocc;
431 nsocc_=newsocc;
432 } else if (most_recent_pg_.nonnull()
433 && most_recent_pg_->equiv(molecule()->point_group())) {
434 // test to see if newocc is different from ndocc_
435 for (i=0; i < nirrep_; i++) {
436 if (ndocc_[i] != newdocc[i]) {
437 ExEnv::err0() << indent << "HSOSSCF::set_occupations: WARNING!!!!\n"
438 << incindent << indent
439 << scprintf("occupations for irrep %d have changed\n",i+1)
440 << indent
441 << scprintf("ndocc was %d, changed to %d", ndocc_[i], newdocc[i])
442 << endl << decindent;
443 }
444 if (nsocc_[i] != newsocc[i]) {
445 ExEnv::err0() << indent << "HSOSSCF::set_occupations: WARNING!!!!\n"
446 << incindent << indent
447 << scprintf("occupations for irrep %d have changed\n",i+1)
448 << indent
449 << scprintf("nsocc was %d, changed to %d", nsocc_[i], newsocc[i])
450 << endl << decindent;
451 }
452 }
453
454 memcpy(ndocc_,newdocc,sizeof(int)*nirrep_);
455 memcpy(nsocc_,newsocc,sizeof(int)*nirrep_);
456 delete[] newdocc;
457 delete[] newsocc;
458 }
459
460 if (!initial_ndocc_
461 || initial_pg_->equiv(molecule()->point_group())) {
462 delete[] initial_ndocc_;
463 initial_ndocc_ = new int[nirrep_];
464 memcpy(initial_ndocc_,ndocc_,sizeof(int)*nirrep_);
465 }
466
467 if (!initial_nsocc_
468 || initial_pg_->equiv(molecule()->point_group())) {
469 delete[] initial_nsocc_;
470 initial_nsocc_ = new int[nirrep_];
471 memcpy(initial_nsocc_,nsocc_,sizeof(int)*nirrep_);
472 }
473
474 most_recent_pg_ = new PointGroup(molecule()->point_group());
475 }
476
477 void
symmetry_changed()478 HSOSSCF::symmetry_changed()
479 {
480 SCF::symmetry_changed();
481 cl_fock_.result_noupdate()=0;
482 op_fock_.result_noupdate()=0;
483 nirrep_ = molecule()->point_group()->char_table().ncomp();
484 set_occupations(0);
485 }
486
487 //////////////////////////////////////////////////////////////////////////////
488 //
489 // scf things
490 //
491
492 void
init_vector()493 HSOSSCF::init_vector()
494 {
495 init_threads();
496
497 // allocate storage for other temp matrices
498 cl_dens_ = hcore_.clone();
499 cl_dens_.assign(0.0);
500
501 cl_dens_diff_ = hcore_.clone();
502 cl_dens_diff_.assign(0.0);
503
504 op_dens_ = hcore_.clone();
505 op_dens_.assign(0.0);
506
507 op_dens_diff_ = hcore_.clone();
508 op_dens_diff_.assign(0.0);
509
510 // gmat is in AO basis
511 cl_gmat_ = basis()->matrixkit()->symmmatrix(basis()->basisdim());
512 cl_gmat_.assign(0.0);
513
514 op_gmat_ = cl_gmat_.clone();
515 op_gmat_.assign(0.0);
516
517 if (cl_fock_.result_noupdate().null()) {
518 cl_fock_ = hcore_.clone();
519 cl_fock_.result_noupdate().assign(0.0);
520 op_fock_ = hcore_.clone();
521 op_fock_.result_noupdate().assign(0.0);
522 }
523
524 // set up trial vector
525 initial_vector(1);
526
527 oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
528 }
529
530 void
done_vector()531 HSOSSCF::done_vector()
532 {
533 done_threads();
534
535 cl_gmat_ = 0;
536 cl_dens_ = 0;
537 cl_dens_diff_ = 0;
538 op_gmat_ = 0;
539 op_dens_ = 0;
540 op_dens_diff_ = 0;
541
542 oso_scf_vector_ = 0;
543 }
544
545 RefSymmSCMatrix
alpha_density()546 HSOSSCF::alpha_density()
547 {
548 RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());
549 RefSymmSCMatrix dens2(so_dimension(), basis_matrixkit());
550
551 so_density(dens1, 2.0);
552 so_density(dens2, 1.0);
553 dens1.accumulate(dens2);
554 dens2=0;
555
556 return dens1;
557 }
558
559 RefSymmSCMatrix
beta_density()560 HSOSSCF::beta_density()
561 {
562 RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
563 so_density(dens, 2.0);
564 return dens;
565 }
566
567 void
reset_density()568 HSOSSCF::reset_density()
569 {
570 cl_gmat_.assign(0.0);
571 cl_dens_diff_.assign(cl_dens_);
572
573 op_gmat_.assign(0.0);
574 op_dens_diff_.assign(op_dens_);
575 }
576
577 double
new_density()578 HSOSSCF::new_density()
579 {
580 // copy current density into density diff and scale by -1. later we'll
581 // add the new density to this to get the density difference.
582 cl_dens_diff_.assign(cl_dens_);
583 cl_dens_diff_.scale(-1.0);
584
585 op_dens_diff_.assign(op_dens_);
586 op_dens_diff_.scale(-1.0);
587
588 so_density(cl_dens_, 2.0);
589 cl_dens_.scale(2.0);
590
591 so_density(op_dens_, 1.0);
592
593 cl_dens_.accumulate(op_dens_);
594
595 cl_dens_diff_.accumulate(cl_dens_);
596 op_dens_diff_.accumulate(op_dens_);
597
598 Ref<SCElementScalarProduct> sp(new SCElementScalarProduct);
599 cl_dens_diff_.element_op(sp.pointer(), cl_dens_diff_);
600
601 double delta = sp->result();
602 delta = sqrt(delta/i_offset(cl_dens_diff_.n()));
603
604 return delta;
605 }
606
607 RefSymmSCMatrix
density()608 HSOSSCF::density()
609 {
610 if (!density_.computed()) {
611 RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
612 RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());
613 so_density(dens, 2.0);
614 dens.scale(2.0);
615
616 so_density(dens1, 1.0);
617 dens.accumulate(dens1);
618 dens1=0;
619
620 density_ = dens;
621 // only flag the density as computed if the calc is converged
622 if (!value_needed()) density_.computed() = 1;
623 }
624
625 return density_.result_noupdate();
626 }
627
628 double
scf_energy()629 HSOSSCF::scf_energy()
630 {
631 RefSymmSCMatrix t = cl_fock_.result_noupdate().copy();
632 t.accumulate(hcore_);
633
634 RefSymmSCMatrix go = op_fock_.result_noupdate().copy();
635 go.scale(-1.0);
636 go.accumulate(cl_fock_.result_noupdate());
637
638 SCFEnergy *eop = new SCFEnergy;
639 eop->reference();
640 Ref<SCElementOp2> op = eop;
641 t.element_op(op, cl_dens_);
642
643 double cl_e = eop->result();
644
645 eop->reset();
646 go.element_op(op, op_dens_);
647 double op_e = eop->result();
648
649 op=0;
650 eop->dereference();
651 delete eop;
652
653 return cl_e-op_e;
654 }
655
656 Ref<SCExtrapData>
extrap_data()657 HSOSSCF::extrap_data()
658 {
659 Ref<SCExtrapData> data =
660 new SymmSCMatrix2SCExtrapData(cl_fock_.result_noupdate(),
661 op_fock_.result_noupdate());
662 return data;
663 }
664
665 RefSymmSCMatrix
effective_fock()666 HSOSSCF::effective_fock()
667 {
668 // use fock() instead of cl_fock_ just in case this is called from
669 // someplace outside SCF::compute_vector()
670 RefSymmSCMatrix mofock(oso_dimension(), basis_matrixkit());
671 mofock.assign(0.0);
672
673 RefSymmSCMatrix mofocko(oso_dimension(), basis_matrixkit());
674 mofocko.assign(0.0);
675
676 // use eigenvectors if oso_scf_vector_ is null
677 if (oso_scf_vector_.null()) {
678 mofock.accumulate_transform(eigenvectors(), fock(0),
679 SCMatrix::TransposeTransform);
680 mofocko.accumulate_transform(eigenvectors(), fock(1),
681 SCMatrix::TransposeTransform);
682 } else {
683 RefSCMatrix so_to_oso_tr = so_to_orthog_so().t();
684 mofock.accumulate_transform(so_to_oso_tr * oso_scf_vector_, fock(0),
685 SCMatrix::TransposeTransform);
686 mofocko.accumulate_transform(so_to_oso_tr * oso_scf_vector_, fock(1),
687 SCMatrix::TransposeTransform);
688 }
689
690 Ref<SCElementOp2> op = new GSGeneralEffH(this);
691 mofock.element_op(op, mofocko);
692
693 return mofock;
694 }
695
696 /////////////////////////////////////////////////////////////////////////////
697
698 void
init_gradient()699 HSOSSCF::init_gradient()
700 {
701 // presumably the eigenvectors have already been computed by the time
702 // we get here
703 oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
704 }
705
706 void
done_gradient()707 HSOSSCF::done_gradient()
708 {
709 cl_dens_=0;
710 op_dens_=0;
711 oso_scf_vector_ = 0;
712 }
713
714 /////////////////////////////////////////////////////////////////////////////
715
716 // MO lagrangian
717 // c o v
718 // c |2*FC|2*FC|0|
719 // -------------
720 // o |2*FC| FO |0|
721 // -------------
722 // v | 0 | 0 |0|
723 //
724 RefSymmSCMatrix
lagrangian()725 HSOSSCF::lagrangian()
726 {
727 RefSCMatrix so_to_oso_tr = so_to_orthog_so().t();
728
729 RefSymmSCMatrix mofock(oso_dimension(), basis_matrixkit());
730 mofock.assign(0.0);
731 mofock.accumulate_transform(so_to_oso_tr * oso_scf_vector_,
732 cl_fock_.result_noupdate(),
733 SCMatrix::TransposeTransform);
734
735 RefSymmSCMatrix mofocko(oso_dimension(), basis_matrixkit());
736 mofocko.assign(0.0);
737 mofocko.accumulate_transform(so_to_oso_tr * oso_scf_vector_,
738 op_fock_.result_noupdate(),
739 SCMatrix::TransposeTransform);
740
741 mofock.scale(2.0);
742
743 Ref<SCElementOp2> op = new MOLagrangian(this);
744 mofock.element_op(op, mofocko);
745 mofocko=0;
746
747 // transform MO lagrangian to SO basis
748 RefSymmSCMatrix so_lag(so_dimension(), basis_matrixkit());
749 so_lag.assign(0.0);
750 so_lag.accumulate_transform(so_to_oso_tr * oso_scf_vector_, mofock);
751
752 // and then from SO to AO
753 Ref<PetiteList> pl = integral()->petite_list();
754 RefSymmSCMatrix ao_lag = pl->to_AO_basis(so_lag);
755
756 ao_lag.scale(-1.0);
757
758 return ao_lag;
759 }
760
761 RefSymmSCMatrix
gradient_density()762 HSOSSCF::gradient_density()
763 {
764 cl_dens_ = basis_matrixkit()->symmmatrix(so_dimension());
765 op_dens_ = cl_dens_.clone();
766
767 so_density(cl_dens_, 2.0);
768 cl_dens_.scale(2.0);
769
770 so_density(op_dens_, 1.0);
771
772 Ref<PetiteList> pl = integral()->petite_list(basis());
773
774 cl_dens_ = pl->to_AO_basis(cl_dens_);
775 op_dens_ = pl->to_AO_basis(op_dens_);
776
777 RefSymmSCMatrix tdens = cl_dens_.copy();
778 tdens.accumulate(op_dens_);
779
780 op_dens_.scale(2.0);
781
782 return tdens;
783 }
784
785 /////////////////////////////////////////////////////////////////////////////
786
787 void
init_hessian()788 HSOSSCF::init_hessian()
789 {
790 }
791
792 void
done_hessian()793 HSOSSCF::done_hessian()
794 {
795 }
796
797 /////////////////////////////////////////////////////////////////////////////
798
799 void
two_body_deriv_hf(double * tbgrad,double exchange_fraction)800 HSOSSCF::two_body_deriv_hf(double * tbgrad, double exchange_fraction)
801 {
802 Ref<SCElementMaxAbs> m = new SCElementMaxAbs;
803 cl_dens_.element_op(m.pointer());
804 op_dens_.element_op(m.pointer());
805 double pmax = m->result();
806 m=0;
807
808 // now try to figure out the matrix specialization we're dealing with.
809 // if we're using Local matrices, then there's just one subblock, or
810 // see if we can convert P to local matrices
811
812 if (local_ || local_dens_) {
813 // grab the data pointers from the P matrices
814 double *pmat, *pmato;
815 RefSymmSCMatrix ptmp = get_local_data(cl_dens_, pmat, SCF::Read);
816 RefSymmSCMatrix potmp = get_local_data(op_dens_, pmato, SCF::Read);
817
818 Ref<PetiteList> pl = integral()->petite_list();
819 LocalHSOSGradContribution l(pmat,pmato);
820
821 int i;
822 int na3 = molecule()->natom()*3;
823 int nthread = threadgrp_->nthread();
824 double **grads = new double*[nthread];
825 Ref<TwoBodyDerivInt> *tbis = new Ref<TwoBodyDerivInt>[nthread];
826 for (i=0; i < nthread; i++) {
827 tbis[i] = integral()->electron_repulsion_deriv();
828 grads[i] = new double[na3];
829 memset(grads[i], 0, sizeof(double)*na3);
830 }
831
832 LocalTBGrad<LocalHSOSGradContribution> **tblds =
833 new LocalTBGrad<LocalHSOSGradContribution>*[nthread];
834
835 for (i=0; i < nthread; i++) {
836 tblds[i] = new LocalTBGrad<LocalHSOSGradContribution>(
837 l, tbis[i], pl, basis(), scf_grp_, grads[i], pmax,
838 desired_gradient_accuracy(), nthread, i, exchange_fraction);
839 threadgrp_->add_thread(i, tblds[i]);
840 }
841
842 if (threadgrp_->start_threads() < 0
843 ||threadgrp_->wait_threads() < 0) {
844 ExEnv::err0() << indent
845 << "HSOSSCF: error running threads" << endl;
846 abort();
847 }
848
849 for (i=0; i < nthread; i++) {
850 for (int j=0; j < na3; j++)
851 tbgrad[j] += grads[i][j];
852
853 delete[] grads[i];
854 delete tblds[i];
855 }
856
857 scf_grp_->sum(tbgrad, na3);
858 }
859
860 // for now quit
861 else {
862 ExEnv::err0() << indent
863 << "HSOSSCF::two_body_deriv: can't do gradient yet\n";
864 abort();
865 }
866 }
867
868 /////////////////////////////////////////////////////////////////////////////
869
870 // Local Variables:
871 // mode: c++
872 // c-file-style: "ETS"
873 // End:
874