1 //
2 // wfn.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.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 <stdlib.h>
33 #include <math.h>
34
35 #include <stdexcept>
36 #include <iostream>
37 #include <iterator>
38
39 #include <util/keyval/keyval.h>
40 #include <util/misc/timer.h>
41 #include <util/misc/formio.h>
42 #include <util/misc/autovec.h>
43 #include <util/state/stateio.h>
44 #include <chemistry/qc/basis/obint.h>
45 #include <chemistry/qc/basis/symmint.h>
46 #include <chemistry/qc/intv3/intv3.h>
47
48 #include <chemistry/qc/wfn/wfn.h>
49
50 using namespace std;
51 using namespace sc;
52
53 #define CHECK_SYMMETRIZED_INTEGRALS 0
54
55 /////////////////////////////////////////////////////////////////////////
56
57 // This maps a TwoBodyThreeCenterInt into a OneBodyInt
58 namespace sc {
59 class ChargeDistInt: public OneBodyInt {
60 Ref<TwoBodyThreeCenterInt> tbint_;
61 Ref<Molecule> mol_;
62 Ref<GaussianBasisSet> ebasis0_;
63 Ref<GaussianBasisSet> ebasis1_;
64 Ref<GaussianBasisSet> atom_basis_;
65 std::vector<int> i_cd_;
66 const double *coef_;
67 public:
68 // The electronic basis are bs1 and bs2 in tbint and the
69 // nuclear basis is bs3.
70 ChargeDistInt(const Ref<TwoBodyThreeCenterInt>& tbint,
71 const double *coef);
72
73 void compute_shell(int,int);
74
75 bool cloneable();
76 };
77
ChargeDistInt(const Ref<TwoBodyThreeCenterInt> & tbint,const double * coef)78 ChargeDistInt::ChargeDistInt(const Ref<TwoBodyThreeCenterInt>& tbint,
79 const double *coef):
80 OneBodyInt(tbint->integral(),tbint->basis1(),tbint->basis2()),
81 tbint_(tbint),
82 coef_(coef)
83 {
84 ebasis0_ = tbint->basis1();
85 ebasis1_ = tbint->basis2();
86 atom_basis_ = tbint->basis3();
87 mol_ = atom_basis_->molecule();
88 buffer_ = new double[tbint->basis1()->max_nfunction_in_shell()
89 *tbint->basis2()->max_nfunction_in_shell()];
90
91 for (int i=0; i<atom_basis_->ncenter(); i++) {
92 if (atom_basis_->nshell_on_center(i) > 0) i_cd_.push_back(i);
93 }
94 }
95
96 void
compute_shell(int ish,int jsh)97 ChargeDistInt::compute_shell(int ish,int jsh)
98 {
99 // std::cout << "starting " << ish << " " << jsh << std::endl;
100 int nijbf
101 = ebasis0_->shell(ish).nfunction()
102 * ebasis1_->shell(jsh).nfunction();
103 memset(buffer_,0,nijbf*sizeof(buffer_[0]));
104 const double *tbbuffer = tbint_->buffer();
105 int ksh = 0;
106 int coef_offset = 0;
107 for (int ii=0; ii<i_cd_.size(); ii++) {
108 int i = i_cd_[ii];
109 int nshell = atom_basis_->nshell_on_center(i);
110 for (int j=0; j<nshell; j++, ksh++) {
111 std::cout << "computing " << ish << " " << jsh << " " << ksh << std::endl;
112 tbint_->compute_shell(ish,jsh,ksh);
113 int nbfk = atom_basis_->shell(i).nfunction();
114 for (int ijbf=0; ijbf<nijbf; ijbf++) {
115 for (int kbf=0; kbf<nbfk; kbf++) {
116 buffer_[ijbf]
117 -= coef_[coef_offset+kbf] * tbbuffer[ijbf*nbfk + kbf];
118 // std::cout << " adding "
119 // << coef_[coef_offset+kbf] * tbbuffer[ijbf*nbfk + kbf]
120 // << " = "
121 // << coef_[coef_offset+kbf]
122 // << " * "
123 // << tbbuffer[ijbf*nbfk + kbf]
124 // << " at location "
125 // << ijbf
126 // << std::endl;
127 }
128 }
129 }
130 coef_offset += nshell;
131 }
132 // std::cout << "done with " << ish << " " << jsh << std::endl;
133 }
134
135 bool
cloneable()136 ChargeDistInt::cloneable()
137 {
138 // not cloneable because tbint is not cloneable
139 return false;
140 }
141
142 }
143
144 /////////////////////////////////////////////////////////////////////////
145
146 static ClassDesc Wavefunction_cd(
147 typeid(Wavefunction),"Wavefunction",7,"public MolecularEnergy",
148 0, 0, 0);
149
Wavefunction(const Ref<KeyVal> & keyval)150 Wavefunction::Wavefunction(const Ref<KeyVal>&keyval):
151 // this will let molecule be retrieved from basis
152 // MolecularEnergy(new AggregateKeyVal(keyval,
153 // new PrefixKeyVal("basis", keyval))),
154 MolecularEnergy(keyval),
155 overlap_(this),
156 hcore_(this),
157 natural_orbitals_(this),
158 natural_density_(this),
159 bs_values(0),
160 bsg_values(0)
161 {
162 overlap_.compute() = 0;
163 hcore_.compute() = 0;
164 natural_orbitals_.compute() = 0;
165 natural_density_.compute() = 0;
166
167 overlap_.computed() = 0;
168 hcore_.computed() = 0;
169 natural_orbitals_.computed() = 0;
170 natural_density_.computed() = 0;
171
172 print_nao_ = keyval->booleanvalue("print_nao");
173 print_npa_ = keyval->booleanvalue("print_npa");
174 KeyValValuedouble lindep_tol_def(1.e-8);
175 lindep_tol_ = keyval->doublevalue("lindep_tol", lindep_tol_def);
176 if (keyval->exists("symm_orthog")) {
177 ExEnv::out0() << indent
178 << "WARNING: using obsolete \"symm_orthog\" keyword"
179 << endl;
180 if (keyval->booleanvalue("symm_orthog")) {
181 orthog_method_ = OverlapOrthog::Symmetric;
182 }
183 else {
184 orthog_method_ = OverlapOrthog::Canonical;
185 }
186 }
187 else {
188 char *orthog_name = keyval->pcharvalue("orthog_method");
189 if (!orthog_name) {
190 orthog_method_ = OverlapOrthog::Symmetric;
191 }
192 else if (::strcmp(orthog_name, "canonical") == 0) {
193 orthog_method_ = OverlapOrthog::Canonical;
194 }
195 else if (::strcmp(orthog_name, "symmetric") == 0) {
196 orthog_method_ = OverlapOrthog::Symmetric;
197 }
198 else if (::strcmp(orthog_name, "gramschmidt") == 0) {
199 orthog_method_ = OverlapOrthog::GramSchmidt;
200 }
201 else {
202 ExEnv::errn() << "ERROR: bad orthog_method: \""
203 << orthog_name << "\"" << endl;
204 abort();
205 }
206 delete[] orthog_name;
207 }
208
209 debug_ = keyval->intvalue("debug");
210
211 gbs_ = require_dynamic_cast<GaussianBasisSet*>(
212 keyval->describedclassvalue("basis").pointer(),
213 "Wavefunction::Wavefunction\n"
214 );
215
216 atom_basis_ << keyval->describedclassvalue("atom_basis");
217 if (atom_basis_.nonnull()) {
218 atom_basis_coef_ = new double[atom_basis_->nbasis()];
219 for (int i=0; i<atom_basis_->nbasis(); i++) {
220 atom_basis_coef_[i] = keyval->doublevalue("atom_basis_coef",i);
221 }
222 scale_atom_basis_coef();
223
224 }
225 else {
226 atom_basis_coef_ = 0;
227 }
228
229 integral_ << keyval->describedclassvalue("integrals");
230 if (integral_.null()) {
231 Integral* default_intf = Integral::get_default_integral();
232 integral_ = default_intf->clone();
233 }
234
235 integral_->set_basis(gbs_);
236 Ref<PetiteList> pl = integral_->petite_list();
237
238 sodim_ = pl->SO_basisdim();
239 aodim_ = pl->AO_basisdim();
240 basiskit_ = gbs_->so_matrixkit();
241 }
242
Wavefunction(StateIn & s)243 Wavefunction::Wavefunction(StateIn&s):
244 SavableState(s),
245 MolecularEnergy(s),
246 overlap_(this),
247 hcore_(this),
248 natural_orbitals_(this),
249 natural_density_(this),
250 bs_values(0),
251 bsg_values(0)
252 {
253 debug_ = 0;
254
255 overlap_.compute() = 0;
256 hcore_.compute() = 0;
257 natural_orbitals_.compute() = 0;
258 natural_density_.compute() = 0;
259
260 overlap_.computed() = 0;
261 hcore_.computed() = 0;
262 natural_orbitals_.computed() = 0;
263 natural_density_.computed() = 0;
264
265 if (s.version(::class_desc<Wavefunction>()) >= 2) {
266 s.get(print_nao_);
267 s.get(print_npa_);
268 }
269 else {
270 print_nao_ = 0;
271 print_npa_ = 0;
272 }
273
274 if (s.version(::class_desc<Wavefunction>()) >= 5) {
275 int orthog_enum;
276 s.get(orthog_enum);
277 orthog_method_ = (OverlapOrthog::OrthogMethod) orthog_enum;
278 }
279 else if (s.version(::class_desc<Wavefunction>()) >= 3) {
280 int symm_orthog;
281 s.get(symm_orthog);
282 if (symm_orthog) orthog_method_ = OverlapOrthog::Symmetric;
283 else orthog_method_ = OverlapOrthog::Canonical;
284 }
285 else {
286 orthog_method_ = OverlapOrthog::Symmetric;
287 }
288
289 if (s.version(::class_desc<Wavefunction>()) >= 4) {
290 s.get(lindep_tol_);
291 }
292 else {
293 lindep_tol_ = 1.e-6;
294 }
295
296 gbs_ << SavableState::restore_state(s);
297 integral_ << SavableState::restore_state(s);
298
299 if (s.version(::class_desc<Wavefunction>()) >= 6) {
300 Ref<KeyVal> original_override = s.override();
301 Ref<AssignedKeyVal> matrixkit_override = new AssignedKeyVal;
302 matrixkit_override->assign("matrixkit", gbs_->so_matrixkit().pointer());
303 Ref<KeyVal> new_override
304 = new AggregateKeyVal(matrixkit_override.pointer(),
305 original_override.pointer());
306 s.set_override(new_override);
307 orthog_ << SavableState::restore_state(s);
308 s.set_override(original_override);
309 }
310
311 if (s.version(::class_desc<Wavefunction>()) >= 7) {
312 atom_basis_ << SavableState::restore_state(s);
313 if (atom_basis_.nonnull()) {
314 s.get_array_double(atom_basis_coef_, atom_basis_->nbasis());
315 }
316 }
317 else {
318 atom_basis_coef_ = 0;
319 }
320
321 integral_->set_basis(gbs_);
322 Ref<PetiteList> pl = integral_->petite_list();
323
324 sodim_ = pl->SO_basisdim();
325 aodim_ = pl->AO_basisdim();
326 basiskit_ = gbs_->so_matrixkit();
327 }
328
329 void
symmetry_changed()330 Wavefunction::symmetry_changed()
331 {
332 MolecularEnergy::symmetry_changed();
333
334 Ref<PetiteList> pl = integral_->petite_list();
335 sodim_ = pl->SO_basisdim();
336 aodim_ = pl->AO_basisdim();
337 overlap_.result_noupdate() = 0;
338 basiskit_ = gbs_->so_matrixkit();
339
340 orthog_ = 0;
341 }
342
~Wavefunction()343 Wavefunction::~Wavefunction()
344 {
345 if (bs_values) {
346 delete[] bs_values;
347 bs_values=0;
348 }
349 if (bsg_values) {
350 delete[] bsg_values;
351 bsg_values=0;
352 }
353 }
354
355 void
save_data_state(StateOut & s)356 Wavefunction::save_data_state(StateOut&s)
357 {
358 MolecularEnergy::save_data_state(s);
359
360 // overlap and hcore integrals are cheap so don't store them.
361 // same goes for natural orbitals
362
363 s.put(print_nao_);
364 s.put(print_npa_);
365 s.put((int)orthog_method_);
366 s.put(lindep_tol_);
367
368 SavableState::save_state(gbs_.pointer(),s);
369 SavableState::save_state(integral_.pointer(),s);
370 SavableState::save_state(orthog_.pointer(),s);
371 SavableState::save_state(atom_basis_.pointer(), s);
372 if (atom_basis_.nonnull()) {
373 s.put_array_double(atom_basis_coef_,atom_basis_->nbasis());
374 }
375 }
376
377 double
charge()378 Wavefunction::charge()
379 {
380 return molecule()->nuclear_charge() - nelectron();
381 }
382
383 RefSymmSCMatrix
ao_density()384 Wavefunction::ao_density()
385 {
386 return integral()->petite_list()->to_AO_basis(density());
387 }
388
389 RefSCMatrix
natural_orbitals()390 Wavefunction::natural_orbitals()
391 {
392 if (!natural_orbitals_.computed()) {
393 RefSymmSCMatrix dens = density();
394
395 // transform the density into an orthogonal basis
396 RefSCMatrix ortho = so_to_orthog_so();
397
398 RefSymmSCMatrix densortho(oso_dimension(), basis_matrixkit());
399 densortho.assign(0.0);
400 densortho.accumulate_transform(so_to_orthog_so_inverse().t(),dens);
401
402 RefSCMatrix natorb(oso_dimension(), oso_dimension(),
403 basis_matrixkit());
404 RefDiagSCMatrix natden(oso_dimension(), basis_matrixkit());
405
406 densortho.diagonalize(natden, natorb);
407
408 // natorb is the ortho SO to NO basis transform
409 // make natural_orbitals_ the SO to the NO basis transform
410 natural_orbitals_ = so_to_orthog_so().t() * natorb;
411 natural_density_ = natden;
412
413 natural_orbitals_.computed() = 1;
414 natural_density_.computed() = 1;
415 }
416
417 return natural_orbitals_.result_noupdate();
418 }
419
420 RefDiagSCMatrix
natural_density()421 Wavefunction::natural_density()
422 {
423 if (!natural_density_.computed()) {
424 natural_orbitals();
425 }
426
427 return natural_density_.result_noupdate();
428 }
429
430 RefSymmSCMatrix
overlap()431 Wavefunction::overlap()
432 {
433 if (!overlap_.computed()) {
434 integral()->set_basis(gbs_);
435 Ref<PetiteList> pl = integral()->petite_list();
436
437 #if ! CHECK_SYMMETRIZED_INTEGRALS
438 // first form skeleton s matrix
439 RefSymmSCMatrix s(basis()->basisdim(), basis()->matrixkit());
440 Ref<SCElementOp> ov =
441 new OneBodyIntOp(new SymmOneBodyIntIter(integral()->overlap(), pl));
442
443 s.assign(0.0);
444 s.element_op(ov);
445 ov=0;
446
447 if (debug_ > 1) s.print("AO skeleton overlap");
448
449 // then symmetrize it
450 RefSymmSCMatrix sb(so_dimension(), basis_matrixkit());
451 pl->symmetrize(s,sb);
452
453 overlap_ = sb;
454 #else
455 ExEnv::out0() << "Checking symmetrized overlap" << endl;
456
457 RefSymmSCMatrix s(basis()->basisdim(), basis()->matrixkit());
458 Ref<SCElementOp> ov =
459 new OneBodyIntOp(new OneBodyIntIter(integral()->overlap()));
460
461 s.assign(0.0);
462 s.element_op(ov);
463 ov=0;
464
465 overlap_ = pl->to_SO_basis(s);
466
467 //// use petite list to form saopl
468
469 // form skeleton Hcore in AO basis
470 RefSymmSCMatrix saopl(basis()->basisdim(), basis()->matrixkit());
471 saopl.assign(0.0);
472
473 ov = new OneBodyIntOp(new SymmOneBodyIntIter(integral_->overlap(), pl));
474 saopl.element_op(ov);
475 ov=0;
476
477 // also symmetrize using pl->symmetrize():
478 RefSymmSCMatrix spl(so_dimension(), basis_matrixkit());
479 pl->symmetrize(saopl,spl);
480
481 //// compare the answers
482
483 int n = overlap_.result_noupdate().dim().n();
484 int me = MessageGrp::get_default_messagegrp()->me();
485 for (int i=0; i<n; i++) {
486 for (int j=0; j<=i; j++) {
487 double val1 = overlap_.result_noupdate().get_element(i,j);
488 double val2 = spl.get_element(i,j);
489 if (me == 0) {
490 if (fabs(val1-val2) > 1.0e-6) {
491 ExEnv::out0() << "bad overlap vals for " << i << " " << j
492 << ": " << val1 << " " << val2 << endl;
493 }
494 }
495 }
496 }
497 #endif
498
499 overlap_.computed() = 1;
500 }
501
502 return overlap_.result_noupdate();
503 }
504
505 RefSymmSCMatrix
core_hamiltonian()506 Wavefunction::core_hamiltonian()
507 {
508 if (!hcore_.computed()) {
509 integral()->set_basis(gbs_);
510 Ref<PetiteList> pl = integral()->petite_list();
511
512 #if ! CHECK_SYMMETRIZED_INTEGRALS
513 // form skeleton Hcore in AO basis
514 RefSymmSCMatrix hao(basis()->basisdim(), basis()->matrixkit());
515 hao.assign(0.0);
516
517 Ref<SCElementOp> hc =
518 new OneBodyIntOp(new SymmOneBodyIntIter(integral_->kinetic(), pl));
519 hao.element_op(hc);
520 hc=0;
521
522 if (atom_basis_.null()) {
523 Ref<OneBodyInt> nuc = integral_->nuclear();
524 nuc->reinitialize();
525 hc = new OneBodyIntOp(new SymmOneBodyIntIter(nuc, pl));
526 hao.element_op(hc);
527 hc=0;
528 }
529 else {
530 // we have an atom_basis, so some nuclear charges will be computed
531 // with the two electron integral code and some will be computed
532 // with the point charge code
533
534 // find out which atoms are point charges and which ones are charge
535 // distributions
536 std::vector<int> i_pc; // point charge centers
537 for (int i=0; i<atom_basis_->ncenter(); i++) {
538 if (atom_basis_->nshell_on_center(i) == 0) i_pc.push_back(i);
539 }
540 int n_pc = i_pc.size();
541
542 // initialize the point charge data
543 auto_vec<double> pc_charges(new double[n_pc]);
544 auto_vec<double*> pc_xyz(new double*[n_pc]);
545 auto_vec<double> pc_xyz0(new double[n_pc*3]);
546 pc_xyz[0] = pc_xyz0.get();
547 for (int i=1; i<n_pc; i++) pc_xyz[i] = pc_xyz[i-1] + 3;
548 for (int i=0; i<n_pc; i++) {
549 pc_charges[i] = molecule()->charge(i_pc[i]);
550 for (int j=0; j<3; j++) pc_xyz[i][j] = molecule()->r(i_pc[i],j);
551 }
552 Ref<PointChargeData> pc_data
553 = new PointChargeData(n_pc, pc_xyz.get(), pc_charges.get());
554
555 // compute the point charge contributions
556 Ref<OneBodyInt> pc_int = integral_->point_charge(pc_data);
557 hc = new OneBodyIntOp(new SymmOneBodyIntIter(pc_int,pl));
558 hao.element_op(hc);
559 hc=0;
560 pc_int=0;
561
562 // compute the charge distribution contributions
563 // H_ij += sum_A -q_A sum_k N_A_k (ij|A_k)
564 integral()->set_basis(gbs_,gbs_,atom_basis_);
565 Ref<TwoBodyThreeCenterInt> cd_tbint
566 = integral_->electron_repulsion3();
567 Ref<OneBodyInt> cd_int = new ChargeDistInt(cd_tbint, atom_basis_coef_);
568 hc = new OneBodyIntOp(new SymmOneBodyIntIter(cd_int,pl));
569 hao.element_op(hc);
570 integral()->set_basis(gbs_);
571 hc=0;
572 cd_int=0;
573 cd_tbint=0;
574 }
575
576 // now symmetrize Hso
577 RefSymmSCMatrix h(so_dimension(), basis_matrixkit());
578 pl->symmetrize(hao,h);
579
580 hcore_ = h;
581 #else
582 ExEnv::out0() << "Checking symmetrized hcore" << endl;
583
584 RefSymmSCMatrix hao(basis()->basisdim(), basis()->matrixkit());
585 hao.assign(0.0);
586
587 Ref<SCElementOp> hc =
588 new OneBodyIntOp(new OneBodyIntIter(integral_->kinetic()));
589 hao.element_op(hc);
590 hc=0;
591
592 // std::cout << "null atom_basis" << std::endl;
593 Ref<OneBodyInt> nuc = integral_->nuclear();
594 nuc->reinitialize();
595 hc = new OneBodyIntOp(new OneBodyIntIter(nuc));
596 hao.element_op(hc);
597 hc=0;
598
599 hcore_ = pl->to_SO_basis(hao);
600
601 //// use petite list to form haopl
602
603 // form skeleton Hcore in AO basis
604 RefSymmSCMatrix haopl(basis()->basisdim(), basis()->matrixkit());
605 haopl.assign(0.0);
606
607 hc = new OneBodyIntOp(new SymmOneBodyIntIter(integral_->kinetic(), pl));
608 haopl.element_op(hc);
609 hc=0;
610
611 nuc = integral_->nuclear();
612 nuc->reinitialize();
613 hc = new OneBodyIntOp(new SymmOneBodyIntIter(nuc, pl));
614 haopl.element_op(hc);
615 hc=0;
616
617 // also symmetrize using pl->symmetrize():
618 RefSymmSCMatrix h(so_dimension(), basis_matrixkit());
619 pl->symmetrize(haopl,h);
620
621 //// compare the answers
622
623 int n = hcore_.result_noupdate().dim().n();
624 int me = MessageGrp::get_default_messagegrp()->me();
625 for (int i=0; i<n; i++) {
626 for (int j=0; j<=i; j++) {
627 double val1 = hcore_.result_noupdate().get_element(i,j);
628 double val2 = h.get_element(i,j);
629 if (me == 0) {
630 if (fabs(val1-val2) > 1.0e-6) {
631 ExEnv::outn() << "bad hcore vals for " << i << " " << j
632 << ": " << val1 << " " << val2 << endl;
633 }
634 }
635 }
636 }
637 #endif
638
639 hcore_.computed() = 1;
640 }
641
642 return hcore_.result_noupdate();
643 }
644
645 // Compute lists of centers that are point charges and lists that
646 // are charge distributions.
647 void
set_up_charge_types(std::vector<int> & q_pc,std::vector<int> & q_cd,std::vector<int> & n_pc,std::vector<int> & n_cd)648 Wavefunction::set_up_charge_types(
649 std::vector<int> &q_pc,
650 std::vector<int> &q_cd,
651 std::vector<int> &n_pc,
652 std::vector<int> &n_cd)
653 {
654 q_pc.clear();
655 q_cd.clear();
656 n_pc.clear();
657 n_cd.clear();
658
659 for (int i=0; i<atom_basis_->ncenter(); i++) {
660 bool is_Q = (molecule()->atom_symbol(i) == "Q");
661 if (atom_basis_->nshell_on_center(i) > 0) {
662 if (is_Q) q_cd.push_back(i);
663 else n_cd.push_back(i);
664 }
665 else {
666 if (is_Q) q_pc.push_back(i);
667 else n_pc.push_back(i);
668 }
669 }
670 }
671
672 double
nuclear_repulsion_energy()673 Wavefunction::nuclear_repulsion_energy()
674 {
675 if (atom_basis_.null()) return molecule()->nuclear_repulsion_energy();
676
677 double nucrep = 0.0;
678
679 std::vector<int> q_pc, q_cd, n_pc, n_cd;
680 set_up_charge_types(q_pc,q_cd,n_pc,n_cd);
681
682 // compute point charge - point charge contribution
683 nucrep += nuc_rep_pc_pc(n_pc, n_pc, true /* i > j */);
684 nucrep += nuc_rep_pc_pc(q_pc, n_pc, false /* all i j */);
685 if (molecule()->include_qq()) {
686 nucrep += nuc_rep_pc_pc(q_pc, q_pc, true /* i > j */);
687 }
688
689 // compute point charge - charge distribution contribution
690 nucrep += nuc_rep_pc_cd(n_pc, n_cd);
691 nucrep += nuc_rep_pc_cd(q_pc, n_cd);
692 nucrep += nuc_rep_pc_cd(n_pc, q_cd);
693 if (molecule()->include_qq()) {
694 nucrep += nuc_rep_pc_cd(q_pc, q_cd);
695 }
696
697 // compute the charge distribution - charge distribution contribution
698 nucrep += nuc_rep_cd_cd(n_cd, n_cd, true /* i > j */);
699 nucrep += nuc_rep_cd_cd(q_cd, n_cd, false /* all i j */);
700 if (molecule()->include_qq()) {
701 nucrep += nuc_rep_cd_cd(q_cd, q_cd, true /* i > j */);
702 }
703
704 return nucrep;
705 }
706
707 double
nuc_rep_pc_pc(const std::vector<int> & c1,const std::vector<int> & c2,bool uniq)708 Wavefunction::nuc_rep_pc_pc(const std::vector<int>&c1,
709 const std::vector<int>&c2,
710 bool uniq)
711 {
712 double e = 0.0;
713
714 if (c1.size() == 0 || c2.size() == 0) return e;
715
716 for (int ii=0; ii<c1.size(); ii++) {
717 int i = c1[ii];
718 SCVector3 ai(molecule()->r(i));
719 double Zi = molecule()->charge(i);
720 int jfence = (uniq?ii:c2.size());
721 for (int jj=0; jj<jfence; jj++) {
722 int j = c2[jj];
723 SCVector3 aj(molecule()->r(j));
724 e += Zi * molecule()->charge(j) / ai.dist(aj);
725 }
726 }
727
728 return e;
729 }
730
731 double
nuc_rep_pc_cd(const std::vector<int> & pc,const std::vector<int> & cd)732 Wavefunction::nuc_rep_pc_cd(const std::vector<int>&pc,
733 const std::vector<int>&cd)
734 {
735 double e = 0.0;
736
737 if (pc.size() == 0 || cd.size() == 0) return e;
738
739 integral()->set_basis(atom_basis());
740
741 sc::auto_vec<double> charges(new double[pc.size()]);
742 sc::auto_vec<double*> positions(new double*[pc.size()]);
743 sc::auto_vec<double> xyz(new double[pc.size()*3]);
744 for (int i=0; i<pc.size(); i++) {
745 positions.get()[i] = &xyz.get()[i*3];
746 }
747
748 for (int ii=0; ii<pc.size(); ii++) {
749 int i=pc[ii];
750 charges.get()[ii] = molecule()->charge(i);
751 for (int j=0; j<3; j++) positions.get()[ii][j] = molecule()->r(i,j);
752 }
753
754 Ref<PointChargeData> pcdata = new PointChargeData(pc.size(),
755 positions.get(),
756 charges.get());
757 Ref<OneBodyOneCenterInt> ob = integral()->point_charge1(pcdata);
758
759 const double *buffer = ob->buffer();
760
761 for (int ii=0,icoef=0; ii<cd.size(); ii++) {
762 int icenter = cd[ii];
763 int joff = atom_basis()->shell_on_center(icenter,0);
764 for (int j=0; j<atom_basis()->nshell_on_center(icenter); j++) {
765 int jsh = j + joff;
766 ob->compute_shell(jsh);
767 int nfunc = atom_basis()->shell(jsh).nfunction();
768 for (int k=0; k<nfunc; k++,icoef++) {
769 e -= atom_basis_coef_[icoef] * buffer[k];
770 }
771 }
772 }
773
774 integral()->set_basis(basis());
775
776 return e;
777 }
778
779 double
nuc_rep_cd_cd(const std::vector<int> & c1,const std::vector<int> & c2,bool uniq)780 Wavefunction::nuc_rep_cd_cd(const std::vector<int>&c1,
781 const std::vector<int>&c2,
782 bool uniq)
783 {
784 double e = 0.0;
785
786 if (c1.size() == 0 || c2.size() == 0) return e;
787
788 integral()->set_basis(atom_basis());
789
790 Ref<TwoBodyTwoCenterInt> tb = integral()->electron_repulsion2();
791
792 const double *buffer = tb->buffer();
793
794 for (int ii=0; ii<c1.size(); ii++) {
795 int icenter = c1[ii];
796 int inshell = atom_basis()->nshell_on_center(icenter);
797 int ish = atom_basis()->shell_on_center(icenter,0);
798 for (int iish=0; iish<inshell; iish++,ish++) {
799 int infunc = atom_basis()->shell(ish).nfunction();
800 int ifuncoff = atom_basis()->shell_to_function(ish);
801 int jjfence = (uniq?ii:c2.size());
802 for (int jj=0; jj<jjfence; jj++) {
803 int jcenter = c2[jj];
804 int jnshell = atom_basis()->nshell_on_center(jcenter);
805 int jsh = atom_basis()->shell_on_center(jcenter,0);
806 for (int jjsh=0; jjsh<jnshell; jjsh++,jsh++) {
807 int jnfunc = atom_basis()->shell(jsh).nfunction();
808 int jfuncoff = atom_basis()->shell_to_function(jsh);
809 tb->compute_shell(ish,jsh);
810 for (int ifunc=0, ijfunc=0; ifunc<infunc; ifunc++) {
811 for (int jfunc=0; jfunc<jnfunc; jfunc++, ijfunc++) {
812 e += atom_basis_coef_[ifuncoff+ifunc]
813 * atom_basis_coef_[jfuncoff+jfunc]
814 * buffer[ijfunc];
815 }
816 }
817 }
818 }
819 }
820 }
821
822 integral()->set_basis(basis());
823
824 return e;
825 }
826
827 void
nuclear_repulsion_energy_gradient(double * g)828 Wavefunction::nuclear_repulsion_energy_gradient(double *g)
829 {
830 int natom = molecule()->natom();
831 sc::auto_vec<double*> gp(new double*[natom]);
832 for (int i=0; i<natom; i++) {
833 gp.get()[i] = &g[i*3];
834 }
835 nuclear_repulsion_energy_gradient(gp.get());
836 }
837
838 void
nuclear_repulsion_energy_gradient(double ** g)839 Wavefunction::nuclear_repulsion_energy_gradient(double **g)
840 {
841 if (atom_basis_.null()) {
842 int natom = molecule()->natom();
843 for (int i=0; i<natom; i++) {
844 molecule()->nuclear_repulsion_1der(i,g[i]);
845 }
846 return;
847 }
848
849 // zero the gradient
850 for (int i=0; i<molecule()->natom(); i++) {
851 for (int j=0; j<3; j++) g[i][j] = 0.0;
852 }
853
854 // compute charge types
855 std::vector<int> q_pc, q_cd, n_pc, n_cd;
856 set_up_charge_types(q_pc,q_cd,n_pc,n_cd);
857
858 // compute point charge - point charge contribution
859 nuc_rep_grad_pc_pc(g, n_pc, n_pc, true /* i > j */);
860 nuc_rep_grad_pc_pc(g, q_pc, n_pc, false /* all i j */);
861 if (molecule()->include_qq()) {
862 nuc_rep_grad_pc_pc(g, q_pc, q_pc, true /* i > j */);
863 }
864
865 // compute point charge - charge distribution contribution
866 nuc_rep_grad_pc_cd(g, n_pc, n_cd);
867 nuc_rep_grad_pc_cd(g, q_pc, n_cd);
868 nuc_rep_grad_pc_cd(g, n_pc, q_cd);
869 if (molecule()->include_qq()) {
870 nuc_rep_grad_pc_cd(g, q_pc, q_cd);
871 }
872
873 // compute the charge distribution - charge distribution contribution
874 nuc_rep_grad_cd_cd(g, n_cd, n_cd, true /* i > j */);
875 nuc_rep_grad_cd_cd(g, q_cd, n_cd, false /* all i j */);
876 if (molecule()->include_qq()) {
877 nuc_rep_grad_cd_cd(g, q_cd, q_cd, true /* i > j */);
878 }
879
880 // note: the electronic terms still need to be done in
881 // a new hcore_deriv implemented in Wavefunction.
882 throw std::runtime_error("Wavefunction::nuclear_repulsion_energy_gradient: not done");
883
884 }
885
886 void
nuc_rep_grad_pc_pc(double ** grad,const std::vector<int> & c1,const std::vector<int> & c2,bool uniq)887 Wavefunction::nuc_rep_grad_pc_pc(double **grad,
888 const std::vector<int>&c1,
889 const std::vector<int>&c2,
890 bool uniq)
891 {
892 if (c1.size() == 0 || c2.size() == 0) return;
893
894 for (int ii=0; ii<c1.size(); ii++) {
895 int i = c1[ii];
896 SCVector3 ai(molecule()->r(i));
897 double Zi = molecule()->charge(i);
898 int jfence = (uniq?ii:c2.size());
899 for (int jj=0; jj<jfence; jj++) {
900 int j = c2[jj];
901 SCVector3 aj(molecule()->r(j));
902 double Zj = molecule()->charge(j);
903 SCVector3 rdiff = ai - aj;
904 double r2 = rdiff.dot(rdiff);
905 double factor = - Zi * Zj / (r2*sqrt(r2));
906 for (int k=0; k<3; k++) {
907 grad[i][k] += factor * rdiff[k];
908 grad[j][k] -= factor * rdiff[k];
909 }
910 }
911 }
912 }
913
914 void
nuc_rep_grad_pc_cd(double ** grad,const std::vector<int> & pc,const std::vector<int> & cd)915 Wavefunction::nuc_rep_grad_pc_cd(double **grad,
916 const std::vector<int>&pc,
917 const std::vector<int>&cd)
918 {
919 if (pc.size() == 0 || cd.size() == 0) return;
920
921 throw std::runtime_error("Wavefunction::nuclear_repulsion_energy_gradient: not done");
922 }
923
924 void
nuc_rep_grad_cd_cd(double ** grad,const std::vector<int> & c1,const std::vector<int> & c2,bool uniq)925 Wavefunction::nuc_rep_grad_cd_cd(double **grad,
926 const std::vector<int>&c1,
927 const std::vector<int>&c2,
928 bool uniq)
929 {
930 if (c1.size() == 0 || c2.size() == 0) return;
931
932 throw std::runtime_error("Wavefunction::nuclear_repulsion_energy_gradient: not done");
933 }
934
935 // returns the orthogonalization matrix
936 RefSCMatrix
so_to_orthog_so()937 Wavefunction::so_to_orthog_so()
938 {
939 if (orthog_.null()) init_orthog();
940 return orthog_->basis_to_orthog_basis();
941 }
942
943 RefSCMatrix
so_to_orthog_so_inverse()944 Wavefunction::so_to_orthog_so_inverse()
945 {
946 if (orthog_.null()) init_orthog();
947 return orthog_->basis_to_orthog_basis_inverse();
948 }
949
950 Ref<GaussianBasisSet>
basis() const951 Wavefunction::basis() const
952 {
953 return gbs_;
954 }
955
956 Ref<Molecule>
molecule() const957 Wavefunction::molecule() const
958 {
959 return basis()->molecule();
960 }
961
962 Ref<GaussianBasisSet>
atom_basis() const963 Wavefunction::atom_basis() const
964 {
965 return atom_basis_;
966 }
967
968 const double *
atom_basis_coef() const969 Wavefunction::atom_basis_coef() const
970 {
971 return atom_basis_coef_;
972 }
973
974 Ref<Integral>
integral()975 Wavefunction::integral()
976 {
977 return integral_;
978 }
979
980 RefSCDimension
so_dimension()981 Wavefunction::so_dimension()
982 {
983 return sodim_;
984 }
985
986 RefSCDimension
ao_dimension()987 Wavefunction::ao_dimension()
988 {
989 return aodim_;
990 }
991
992 RefSCDimension
oso_dimension()993 Wavefunction::oso_dimension()
994 {
995 if (orthog_.null()) init_orthog();
996 return orthog_->orthog_dim();
997 }
998
999 Ref<SCMatrixKit>
basis_matrixkit()1000 Wavefunction::basis_matrixkit()
1001 {
1002 return basiskit_;
1003 }
1004
1005 void
print(ostream & o) const1006 Wavefunction::print(ostream&o) const
1007 {
1008 MolecularEnergy::print(o);
1009 ExEnv::out0() << indent << "Electronic basis:" << std::endl;
1010 ExEnv::out0() << incindent;
1011 basis()->print_brief(o);
1012 ExEnv::out0() << decindent;
1013 if (atom_basis_.nonnull()) {
1014 ExEnv::out0() << indent << "Nuclear basis:" << std::endl;
1015 ExEnv::out0() << incindent;
1016 atom_basis_->print_brief(o);
1017 ExEnv::out0() << decindent;
1018 }
1019 // the other stuff is a wee bit too big to print
1020 if (print_nao_ || print_npa_) {
1021 tim_enter("NAO");
1022 RefSCMatrix naos = ((Wavefunction*)this)->nao();
1023 tim_exit("NAO");
1024 if (print_nao_) naos.print("NAO");
1025 }
1026 }
1027
1028 RefSymmSCMatrix
alpha_density()1029 Wavefunction::alpha_density()
1030 {
1031 if (!spin_polarized()) {
1032 RefSymmSCMatrix result = density().copy();
1033 result.scale(0.5);
1034 return result;
1035 }
1036 ExEnv::errn() << class_name() << "::alpha_density not implemented" << endl;
1037 abort();
1038 return 0;
1039 }
1040
1041 RefSymmSCMatrix
beta_density()1042 Wavefunction::beta_density()
1043 {
1044 if (!spin_polarized()) {
1045 RefSymmSCMatrix result = density().copy();
1046 result.scale(0.5);
1047 return result;
1048 }
1049 ExEnv::errn() << class_name() << "::beta_density not implemented" << endl;
1050 abort();
1051 return 0;
1052 }
1053
1054 RefSymmSCMatrix
alpha_ao_density()1055 Wavefunction::alpha_ao_density()
1056 {
1057 return integral()->petite_list()->to_AO_basis(alpha_density());
1058 }
1059
1060 RefSymmSCMatrix
beta_ao_density()1061 Wavefunction::beta_ao_density()
1062 {
1063 return integral()->petite_list()->to_AO_basis(beta_density());
1064 }
1065
1066 void
obsolete()1067 Wavefunction::obsolete()
1068 {
1069 orthog_ = 0;
1070
1071 MolecularEnergy::obsolete();
1072 }
1073
1074 void
copy_orthog_info(const Ref<Wavefunction> & wfn)1075 Wavefunction::copy_orthog_info(const Ref<Wavefunction>&wfn)
1076 {
1077 if (orthog_.nonnull()) {
1078 ExEnv::errn() << "WARNING: Wavefunction: orthogonalization info changing"
1079 << endl;
1080 }
1081 if (wfn->orthog_.null())
1082 wfn->init_orthog();
1083 orthog_ = wfn->orthog_->copy();
1084 }
1085
1086 void
init_orthog()1087 Wavefunction::init_orthog()
1088 {
1089 orthog_ = new OverlapOrthog(
1090 orthog_method_,
1091 overlap(),
1092 basis_matrixkit(),
1093 lindep_tol_,
1094 debug_
1095 );
1096 }
1097
1098 double
min_orthog_res()1099 Wavefunction::min_orthog_res()
1100 {
1101 return orthog_->min_orthog_res();
1102 }
1103
1104 double
max_orthog_res()1105 Wavefunction::max_orthog_res()
1106 {
1107 if (orthog_.null()) init_orthog();
1108 return orthog_->max_orthog_res();
1109 }
1110
1111 OverlapOrthog::OrthogMethod
orthog_method() const1112 Wavefunction::orthog_method() const
1113 {
1114 return orthog_method_;
1115 }
1116
1117 void
set_orthog_method(const OverlapOrthog::OrthogMethod & omethod)1118 Wavefunction::set_orthog_method(const OverlapOrthog::OrthogMethod& omethod)
1119 {
1120 if (orthog_method_ != omethod) {
1121 orthog_method_ = omethod;
1122 init_orthog();
1123 obsolete();
1124 }
1125 }
1126
1127 double
lindep_tol() const1128 Wavefunction::lindep_tol() const
1129 {
1130 return lindep_tol_;
1131 }
1132
1133 void
set_lindep_tol(double lindep_tol)1134 Wavefunction::set_lindep_tol(double lindep_tol)
1135 {
1136 if (lindep_tol_ != lindep_tol) {
1137 lindep_tol_ = lindep_tol;
1138 init_orthog();
1139 obsolete();
1140 }
1141 }
1142
1143 void
scale_atom_basis_coef()1144 Wavefunction::scale_atom_basis_coef()
1145 {
1146 std::vector<int> i_cd;
1147 for (int i=0; i<atom_basis_->ncenter(); i++) {
1148 if (atom_basis_->nshell_on_center(i) > 0) i_cd.push_back(i);
1149 }
1150
1151 if (atom_basis_->max_angular_momentum() > 0) {
1152 // Only s functions will preserve the full symmetry.
1153 // Can only normalize functions that don't integrate to zero.
1154 throw std::runtime_error("ChargeDistInt: max am larger than 0");
1155 }
1156
1157 int coef_offset = 0;
1158 int icoef = 0;
1159 for (int ii=0; ii<i_cd.size(); ii++) {
1160 int i = i_cd[ii];
1161 int nshell = atom_basis_->nshell_on_center(i);
1162 int nfunc = 0;
1163 if (nshell > 0) {
1164 double raw_charge = 0.0;
1165 for (int jj=0, icoef=coef_offset; jj<nshell; jj++) {
1166 int j = atom_basis_->shell_on_center(i,jj);
1167 const GaussianShell &shell = atom_basis_->shell(j);
1168 // loop over the contractions
1169 // the number of functions in each contraction is 1
1170 // since only s functions are allowed
1171 for (int k=0; k<shell.ncontraction(); k++, icoef++) {
1172 for (int l=0; l<shell.nprimitive(); l++) {
1173 double alpha = shell.exponent(l);
1174 double con_coef = shell.coefficient_unnorm(k,l);
1175 // The exponent is halved because the normalization
1176 // formula is for the s function squared.
1177 double integral = ::pow(M_PI/alpha,1.5);
1178 raw_charge += atom_basis_coef_[icoef] * con_coef * integral;
1179 // std::cout << "con_coef = " << con_coef
1180 // << " integral = " << integral
1181 // << std::endl;
1182 }
1183 }
1184 nfunc += shell.ncontraction();
1185 }
1186 double charge = atom_basis_->molecule()->charge(i);
1187 double correction = charge/raw_charge;
1188 for (int icoef=coef_offset; icoef<coef_offset+nfunc; icoef++) {
1189 atom_basis_coef_[icoef] = correction * atom_basis_coef_[icoef];
1190 }
1191 }
1192 coef_offset += nshell;
1193 }
1194 }
1195
1196 /////////////////////////////////////////////////////////////////////////////
1197
1198 // Local Variables:
1199 // mode: c++
1200 // c-file-style: "ETS"
1201 // End:
1202