1 //
2 // eht.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 #include <util/misc/formio.h>
29 #include <util/state/stateio.h>
30
31 #include <math/scmat/blocked.h>
32
33 #include <chemistry/qc/wfn/eht.h>
34 #include <chemistry/qc/basis/integral.h>
35 #include <chemistry/qc/basis/petite.h>
36
37 using namespace std;
38 using namespace sc;
39
40 /////////////////////////////////////////////////////////////////////////
41
42 static ClassDesc ExtendedHuckelWfn_cd(
43 typeid(ExtendedHuckelWfn),"ExtendedHuckelWfn",1,"public OneBodyWavefunction",
44 0, create<ExtendedHuckelWfn>, create<ExtendedHuckelWfn>);
45
ExtendedHuckelWfn(StateIn & s)46 ExtendedHuckelWfn::ExtendedHuckelWfn(StateIn& s) :
47 SavableState(s),
48 OneBodyWavefunction(s)
49 {
50 s.get(nirrep_);
51 s.get(docc_);
52 s.get(socc_);
53 s.get(user_occ_);
54 s.get(total_charge_);
55 }
56
ExtendedHuckelWfn(const Ref<KeyVal> & keyval)57 ExtendedHuckelWfn::ExtendedHuckelWfn(const Ref<KeyVal>&keyval):
58 OneBodyWavefunction(keyval)
59 {
60 CharacterTable ct = molecule()->point_group()->char_table();
61
62 nirrep_ = ct.ncomp();
63 docc_ = new int[nirrep_];
64 socc_ = new int[nirrep_];
65
66 user_occ_ = 0;
67 total_charge_ = keyval->intvalue("total_charge");
68
69 int nuclear_charge = int(molecule()->nuclear_charge());
70 int computed_charge = nuclear_charge;
71
72 for (int i=0; i < nirrep_; i++) {
73 docc_[i]=0;
74 socc_[i]=0;
75
76 if (keyval->exists("docc",i)) {
77 docc_[i] = keyval->intvalue("docc",i);
78 computed_charge -= 2*docc_[i];
79 user_occ_ = 1;
80 }
81 if (keyval->exists("socc",i)) {
82 socc_[i] = keyval->intvalue("socc",i);
83 computed_charge -= 1*socc_[i];
84 user_occ_ = 1;
85 }
86 }
87 if (!keyval->exists("total_charge")) {
88 if (user_occ_) total_charge_ = computed_charge;
89 else total_charge_ = 0;
90 }
91 else if (total_charge_ != computed_charge && user_occ_) {
92 ExEnv::err0() << indent
93 << "ERROR: ExtendedHuckelWfn: total_charge != computed_charge"
94 << endl;
95 abort();
96 }
97 if (total_charge_ > nuclear_charge) {
98 ExEnv::err0() << indent
99 << "ERROR: ExtendedHuckelWfn: total_charge > nuclear_charge"
100 << endl;
101 abort();
102 }
103 }
104
~ExtendedHuckelWfn()105 ExtendedHuckelWfn::~ExtendedHuckelWfn()
106 {
107 delete[] docc_;
108 delete[] socc_;
109 }
110
111 void
save_data_state(StateOut & s)112 ExtendedHuckelWfn::save_data_state(StateOut&s)
113 {
114 OneBodyWavefunction::save_data_state(s);
115 s.put(nirrep_);
116 s.put(docc_,nirrep_);
117 s.put(socc_,nirrep_);
118 s.put(user_occ_);
119 s.put(total_charge_);
120 }
121
122 RefSymmSCMatrix
h_eht_oso()123 ExtendedHuckelWfn::h_eht_oso()
124 {
125 Ref<PetiteList> pl = integral()->petite_list();
126
127 // Compute H in the AO basis
128 double K = 1.75;
129 Ref<AtomInfo> atominfo = molecule()->atominfo();
130 RefSymmSCMatrix h_ao = pl->to_AO_basis(overlap());
131 int natom = basis()->ncenter();
132 int funcoff1 = 0;
133 for (int atom1=0; atom1<natom; atom1++) {
134 int nbasis1 = basis()->nbasis_on_center(atom1);
135 double I1 = atominfo->ip(molecule()->Z(atom1));
136 int funcoff2 = 0;
137 for (int atom2=0; atom2<=atom1; atom2++) {
138 int nbasis2 = basis()->nbasis_on_center(atom2);
139 double I2 = atominfo->ip(molecule()->Z(atom2));
140 for (int func1=0; func1<nbasis1; func1++) {
141 if (atom1 == atom2) nbasis2 = func1 + 1;
142 for (int func2=0; func2<nbasis2; func2++) {
143 int i1 = funcoff1+func1;
144 int i2 = funcoff2+func2;
145 double val = h_ao(i1,i2);
146 if (atom1 == atom2 && func1 == func2) {
147 // The overlap integral is not a part of the diagonal
148 // element in standard EHT formulae. It is here though,
149 // since basis functions are not always normalized (some
150 // d shell components for example).
151 val *= -I1;
152 }
153 else {
154 val *= -0.5*K*(I1+I2);
155 }
156 h_ao(i1,i2) = val;
157 }
158 }
159 funcoff2 += nbasis2;
160 }
161 funcoff1 += nbasis1;
162 }
163
164 if (debug_ > 1) h_ao.print("h in the AO basis");
165
166 // Compute H in the SO basis
167 RefSymmSCMatrix h_so = pl->to_SO_basis(h_ao);
168
169 if (debug_ > 1) {
170 pl->to_AO_basis(overlap()).print("S in AO basis");
171 overlap().print("S in SO basis");
172 pl->aotoso().print("AO to SO transform");
173 h_so.print("h in the SO basis");
174 }
175
176 // Compute H in the OSO basis
177 RefSymmSCMatrix h_oso(oso_dimension(), basis_matrixkit());
178 h_oso->assign(0.0);
179 h_oso->accumulate_transform(so_to_orthog_so(),h_so);
180
181 return h_oso;
182 }
183
184 RefSCMatrix
oso_eigenvectors()185 ExtendedHuckelWfn::oso_eigenvectors()
186 {
187 if (!oso_eigenvectors_.computed() || !eigenvalues_.computed()) {
188 RefSymmSCMatrix h_oso = h_eht_oso();
189
190 if (debug_ > 1) {
191 h_oso.print("h in ortho SO basis");
192 }
193
194 RefSCMatrix vec(oso_dimension(), oso_dimension(), basis_matrixkit());
195 RefDiagSCMatrix val(oso_dimension(), basis_matrixkit());
196
197 h_oso.diagonalize(val,vec);
198
199 if (debug_ > 1) {
200 val.print("h eigenvalues in ortho SO basis");
201 vec.print("h eigenvectors in ortho SO basis");
202 }
203 oso_eigenvectors_=vec;
204 oso_eigenvectors_.computed() = 1;
205
206 eigenvalues_ = val;
207 eigenvalues_.computed() = 1;
208
209 if (!user_occ_) {
210 int nelectron = int(molecule()->nuclear_charge()) - total_charge_;
211 int ndocc = nelectron/2;
212 int nsocc = nelectron%2;
213 fill_occ(val, ndocc, docc_, nsocc, socc_);
214
215 ExEnv::out0() << indent << "docc = [";
216 for (int i=0; i<nirrep_; i++) ExEnv::out0() << " " << docc_[i];
217 ExEnv::out0() << "]" << endl;
218
219 ExEnv::out0() << indent << "socc = [";
220 for (int i=0; i<nirrep_; i++) ExEnv::out0() << " " << socc_[i];
221 ExEnv::out0() << "]" << endl;
222 }
223 }
224
225 return oso_eigenvectors_.result_noupdate();
226 }
227
228 RefDiagSCMatrix
eigenvalues()229 ExtendedHuckelWfn::eigenvalues()
230 {
231 if (!eigenvalues_.computed()) {
232 oso_eigenvectors();
233 }
234
235 return eigenvalues_.result_noupdate();
236 }
237
238 RefSymmSCMatrix
density()239 ExtendedHuckelWfn::density()
240 {
241 if (!density_.computed()) {
242 // Make sure the eigenvectors are computed before
243 // the MO density is computed, otherwise, docc and
244 // socc might not have been initialized.
245 oso_eigenvectors();
246
247 RefDiagSCMatrix mo_density(oso_dimension(), basis_matrixkit());
248 BlockedDiagSCMatrix *modens
249 = dynamic_cast<BlockedDiagSCMatrix*>(mo_density.pointer());
250 if (!modens) {
251 ExEnv::err0() << indent
252 << "ExtendedHuckelWfn::density: wrong MO matrix kit" << endl;
253 abort();
254 }
255
256 modens->assign(0.0);
257 for (int iblock=0; iblock < modens->nblocks(); iblock++) {
258 RefDiagSCMatrix modens_ib = modens->block(iblock);
259 int i;
260 for (i=0; i < docc_[iblock]; i++)
261 modens_ib->set_element(i, 2.0);
262 for ( ; i < docc_[iblock]+socc_[iblock]; i++)
263 modens_ib->set_element(i, 1.0);
264 }
265
266 if (debug_ > 1) mo_density.print("MO Density");
267
268 RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
269 dens->assign(0.0);
270 dens->accumulate_transform(so_to_orthog_so().t() * mo_to_orthog_so(),
271 mo_density);
272
273 if (debug_ > 1) {
274 mo_density.print("MO Density");
275 dens.print("SO Density");
276 ExEnv::out0()
277 << indent << "Nelectron(MO) = " << mo_density.trace()
278 << endl
279 << indent << "Nelectron(SO) = "
280 << (overlap()*dens).trace()
281 << endl;
282 }
283
284 density_ = dens;
285 density_.computed() = 1;
286 }
287
288 return density_.result_noupdate();
289 }
290
291 double
occupation(int ir,int i)292 ExtendedHuckelWfn::occupation(int ir, int i)
293 {
294 if (i < docc_[ir])
295 return 2.0;
296 else if (i < docc_[ir]+socc_[ir])
297 return 1.0;
298 else
299 return 0.0;
300 }
301
302 int
spin_polarized()303 ExtendedHuckelWfn::spin_polarized()
304 {
305 return 0;
306 }
307
308 int
spin_unrestricted()309 ExtendedHuckelWfn::spin_unrestricted()
310 {
311 return 0;
312 }
313
314 void
compute()315 ExtendedHuckelWfn::compute()
316 {
317 double e = (density()*core_hamiltonian()).trace();
318 set_energy(e);
319 set_actual_value_accuracy(desired_value_accuracy());
320 return;
321 }
322
323 int
value_implemented() const324 ExtendedHuckelWfn::value_implemented() const
325 {
326 return 1;
327 }
328
329 void
fill_occ(const RefDiagSCMatrix & evals,int ndocc,int * docc,int nsocc,int * socc)330 ExtendedHuckelWfn::fill_occ(const RefDiagSCMatrix &evals,int ndocc,int *docc,
331 int nsocc, int *socc)
332 {
333 BlockedDiagSCMatrix *bval
334 = require_dynamic_cast<BlockedDiagSCMatrix*>(evals.pointer(),
335 "ExtendedHuckelWfn: getting occupations");
336 int nblock = bval->nblocks();
337 if (nblock != nirrep_) {
338 ExEnv::errn() << "ERROR: ExtendedHuckelWfn: fill_occ: nblock != nirrep" << endl
339 << " nblock = " << nblock << endl
340 << " nirrep = " << nirrep_ << endl;
341 abort();
342 }
343 memset(docc,0,sizeof(docc[0])*nblock);
344 memset(socc,0,sizeof(socc[0])*nblock);
345 for (int i=0; i<ndocc; i++) {
346 double lowest;
347 int lowest_j = -1;
348 for (int j=0; j<nblock; j++) {
349 RefDiagSCMatrix block = bval->block(j);
350 if (block.null()) continue;
351 double current = block->get_element(docc[j]);
352 if (lowest_j < 0 || lowest > current) {
353 lowest = current;
354 lowest_j = j;
355 }
356 }
357 docc[lowest_j]++;
358 }
359 for (int i=0; i<nsocc; i++) {
360 double lowest;
361 int lowest_j = -1;
362 for (int j=0; j<nblock; j++) {
363 double current = bval->block(j)->get_element(docc[j]+socc[j]);
364 if (lowest_j < 0 || lowest > current) {
365 lowest = current;
366 lowest_j = j;
367 }
368 }
369 socc[lowest_j]++;
370 }
371 }
372
373 /////////////////////////////////////////////////////////////////////////////
374
375 // Local Variables:
376 // mode: c++
377 // c-file-style: "ETS"
378 // End:
379