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 <iostream>
25 #include <cmath>
26 #include <boost/format.hpp>
27 using boost::format;
28 #include "Ir.h"
29 #include "CxDiis.h"
30 #include "CtRhf.h"
31 #include "CtIo.h"
32 #include "CtTiming.h"
33 #include "CxPodArray.h"
34 #if 0
35 #include "RDTSC.h"
36 #else
37 #define RESET_CLOCKS
38 #define RESUME_CLOCK(x)
39 #define PAUSE_CLOCK(x)
40 #endif
41
42 #include "CtConstants.h"
43 #include "CtDft.h"
44
45
46 #include "CtDfti.h"
47
48 namespace ct {
49
50 double g_fThrOrb = 1e-12;
51 // double g_fThrOrb = 1e-10;
52
MakeLogThrOrb(double ThrOrb)53 double MakeLogThrOrb(double ThrOrb) {
54 // hm... this just used a constant 40 before. Seems a bit much.
55 // set it to ThrOrb * 1e-3 or something?
56 return -std::log(.1 * ThrOrb);
57 // return 40.;
58 }
59
~FFockComponentBuilder()60 FFockComponentBuilder::~FFockComponentBuilder()
61 {}
62
AccFock(FMatrixView & FockC,FMatrixView & FockO,FBasisSet * pBasis,FMatrixView const & COccC,FMatrixView const & COccO,uint Flags,FMemoryStack & Mem)63 void FFockComponentBuilder::AccFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FMemoryStack &Mem)
64 {
65 size_t
66 nBf = pBasis->nFn();
67 if (!((FockC.nRows == nBf && FockC.nCols == nBf) &&
68 (FockO.pData == 0 || (FockO.nRows == nBf && FockO.nCols == nBf)) &&
69 (COccC.nRows == nBf) &&
70 (COccO.pData == 0 || COccO.nRows == nBf)))
71 throw std::runtime_error("FFockComponentBuilder: Input orbitals not consistent with orbital basis sets.");
72 IR_SUPPRESS_UNUSED_WARNING(Flags);
73 IR_SUPPRESS_UNUSED_WARNING(Mem);
74 }
75
PrintEnergyContribs()76 void FFockComponentBuilder::PrintEnergyContribs()
77 {
78 m_Log.Write("FFockComponentBuilder::PrintEnergyContribs: not implemented for current Fock builder. Fix this.");
79 }
80
AccGradient(FMatrixView Gradient,FMatrixView COccC,FMatrixView COccO,FMemoryStack & Mem)81 void FFockComponentBuilder::AccGradient(FMatrixView Gradient, FMatrixView COccC, FMatrixView COccO, FMemoryStack &Mem)
82 {
83 m_Log.Write("FFockComponentBuilder::AccGradient: not implemented for current Fock builder. Fix this.");
84 }
85
SwitchToRefineGrid(FDftGridParams const & NewGridParams)86 bool FFockComponentBuilder::SwitchToRefineGrid(FDftGridParams const &NewGridParams)
87 {
88 return false;
89 IR_SUPPRESS_UNUSED_WARNING(NewGridParams);
90 }
91
92
93
94 // forms the integral matrix set (\mu\nu|F) for a single set of F shells.
95 // \mu comes from OrbBasis1, \nu comes from OrbBasis2. If both pointers are identical,
96 // symmetric integrals will be assumed [(\mu\nu|F) = (\nu\mu|F)] and this symmetry will be used.
FormIntMNF(ir::FRawShell const & ShF,ir::FIntegralKernel * pIntKernel,FRawBasis const * pOrbBasisA,FRawBasis const * pOrbBasisB,FRawBasis const * pFitBasis,FMemoryStack & Mem,FMatrixView ScrDen,double fThr)97 double *FormIntMNF(ir::FRawShell const &ShF,
98 ir::FIntegralKernel *pIntKernel, FRawBasis const *pOrbBasisA, FRawBasis const *pOrbBasisB, FRawBasis const *pFitBasis, FMemoryStack &Mem, FMatrixView ScrDen, double fThr)
99 {
100 size_t
101 nAoA = pOrbBasisA->nFn(),
102 nAoB = pOrbBasisB->nFn(),
103 nFnF = ShF.nFn();
104 double
105 *pMNF;
106 bool
107 Symmetric = (pOrbBasisA == pOrbBasisB);
108 Mem.ClearAlloc(pMNF, nAoA * nAoB * nFnF);
109
110 std::size_t nIntTotal = 0, nIntRetained = 0;
111
112 for ( size_t iShB = 0; iShB < pOrbBasisB->Shells.size(); ++ iShB ){
113 ir::FRawShell const &ShB = pOrbBasisB->Shells[iShB];
114 size_t nFnB = ShB.nFn();
115 size_t iShA_First = Symmetric ? iShB : 0;
116 for ( size_t iShA = iShA_First; iShA < pOrbBasisA->Shells.size(); ++ iShA ) {
117 ir::FRawShell const &ShA = pOrbBasisA->Shells[iShA];
118 size_t nFnA = ShA.nFn();
119 nIntTotal += nFnA * nFnB;
120
121 // if (ScrDen(iShA, iShB) < fThr)
122 // continue;
123 double
124 fDistSqAB = DistSq(FVector3(ShA.vCen), FVector3(ShB.vCen));
125 if (ShA.pRange && ShB.pRange && sqr(ShA.MaxCoRange() + ShB.MaxCoRange()) < fDistSqAB)
126 continue;
127 size_t
128 Strides[3] = {1, nFnA, nFnA * nFnB};
129
130 // double fElecInPair = 0.;
131 // for (size_t iA = 0; iA < nFnA; ++ iA)
132 // for (size_t iB = 0; iB < nFnB; ++ iB)
133 // // fElecInPair += ScrDen((pOrbBasisA->iFn(iShA) + iA), (pOrbBasisB->iFn(iShB) + iB));
134 // fElecInPair += sqr(ScrDen((pOrbBasisA->iFn(iShA) + iA), (pOrbBasisB->iFn(iShB) + iB)));
135 // if (std::abs(fElecInPair) < sqr(fThr))
136 // continue;
137 nIntRetained += nFnA * nFnB;
138
139 double
140 *pIntData;
141 Mem.Alloc(pIntData, nFnA * nFnB * nFnF );
142
143 RESUME_CLOCK(100)
144 ir::EvalInt2e3c(pIntData, Strides, &ShA, &ShB, &ShF,1, 1.0, pIntKernel, Mem);
145 PAUSE_CLOCK(100)
146
147 if (Symmetric) {
148 assert(pOrbBasisA == pOrbBasisB && nAoA == nAoB);
149 for ( size_t iF = 0; iF < nFnF; ++ iF )
150 for ( size_t iB = 0; iB < nFnB; ++ iB )
151 for ( size_t iA = 0; iA < nFnA; ++ iA ) {
152 double
153 f = pIntData[iA + nFnA * (iB + nFnB * iF)];
154 // assign to (\mu\nu| and (\nu\mu|. (int has perm symmetry).
155 pMNF[(pOrbBasisA->iFn(iShA) + iA) + nAoA*(pOrbBasisA->iFn(iShB) + iB) + nAoA*nAoA*iF] = f;
156 pMNF[(pOrbBasisA->iFn(iShB) + iB) + nAoA*(pOrbBasisA->iFn(iShA) + iA) + nAoA*nAoA*iF] = f;
157 }
158 } else {
159 for ( size_t iF = 0; iF < nFnF; ++ iF )
160 for ( size_t iB = 0; iB < nFnB; ++ iB )
161 for ( size_t iA = 0; iA < nFnA; ++ iA ) {
162 double
163 f = pIntData[iA + nFnA * (iB + nFnB * iF)];
164 pMNF[(pOrbBasisA->iFn(iShA) + iA) + nAoA*(pOrbBasisB->iFn(iShB) + iB) + nAoA*nAoB*iF] = f;
165 }
166 }
167
168 Mem.Free(pIntData);
169 }
170 }
171
172 // if (&ShF == &pFitBasis->Shells[0] && fThr != 0.) {
173 // if (&ShF == &pFitBasis->Shells[0]) {
174 // double f = 0, g = 0;
175 // for (std::size_t i = 0; i < ScrDen.GetStridedSize(); ++ i)
176 // f += ScrDen.pData[i];
177 // for (std::size_t i = 0; i < nAo; ++ i)
178 // g += ScrDen.pData[i*(nAo+1)];
179 // xout << format("ScrAB: %i of %i (\\mu\\nu| integrals evaluated (%.2f%% screened, %.6f elec total / %.6f diag).\n")
180 // % nIntRetained % nIntTotal % (100.*(1. - nIntRetained/static_cast<double>(nIntTotal))) % f % g;
181 // }
182 // if (&ShF == &pFitBasis->Shells[0]) {
183 // xout << format("ScrAB: %i of %i (\\mu\\nu| integrals evaluated (%.2f%% screened).\n")
184 // % nIntRetained % nIntTotal % (100.*(1. - nIntRetained/static_cast<double>(nIntTotal)));
185 // }
186
187 return pMNF;
188 IR_SUPPRESS_UNUSED_WARNING(pFitBasis);
189 IR_SUPPRESS_UNUSED_WARNING(ScrDen);
190 IR_SUPPRESS_UNUSED_WARNING(fThr);
191 }
192
193
FormIntMNF(ir::FRawShell const & ShF,ir::FIntegralKernel * pIntKernel,FRawBasis const * pOrbBasis,FRawBasis const * pFitBasis,FMemoryStack & Mem,FMatrixView ScrDen,double fThr)194 double *FormIntMNF(ir::FRawShell const &ShF,
195 ir::FIntegralKernel *pIntKernel, FRawBasis const *pOrbBasis, FRawBasis const *pFitBasis, FMemoryStack &Mem, FMatrixView ScrDen, double fThr)
196 {
197 return FormIntMNF(ShF, pIntKernel, pOrbBasis, pOrbBasis, pFitBasis, Mem, ScrDen, fThr);
198 }
199
200
201
202
203
204
205
FFockComponentBuilderDfCoulXcCached(FDfJkOptions const & JkOptions_,FDftGridParams const & GridParams_,std::string const & XcFunctionalName_,FLog & Log_,FTimerSet * pTimers_)206 FFockComponentBuilderDfCoulXcCached::FFockComponentBuilderDfCoulXcCached(FDfJkOptions const &JkOptions_, FDftGridParams const &GridParams_, std::string const &XcFunctionalName_, FLog &Log_, FTimerSet *pTimers_)
207 : FFockComponentBuilder(Log_, pTimers_), XcFunctionalName(XcFunctionalName_), JkOptions(JkOptions_), GridParams(GridParams_)
208 {
209 if ( !XcFunctionalName.empty() ) {
210 g_pXcFunctional = new FXcFunctional(XcFunctionalName);
211 pXcFn = g_pXcFunctional;
212 m_Log.Write("\n"+pXcFn->Desc());
213 }
214 EnergyCoulomb = 0.;
215 EnergyXc = 0.;
216 fElecTotal = 0.;
217 AuxExpandXc = true;
218 // AuxExpandXc = false;
219 // m_Log.Write(" !!xc/j builder initialized with grid level: {}/{}", GridParams.nLevel, GridParams_.nLevel);
220 }
221
~FFockComponentBuilderDfCoulXcCached()222 FFockComponentBuilderDfCoulXcCached::~FFockComponentBuilderDfCoulXcCached()
223 {
224 }
225
Init(FWfDecl const & WfDecl_,FBasisSet * pOrbBasis_,FAtomSet const & Atoms_,FHfOptions const & Options_,FMemoryStack & Mem)226 void FFockComponentBuilderDfCoulXcCached::Init(FWfDecl const &WfDecl_, FBasisSet *pOrbBasis_, FAtomSet const &Atoms_, FHfOptions const &Options_, FMemoryStack &Mem)
227 {
228 m_WfDecl = WfDecl_;
229 fElecTotalAnalytic = double(WfDecl_.nElec);
230
231 pAtoms = &Atoms_;
232 pOrbBasis = pOrbBasis_;
233 nAo = pOrbBasis->nFn();
234
235 if (1) {
236 pFitBasis = new FBasisSet(*pAtoms, BASIS_JFit);
237 } else {
238 pFitBasis = new FBasisSet(*pAtoms, BASIS_JkFit);
239 m_Log.Write(" NOTE: Using JkFit basis as JFit. Fix this!");
240 }
241
242 SwitchToRefineGrid(GridParams);
243 // if (pXcFn) {
244 // // find minimum L which should be integrated exactly. Depends on the basis to some degree.
245 // // (this is further adjusted inside the grid generator based on atom types)
246 // GridParams.iMinL = AuxExpandXc? (2*pFitBasis->nMaxL()) : (2*pOrbBasis->nMaxL());
247 // // FTimer
248 // // tDftGrid;
249 // pDftGrid = new FDftGrid(Atoms_, GridParams, &m_Log);
250 // // m_Log.Write(" Generated DFT grid with {} points for {} atoms in {:.2} sec.\n", pDftGrid->Points.size(), Atoms_.size(), (double)tDftGrid);
251 // // xout << "" << format(pTimingFmt) % "DFT integration grid" % (double)tDftGrid; xout.flush();
252 // }
253
254 // Make fitting coefficients J^{-1/2}.
255 FTimer TimerJcd;
256 // xout << *pFitBasis;
257 pFitBasisRaw = pFitBasis->pRawBasis.get();
258 pOrbBasisRaw = pOrbBasis->pRawBasis.get();
259
260 nFit = pFitBasis->nFn();
261 Jcd = MakeStackMatrix(nFit, nFit, Mem);
262 MakeIntMatrix(Jcd, *pFitBasis, *pFitBasis, FKrn2i_Direct(&ir::g_IrCoulombKernel), Mem);
263 CalcCholeskyFactors(Jcd);
264 m_Log.WriteTiming("fitting metric (coul)", (double)TimerJcd);
265
266 // make the cached integrals themselves.
267 size_t
268 nAoTr = nAo * (nAo+1)/2;
269 Int3ixStorage.resize(nAoTr * nFit);
270 Int3ix = FMatrixView(&Int3ixStorage[0], nAoTr, nFit);
271
272 FTimer Timer3ix;
273 FMatrixView ScrDen;
274 {
275 FMemoryStackArray MemStacks(Mem);
276 #pragma omp parallel for schedule(dynamic)
277 for ( int iShF__ = 0; iShF__ < int(pFitBasisRaw->Shells.size()); ++ iShF__ ) {
278 size_t iShF = size_t(iShF__); // that one's for OpenMP.
279 //for ( size_t iShF = 0; iShF < pFitBasisRaw->Shells.size(); ++ iShF ) {
280 if (m_Log.StatusOkay()) {
281 FMemoryStack &Mem1 = MemStacks.GetStackOfThread();
282 ir::FRawShell const &ShF = pFitBasisRaw->Shells[iShF];
283 size_t nFnF = ShF.nFn();
284 // FMemoryStack2
285 // // FIXME: this is a hack around the workspace memory size computation...
286 // // the *much* better alternative would be to do this op with less space.
287 // Mem1(nAo * nAo * nFnF * sizeof(double) + 100000);
288
289 double
290 // (\mu\nu|F): nAo x nAo x nFnF
291 *pMNF = FormIntMNF(ShF, &ir::g_IrCoulombKernel, pOrbBasisRaw, pFitBasisRaw, Mem1, ScrDen, 0.);
292
293 for (size_t iF = 0; iF < nFnF; ++ iF) {
294 FMatrixView
295 MN(pMNF + nAo*nAo*iF, nAo, nAo);
296 MN.TriangularReduce(1, &Int3ix(0, iF + pFitBasisRaw->iFn(iShF)));
297
298 assert(MN.TriangularStorageSize1() == nAoTr);
299 }
300
301 Mem1.Free(pMNF);
302 }
303 }
304 }
305 m_Log.CheckStatus(); // may raise exception.
306
307 double
308 fIntSizeMb = double(Int3ix.GetStridedSize()) * double(sizeof(Int3ix[0])) / double(1<<20);
309 // xout << format(pTimingFmt) % str(format("3-index integrals (%.2f MB)") % fIntSizeMb) % (double)Timer3ix; xout.flush();
310 // xout << format(pTimingFmt) % str(format("3-index integrals (%i MB)") % (size_t)fIntSizeMb) % (double)Timer3ix; xout.flush();
311 m_Log.WriteTiming(fmt::format("3-index integrals ({} MB)", (size_t)fIntSizeMb), (double)Timer3ix);
312 IR_SUPPRESS_UNUSED_WARNING(Options_);
313 IR_SUPPRESS_UNUSED_WARNING(WfDecl_);
314 }
315
316
317
318
319
320
321
322 // fixme: this should probably make separate open- and closed-shell densities. open-shell needed for xc (but not for coulomb).
Make1ixDensity(FMatrixView jgamma,FMatrixView AuxDen,FMatrixView COccC,FMatrixView COccO,FMemoryStack & Mem)323 void FFockComponentBuilderDfCoulXcCached::Make1ixDensity(FMatrixView jgamma, FMatrixView AuxDen, FMatrixView COccC, FMatrixView COccO, FMemoryStack &Mem)
324 {
325 FStackMatrix
326 Density(nAo, nAo, &Mem); // OrbC x OrbC.T + OrbO x OrbO.T
327 // jgamma(nFit, 1, &Mem), // (A|rs) gamma[r,s].
328 // AuxDen(nFit, 1, &Mem); // (A|B)^{-1} jgamma[B]
329 m_pTimers->Resume(0x201, "DF-J/XC (RDM)");
330 Mxm(Density, COccC, Transpose(COccC));
331 if (COccO.nCols != 0)
332 Mxm(Density, COccO, Transpose(COccO), MXM_Add);
333
334 Density.TriangularReduce(1, 0, 2.); // reduce to triangular form; multiply off diagonal elements by 2.0
335 m_pTimers->Pause(0x201, "DF-J/XC (RDM)");
336
337 m_pTimers->Resume(0x202, "DF-J/XC (1x density fit)");
338 Mxva(jgamma.pData, Transpose(Int3ix), Density.pData);
339
340 Move(AuxDen, jgamma);
341
342 // Solve[Jcd[AB] D[\nu B i]] -> D[\nu A I] to get density coefficients
343 CholeskySolve(AuxDen, Jcd);
344 m_pTimers->Pause(0x202, "DF-J/XC (Density fit)");
345 }
346
AccFock(FMatrixView & FockC,FMatrixView & FockO,FBasisSet * pBasis,FMatrixView const & OrbC,FMatrixView const & OrbO,uint Flags,FMemoryStack & Mem)347 void FFockComponentBuilderDfCoulXcCached::AccFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &OrbC, FMatrixView const &OrbO, uint Flags, FMemoryStack &Mem)
348 {
349 FFockComponentBuilder::AccFock(FockC, FockO, pBasis, OrbC, OrbO, Flags, Mem); // consistency checks.
350
351 // if (1) {
352 // xout << "CALLED GEMMV-INIT! -- CALL" << std::endl;
353 // double a[2] = {0}, b[4] = {0}, c[2] = {0};
354 // // a little hack around a thread-init problem in eigen... to prevent that the first call of DGEMV
355 // // happens inside a OpenMP block. This collides with the static initializations done in eigen's version
356 // // of those routines, which are not quite threadsafe (crashes in -O3 mode due to init flag re-arrangements by the compiler).
357 // Mxv(&a[0], 1, &b[0],1,2, &c[0],1, 2,2);
358 // xout << "CALLED GEMMV-INIT! -- DONE" << std::endl;
359 // }
360
361 if (pBasis != pOrbBasis)
362 throw std::runtime_error("FFockComponentBuilderDfCoulXcCached must now be used with orb-projected initial guess (not Fock-projected guess!)");
363 // if (pBasis != pOrbBasis) {
364 // m_pTimers->Enter(0x201, "DF-J/XC (Guess)");
365 // // hm hm. This might be a problem since we project Fock matrices...
366 // // (note: not sure if this actually *accumulates* matrices, instead of overwriting them)
367 // FockC.Clear();
368 // MakeCoul(FockC, &ir::g_IrCoulombKernel, &*pBasis->pRawBasis, &*pFitBasis->pRawBasis, OrbC, OrbO, Jcd, Mem);
369 //
370 // if (pXcFn) {
371 // using namespace dfti;
372 //
373 // size_t nGu = pBasis->nFn();
374 //
375 // FStackMatrix
376 // DenC(nGu, nGu, &Mem);
377 // Mxm(DenC, OrbC, Transpose(OrbC));
378 // DenC.TriangularReduce(1, 0, 1.);
379 //
380 // FStackMatrix
381 // FockTr(nGu, nGu, &Mem);
382 // FockTr.Clear();
383 //
384 // EnergyXc = 0.;
385 // FDftiArgs DftiArgs = {
386 // DFTI_MakeXc,
387 // &EnergyXc, 0, FockTr.pData, 0,
388 // DenC.pData, 0, 0, 0, 0, 0, // no orbitals and no open-shell density provided.
389 // pBasis->pRawBasis.get(),
390 // pXcFn.get(),
391 // pDftGrid.get(),
392 // g_fThrOrb, MakeLogThrOrb(g_fThrOrb), // 1e-1 * ThrDen?
393 // };
394 // AccXc(DftiArgs, m_Log, Mem);
395 // FockTr.TriangularExpand();
396 // Add(FockC, FockTr, 1.0);
397 // xout << format(pResultFmt) % "Density functional" % EnergyXc;
398 // }
399 //
400 // m_pTimers->Leave(0x201, "DF-J/XC (Guess)");
401 // return;
402 // }
403
404 if ( nAo != OrbC.nRows || nAo != OrbO.nRows || nFit != Jcd.nRows )
405 throw std::runtime_error("AccFock: Input orbitals not consistent with orbital basis set.");
406
407 m_pTimers->Enter(0x200, "DF-J/XC");
408
409 FStackMatrix
410 jgamma(nFit, 1, &Mem), // (A|rs) gamma[r,s].
411 AuxDen(nFit, 1, &Mem); // (A|B)^{-1} jgamma[B]
412 Make1ixDensity(jgamma, AuxDen, OrbC, OrbO, Mem);
413
414 if (1 && AuxExpandXc) {
415 // compute number of electrons in the auxiliary basis. Since we fit densities,
416 // this number may not be exact. Note that only s functions carry electrons.
417 // Mathematica says: Assuming[\[Alpha] > 0, Integrate[Exp[-\[Alpha]*r^2]*r^2*4*Pi, {r, 0, Infinity}]] = (pi/alpha)^(3/2)
418 fElecTotalAnalytic = 0.;
419 for (size_t iSh = 0; iSh < pFitBasisRaw->Shells.size(); ++ iSh) {
420 ir::FRawShell const &Sh = pFitBasisRaw->Shells[iSh];
421 if (Sh.l != 0) continue;
422 for (size_t iExp = 0; iExp < Sh.nExp; ++ iExp) {
423 double fDen = std::pow(M_PI/Sh.pExp[iExp], 1.5);
424 double *pCoeff = &AuxDen[pFitBasisRaw->iFn(iSh)];
425 for (size_t iCo = 0; iCo < Sh.nCo; ++ iCo)
426 fElecTotalAnalytic += pCoeff[iCo] * Sh.pCo[iExp + Sh.nExp*iCo] * fDen;
427 }
428 }
429 }
430
431
432 if (1) {
433 m_LastDensity.resize(nFit);
434 assert(sizeof(AuxDen[0]) == sizeof(m_LastDensity[0]));
435 memcpy(&m_LastDensity[0], &AuxDen[0], sizeof(AuxDen[0]) * nFit);
436 }
437
438
439
440 // Transpose(AuxDen).Print(xout, "aux density (iters)");
441 // FStackMatrix
442 // Density(nAo, nAo, &Mem), // OrbC x OrbC.T + OrbO x OrbO.T
443 // jgamma(nFit, 1, &Mem), // (A|rs) gamma[r,s].
444 // AuxDen(nFit, 1, &Mem); // (A|B)^{-1} jgamma[B]
445 // m_pTimers->Resume(0x201, "DF-J/XC (RDM)");
446 // Mxm(Density, OrbC, Transpose(OrbC));
447 // Mxm(Density, OrbO, Transpose(OrbO), MXM_Add);
448 //
449 // Density.TriangularReduce(1, 0, 2.); // reduce to triangular form; multiply off diagonal elements by 2.0
450 // m_pTimers->Pause(0x201, "DF-J/XC (RDM)");
451 //
452 // m_pTimers->Resume(0x202, "DF-J/XC (1x density fit)");
453 // Mxva(jgamma.pData, Transpose(Int3ix), Density.pData);
454 //
455 // Move(AuxDen, jgamma);
456 //
457 // // Solve[Jcd[AB] D[\nu B i]] -> D[\nu A I] to get density coefficients
458 // CholeskySolve(AuxDen, Jcd);
459 // m_pTimers->Pause(0x202, "DF-J/XC (Density fit)");
460
461 FStackMatrix
462 AuxDen1(nFit, 1, &Mem);
463 Move(AuxDen1, AuxDen);
464
465 if (1) {
466 m_pTimers->Resume(0x203, "DF-J/XC (coul. energy)");
467 // make coulomb energy from 2ix integrals:
468 // Ecoulomb = .5 j[A] (A|B) j[B]
469 // (since we use a robust fit this is equal to the density matrix contraction.
470 // the advantage here is that we can do it also in the xc auxiliary expansion case,
471 // where we will not get a pure j matrix).
472 // FStackMatrix
473 // cgamma(nFit, 1, &Mem);
474 // Move(cgamma, jgamma);
475 // TriangularMxm(Transpose(cgamma), Jcd, 'R');
476 // EnergyCoulomb = .5 * Dot2(cgamma.pData, cgamma.pData, nFit);
477 EnergyCoulomb = .5 * Dot(AuxDen.pData, jgamma.pData, nFit);
478 m_pTimers->Pause(0x203, "DF-J/XC (coul. energy)");
479 }
480
481 if (pXcFn && AuxExpandXc) {
482 m_pTimers->Resume(0x204, "DF-J/XC (xc contrib.)");
483 using namespace dfti;
484 FStackMatrix
485 vxc(nFit, 1, &Mem),
486 AuxDenRenorm(nFit, 1, &Mem);
487
488 vxc.Clear();
489 Move(AuxDenRenorm, AuxDen);
490
491
492 FDftiArgs DftiArgs = {
493 DFTI_MakeXc | DFTI_AuxiliaryExpand,
494 &EnergyXc, 0, vxc.pData, 0,
495 AuxDenRenorm.pData, 0, 0, 0, 0, 0, // no orbitals and no open-shell density provided.
496 pFitBasisRaw,
497 pXcFn.get(),
498 pDftGrid.get(),
499 g_fThrOrb, MakeLogThrOrb(g_fThrOrb), 0, &fElecTotal // 1e-1 * ThrDen?
500 };
501 // xout << "ACC-XC: ENTER" << std::endl;
502 AccXc(DftiArgs, m_Log, Mem);
503
504
505 // xout << "ACC-CX: LEAVE" << std::endl;
506 // vxc.Print(xout, "AUXILIARY XC POTENTIAL (pure, pre-J^1).");
507 if (0) {
508 FStackMatrix
509 vxcd(nFit, 1, &Mem),
510 Scd(nFit, nFit, &Mem);
511 MakeIntMatrix(Scd, *pFitBasis, *pFitBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
512 // CalcCholeskyFactors(Scd);
513 Move(vxcd, vxc);
514 // CholeskySolve(vxcd, Scd);
515 // Mxm(vxcd, Scd, vxc);
516 vxcd.Print(xout, "AUXILIARY XC POTENTIAL (S^-1).");
517 xout << format(" <vxc,auxden> = %18.12f\n") % Dot(vxcd.pData, AuxDen1.pData, nFit);
518 }
519 CholeskySolve(vxc, Jcd);
520 Add(AuxDen, vxc);
521 // AuxDen1.Print(xout, "AUXILIARY DENSITY.");
522 // vxc.Print(xout, "DENSITY-LIKE XC POTENTIAL.");
523 // xout << format(pResultFmt) % "Density functional" % DfuEnergy;
524 // xout << format(pResultFmt) % "Functional energy" % DfuEnergy;
525 // EnergyXc
526 m_pTimers->Pause(0x204, "DF-J/XC (xc contrib.)");
527 }
528
529
530
531 // recalculate integrals to form the coulomb matrix.
532 // FormIntMNF_ContractF(Coul, pIntKernel, pOrbBasis, pFitBasis, AuxDen.pData, Mem);
533 FStackMatrix
534 FockTr(nAo, nAo, &Mem);
535
536 m_pTimers->Resume(0x205, "DF-J/XC (j/vxc matrix)");
537 Mxva(FockTr.pData, Int3ix, AuxDen.pData);
538 m_pTimers->Pause(0x205, "DF-J/XC (j/vxc matrix)");
539
540 // coulomb energy.
541 // EnergyCoulomb = .5 * Dot(Density.pData, FockTr.pData, FockTr.TriangularStorageSize1());
542
543
544 // FockTr.Clear(); // FIXME: REMOVE THIS.
545 if (pXcFn && !AuxExpandXc) {
546 m_pTimers->Resume(0x204, "DF-J/XC (xc contrib.)");
547 using namespace dfti;
548 FStackMatrix
549 DenC(nAo, nAo, &Mem);
550 Mxm(DenC, OrbC, Transpose(OrbC));
551 DenC.TriangularReduce(1, 0, 1.);
552 FMatrixView
553 OccOrbC(0,0,0); // note: nOcc x nAo with absorbed occupation numbers.
554 OccOrbC = MakeStackMatrix(OrbC.nCols, OrbC.nRows, Mem);
555 Move(OccOrbC, Transpose(OrbC));
556
557 FDftiArgs DftiArgs = {
558 DFTI_MakeXc,
559 &EnergyXc, 0, FockTr.pData, 0,
560 DenC.pData, 0, 0, 0, 0, 0, // no orbitals and no open-shell density provided.
561 pOrbBasisRaw,
562 pXcFn.get(),
563 pDftGrid.get(),
564 g_fThrOrb, MakeLogThrOrb(g_fThrOrb), 0, &fElecTotal // 1e-1 * ThrDen?
565 };
566 AccXc(DftiArgs, m_Log, Mem);
567 // xout << format(pResultFmt) % "Density functional" % DfuEnergy;
568 // xout << format(pResultFmt) % "Functional energy" % DfuEnergy;
569 m_pTimers->Pause(0x204, "DF-J/XC (xc contrib.)");
570 }
571
572 Energy = EnergyCoulomb + EnergyXc;
573
574 FockTr.TriangularExpand();
575 // FockTr.Print(xout, "FOCK + XC");
576 // throw std::runtime_error("absicht.");
577 Add(FockC, FockTr, 1.0);
578
579 m_pTimers->Leave(0x200, "DF-J/XC");
580 }
581
PrintEnergyContribs()582 void FFockComponentBuilderDfCoulXcCached::PrintEnergyContribs()
583 {
584 // m_Log.WriteResult("Number of electrons", fElecTotal);
585 // m_Log.WriteResult("Number of electrons", fElecTotal - double(m_WfDecl.nElec()));
586 // m_Log.WriteResult("Lost electrons", fElecTotal - double(m_WfDecl.nElec));
587 m_Log.Write(" {:<32}{:18.4e}\n", "Lost electrons", fElecTotal - fElecTotalAnalytic);
588
589 m_Log.WriteResult("Coulomb energy", EnergyCoulomb);
590 m_Log.WriteResult("Density functional energy", EnergyXc);
591 }
592
AccGradient(FMatrixView Gradient,FMatrixView COccC,FMatrixView COccO,FMemoryStack & Mem)593 void FFockComponentBuilderDfCoulXcCached::AccGradient(FMatrixView Gradient, FMatrixView COccC, FMatrixView COccO, FMemoryStack &Mem)
594 {
595 // Gradient.Clear(); // FIXME: REMOVE THIS.
596 }
597
598
SwitchToRefineGrid(FDftGridParams const & NewGridParams)599 bool FFockComponentBuilderDfCoulXcCached::SwitchToRefineGrid(FDftGridParams const &NewGridParams)
600 {
601 GridParams = NewGridParams;
602 if (pXcFn) {
603 // find minimum L which should be integrated exactly. Depends on the basis to some degree.
604 // (this is further adjusted inside the grid generator based on atom types)
605 // GridParams.iMinL = AuxExpandXc? (2*pFitBasis->nMaxL()) : (2*pOrbBasis->nMaxL());
606 GridParams.iMinL = AuxExpandXc? (pFitBasis->nMaxL()) : (2*pOrbBasis->nMaxL());
607 // FTimer
608 // tDftGrid;
609 pDftGrid = new FDftGrid(*pAtoms, GridParams, &m_Log);
610 // m_Log.Write(" Generated DFT grid with {} points for {} atoms in {:.2} sec.\n", pDftGrid->Points.size(), Atoms_.size(), (double)tDftGrid);
611 // xout << "" << format(pTimingFmt) % "DFT integration grid" % (double)tDftGrid; xout.flush();
612 }
613 return true;
614 }
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629 } // namespace ct
630