1 /* Copyright (c) 2015 Gerald Knizia
2 *
3 * This file is part of the IboView program (see: http://www.iboview.org)
4 *
5 * IboView is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * IboView is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
16 *
17 * Please see IboView documentation in README.txt for:
18 * -- A list of included external software and their licenses. The included
19 * external software's copyright is not touched by this agreement.
20 * -- Notes on re-distribution and contributions to/further development of
21 * the IboView software
22 */
23
24 #include "Iv.h"
25
26 #include <iostream>
27 #include <boost/format.hpp>
28 #include <algorithm>
29 #include <set>
30
31 #include "CxAlgebra.h"
32 #include "CtMatrix.h"
33 #include "CtTiming.h"
34 #include "IvIao.h"
35 #include "IvDocument.h"
36
37 using namespace ct;
38
39 enum FIaoBasisFlags {
40 IAO_OrthType = 0x03,
41 IAO_OrthNone = 0x00,
42 IAO_OrthSym = 0x01,
43 IAO_OrthZbd = 0x02,
44 IAO_NormalizeInput = 0x04
45 };
46
47
48 // // symmetrically orthogonalize a set of vectors with respect to overlap matrix S.
49 // void SymOrth(FMatrixView Orbs, FMatrixView const S, FMemoryStack &Mem)
50 // {
51 // FStackMatrix
52 // SmhOrb(Orbs.nCols, Orbs.nCols, &Mem),
53 // T1(Orbs.nRows, Orbs.nCols, &Mem);
54 // ChainMxm(SmhOrb, Transpose(Orbs), S, Orbs, Mem);
55 // CalcSmhMatrix(SmhOrb, Mem, FSmhOptions(1e-15,1e-15,0,0));
56 // Mxm(T1, Orbs, SmhOrb);
57 // Move(Orbs, T1);
58 // }
59
60
61
62 // note: CIb must be provided from outside.
MakeIaoBasisNew(FMatrixView CIb,FBasisSet * pMinBasis,FAtomSet const &,FBasisSet const * pOrbBasis,FMatrixView const & COcc_,FMatrixView const & S1,FMatrixView const & S1cd,uint Flags,FMemoryStack & Mem)63 void MakeIaoBasisNew(FMatrixView CIb, FBasisSet *pMinBasis,
64 FAtomSet const &/*Atoms*/, FBasisSet const *pOrbBasis,
65 FMatrixView const &COcc_, FMatrixView const &S1, FMatrixView const &S1cd, uint Flags,
66 FMemoryStack &Mem)
67 {
68 uint
69 nBf = pOrbBasis->nFn(), // number of 'real basis' functions
70 nOcc = COcc_.nCols,
71 nIb = pMinBasis->nFn(); // actual AOs in minimal basis
72 // CIb = MakeStackMatrix(nBf, nIb, Mem);
73 assert(CIb.nRows == nBf && CIb.nCols == nIb);
74
75 FMatrixView
76 COcc;
77 if (0 == (Flags&IAO_NormalizeInput)) {
78 COcc = COcc_;
79 } else {
80 // input has absorbed occupation numbers. Remove this.
81 // (todo: find a better way...)
82 assert(COcc_.nRows == S1.nRows);
83 COcc = MakeStackMatrix(COcc_.nRows, nOcc, Mem);
84 Move(COcc, COcc_);
85 SymOrth(COcc, S1, Mem);
86 // Mem.Alloc(pOccScales, nOcc);
87 // FStackMatrix
88 // SOcc(nOcc, nOcc, &Mem);
89 // ChainMxm(SOcc, Transpose(COcc_), S1, COcc_, Mem);
90 // SOcc.Print(xout, "Occ self overlap");
91 // for (uint iOcc = 0; iOcc < nOcc; ++ iOcc) {
92 // pOccScales[iOcc] = std::sqrt(SOcc(iOcc,iOcc));
93 // double f = 1./pOccScales[iOcc];
94 // for (uint i = 0; i < COcc_.nRows; ++ i)
95 // COcc(i,iOcc) = f * COcc_(i,iOcc);
96 // }
97 }
98
99 {
100 FStackMatrix
101 S2cd(nIb, nIb, &Mem),
102 P12(nBf, nIb, &Mem),
103 COcc2(nIb, nOcc, &Mem),
104 CTil(nIb, nOcc, &Mem),
105 STilCd(nIb, nIb, &Mem),
106 CTil2BarT(nOcc, nIb, &Mem),
107 T4(nBf, nOcc, &Mem);
108 MakeIntMatrix(P12, *pOrbBasis, *pMinBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
109 MakeIntMatrix(S2cd, *pMinBasis, *pMinBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
110
111 // COcc2 = mdot(S21, COcc) # O(N m^2)
112 Mxm(COcc2, Transpose(P12), COcc);
113
114 // P12 = la.solve(S1, S12) # P12 = S1^{-1} S12
115 CholeskySolve(P12, S1cd);
116
117 // CTil = la.solve(S2, COcc2) # O(m^3)
118 CalcCholeskyFactors(S2cd);
119 Move(CTil, COcc2);
120 CholeskySolve(CTil, S2cd);
121
122 // STil = mdot(COcc2.T, CTil) # O(m^3)
123 Mxm(STilCd, Transpose(COcc2), CTil);
124
125 // CTil2Bar = la.solve(STil, CTil.T).T # O(m^3)
126 CalcCholeskyFactors(STilCd);
127 Move(CTil2BarT, Transpose(CTil));
128 CholeskySolve(CTil2BarT, STilCd);
129
130 // T4 = COcc - mdot(P12, CTil2Bar) # O(N m^2)
131 Mxm(T4, P12, Transpose(CTil2BarT), -1.0);
132 Add(T4, COcc);
133
134 // CIb = P12 + mdot(T4, COcc2.T) # O(N m^2)
135 Mxm(CIb, T4, Transpose(COcc2));
136 Add(CIb, P12);
137
138 if ((Flags & IAO_OrthType) == IAO_OrthSym)
139 SymOrth(CIb, S1, Mem);
140 }
141
142 if (COcc.pData != COcc_.pData)
143 Mem.Free(COcc.pData);
144
145 // still on stack: CIb.
146 }
147
148
149
150
151
152
153 namespace ct {
154
155 struct FLocalizeOptions
156 {
157 uint
158 nLocExp; // 2 or 4, exponent p in localization functional: L = sum[iA] n(i,A)^p
159 uint
160 nMaxIt;
161 double
162 ThrLoc;
163 std::ostream
164 *pxout; // if non-0, write output here.
165 uint
166 Verbosity;
167
FLocalizeOptionsct::FLocalizeOptions168 FLocalizeOptions()
169 : nLocExp(4), nMaxIt(2048), ThrLoc(1e-8), pxout(0), Verbosity(1)
170 {}
171 };
172
pow2(double x)173 inline double pow2(double x) { return x*x; }
pow3(double x)174 inline double pow3(double x) { return x*x*x; }
pow4(double x)175 inline double pow4(double x) { double x2 = x*x; return x2*x2; }
176
177 // execute a 2x2 rotation of the vectors pA and pB (both with length nSize)
178 // of angle phi. In-place.
Rot2x2(double * RESTRICT pA,double * RESTRICT pB,size_t nSize,double phi)179 void Rot2x2(double *RESTRICT pA, double *RESTRICT pB, size_t nSize, double phi)
180 {
181 double
182 cs = std::cos(phi),
183 ss = std::sin(phi);
184 for (std::size_t i = 0; i < nSize; ++ i) {
185 double
186 tempA = cs * pA[i] + ss * pB[i],
187 tempB = -ss * pA[i] + cs * pB[i];
188 pA[i] = tempA;
189 pB[i] = tempB;
190 }
191 }
192
193 // localize a set of vectors onto
194 // L = \sum[i,A] \sum[mu \in A] CVec(mu,i)^p -> min.
195 // in a PM-like fashion.
196 // Parameters:
197 // - pGroupOffsets: Array of length nGroups+1.
198 // [pGroupOffsets[A], pGroupOffsets[A+1]) gives the association of
199 // vector components to groups.
200 // Note: This function can be extended to localize on other criteria
201 // than charge, provided a replacement of the calculation of the Cii/Cjj/Cij
202 // matrix elements is provided. Most useful in the local exchange case is likely
203 // the optimization of <i|(r-A)^4|i>. It can also be used as a replacement for
204 // ZBD orthogonalization then (which is equivalent to optimizing on (r-A)^2).
LocalizeVectors(FMatrixView CVec,size_t const * pGroupOffsets,size_t nGroups,FLocalizeOptions const & Opt,FMemoryStack const & Mem)205 void LocalizeVectors(FMatrixView CVec, size_t const *pGroupOffsets, size_t nGroups, FLocalizeOptions const &Opt, FMemoryStack const &Mem)
206 {
207 IR_SUPPRESS_UNUSED_WARNING(Mem);
208
209 size_t
210 iIt;
211 double
212 L = -1., fVar2 = -1.;
213 if (Opt.pxout && Opt.Verbosity >= 2) {
214 *Opt.pxout << " ITER. LOC(Orb) GRADIENT" << std::endl;
215 }
216 for (iIt = 0; iIt < Opt.nMaxIt; ++ iIt) {
217 // calculate the value of the functional.
218 L = 0.;
219 // loop over groups (normally a ``group'' is all functions on an atom)
220 for (size_t iGrp = 0; iGrp < nGroups; ++ iGrp)
221 for (size_t iVec = 0; iVec < CVec.nCols; ++ iVec) {
222 double nA = 0.; // population on the atom.
223 for (size_t iFn = pGroupOffsets[iGrp]; iFn != pGroupOffsets[iGrp+1]; ++ iFn)
224 nA += pow2(CVec(iFn,iVec));
225 L += std::pow(nA, (int)Opt.nLocExp);
226 }
227 L = std::pow(L, 1./double(Opt.nLocExp)); // <- easier to compare different powers that way.
228
229 fVar2 = 0.;
230 // loop over vector pairs.
231 for (size_t iVec = 0; iVec < CVec.nCols; ++ iVec)
232 for (size_t jVec = 0; jVec < iVec; ++ jVec) {
233 double
234 Aij = 0.0, // hessian for 2x2 rotation (with variable atan(phi/4) substituted)
235 Bij = 0.0; // gradient for 2x2 rotation (with variable atan(phi/4) substituted)
236 // loop over groups (normally atoms)
237 for (size_t iGrp = 0; iGrp < nGroups; ++ iGrp) {
238 // calculate the charge matrix elements
239 // Cii = <i|A|i>,
240 // Cjj = <j|A|j>,
241 // and Cij = <i|A|j>.
242 double
243 Cii = 0., Cij = 0., Cjj = 0.;
244 for (size_t iFn = pGroupOffsets[iGrp]; iFn != pGroupOffsets[iGrp+1]; ++ iFn) {
245 Cii += CVec(iFn,iVec) * CVec(iFn,iVec);
246 Cjj += CVec(iFn,jVec) * CVec(iFn,jVec);
247 Cij += CVec(iFn,iVec) * CVec(iFn,jVec);
248 }
249
250 // see Supporting Information for `` Intrinsic atomic orbitals:
251 // An unbiased bridge between quantum theory and chemical concepts''
252 // for a description of what this does.
253 if (Opt.nLocExp == 2 || Opt.nLocExp == 3) {
254 Aij += sqr(Cij) - .25*sqr(Cii - Cjj);
255 Bij += Cij*(Cii - Cjj);
256 } else if (Opt.nLocExp == 4) {
257 Aij += -pow4(Cii) - pow4(Cjj) + 6.*(pow2(Cii) + pow2(Cjj))*pow2(Cij) + pow3(Cii)*Cjj + Cii*pow3(Cjj);
258 Bij += 4.*Cij*(pow3(Cii) - pow3(Cjj));
259 }
260 }
261
262 // non-degenerate rotation? This condition is not quite right,
263 // and in Python it was not required. However, the Fortran atan2
264 // sometimes did mysterious things without it.
265 if (pow2(Aij) + pow2(Bij) > Opt.ThrLoc * 1e-2) {
266 double
267 phi = .25 * std::atan2(Bij,-Aij);
268 // 2x2 rotate vectors iVec and jVec.
269 Rot2x2(&CVec(0,iVec), &CVec(0,jVec), CVec.nRows, phi);
270 fVar2 += pow2(phi);
271 }
272 }
273 fVar2 = std::sqrt(fVar2 / pow2(CVec.nCols));
274
275 if (Opt.pxout && Opt.Verbosity >= 2) {
276 *Opt.pxout << boost::format("%6i %12.6f %8.2e\n") % (1+iIt) % L % fVar2;
277 }
278
279 if (fVar2 < Opt.ThrLoc)
280 break;
281 }
282
283 if (Opt.pxout && Opt.Verbosity >= 1) {
284 if (Opt.Verbosity >= 2)
285 *Opt.pxout << "\n";
286 // *Opt.pxout << boost::format(" Iterative localization: IB/PM, %i iter; Final gradient %8.2e") % iIt % fVar2 << std::endl;
287 *Opt.pxout << boost::format(" Iterative localization: IB/PM, %i iter; Final gradient %8.2e") % iIt % fVar2;
288 Opt.pxout->flush();
289 }
290 }
291
292
293
294 // void LocalizeOrbitals(FBasisSetPtr pBasis, FMatrixView COrb, FLocalizeOptions const &Opt);
295
296
297 } // namespace ct
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314 // // note: CIb is allocated on Mem.
315 // void MakeIaoBasis(FMatrixView &CIb, uint &nIb, FBasisSetPtr &pMinBasis,
316 // FAtomSet const &Atoms, FBasisSet const *pOrbBasis, FMatrixView &COcc, FMemoryStack &Mem)
317 // {
318 // pMinBasis = new FBasisSet(Atoms, BASIS_Guess);
319 // uint
320 // nBf = pOrbBasis->nFn(), // number of 'real basis' functions
321 // nOcc = COcc.nCols;
322 // nIb = pMinBasis->nFn(); // actual AOs in minimal basis
323 // CIb = MakeStackMatrix(nBf, nIb, Mem);
324 //
325 // FStackMatrix
326 // S1(nBf, nBf, &Mem),
327 // S1cd(nBf, nBf, &Mem);
328 // MakeIntMatrix(S1, *pOrbBasis, *pOrbBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
329 // Move(S1cd,S1);
330 // CalcCholeskyFactors(S1cd);
331 //
332 // MakeIaoBasisNew(CIb, &*pMinBasis, Atoms, pOrbBasis, COcc, S1, S1cd, IAO_OrthSym|IAO_NormalizeInput, Mem);
333 // // MakeIaoBasisNew(CIb, nIb, &*pMinBasis, Atoms, pOrbBasis, COcc, S1, S1cd, IAO_OrthSym, Mem);
334 // // SymOrth(CIb, S1, Mem);
335 // }
336
337
CalcMatrixRmsNorm(FMatrixView M)338 double CalcMatrixRmsNorm(FMatrixView M)
339 {
340 double d = 0;
341 for (uint iCol = 0; iCol < M.nCols; ++ iCol)
342 for (uint iRow = 0; iRow < M.nRows; ++ iRow)
343 d += sqr(M(iRow,iCol));
344 return std::sqrt(d);
345 }
346
347 // TODO: merge these back into IvDocument.cpp. It has a bunch of very similar
348 // functions in FOrbital::MakeFullDesc and FOrbital::MakeIaoCharges.
MakeIaoChargesRaw(double * pOut,double fOcc,double const * pIaoCoeffs,ct::FBasisSet * pMinBasis,ct::FAtomSet * pAtoms)349 void MakeIaoChargesRaw(double *pOut, double fOcc, double const *pIaoCoeffs, ct::FBasisSet *pMinBasis, ct::FAtomSet *pAtoms)
350 {
351 size_t
352 iShOf = 0;
353 memset(pOut, 0, sizeof(pOut[0]) * pAtoms->size());
354 for (size_t iSh = 0; iSh < pMinBasis->Shells.size(); ++ iSh) {
355 ct::FBasisShell
356 &Sh = pMinBasis->Shells[iSh];
357 uint
358 nShFn = Sh.nFn();
359 double
360 fChg = 0;
361 for (uint iFn = 0; iFn < nShFn; ++ iFn)
362 fChg += sqr(pIaoCoeffs[iFn + iShOf]);
363 assert(size_t(Sh.iCenter) < pAtoms->size());
364 pOut[Sh.iCenter] += fOcc * fChg;
365 iShOf += nShFn;
366 }
367 assert(iShOf == pMinBasis->nFn());
368 }
369
370 struct FPrintOrbitalOptions
371 {
372 double
373 ThrPrint;
374 size_t
375 nMaxAtoms;
376 bool
377 PrintHeader;
378 size_t
379 iOrbOffset;
FPrintOrbitalOptionsFPrintOrbitalOptions380 FPrintOrbitalOptions()
381 : ThrPrint(0.02), nMaxAtoms(20), PrintHeader(true), iOrbOffset(0)
382 {}
383 };
384
385
PrintOrbitalChargeComposition(FLog & Log,FMatrixView CIbOcc,FMatrixView Ew,ct::FAtomSet * pAtoms,FBasisSet * pMinBasis,std::string const sGrp,double fOcc1,FPrintOrbitalOptions & Opt,ct::FMemoryStack & Mem)386 static void PrintOrbitalChargeComposition(FLog &Log, FMatrixView CIbOcc, FMatrixView Ew, ct::FAtomSet *pAtoms,
387 FBasisSet *pMinBasis, std::string const sGrp, double fOcc1, FPrintOrbitalOptions &Opt,
388 ct::FMemoryStack &Mem)
389 {
390 size_t const
391 nAt = pAtoms->size();
392 // *pAtOffsets = &pMinBasis->pRawBasis->CenterOffsets[0];
393 if (Opt.PrintHeader) {
394 Log.Write(" Summary of localized orbital composition: [THRPRINT={:6.3f}]", Opt.ThrPrint);
395 Log.Write("\n ORB GRP ORB.ENERGY CENTERS/CHARGES");
396 Opt.PrintHeader = false;
397 // ^- that's for appending lines using multiple calls.
398 }
399 for (size_t iOrb = 0; iOrb < CIbOcc.nCols; ++ iOrb) {
400 TMemoryLock<double>
401 Charges(nAt, &Mem);
402 MakeIaoChargesRaw(Charges, fOcc1, &CIbOcc(0,iOrb), pMinBasis, pAtoms);
403 TMemoryLock<size_t>
404 Indices(nAt, &Mem);
405 ct::ArgSort1(&Indices[0], &Charges[0], 1, nAt, true); // last: reverse. Largest first.
406
407 double
408 fChgRest = 0.;
409 size_t
410 nAtomsEmitted = 0;
411 fmt::MemoryWriter
412 out;
413 out.write("{:6} {:>3} {:12.6f} ", iOrb+1, sGrp, Ew[iOrb]);
414 for (size_t ii = 0; ii < nAt; ++ ii) {
415 size_t
416 iAt = Indices[ii];
417 if (iAt >= nAt)
418 continue; // not supposed to happen.
419 double
420 fAtChg = Charges[iAt];
421 if (fAtChg < Opt.ThrPrint || nAtomsEmitted >= Opt.nMaxAtoms) {
422 fChgRest += fAtChg;
423 continue;
424 }
425 ct::FAtom
426 &At = (*pAtoms)[iAt];
427 out.write(" {:>2}{:3} {:6.3f}", At.GetElementName(), (1+iAt), Charges[iAt]);
428 nAtomsEmitted += 1;
429 }
430 if (fChgRest > Opt.ThrPrint)
431 out.write(" (other: {:6.3f})", fChgRest);
432 Log.Write(out.str());
433 }
434 Opt.iOrbOffset += CIbOcc.nCols;
435 }
436
437 // predicate: sorts orbitals first by type, then by occupation, and finally
438 // by energy.
439 struct FSortOrbitalPred
440 {
operator ()FSortOrbitalPred441 bool operator() (FDataSetPtr pA, FDataSetPtr pB) {
442 FOrbital
443 *pOrbA = dynamic_cast<FOrbital*>(pA.get()),
444 *pOrbB = dynamic_cast<FOrbital*>(pB.get());
445 if (pOrbA == 0 && pOrbB == 0)
446 return false;
447 if (pOrbA == 0)
448 return true;
449 if (pOrbB == 0)
450 return false;
451 // so... both are orbitals.
452 if (pOrbA->info.Spin < pOrbB->info.Spin)
453 return true;
454 if (pOrbB->info.Spin < pOrbA->info.Spin)
455 return false;
456 // if (pOrbA->info.fOcc < pOrbB->info.fOcc)
457 // return true;
458 // if (pOrbB->info.fOcc < pOrbA->info.fOcc)
459 // return false;
460 return pOrbA->info.fEnergy < pOrbB->info.fEnergy;
461 }
462 };
463
464
465
RunIaoAnalysis(FLog & Log,FWfOptions * pWfOptions,bool AllowLocalize,ct::FMemoryStack & Mem)466 void FFrame::RunIaoAnalysis(FLog &Log, FWfOptions *pWfOptions, bool AllowLocalize, ct::FMemoryStack &Mem)
467 {
468 if (!HaveOrbitals())
469 // nothing to make coefficients of...
470 return;
471
472 Log.WriteProgramIntro("INTRINSIC BASIS BONDING ANALYSIS (IBBA)");
473 FTimer
474 tIbbaTotal;
475
476 FWfType
477 WfType = GetWfType();
478 if (WfType != WFTYPE_Rhf && WfType != WFTYPE_Uhf) {
479 Log.Write("Sorry, wave function type not supported. Can only do RHF/RKS and UHF/UKS at this moment.");
480 return;
481 }
482
483 using namespace ct;
484 FAtomSet
485 *pAtoms = pGetAtoms();
486
487 TMemoryLock<char>
488 pBaseOfMemory(0, &Mem);
489
490 FBasisSetPtr
491 pBasis,
492 pMinBasis = new FBasisSet(*pAtoms, BASIS_Guess);
493 FMatrixView
494 COccAC, // alpha and closed
495 COccBC; // beta and closed
496 {
497 TArray<FOrbital*>
498 RefOrbitalsAC,
499 RefOrbitalsBC;
500 MakeOrbitalMatrix(COccAC, pBasis, &RefOrbitalsAC, ORBMAT_OccupiedOnly | ORBMAT_AlphaAndClosedOnly, Mem);
501 assert(RefOrbitalsAC.size() == COccAC.nCols);
502 if (WfType == WFTYPE_Uhf) {
503 MakeOrbitalMatrix(COccBC, pBasis, &RefOrbitalsBC, ORBMAT_OccupiedOnly | ORBMAT_BetaAndClosedOnly, Mem);
504 assert(RefOrbitalsBC.size() == COccBC.nCols);
505 } else
506 COccBC = COccAC;
507 }
508
509
510 uint
511 nBf = pBasis->nFn(),
512 nIb = pMinBasis->nFn(),
513 nOccAC = COccAC.nCols,
514 nOccBC = COccBC.nCols;
515
516 if (nOccAC > nIb || nOccBC > nIb) {
517 Log.Write(" More occupied orbitals than valence orbitals. IAO analysis aborted.");
518 return;
519 }
520
521 // form orthogonal IAO vectors.
522 uint
523 IaoFlags = IAO_NormalizeInput;
524 char const
525 *pOrthMethodName = "Unknown";
526 if (pWfOptions->GetAoType() == "IAO (Sym Orth.)") {
527 IaoFlags |= IAO_OrthSym;
528 pOrthMethodName = "Symmetric (Loewdin)";
529 } else if (pWfOptions->GetAoType() == "IAO (ZBD Orth.)") {
530 IaoFlags |= IAO_OrthZbd;
531 pOrthMethodName = "Zero-Bond-Dipole (Laikov)";
532 } else {
533 IaoFlags |= IAO_OrthSym;
534 IvNotify(NOTIFY_Error, IvFmt("Intrinsic basis type '%1' not recognized. Using IAO (Sym Orth.) instead.", pWfOptions->GetAoType()));
535 pOrthMethodName = "Symmetric (Loewdin)";
536 }
537
538 Log.Write(" {:<34}{}", "Polarized AO type:", "IAO/2014");
539 Log.Write(" {:<34}{}", "Orthogonalization method:", pOrthMethodName);
540
541 FMatrixView
542 // intrinsic basis spanning the closed & alpha orbitals. In RHF, this includes the span of the beta orbitals.
543 CIbAC = MakeStackMatrix(nBf, nIb, Mem),
544 // intrinsic basis spanning closed & beta orbitals
545 CIbBC;
546
547 {
548 FStackMatrix
549 S(nBf, nBf, &Mem),
550 Scd(nBf, nBf, &Mem);
551 MakeIntMatrix(S, *pBasis, *pBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
552 // S.Print(std::cout, "MAIN BASIS OVERLAP MATRIX #24");
553 if (1) {
554 // apply a sanity check to the orbitals.
555 FStackMatrix
556 SMo(COccAC.nCols, COccAC.nCols, &Mem);
557 ChainMxm(SMo, Transpose(COccAC), S, COccAC, Mem);
558 double
559 fRmsd = SMo.fRmsdFromIdentity();
560 if (fRmsd > 1e-5) {
561 fmt::MemoryWriter ss;
562 ss.write("Sanity check for input orbitals failed! Rmsd from orthogonality: {:8.2e}", fRmsd);
563 Log.EmitWarning(ss.str());
564 Log.Write(" If the input orbitals came from a .molden file: Beware of differnt normalizations for"
565 "\n different programs. Options: (1) Check for a newer version of IboView, (2) If it"
566 "\n still happens with the current version, contact Gerald Knizia.");
567 }
568 // SMo.Print(std::cout, "MO BASIS OVERLAP MATRIX #22");
569 }
570 Move(Scd, S);
571 CalcCholeskyFactors(Scd);
572 MakeIaoBasisNew(CIbAC, &*pMinBasis, *pAtoms, &*pBasis, COccAC, S, Scd, IaoFlags, Mem);
573 if (WfType == WFTYPE_Uhf) {
574 CIbBC = MakeStackMatrix(nBf, nIb, Mem);
575 MakeIaoBasisNew(CIbBC, &*pMinBasis, *pAtoms, &*pBasis, COccBC, S, Scd, IaoFlags, Mem);
576 } else {
577 CIbBC = CIbAC;
578 }
579 }
580
581
582 if (0) {
583 FStackMatrix
584 Rdm(nBf, nBf, &Mem),
585 S(nBf, nBf, &Mem),
586 Op(nBf, nBf, &Mem);
587 Mxm(Rdm, COccAC, Transpose(COccAC));
588 Op.Clear();
589 pAtoms->AddKineticMatrix(Op, *pBasis, *pBasis, Mem);
590 MakeIntMatrix(S, *pBasis, *pBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
591 // MakeIntMatrix(Op, *pBasis, *pBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
592 Log.WriteResult("Nuclear repulsion energy", pAtoms->NuclearRepulsionEnergy());
593 Log.WriteResult("Kinetic energy", 2*Dot(Rdm,Op));
594 Log.WriteResult("Number of electrons", Dot(Rdm,S));
595 }
596
597
598 // if (1) {
599 // // keep a copy of the IAO coefficients and the basis sets. May be needed to form hybrids.
600 // // Not sure if this is a good idea; could be better to just make a new one whenever needed.
601 // m_pOrbBasis = pBasis;
602 // m_pMinBasis = pMinBasis;
603 // m_CIaoData.resize(CIb.GetStridedSize());
604 // memcpy(&m_CIaoData[0], &CIb[0], sizeof(CIb[0]) * CIb.GetStridedSize());
605 // m_CIb = CIb;
606 // m_CIb.pData = &m_CIaoData[0];
607 // }
608
609 // setup the localization options.
610 FLocalizeOptions
611 LocOpt;
612 QString
613 LocMethod = pWfOptions->GetLocMethod().toLower();
614 if (AllowLocalize) {
615 LocOpt.nLocExp = 2;
616 if (LocMethod == "ibo (exponent 2)")
617 LocOpt.nLocExp = 2;
618 else if (LocMethod == "ibo (exponent 4)")
619 LocOpt.nLocExp = 4;
620 else if (LocMethod == "none")
621 LocOpt.nLocExp = 0;
622 else {
623 LocOpt.nLocExp = 2;
624 Log.EmitWarning(q2s(IvFmt("Localization type '%1' not recognized. Defaulting to IBO (exponent 2).", LocMethod)));
625 }
626 Log.Write(" {:<34}L = sum[A,i] (i|n_A|i)^{} -> max", "Localization functional:", LocOpt.nLocExp);
627 Log.WriteLine();
628 } else {
629 LocOpt.nLocExp = 0;
630 LocMethod = "none";
631 }
632
633 LocOpt.ThrLoc = std::min(double(1e-6), double(1e-2 * pWfOptions->GetThrGrad()));
634 // ^- WfOptions::ThrGrad is actually the SCF gradient threshold.
635
636 typedef std::set<FOrbital*>
637 FOrbitalRefSet;
638 FOrbitalRefSet
639 ProcessedOrbitals;
640
641 bool
642 InputOrbitalsChanged = false;
643 for (int iCase = 0; iCase < 2; ++ iCase)
644 {
645 // need to define:
646 // - RefOrbitals set and COcc
647 // - Coeffs of IB vectors which span this space
648 // - factors for spin and total densities?
649
650 QString
651 SpaceDesc = "";
652 TArray<FOrbital*>
653 RefOrbitals; // orbitals in current subspace which are localized together (i.e., freely mixable w/o changing the wave function)
654 FMatrixView
655 COcc; // corresponding orbital matrix
656 double
657 fChargeFactor,
658 fSpinFactor;
659 FMatrixView
660 CIb; // IB which spanns the orbitals of the current subspace.
661 if (WfType == WFTYPE_Rhf) {
662 CIb = CIbAC; // <- that's all the occupied orbtials in RHF.
663 if (iCase == 0) {
664 // closed-shell orbitals.
665 SpaceDesc = "AB";
666 fChargeFactor = 2.;
667 fSpinFactor = 0.;
668 MakeOrbitalMatrix(COcc, pBasis, &RefOrbitals, ORBMAT_OccupiedOnly | ORBMAT_ClosedOnly, Mem);
669 } else if (iCase == 1) {
670 // open-shell orbitals (alpha occupied).
671 SpaceDesc = "A_";
672 fChargeFactor = 1.;
673 fSpinFactor = 1.;
674 MakeOrbitalMatrix(COcc, pBasis, &RefOrbitals, ORBMAT_OccupiedOnly | ORBMAT_AlphaOnly, Mem);
675 } else if (iCase == 2) {
676 // external orbitals (unoccupied)
677 SpaceDesc = "__";
678 fChargeFactor = 1.;
679 fSpinFactor = 1.;
680 MakeOrbitalMatrix(COcc, pBasis, &RefOrbitals, ORBMAT_VirtualOnly, Mem);
681 // ^- TODO: if I actually had all the virtual orbitals, would the LeastSquaresSolve be sufficient
682 // to get the span of the anti-bonds? (w/o explicitly constructing them).
683 // Would have to update the nOcc number, but if it is ONLY that...
684 } else {
685 continue;
686 }
687 } else if (WfType == WFTYPE_Uhf) {
688 if (iCase == 0) {
689 // alpha orbitals.
690 CIb = CIbAC;
691 SpaceDesc = "A_";
692 fChargeFactor = 1.;
693 fSpinFactor = 1.;
694 MakeOrbitalMatrix(COcc, pBasis, &RefOrbitals, ORBMAT_OccupiedOnly | ORBMAT_AlphaOnly, Mem);
695 } else {
696 // beta orbitals.
697 CIb = CIbBC;
698 SpaceDesc = "_B";
699 fChargeFactor = 1.;
700 fSpinFactor = -1.;
701 MakeOrbitalMatrix(COcc, pBasis, &RefOrbitals, ORBMAT_OccupiedOnly | ORBMAT_BetaOnly, Mem);
702 }
703 } else {
704 Log.EmitError("Unrecognized wave function type in RunIaoAnalysis (inner).");
705 continue;
706 }
707 IR_SUPPRESS_UNUSED_WARNING2(fSpinFactor, fChargeFactor);
708 uint
709 nOcc = COcc.nCols;
710 if (nOcc == 0)
711 continue; // no orbitals in this class.
712
713 // re-express occupied vectors in terms of the IAO vectors. This is an exact
714 // transformation.
715 // (note: this probably can be done in a different way than via SVD... probably via
716 // Sib^{-1} CIb.T S1 COcc, where Sib = CIb.T S1 CIb
717 // I just don't want to think about it for the moment)
718 // Note: BEWARE OF THE ABSORBED OCCUPATION NUMBERS!!
719
720 // solve CIb CIbOcc = COcc for CIbOcc
721 FStackMatrix
722 CIbOcc(nIb, nOcc, &Mem);
723 LeastSquaresSolveSafe(CIb, CIbOcc, COcc, 1e-9, Mem);
724
725 if (0) {
726 // test if IB works.
727 uint
728 nBf = pBasis->nFn();
729 FStackMatrix
730 // COcc reconstructed from IB and CIbOcc.
731 COccRec(nBf, nIb, &Mem);
732 Mxm(COccRec, CIb, CIbOcc);
733 Add(COccRec, COcc, -1.);
734 std::cout << boost::format("Norm(COccRec - COcc) = %8.2e [nBf=%i, nOcc=%i]\n") % CalcMatrixRmsNorm(COccRec) % nBf % COcc.nCols;
735 }
736
737
738 // make a fock matrix by assuming that the input orbitals are canonical.
739 // (FIXME: need all orbitals for this. at least when making anti-bonds.)
740 FStackMatrix
741 Ew(nOcc, 1, &Mem);
742 for (size_t iOcc = 0; iOcc < nOcc; ++ iOcc)
743 Ew[iOcc] = RefOrbitals[iOcc]->info.fEnergy;
744
745
746 // {
747 // Log.Write("!INITIAL ORBITALS:");
748 // FPrintOrbitalOptions PrintOpt;
749 // PrintOrbitalChargeComposition(Log, CIbOcc, Ew, pAtoms, &*pMinBasis, "AB", 2., PrintOpt, Mem);
750 // }
751
752
753 // Log.Write("Run localize? {} AllowLocalize? {}", LocOpt.nLocExp, (int)AllowLocalize);
754 if (LocOpt.nLocExp != 0) {
755 std::stringstream
756 ss;
757 LocOpt.Verbosity = 1;
758 LocOpt.pxout = &ss;
759 size_t
760 nAt = pAtoms->size();
761 TMemoryLock<size_t>
762 pAtOffsets(nAt+1, &Mem);
763 // ^- that is NOT the same as pRawBasis->CenterOffsets! The latter is in terms of shells!!!
764 for (size_t iAt = 0; iAt <= nAt; ++ iAt)
765 pAtOffsets[iAt] = pMinBasis->pRawBasis->iCenFn(iAt);
766 if (nAt + 1 != pMinBasis->pRawBasis->CenterOffsets.size())
767 Log.EmitError("Number of atoms in minimal basis inconsistent with frame main basis (?).");
768 if (CIbOcc.nRows != pMinBasis->pRawBasis->nFn())
769 Log.EmitError("Number of basis functions in minimal basis inconsistent with frame main basis (?).");
770
771 if (0) {
772 std::stringstream str;
773 pMinBasis->Print(str);
774 Log.Write(str.str());
775 Log.Write("\ncenter offsets:");
776 for (size_t i = 0; i < nAt+1; ++ i)
777 Log.Write(" {:4}. {:6}", i, pAtOffsets[i]);
778 }
779
780 FStackMatrix
781 PrevCIbOcc(nIb, nOcc, &Mem),
782 PrevEw(nOcc, 1, &Mem);
783 Move(PrevCIbOcc, CIbOcc);
784 Move(PrevEw, Ew);
785
786 if (iCase == 0)
787 Log.WriteLine();
788 LocalizeVectors(CIbOcc, pAtOffsets, nAt, LocOpt, Mem);
789 Log.Write("{}; Class: {}", ss.str(), q2s(SpaceDesc));
790
791 // make new set of vectors by transforming CIbOcc back to the
792 // main basis.
793 Mxm(COcc, CIb, CIbOcc);
794
795 if (1) {
796 // compute new orbital energies---assuming old orbitals were canonical.
797 FStackMatrix
798 U(nOcc, nOcc, &Mem);
799 Mxm(U, Transpose(PrevCIbOcc), CIbOcc);
800 for (size_t iOcc = 0; iOcc < nOcc; ++ iOcc) {
801 Ew[iOcc] = 0;
802 for (size_t iPrev = 0; iPrev < nOcc; ++ iPrev)
803 Ew[iOcc] += PrevEw[iPrev] * sqr(U(iPrev, iOcc));
804 }
805 }
806 // fixme: fock matrix? occupation numbers?
807 InputOrbitalsChanged = true;
808 }
809
810 // {
811 // // Log.Write("!FINAL ORBITALS:");
812 // FPrintOrbitalOptions PrintOpt;
813 // PrintOpt.PrintHeader = iCase == 0;
814 // PrintOrbitalChargeComposition(Log, CIbOcc, Ew, pAtoms, &*pMinBasis, q2s(SpaceDesc), fChargeFactor, PrintOpt, Mem);
815 // }
816
817 // store IAO coefficients of orbitals in the orbital data sets,
818 // and also a link to the minimal basis (to distinguish what is what in the basis)
819 for (uint iOcc = 0; iOcc < COcc.nCols; ++ iOcc) {
820 FOrbital *pOrb = RefOrbitals[iOcc];
821 if (ProcessedOrbitals.find(pOrb) != ProcessedOrbitals.end())
822 IvNotify(NOTIFY_Warning, "An orbital occured in more than one localization group!");
823
824 ProcessedOrbitals.insert(pOrb);
825 pOrb->pMinBasis = pMinBasis;
826 pOrb->pIaoCoeffs.clear();
827 pOrb->pIaoCoeffs.insert(pOrb->pIaoCoeffs.begin(), &CIbOcc(0,iOcc), &CIbOcc(0,iOcc)+nIb);
828 if (InputOrbitalsChanged) {
829 RefOrbitals[iOcc]->info.fEnergy = Ew[iOcc];
830 pOrb->pCoeffs.clear();
831 pOrb->pCoeffs.insert(pOrb->pCoeffs.begin(), &COcc(0,iOcc), &COcc(0,iOcc)+nBf);
832 }
833 }
834 }
835 if (InputOrbitalsChanged) {
836 // remove all orbitals which were not processed. This is used to delete either
837 // all or non-valence virtual orbitals after the localization procedure.
838 if (1) {
839 FDataSetList
840 NewData;
841 for (size_t i = 0; i < m_Data.size(); ++ i) {
842 FOrbital
843 *pOrb = dynamic_cast<FOrbital*>(m_Data[i].get());
844 if (pOrb == 0 || ProcessedOrbitals.find(pOrb) != ProcessedOrbitals.end())
845 // copy pointers to all data sets which either are not orbitals,
846 // or are orbitals which were processed.
847 NewData.push_back(m_Data[i]);
848 }
849 m_Data.swap(NewData);
850 }
851
852 // sort new orbitals (by type and energy) and update their descriptions
853 std::stable_sort(m_Data.begin(), m_Data.end(), FSortOrbitalPred());
854
855 size_t iOrb = 0;
856 for (size_t i = 0; i < m_Data.size(); ++ i) {
857 FOrbital *pOrb = dynamic_cast<FOrbital*>(m_Data[i].get());
858 if (pOrb) {
859 pOrb->info.iOriginalIndex = i;
860 pOrb->UpdateDescFromInfo(iOrb);
861 // Log.Write(" pOrb[{:4}] E = {:12.5f} NewDesc = '{}'", i, pOrb->info.fEnergy, q2s(pOrb->GetDesc()));
862 iOrb += 1;
863 }
864 }
865
866 // print orbital composition in terms of atomic partial charges (with new order)
867 double ThrPrint = 0.02;
868 iOrb = 0;
869 Log.Write("\n Summary of localized orbital composition: [THRPRINT={:6.3f}]", ThrPrint);
870 Log.Write("\n ORB GRP ORB.ENERGY CENTERS/CHARGES");
871 for (size_t i = 0; i < m_Data.size(); ++ i) {
872 FOrbital *pOrb = dynamic_cast<FOrbital*>(m_Data[i].get());
873 if (pOrb) {
874 Log.Write("{:6} {:>3} {:12.6f} {}", iOrb+1, pOrbDescFromType(pOrb->info.Spin, pOrb->info.fOcc),
875 pOrb->info.fEnergy, q2s(pOrb->MakeFullDesc(ThrPrint, FOrbital::ORBDESC_ChargesOnly)));
876 iOrb += 1;
877 }
878 }
879 }
880
881 if (1) {
882 Log.WriteLine();
883 RunChargeAnalysis(Log, 0);
884 }
885
886 Log.WriteLine();
887 Log.WriteTiming("IBBA", (double)tIbbaTotal);
888 Log.WriteLine();
889
890 // GetBondOrderType
891 // GetThrGrad
892 // GetOrbDisplay
893 // GetOrbDivision
894 // GetLocMethod
895 }
896
897
898
899
MakeIaoBasis(ct::FAtomSet * pAtoms,ct::FMemoryStack & Mem)900 ct::FMatrixView FFrame::MakeIaoBasis(ct::FAtomSet *pAtoms, ct::FMemoryStack &Mem)
901 {
902 IR_SUPPRESS_UNUSED_WARNING2(pAtoms, Mem);
903 return m_CIb;
904 }
905
906
MakeHybridsForSelectedAtomGroup()907 void FDocument::MakeHybridsForSelectedAtomGroup()
908 {
909 }
910