1 //
2 // nao.cc
3 //
4 // Copyright (C) 1997 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 <chemistry/qc/wfn/wfn.h>
30 #include <chemistry/qc/basis/petite.h>
31 #include <chemistry/qc/basis/transform.h>
32
33 using namespace std;
34 using namespace sc;
35
36 #undef DEBUG
37
38 namespace sc {
39 static RefSCMatrix
operator *(const RefDiagSCMatrix & d,const RefSymmSCMatrix & s)40 operator *(const RefDiagSCMatrix &d, const RefSymmSCMatrix &s)
41 {
42 RefSCMatrix ret(s.dim(), s.dim(), s.kit());
43 int n = s.dim()->n();
44 for (int i=0; i<n; i++) {
45 for (int j=0; j<n; j++) {
46 ret.set_element(i,j, d.get_element(i)*s.get_element(i,j));
47 }
48 }
49 return ret;
50 }
51 }
52
53 static RefSymmSCMatrix
weight_matrix(const RefDiagSCMatrix & d,const RefSymmSCMatrix & s)54 weight_matrix(const RefDiagSCMatrix &d, const RefSymmSCMatrix &s)
55 {
56 RefSymmSCMatrix ret = s.clone();
57 int n = s.dim()->n();
58 for (int i=0; i<n; i++) {
59 for (int j=0; j<=i; j++) {
60 ret.set_element(i,j, s.get_element(i,j)
61 *d.get_element(i)*d.get_element(j));
62 }
63 }
64 return ret;
65 }
66
67 static int
nnmb_atom(int z,int l)68 nnmb_atom(int z, int l)
69 {
70 if (l==0) {
71 if (z <= 2) return 1;
72 else if (z <= 10) return 2;
73 else if (z <= 18) return 3;
74 }
75 else if (l==1) {
76 if (z <= 4) return 0;
77 else if (z <= 12) return 1;
78 else if (z <= 20) return 2;
79 }
80 else if (l==2) {
81 if (z <= 20) return 0;
82 }
83 else if (l==3) {
84 if (z <= 56) return 0;
85 }
86 else {
87 return 0;
88 }
89 ExEnv::errn() << "NAO: z too big" << endl;
90 abort();
91 return 0;
92 }
93
94 static int
nnmb_all_atom(int z,int maxl)95 nnmb_all_atom(int z, int maxl)
96 {
97 int ret = 0;
98 for (int i=0; i<=maxl; i++) {
99 ret += nnmb_atom(z,i) * (2*i+1);
100 }
101 return ret;
102 }
103
104 static RefSymmSCMatrix
mhalf(const RefSymmSCMatrix & S)105 mhalf(const RefSymmSCMatrix &S)
106 {
107 RefSCDimension tdim = S.dim();
108 Ref<SCMatrixKit> kit = S.kit();
109
110 // find a symmetric orthogonalization transform
111 RefSCMatrix trans(tdim,tdim,kit);
112 RefDiagSCMatrix eigval(tdim,kit);
113
114 S.diagonalize(eigval,trans);
115
116 Ref<SCElementOp> squareroot = new SCElementSquareRoot;
117 eigval.element_op(squareroot);
118
119 Ref<SCElementOp> invert = new SCElementInvert(1.0e-12);
120 eigval.element_op(invert);
121
122 RefSymmSCMatrix OL(tdim,kit);
123 OL.assign(0.0);
124 // OL = trans * eigval * trans.t();
125 OL.accumulate_transform(trans, eigval);
126 return OL;
127 }
128
129 static void
delete_partition_info(int natom,int * maxam_on_atom,int ** nam_on_atom,int *** amoff_on_atom)130 delete_partition_info(int natom, int *maxam_on_atom,
131 int **nam_on_atom, int ***amoff_on_atom)
132 {
133 int i, j;
134 for (i=0; i<natom; i++) {
135 for (j=0; j<=maxam_on_atom[i]; j++) {
136 delete[] amoff_on_atom[i][j];
137 }
138 delete[] nam_on_atom[i];
139 delete[] amoff_on_atom[i];
140 }
141 delete[] maxam_on_atom;
142 delete[] nam_on_atom;
143 delete[] amoff_on_atom;
144 }
145
146 #ifdef DEBUG
147 static double
ttrace(const RefSCMatrix & N,const RefSymmSCMatrix & P,const RefSymmSCMatrix & S)148 ttrace(const RefSCMatrix &N,
149 const RefSymmSCMatrix &P,
150 const RefSymmSCMatrix &S)
151 {
152 RefSCMatrix Nt = N.t();
153 RefSymmSCMatrix Pt = P.clone();
154 Pt.assign(0.0);
155 Pt.accumulate_transform(Nt, P);
156 RefSymmSCMatrix St = S.clone();
157 St.assign(0.0);
158 St.accumulate_transform(Nt, S);
159 return (mhalf(St)*Pt*mhalf(St)).trace();
160 }
161
162 // for N giving an orthonormal basis
163 static double
ttrace(const RefSCMatrix & N,const RefSymmSCMatrix & P)164 ttrace(const RefSCMatrix &N,
165 const RefSymmSCMatrix &P)
166 {
167 RefSCMatrix Nt = N.t();
168 RefSymmSCMatrix Pt = P.clone();
169 Pt.assign(0.0);
170 Pt.accumulate_transform(Nt, P);
171 return Pt.trace();
172 }
173 #endif
174
175 static RefSCMatrix
assemble(const RefSCDimension dim,const RefSCMatrix & Nm,int * Nm_map,const RefSCMatrix & Nr1,int * Nr1_map,const RefSCMatrix & Nr2=0,int * Nr2_map=0)176 assemble(const RefSCDimension dim,
177 const RefSCMatrix &Nm, int *Nm_map,
178 const RefSCMatrix &Nr1, int *Nr1_map,
179 const RefSCMatrix &Nr2 = 0, int *Nr2_map = 0)
180 {
181 int nnmb = Nm.ncol();
182 int nr1 = Nr1.ncol();
183 int nr2 = (Nr2.null()?0:Nr2.ncol());
184 int nb = dim.n();
185 if (nb != nnmb + nr1 + nr2) {
186 ExEnv::errn() << "assemble: dim mismatch" << endl;
187 abort();
188 }
189 RefSCMatrix N(Nm.rowdim(), Nm.rowdim(), Nm.kit());
190 // collect Nm, Nr1, and Nr2 back into N
191 int i;
192 for (i=0; i<nnmb; i++) {
193 if (Nm_map[i] < 0 || Nm_map[i] >= nb) {
194 ExEnv::errn() << "assemble: bad Nm_map" << endl;
195 abort();
196 }
197 N.assign_column(Nm.get_column(i), Nm_map[i]);
198 }
199 for (i=0; i<nr1; i++) {
200 if (Nr1_map[i] < 0 || Nr1_map[i] >= nb) {
201 ExEnv::errn() << "assemble: bad Nr1_map" << endl;
202 abort();
203 }
204 N.assign_column(Nr1.get_column(i), Nr1_map[i]);
205 }
206 for (i=0; i<nr2; i++) {
207 if (Nr2_map[i] < 0 || Nr2_map[i] >= nb) {
208 ExEnv::errn() << "assemble: bad Nr2_map" << endl;
209 abort();
210 }
211 N.assign_column(Nr2.get_column(i), Nr2_map[i]);
212 }
213 return N;
214 }
215
216 // form symmetry average NAO for each atom
217 static void
form_nao(const RefSymmSCMatrix & P,const RefSymmSCMatrix & S,const RefSCMatrix & N,const RefDiagSCMatrix & W,int natom,int * maxam_on_atom,int ** nam_on_atom,int *** amoff_on_atom,const Ref<SCMatrixKit> & kit)218 form_nao(const RefSymmSCMatrix &P, const RefSymmSCMatrix &S,
219 const RefSCMatrix &N, const RefDiagSCMatrix &W, int natom,
220 int *maxam_on_atom, int **nam_on_atom, int ***amoff_on_atom,
221 const Ref<SCMatrixKit>& kit)
222 {
223 int i,j,k,l,m;
224
225 N.assign(0.0);
226 W.assign(0.0);
227
228 for (i=0; i<natom; i++) {
229 for (j=0; j<=maxam_on_atom[i]; j++) {
230 int nfunc = 2*j + 1;
231 double oonfunc = 1.0/nfunc;
232 int nt = nam_on_atom[i][j];
233 RefSCDimension tdim(new SCDimension(nt));
234 RefSymmSCMatrix Pt(tdim, kit);
235 RefSymmSCMatrix St(tdim, kit);
236 Pt.assign(0.0);
237 St.assign(0.0);
238 for (k=0; k<nt; k++) {
239 for (l=0; l<nt; l++) {
240 double Stmp = 0.0;
241 double Ptmp = 0.0;
242 for (m=0; m<nfunc; m++) {
243 int ii = amoff_on_atom[i][j][k] + m;
244 int jj = amoff_on_atom[i][j][l] + m;
245 Stmp += S.get_element(ii,jj);
246 Ptmp += P.get_element(ii,jj);
247 }
248 St.set_element(k,l,Stmp*oonfunc);
249 Pt.set_element(k,l,Ptmp*oonfunc);
250 }
251 }
252 // find a symmetric orthogonalization transform
253 RefSymmSCMatrix OL = mhalf(St);
254
255 // transform Pt to the orthogonal basis
256 RefSymmSCMatrix PtL(tdim,kit);
257 PtL.assign(0.0);
258 PtL.accumulate_transform(OL, Pt);
259
260 // diagonalize PtL
261 RefSCMatrix trans(tdim,tdim,kit);
262 RefDiagSCMatrix eigval(tdim,kit);
263 PtL.diagonalize(eigval, trans);
264
265 // transform trans to the nonortho basis
266 trans = OL * trans;
267
268 # ifdef DEBUG
269 eigval.print("eigval");
270 # endif
271 // fill in the elements of W
272 for (k=0; k<nt; k++) {
273 // the eigenvalues come out backwards: reverse them
274 int krev = nt-k-1;
275 double elem = eigval.get_element(krev);
276 for (m=0; m<nfunc; m++) {
277 int ii = amoff_on_atom[i][j][k] + m;
278 # ifdef DEBUG
279 ExEnv::outn().form("W(%2d) = %12.8f\n", ii, elem);
280 # endif
281 W.set_element(ii, elem);
282 }
283 }
284
285 // fill in the elements of N
286 for (k=0; k<nt; k++) {
287 for (l=0; l<nt; l++) {
288 // the eigenvalues come out backwards: reverse them
289 int lrev = nt-l-1;
290 double elem = trans.get_element(k,lrev);
291 for (m=0; m<nfunc; m++) {
292 int ii = amoff_on_atom[i][j][k] + m;
293 int jj = amoff_on_atom[i][j][l] + m;
294 N.set_element(ii,jj, elem);
295 }
296 }
297 }
298 }
299 }
300 }
301
302 // From "Natural Population Analysis", Alan E. Reed, Robert B. Weinstock,
303 // Frank Weinhold, JCP, 83 (1985), p 735.
304 RefSCMatrix
nao(double * atom_charges)305 Wavefunction::nao(double *atom_charges)
306 {
307
308 Ref<GaussianBasisSet> b = basis();
309 Ref<PetiteList> pl = integral()->petite_list();
310
311 // compute S, the ao basis overlap
312 RefSymmSCMatrix blockedS = pl->to_AO_basis(overlap());
313 RefSymmSCMatrix S
314 = dynamic_cast<BlockedSymmSCMatrix*>(blockedS.pointer())->block(0);
315 blockedS = 0;
316 # ifdef DEBUG
317 S.print("S");
318 # endif
319
320 // compute P, the ao basis density
321 RefSymmSCMatrix P
322 = dynamic_cast<BlockedSymmSCMatrix*>(ao_density().pointer())->block(0);
323
324 // why? good question.
325 RefSymmSCMatrix Ptmp = P->clone();
326 Ptmp.assign(0.0);
327 Ptmp->accumulate_transform(S, P);
328 # ifdef DEBUG
329 P.print("P");
330 ExEnv::out0() << "nelec = " << (mhalf(S) * Ptmp * mhalf(S)).trace() << endl;
331 ExEnv::out0() << "nelec(2) = " << (P * S).trace() << endl;
332 # endif
333 P = Ptmp;
334 Ptmp = 0;
335
336 int i,j,k,l;
337 int nb = b->nbasis();
338 int nsh = b->nshell();
339 int natom = molecule()->natom();
340
341 # ifdef DEBUG
342 ExEnv::out0() << "nb = " << nb << endl;
343 ExEnv::out0() << "nsh = " << nsh << endl;
344 ExEnv::out0() << "natom = " << natom << endl;
345 # endif
346
347 // Step 2a. Transform to solid harmonics.
348 // -- for now program will abort if basis does not use only S.H and cart d.
349 RefSCDimension aodim = P.dim();
350 RefSCMatrix Tdfg(aodim, aodim, matrixkit());
351 Tdfg->unit();
352 for (i=0; i<nsh; i++) {
353 const GaussianShell &shell = b->shell(i);
354 int off = b->shell_to_function(i);
355 for (j=0; j<shell.ncontraction(); j++) {
356 if (shell.am(j) == 2 && ! shell.is_pure(j)) {
357 for (k=0; k<6; k++) {
358 for (l=0; l<6; l++) {
359 Tdfg.set_element(off+k,off+l,0.0);
360 }
361 }
362 // this will put the s function first and the d second
363 // first grab the s function
364 SphericalTransformIter *sti;
365 sti = integral()->new_spherical_transform_iter(2,0,0);
366 for (sti->begin(); sti->ready(); sti->next()) {
367 Tdfg->set_element(off + sti->pureindex(),
368 off + sti->cartindex(),
369 sti->coef());
370 }
371 delete sti;
372 // now for the pure d part of the cartesian d shell
373 sti = integral()->new_spherical_transform_iter(2,0,2);
374 for (sti->begin(); sti->ready(); sti->next()) {
375 Tdfg->set_element(off + sti->pureindex() + 1,
376 off + sti->cartindex(),
377 sti->coef());
378 }
379 delete sti;
380 }
381 else if (shell.am(j) > 2 && ! shell.is_pure(j)) {
382 ExEnv::errn() << "NAOs can only be computed for puream if am > 2" << endl;
383 abort();
384 }
385 off += shell.nfunction(j);
386 }
387 }
388
389 // Tdfg should already be orthogonal, normalize them
390 // RefSCMatrix Tdfgo = Tdfg*Tdfg.t();
391 // RefDiagSCMatrix Tdfg_norm(Tdfg.rowdim(), matrixkit());
392 // for (i=0; i<nb; i++) {
393 // double o = Tdfgo.get_element(i,i);
394 // Tdfg_norm.set_element(i,1.0/sqrt(o));
395 // }
396 // Tdfgo = 0;
397 // Tdfg = Tdfg_norm * Tdfg;
398
399 # ifdef DEBUG
400 Tdfg.print("Tdfg");
401 (Tdfg.t() * Tdfg).print("Tdfg.t() * Tdfg");
402 (Tdfg * Tdfg.t()).print("Tdfg * Tdfg.t()");
403 # endif
404
405 RefSymmSCMatrix Pdfg(aodim, matrixkit());
406 // Pdfp = Tdfp.t() * P * Tdfp
407 Pdfg.assign(0.0); Pdfg.accumulate_transform(Tdfg, P);
408 RefSymmSCMatrix Sdfg(aodim, matrixkit());
409 // Sdfp = Tdfp.t() * S * Tdfp
410 Sdfg.assign(0.0); Sdfg.accumulate_transform(Tdfg, S);
411 # ifdef DEBUG
412 ExEnv::out0() << "nelec = " << (mhalf(Sdfg) * Pdfg * mhalf(Sdfg)).trace() << endl;
413 # endif
414
415 // Step 2b. Partitioning and symmetry averaging of P and S
416 // Partitioning:
417 int *maxam_on_atom = new int[natom];
418 int **nam_on_atom = new int*[natom];
419 int ***amoff_on_atom = new int**[natom];
420 int maxam = -1;
421 for (i=0; i<natom; i++) {
422 maxam_on_atom[i] = -1;
423 for (j=0; j<b->nshell_on_center(i); j++) {
424 GaussianShell &shell = b->shell(i,j);
425 int maxam_on_shell = shell.max_angular_momentum();
426 if (maxam_on_atom[i] < maxam_on_shell)
427 maxam_on_atom[i] = maxam_on_shell;
428 }
429 if (maxam_on_atom[i] > maxam) maxam = maxam_on_atom[i];
430 nam_on_atom[i] = new int[maxam_on_atom[i]+1];
431 for (j=0; j<=maxam_on_atom[i]; j++) {
432 nam_on_atom[i][j] = 0;
433 for (k=0; k<b->nshell_on_center(i); k++) {
434 GaussianShell &shell = b->shell(i,k);
435 for (l=0; l<shell.ncontraction(); l++) {
436 if (shell.am(l) == j) nam_on_atom[i][j]++;
437 if (shell.am(l) == 2 && !shell.is_pure(l) && j == 0) {
438 // the s component of a cartesian d
439 nam_on_atom[i][0]++;
440 }
441 }
442 }
443 }
444 amoff_on_atom[i] = new int*[maxam_on_atom[i]+1];
445 for (j=0; j<=maxam_on_atom[i]; j++) {
446 amoff_on_atom[i][j] = new int[nam_on_atom[i][j]];
447 int nam = 0;
448 for (k=0; k<b->nshell_on_center(i); k++) {
449 GaussianShell &shell = b->shell(i,k);
450 int function_offset
451 = b->shell_to_function(b->shell_on_center(i,k));
452 int conoffset = 0;
453 for (l=0; l<shell.ncontraction(); l++) {
454 if (shell.am(l) == j) {
455 amoff_on_atom[i][j][nam]
456 = function_offset + conoffset;
457 if (j == 2 && !shell.is_pure(l)) {
458 // the pure d part of a cartesian d shell is offset
459 amoff_on_atom[i][j][nam]++;
460 }
461 nam++;
462 }
463 if (shell.am(l) == 2 && (!shell.is_pure(l)) && j == 0) {
464 // the s component of a cartesian d
465 amoff_on_atom[i][j][nam]
466 = function_offset + conoffset;
467 nam++;
468 }
469 conoffset += shell.nfunction(l);
470 }
471 }
472 }
473 }
474
475 # ifdef DEBUG
476 ExEnv::out0() << indent << "Basis set partitioning:" << endl;
477 ExEnv::out0() << incindent;
478 for (i=0; i<natom; i++) {
479 ExEnv::out0() << indent << "atom " << i
480 << " maxam = " << maxam_on_atom[i] << endl;
481 ExEnv::out0() << incindent;
482 for (j=0; j<=maxam_on_atom[i]; j++) {
483 ExEnv::out0() << indent << "am = " << j
484 << " n = " << nam_on_atom[i][j] << endl;
485 ExEnv::out0() << incindent;
486 ExEnv::out0() << indent << "offsets =";
487 for (k=0; k<nam_on_atom[i][j]; k++) {
488 ExEnv::out0() << " " << amoff_on_atom[i][j][k];
489 }
490 ExEnv::out0() << endl;
491 ExEnv::out0() << decindent;
492 }
493 ExEnv::out0() << decindent;
494 }
495 ExEnv::out0() << decindent;
496 # endif
497
498 // Symmetry averaging and Step 2c: Formation of pre-NAO's
499 RefSCMatrix N(aodim, aodim, matrixkit());
500 RefDiagSCMatrix W(aodim, matrixkit());
501 form_nao(Pdfg, Sdfg, N, W, natom, maxam_on_atom, nam_on_atom, amoff_on_atom,
502 matrixkit());
503 # ifdef DEBUG
504 N.print("N");
505 W.print("W");
506 ExEnv::out0() << "nelec = " << ttrace(N, Pdfg, Sdfg) << endl;
507 # endif
508
509 // Step 3a: selection of NMB orbitals
510
511 // count the size of nmb
512 int nnmb = 0;
513 for (i=0; i<natom; i++) {
514 nnmb += nnmb_all_atom(molecule()->Z(i),
515 maxam_on_atom[i]);
516 }
517 int nnrb = nb - nnmb;
518
519 # ifdef DEBUG
520 ExEnv::out0() << "nnmb = " << nnmb << endl;
521 ExEnv::out0() << "nnrb = " << nnrb << endl;
522 # endif
523
524 RefSCDimension nmbdim = new SCDimension(nnmb);
525 RefSCDimension nrbdim = new SCDimension(nnrb);
526
527 // split N into the nmb and nrb parts
528 RefSCMatrix Nm(aodim, nmbdim, matrixkit());
529 RefSCMatrix Nr(aodim, nrbdim, matrixkit());
530 RefDiagSCMatrix Wm(nmbdim, matrixkit());
531 RefDiagSCMatrix Wr(nrbdim, matrixkit());
532 int *Nm_map = new int[nnmb];
533 int *Nr_map = new int[nnrb];
534
535 int im = 0;
536 int ir = 0;
537 for (i=0; i<natom; i++) {
538 int z = molecule()->Z(i);
539 for (j=0; j<=maxam_on_atom[i]; j++) {
540 int nnmb_zj = nnmb_atom(z,j);
541 int nt = nam_on_atom[i][j];
542 for (k=0; k<nt; k++) {
543 int iN = amoff_on_atom[i][j][k];
544 if (k<nnmb_zj) {
545 for (l=0; l<(2*j+1); l++) {
546 Nm_map[im] = iN;
547 Wm.set_element(im, W.get_element(iN));
548 Nm.assign_column(N.get_column(iN++),im++);
549 }
550 }
551 else {
552 for (l=0; l<(2*j+1); l++) {
553 Nr_map[ir] = iN;
554 Wr.set_element(ir, W.get_element(iN));
555 Nr.assign_column(N.get_column(iN++),ir++);
556 }
557 }
558 }
559 }
560 }
561 # ifdef DEBUG
562 ExEnv::out0() << "Nmmap:"; for (i=0;i<nnmb;i++) ExEnv::out0()<<" "<<Nm_map[i]; ExEnv::out0()<<endl;
563 ExEnv::out0() << "Nrmap:"; for (i=0;i<nnrb;i++) ExEnv::out0()<<" "<<Nr_map[i]; ExEnv::out0()<<endl;
564 Wm.print("Wm");
565 Wr.print("Wr");
566 Nm.print("Nm");
567 Nr.print("Nr");
568 (Nm.t() * Sdfg * Nr).print("3a Smr");
569 ExEnv::out0() << "nelec = "
570 << ttrace(assemble(aodim,Nm,Nm_map,Nr,Nr_map), Pdfg, Sdfg) << endl;
571 # endif
572
573 // Step 3b: Schmidt interatomic orthogonalization of NRB to NMB orbs
574
575 // orthogonalize the NMB orbs (temporarily, to project them out of NRB)
576 int ii=0;
577 for (i=0; i<nnmb; i++,ii++) {
578 N.assign_column(Nm.get_column(i),ii);
579 }
580 for (i=0; i<nnrb; i++,ii++) {
581 N.assign_column(Nr.get_column(i),ii);
582 }
583 N->schmidt_orthog(Sdfg.pointer(),nnmb);
584
585 RefSCMatrix Nmo = Nm.clone();
586 for (i=0; i<nnmb; i++) {
587 Nmo.assign_column(N.get_column(i),i);
588 }
589 RefSCMatrix OSmr = Nmo.t() * Sdfg * Nr;
590 OSmr.scale(-1.0);
591 Nr.accumulate(Nmo * OSmr);
592 # ifdef DEBUG
593 OSmr.print("OSmr");
594 Nmo.print("Nmo = Nm after temporay orthog");
595 Nr.print("Nr after orthogonalization to NMB");
596 (Nm.t() * Sdfg * Nr).print("3b Smr");
597 # endif
598 Nmo = 0;
599
600 // Step 3c: Restoration of natural character of the NRB
601 // Partitioning:
602 int *r_maxam_on_atom = new int[natom];
603 int **r_nam_on_atom = new int*[natom];
604 int ***r_amoff_on_atom = new int**[natom];
605 int r_offset = 0;
606 for (i=0; i<natom; i++) {
607 int z = molecule()->Z(i);
608 r_maxam_on_atom[i] = maxam_on_atom[i];
609 r_nam_on_atom[i] = new int[r_maxam_on_atom[i]+1];
610 for (j=0; j<=r_maxam_on_atom[i]; j++) {
611 r_nam_on_atom[i][j] = nam_on_atom[i][j] - nnmb_atom(z,j);
612 if (r_nam_on_atom[i][j] < 0) {
613 ExEnv::errn() << "NAO: < 0 rydberg orbitals of a given type" << endl;
614 abort();
615 }
616 }
617 r_amoff_on_atom[i] = new int*[r_maxam_on_atom[i]+1];
618 for (j=0; j<=r_maxam_on_atom[i]; j++) {
619 r_amoff_on_atom[i][j] = new int[r_nam_on_atom[i][j]];
620 for (k=0; k<r_nam_on_atom[i][j]; k++) {
621 r_amoff_on_atom[i][j][k] = r_offset;
622 r_offset += 2*j + 1;
623 }
624 }
625 }
626 RefSymmSCMatrix Pr(nrbdim, matrixkit());
627 // Pr = Nr.t() * Tdfg.t() * P * Tdfg * Nr;
628 Pr.assign(0.0); Pr.accumulate_transform(Nr.t(), Pdfg);
629 RefSymmSCMatrix Sr(nrbdim, matrixkit());
630 // Sr = Nr.t() * Tdfg.t() * S * Tdfg * Nr;
631 Sr.assign(0.0); Sr.accumulate_transform(Nr.t(), Sdfg);
632
633 // Symmetry averaging and restoration of natural character of NRB
634 RefSCMatrix Nrr(nrbdim, nrbdim, matrixkit());
635 form_nao(Pr, Sr, Nrr, Wr,
636 natom, r_maxam_on_atom, r_nam_on_atom, r_amoff_on_atom,
637 matrixkit());
638 Nr = Nr * Nrr;
639 // these are out-of-date
640 Pr = 0; Sr = 0;
641 # ifdef DEBUG
642 Wr.print("Wr after restoring natural character");
643 Nr.print("Nr after restoring natural character");
644 (Nm.t() * Sdfg * Nr).print("3c Smr");
645 # endif
646
647 // Step 4a: Weighted interatomic orthogonalization
648 // nmb part of OW
649 RefSymmSCMatrix Sm(nmbdim, matrixkit());
650 Sm.assign(0.0); Sm.accumulate_transform(Nm.t(), Sdfg);
651 RefSymmSCMatrix SWm = weight_matrix(Wm, Sm);
652 RefSCMatrix OWm = Wm * mhalf(SWm);
653 # ifdef DEBUG
654 Sm.print("Sm before 4a");
655 OWm.print("OWm");
656 (OWm.t() * Sm * OWm).print("Sm after 4a");
657 ExEnv::out0() << "nelec = "
658 << ttrace(assemble(aodim,Nm,Nm_map,Nr,Nr_map), Pdfg, Sdfg) << endl;
659 # endif
660
661 // put OWm into Nm
662 Nm = Nm * OWm;
663
664 # ifdef DEBUG
665 Nm.print("Nm after interatomic orthog");
666 (Nm.t() * Sdfg * Nr).print("4a Smr before r orthog");
667 ExEnv::out0() << "nelec = "
668 << ttrace(assemble(aodim,Nm,Nm_map,Nr,Nr_map), Pdfg, Sdfg)
669 << endl;
670 # endif
671
672 // nrb part of OW
673 // based on Wr, r is split into r1 and r2
674 double tw = 1.0e-4; // the tolerance used for the split
675 int nr1 = 0;
676 int nr2 = 0;
677 for (i=0; i<nnrb; i++) {
678 if (fabs(Wr.get_element(i)) >= tw) nr1++;
679 else nr2++;
680 }
681 RefSCDimension r1dim(new SCDimension(nr1));
682 RefSCDimension r2dim(new SCDimension(nr2));
683 RefSCMatrix Nr1(aodim, r1dim, matrixkit());
684 RefSCMatrix Nr2(aodim, r2dim, matrixkit());
685 RefDiagSCMatrix Wr1(r1dim, matrixkit());
686 int *Nr1_map = new int[nr1];
687 int *Nr2_map = new int[nr2];
688 int ir1 = 0;
689 int ir2 = 0;
690 for (i=0; i<nnrb; i++) {
691 if (fabs(Wr.get_element(i)) >= tw) {
692 Nr1_map[ir1] = Nr_map[i];
693 Wr1.set_element(ir1, Wr.get_element(i));
694 Nr1.assign_column(Nr.get_column(i),ir1++);
695 }
696 else {
697 Nr2_map[ir2] = Nr_map[i];
698 Nr2.assign_column(Nr.get_column(i),ir2++);
699 }
700 }
701 # ifdef DEBUG
702 ExEnv::out0() << "Nr1map:"; for (i=0;i<nr1;i++) ExEnv::out0()<<" "<<Nr1_map[i]; ExEnv::out0()<<endl;
703 ExEnv::out0() << "Nr2map:"; for (i=0;i<nr2;i++) ExEnv::out0()<<" "<<Nr2_map[i]; ExEnv::out0()<<endl;
704 Nr1.print("Nr1");
705 Nr2.print("Nr2");
706 (Nm.t() * Sdfg * Nr1).print("4a Smr1 before r orthog");
707 (Nm.t() * Sdfg * Nr2).print("4a Smr2 before r orthog");
708 ExEnv::out0() << "nelec = "
709 << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
710 << endl;
711 # endif
712
713 // Schmidt orthogonalization of r2 to r1
714 // Collect Nr together again (but in the order: r1, r2)
715 ii=0;
716 for (i=0; i<nr1; i++,ii++) {
717 Nr.assign_column(Nr1.get_column(i),ii);
718 }
719 for (i=0; i<nr2; i++,ii++) {
720 Nr.assign_column(Nr2.get_column(i),ii);
721 }
722 Nr->schmidt_orthog(Sdfg.pointer(),nr1);
723 RefSCMatrix Nr1o = Nr1.copy();
724 for (i=0; i<nr1; i++) {
725 Nr1o.assign_column(Nr.get_column(i), i);
726 }
727 RefSCMatrix Or1r2 = Nr1o.t() * Sdfg * Nr2;
728 Or1r2.scale(-1.0);
729 # ifdef DEBUG
730 (Nm.t() * Sdfg * Nr2).print("4a Smr2 before orthog of r2 to r1");
731 (Nr1.t() * Sdfg * Nr2).print("4a Sr1r2 before orthog of r2 to r1");
732 # endif
733 Nr2.accumulate(Nr1o * Or1r2);
734 # ifdef DEBUG
735 Nr2.print("Nr2 after orthogonalization to r1");
736 (Nm.t() * Sdfg * Nr2).print("4a Smr2 after orthog of r2 to r1");
737 (Nr1.t() * Sdfg * Nr2).print("4a Sr1r2 after orthog of r2 to r1");
738 ExEnv::out0() << "nelec = "
739 << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
740 << endl;
741 # endif
742
743 // weighted symmetric orthog of r1
744 RefSymmSCMatrix Sr1(r1dim, matrixkit());
745 Sr1.assign(0.0); Sr1.accumulate_transform(Nr1.t(), Sdfg);
746 RefSymmSCMatrix SWr1 = weight_matrix(Wr1, Sr1);
747 RefSCMatrix OWr1 = Wr1 * mhalf(SWr1);
748 # ifdef DEBUG
749 OWr1.print("OWr1");
750 (Nr1.t() * Sdfg * Nr1).print("Nr1.t() * Sdfg * Nr1");
751 # endif
752 // Put OWr1 into Nr1
753 Nr1 = Nr1 * OWr1;
754 # ifdef DEBUG
755 Nr1.print("Nr1 after weighted symmetric orthogonalization");
756 ExEnv::out0() << "nelec = "
757 << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
758 << endl;
759 # endif
760
761 // symmetric orthog of r1
762 RefSymmSCMatrix Sr2(r2dim, matrixkit());
763 Sr2.assign(0.0); Sr2.accumulate_transform(Nr2.t(), Sdfg);
764 RefSymmSCMatrix OWr2 = mhalf(Sr2);
765 # ifdef DEBUG
766 OWr2.print("OWr2");
767 # endif
768
769 // Put OWr2 into Nr2
770 Nr2 = Nr2 * OWr2;
771 # ifdef DEBUG
772 Nr2.print("Nr2 after weighted symmetric orthogonalization");
773 ExEnv::out0() << "nelec = "
774 << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
775 << endl;
776 ExEnv::out0() << "nelec(o) = "
777 << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg)
778 << endl;
779 # endif
780
781 // Step 4b. restoration of the natural character of the naos
782
783 // collect Nm, Nr1, and Nr2 back into N
784 N = assemble(aodim, Nm,Nm_map, Nr1,Nr1_map, Nr2,Nr2_map);
785 # ifdef DEBUG
786 N.print("N after 4a");
787 # endif
788
789 // compute the density and overlap in the current basis
790 // N currently has the entire transform, starting from the dfg basis
791 P.assign(0.0);
792 P.accumulate_transform(N.t(), Pdfg);
793 S.assign(0.0);
794 S.accumulate_transform(N.t(), Sdfg);
795 # ifdef DEBUG
796 P.print("P after 4a");
797 S.print("S after 4a");
798 (Nm.t() * Sdfg * Nm).print("4a Sm");
799 (Nr1.t() * Sdfg * Nr1).print("4a Sr1");
800 (Nr2.t() * Sdfg * Nr2).print("4a Sr2");
801 (Nm.t() * Sdfg * Nr1).print("4a Smr1");
802 (Nm.t() * Sdfg * Nr2).print("4a Smr2");
803 (Nr1.t() * Sdfg * Nr2).print("4a Sr1r2");
804 # endif
805
806 RefSCMatrix Nred(aodim, aodim, matrixkit());
807 form_nao(P, S, Nred, W, natom, maxam_on_atom, nam_on_atom, amoff_on_atom,
808 matrixkit());
809 N = N * Nred;
810
811 RefSymmSCMatrix Pfinal(aodim, matrixkit());
812 Pfinal.assign(0.0);
813 Pfinal.accumulate_transform(N.t(), Pdfg);
814 # ifdef DEBUG
815 Nred.print("Nred");
816 N.print("N after 4b");
817 ExEnv::out0() << "nelec = " << ttrace(N, Pdfg, Sdfg) << endl;
818 ExEnv::out0() << "nelec(o) = " << ttrace(N, Pdfg) << endl;
819 Pfinal.print("final P");
820 (N.t() * Sdfg * N).print("final S");
821 ExEnv::out0().form("nelec = trace(final P) = %14.8f", (N.t() * Pdfg * N).trace());
822
823 (mhalf(Sdfg) * Pdfg * mhalf(Sdfg)).print("P in symm orth basis");
824 # endif
825
826 # ifdef DEBUG
827 ExEnv::out0() << "nb = " << nb << endl;
828 ExEnv::out0() << "nnmb = " << nnmb << endl;
829 ExEnv::out0() << "nnrb = " << nnrb << endl;
830 ExEnv::out0() << "nr1 = " << nr1 << endl;
831 ExEnv::out0() << "nr2 = " << nr2 << endl;
832 # endif
833
834 ExEnv::out0() << indent << "Natural Population Analysis:" << endl;
835 ExEnv::out0() << incindent;
836 ExEnv::out0() << indent << " n atom charge ";
837 for (i=0; i<=maxam; i++) {
838 const char *am = "SPDFGH?";
839 int index;
840 if (i>6) index = 6;
841 else index = i;
842 ExEnv::out0() << " ne(" << am[index] << ") ";
843 }
844 ExEnv::out0() << endl;
845 for (i=0; i<natom; i++) {
846 double e = 0.0;
847 for (j=0; j<=maxam_on_atom[i]; j++) {
848 for (k=0; k<nam_on_atom[i][j]; k++) {
849 for (l=0; l<(2*j+1); l++) {
850 e += Pfinal.get_element(amoff_on_atom[i][j][k] + l,
851 amoff_on_atom[i][j][k] + l);
852 }
853 }
854 }
855 std::string symbol(molecule()->atom_symbol(i));
856 ExEnv::out0() << indent
857 << scprintf("%3d %2s % 8.6f",i + 1,
858 symbol.c_str(),
859 double(molecule()->Z(i)) - e);
860 if (atom_charges) {
861 atom_charges[i] = molecule()->Z(i) - e;
862 }
863 for (j=0; j<=maxam_on_atom[i]; j++) {
864 e = 0.0;
865 for (k=0; k<nam_on_atom[i][j]; k++) {
866 for (l=0; l<(2*j+1); l++) {
867 e += Pfinal.get_element(amoff_on_atom[i][j][k] + l,
868 amoff_on_atom[i][j][k] + l);
869 }
870 }
871 ExEnv::out0() << scprintf(" % 8.6f",e);
872 }
873 ExEnv::out0() << endl;
874 }
875 ExEnv::out0() << endl;
876 ExEnv::out0() << decindent;
877
878 delete[] Nm_map;
879 delete[] Nr_map;
880 delete[] Nr1_map;
881 delete[] Nr2_map;
882 delete_partition_info(natom,maxam_on_atom,nam_on_atom,amoff_on_atom);
883 delete_partition_info(natom,r_maxam_on_atom,r_nam_on_atom,r_amoff_on_atom);
884
885 return N;
886 }
887
888 /////////////////////////////////////////////////////////////////////////////
889
890 // Local Variables:
891 // mode: c++
892 // c-file-style: "CLJ"
893 // End:
894