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