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 <boost/format.hpp>
26 #include <cmath>
27
28 #include "CxDefs.h"
29 #include "CxTypes.h"
30
31 #include "CtAtomSet.h"
32 #include "CtDftGrid_ivb.h"
33 #include "CxLebedevGrid.h"
34 #include "CtConstants.h"
35 #include "CtTiming.h"
36
37 using boost::format;
38
39 namespace ct {
40
41
42 // main references:
43 // [1] JCP 102 346 (1995) (Treutler & Ahlrichs)
44 // [2] JCP 108 3226 (1998) (Krack & Koester)
45 // [3] JCP 88 2547 (1988) (Becke)
46 // [4] JCP 104 9848 (1996) (Mura & Knowles)
47 // [5] JCP 101 8894 (1994) (Delley)
48 // [6] CPL 257 213 (1996) (Stratmann, Scuseria, Frisch). This is what Scheffler recommend for partitioning.
49 // (Stratmann, Scuseria, Frisch - Achieving linear scaling in exchange-correlation density functional quadratures.pdf)
50 // [7] See also "Ab initio molecular simulations with numeric atom-centered orbitals" by Scheffler et al,
51 // who give some notes on their grids (mostly referring back to Delly)
52
53
54 // maybe better:
55 // JCP 121 681 (2004) (Koester, Flores-Moreno, Reveles)
56
57 // Further refs which migt be useful:
58 // - Izsak and Neese describe in http://dx.doi.org/10.1063/1.3646921
59 // how they set up grid for COSX computations (Section IV A). Might be worth investigating.
60 // - NWChem's description:
61 // http://www.nwchem-sw.org/index.php/Density_Functional_Theory_for_Molecules#GRID_--_Numerical_Integration_of_the_XC_Potential
62 // - Kakhiani, Tsereteli, Tsereteli -- A program to generate a basis set adaptive radial quadrature grid for density functional theory
63 // 10.1016/j.cpc.2008.10.004
64 // - Lindh, Malmqvist, Gagliardi - Molecular integrals by numerical quadrature. I. Radial integration.pdf
65 // (2001), Theor Chem Acc (2001) 106:17 10.1007/s002140100263
66 // - Weber, Daul, Baltensperger - Radial numerical integrations based on the sinc function.pdf
67 // 10.1016/j.cpc.2004.08.008
68 // - El-Sherbiny, Poirier - An evaluation of the radial part of numerical integration commonly used in DFT.pdf
69 // - The Becke Fuzzy Cells Integration Scheme in the Amsterdam Density Functional Program Suite
70 // http://onlinelibrary.wiley.com/doi/10.1002/jcc.23323/full
71 // - Note: This also has a set of TA Eta scaling parameters for the entire periodic table
72
73
fTargetAccuracy() const74 double FDftGridParams::fTargetAccuracy() const {
75 // would return 1e-5 for level 3 and 1e-6 for level 5
76 return 1e-3 * std::pow(.1, double(this->nLevel+1)/2.);
77 }
78
SetTargetAccuracy(double fAcc)79 void FDftGridParams::SetTargetAccuracy(double fAcc) {
80 // fAcc = 1e-3 * std::pow(.1, double(this->nLevel+1)/2.);
81 // 1e3 * fAcc = std::pow(.1, double(this->nLevel+1)/2.);
82 // std::log(1e3 * fAcc)/std::log(.1) = double(this->nLevel+1)/2.
83 // 2 * std::log(1e3 * fAcc)/std::log(.1) = double(this->nLevel+1)
84 // 2 * std::log(1e3 * fAcc)/std::log(.1) - 1 = this->nLevel
85 nLevel = 2. * (std::log(1e3 * fAcc)/std::log(.1)) - 1.;
86 assert_rt(std::abs(fAcc - fTargetAccuracy()) < 1e-10);
87 }
88
89 // separated from FDftGrid in order to decouple grid generation from grid usage.
90 struct FDftGridGenerator
91 {
92 typedef FDftGrid::FPoint
93 FPoint;
94 typedef FDftGrid::FPointList
95 FPointList;
96 typedef FDftGrid::FGridBlock
97 FGridBlock;
98 typedef FDftGrid::FGridBlockList
99 FGridBlockList;
100
101 FAtomSet const
102 &Atoms;
103 FDftGridParams const
104 &Params;
105 FDftGrid
106 &Grid;
107 FPointList
108 &Points;
109 FGridBlockList
110 &GridBlocks;
111
112
FDftGridGeneratorct::FDftGridGenerator113 FDftGridGenerator( FDftGrid &Grid_, FAtomSet const &Atoms_,
114 FDftGridParams const &Params_ )
115 : Atoms(Atoms_), Params(Params_), Grid(Grid_), Points(Grid_.Points ), GridBlocks(Grid_.GridBlocks)
116 {}
117
118 void Create();
119
120 FScalar GetAtomWeight( FVector3 const &vPos, uint iAtom, double const *pInvDistAt, FMemoryStack &Mem );
121 static FScalar GetAtomPeriodRowFraction( uint ElementNumber );
122 void GetAtomGridParams( uint &nRadialPt, double &AtomicScale, uint &iAngGrid, uint iAtom );
123 void GetAtomRadialGrid( double *r, double *w, uint n, double AtomicScale );
124 void AddAtomGrid( FPointList &Points, uint iAtom, double const *pInvDistAt, FMemoryStack &Mem);
125
126 FScalar GetPairVoronoiR(double Mu, size_t iAtom, size_t iOtherAtom);
127
128 void BlockifyGridR( FGridBlockList &Blocks, uint iFirst, FPoint *pFirst, FPoint *pLast );
129 };
130
131
132 // // JCP 41 3199 (1964). In Angstrom (of the time, strictly)
133 // static FScalar SlaterBraggAtomicRadii[] = {
134 // 0.35, // modified recommendation by Becke in [3]
135 // 0,
136 // 1.45, 1.05, 0.85, 0.70, 0.65, 0.60, 0.50, // Li--F
137 // 0,
138 // 1.8, 1.5, 1.25, 1.1, 1.0, 1.0, 1.0,
139 // 0,
140 // 2.2, 1.8, 1.6, 1.4, 1.35, 1.4, 1.4, 1.4, 1.35, 1.35, 1.35, 1.35, 1.30, 1.25, 1.15, 1.15, 1.15, // K--Br
141 // 0,
142 // 2.35, 2.0, 1.8, 1.55, 1.45, 1.35, 1.35, 1.30, 1.35, 1.40, 1.60, 1.55, 1.55, 1.45, 1.45, 1.40, 1.40, // Rb-I
143 // 2.6, 2.15, 1.95, 1.85, 1.85, 1.85, 1.85, 1.85, 1.85, 1.8, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.55, 1.45, 1.35, 1.35, 1.30, 1.35, 1.35, 1.35, 1.5, 1.9, 1.6, 1.9, // Cs--Po
144 // 0.0, 0.0, 0.0, // At, Rn, Fr
145 // 2.15,1.95,1.8,1.8,1.75,1.75,1.75,1.75 // Ra-Am
146 // };
147
148 // from wikipedia... : http://en.wikipedia.org/wiki/Covalent_radius
149 // Looks better than Slaters values, imo. Original citation:
150 // Beatriz Cordero, Verónica Gómez, Ana E. Platero-Prats, Marc Revés, Jorge Echeverría, Eduard Cremades, Flavia Barragán and Santiago Alvarez. Covalent radii revisited. Dalton Trans., 2008, 2832-2838, doi:10.1039/b801115j
151 // Multiple entries for a atom (e.g., C,sp..sp3) have been averaged.
152 static FScalar CovalentRadii[] = { // H-Cm.
153 0.31, 0.28, 1.28, 0.96, 0.84, 0.73, 0.71, 0.66, 0.57, 0.58, 1.66, 1.41,
154 1.21, 1.11, 1.07, 1.05, 1.02, 1.06, 2.03, 1.76, 1.70, 1.60, 1.53, 1.39,
155 1.50, 1.42, 1.38, 1.24, 1.32, 1.22, 1.22, 1.20, 1.19, 1.20, 1.20, 1.16,
156 2.20, 1.95, 1.90, 1.75, 1.64, 1.54, 1.47, 1.46, 1.42, 1.39, 1.45, 1.44,
157 1.42, 1.39, 1.39, 1.38, 1.39, 1.40, 2.44, 2.15, 2.07, 2.04, 2.03, 2.01,
158 1.99, 1.98, 1.98, 1.96, 1.94, 1.92, 1.92, 1.89, 1.90, 1.87, 1.87, 1.75,
159 1.70, 1.62, 1.51, 1.44, 1.41, 1.36, 1.36, 1.32, 1.45, 1.46, 1.48, 1.40,
160 1.50, 1.50, 2.60, 2.21, 2.15, 2.06, 2.00, 1.96, 1.90, 1.87, 1.80, 1.69
161 };
162
GetPairVoronoiR(double Mu,size_t iAtom,size_t iOtherAtom)163 FScalar FDftGridGenerator::GetPairVoronoiR(double Mu, size_t iAtom, size_t iOtherAtom)
164 {
165 // [1], eq. (4)--(8)
166 double g = Mu;
167 if ( Params.AdjustByAtomicRadii && false ){
168 #ifdef _DEBUG
169 uint const
170 nRadii = sizeof(CovalentRadii)/sizeof(FScalar);
171 assert( Atoms[iAtom].AtomicNumber < nRadii );
172 assert( Atoms[iOtherAtom].AtomicNumber < nRadii );
173 #endif
174 FScalar
175 R1 = CovalentRadii[Atoms[iAtom].AtomicNumber - 1],
176 R2 = CovalentRadii[Atoms[iOtherAtom].AtomicNumber - 1],
177 X = 1.0;
178 if ( R1 != 0.0 && R2 != 0.0 )
179 // ^- data might not be there if Slater values are used.
180 X = std::sqrt( R1/R2 ); // [1], eq. (13)
181 FScalar
182 u = (X - 1.0)/(X + 1.0),
183 a = u/(u*u - 1.0);
184 if ( a < -0.5 ) a = -0.5; // [3], eq. (A3-6)ff
185 if ( a > 0.5 ) a = 0.5;
186
187 g = g + a * (1.0 - g*g); // [3], eq. (A2)
188 }
189
190 if (0) {
191 // Original Becke scheme.
192 // [1], eq. (7) (k=3)
193 g = ( (3./2.)*g - .5*g*g*g );
194 g = ( (3./2.)*g - .5*g*g*g );
195 g = ( (3./2.)*g - .5*g*g*g );
196 } else {
197 // Stratmann's modified Becke scheme [6]
198 double a = 0.64; // comment after eq. 14
199 // [6], eq.11
200 if (g <= -a)
201 g = -1.;
202 else if (g >= a)
203 g = +1;
204 else {
205 // [6], eq.14
206 double ma = g/a; // mu/a
207 g = (1/16.)*(ma*(35 + ma*ma*(-35 + ma*ma*(21 - 5 *ma*ma))));
208 }
209 }
210
211 double s = .5 * ( 1 - g );
212 return s;
213 }
214
215 void ArgSort1(size_t *pOrd, double const *pVals, size_t nValSt, size_t nVals, bool Reverse);
216
217 // this function subdivides the entire space volume into atomic contributions.
218 // Returns weight of atom iAtom at point vPos.
219 // The weight returned by this function is /not normalized/ [to the molecule,
220 // only to the atom].
GetAtomWeight(FVector3 const & vGridPos,uint iAtom,double const * pInvDistAt,FMemoryStack & Mem)221 FScalar FDftGridGenerator::GetAtomWeight( FVector3 const &vGridPos, uint iAtom, double const *pInvDistAt, FMemoryStack &Mem )
222 {
223 // compute distance from grid point to all other atoms.
224 size_t
225 nAt = Atoms.size();
226 double
227 *pDistAg;
228 Mem.Alloc(pDistAg, nAt);
229 for (size_t iAt = 0; iAt < nAt; ++ iAt)
230 pDistAg[iAt] = Dist(vGridPos, Atoms[iAt].vPos);
231
232 // sort atoms by distance from grid point. close ones first.
233 size_t
234 *pAtOrd;
235 Mem.Alloc(pAtOrd, nAt);
236 ArgSort1(pAtOrd, pDistAg, 1, nAt, false);
237
238 FScalar
239 wOut = 0.,
240 wTotal = 0.;
241
242 // [1], eq. (4)--(8)
243 for ( size_t iCenterAtom = 0; iCenterAtom != Atoms.size(); ++ iCenterAtom) {
244 double wCen = 1.;
245
246 for ( size_t iOtherAtom_ = 0; iOtherAtom_ != Atoms.size(); ++ iOtherAtom_ ) {
247 size_t
248 iOtherAtom = pAtOrd[iOtherAtom_];
249
250 if ( iCenterAtom == iOtherAtom )
251 continue;
252 FScalar
253 InvR12 = pInvDistAt[iCenterAtom * nAt + iOtherAtom],
254 r1 = pDistAg[iCenterAtom],
255 r2 = pDistAg[iOtherAtom],
256 // ``confocal elliptical coordinates''. A great word. [3], eq.(9)
257 Mu = (r1-r2)*InvR12;
258 // Mu = (r1-r2)/R12;
259
260 wCen *= GetPairVoronoiR(Mu, iCenterAtom, iOtherAtom);
261 if (wCen < 1e-12) {
262 wCen = 0;
263 break;
264 }
265 }
266 if (iCenterAtom == iAtom)
267 wOut = wCen;
268 wTotal += wCen;
269 }
270
271 Mem.Free(pDistAg);
272 if (wTotal == 0.)
273 return 0.;
274 return wOut/wTotal;
275 }
276
277
278 // 10.1002/jcc.23323, supp info Tab. 1 (replacement for [1], Tab 1.)
279 static FScalar EtaRescalingParameters[] = {
280 0.8, 0.9, 1.8, 1.4, 1.3, 1.1, 0.9, 0.9, 0.9, 0.9, 1.4, 1.3, 1.3, 1.2, 1.1, 1.0, 1.0, 1.0, 1.5, 1.4, 1.3, 1.2, 1.2, 1.2, 1.2, 1.2, // H - Fe
281 1.2, 1.1, 1.1, 1.1, 1.1, 1.0, 0.9, 0.9, 0.9, 0.9, 1.4, 1.4, 1.1, 1.3, 1.0, 1.2, 0.9, 0.9, 0.9, 1.0, 0.9, 1.0, 1.0, 1.3, 1.2, 1.2, // Co-Te
282 0.9, 1.0, 1.7, 1.5, 1.5, 1.3, 1.3, 1.4, 1.8, 1.4, 1.2, 1.3, 1.3, 1.4, 1.1, 1.1, 1.2, 1.6, 1.4, 1.3, 1.2, 1.0, 1.0, 0.9, 1.3, 1.2, // I - Pt
283 1.2, 1.0, 1.2, 1.2, 1.1, 1.2, 1.1, 2.1, 2.2, 1.8, 1.7, 1.3, 1.4, 1.2, 1.2, 1.3, 1.4, 1.4, 1.7, 1.9, 1.9, 2.0, 2.0, 1.6, 2.0 // Au-Lw
284 };
285
286
GetAtomPeriodRowFraction(uint ElementNumber)287 FScalar FDftGridGenerator::GetAtomPeriodRowFraction( uint ElementNumber )
288 {
289 static uint
290 RareGasAtomicNumbers[] = { 0, 2, 10, 18, 36, 54, 86 };
291 static uint const
292 N = sizeof(RareGasAtomicNumbers)/sizeof(uint);
293 if ( ElementNumber == RareGasAtomicNumbers[N-1] )
294 return 1.0;
295 if ( ElementNumber > RareGasAtomicNumbers[N-1] )
296 return 0.5;
297 for ( uint i = 1; i < N; ++ i ){
298 if ( RareGasAtomicNumbers[i] >= ElementNumber ){
299 uint
300 nFirst = RareGasAtomicNumbers[i-1],
301 nLast = RareGasAtomicNumbers[i];
302 return static_cast<FScalar>(ElementNumber - nFirst)/(nLast - nFirst);
303 }
304 }
305 assert_rt(0);
306 return 1.0;
307 }
308
GetAtomPeriod(uint ElementNumber)309 uint GetAtomPeriod(uint ElementNumber)
310 {
311 // returns 1 for H, He; 2 for Li-Ne, 3 for ...
312 static uint
313 RareGasAtomicNumbers[] = { 0, 2, 10, 18, 36, 54, 86 };
314 static uint const
315 N = sizeof(RareGasAtomicNumbers)/sizeof(uint);
316 for ( uint i = 1; i < N; ++ i ){
317 if ( ElementNumber <= RareGasAtomicNumbers[i] )
318 return i;
319 }
320 assert(0);
321 return N;
322 }
323
324
GetAtomGridParams(uint & nRadialPt,double & AtomicScale,uint & iAngGrid,uint iAtom)325 void FDftGridGenerator::GetAtomGridParams(uint &nRadialPt, double &AtomicScale, uint &iAngGrid, uint iAtom)
326 {
327 FAtom const
328 &Atom = Atoms.Atoms[iAtom];
329 if ( 1 ) {
330 if ( Atom.AtomicNumber <= sizeof(EtaRescalingParameters)/sizeof(FScalar) )
331 AtomicScale = EtaRescalingParameters[Atom.AtomicNumber - 1];
332 else {
333 // lerp with 1.5 for alkali metals to 0.9 for rare gases.
334 FScalar
335 f = GetAtomPeriodRowFraction(Atom.AtomicNumber);
336 f *= f;
337 AtomicScale = f*1.5 + (1-f)*0.9;
338 }
339 } else
340 // AtomicScale = std::sqrt(CovalentRadii[Atom.AtomicNumber-1]);
341 AtomicScale = CovalentRadii[Atom.AtomicNumber-1] * (1/ToAng);
342
343 {
344 // if (GetAtomPeriod(Atom.AtomicNumber) != 0)
345 // iL += 1; // use larger angular grid.
346 uint DefaultAngularGridLs[] = {9,11,17,23,29,35,47,59,71,89};
347 uint iRow = GetAtomPeriod(Atom.AtomicNumber);
348 // std::cout << format(" %s is row %i.\n") % Atom.ElementName() % iRow;
349 uint iOffs = Params.nLevel-1 + iRow-1;
350 // std::cout << format(" %s is row %i. grid offs: %i\n") % Atom.ElementName() % iRow % iOffs;
351 uint iMaxL = DefaultAngularGridLs[iOffs];
352 for (iAngGrid = 0; iAngGrid < nAngularGrids; ++ iAngGrid) {
353 if (AngularGridInfo[iAngGrid].MaxL == iMaxL)
354 break;
355 }
356
357 // [5], text after 12c.
358 // def nf(sp,iat): return int(sp*14.*(iat+2)**(1./3.))
359 // >>> [nf(1.,iat) for iat in [1,10,18,36]]
360 // [20, 32, 38, 47]
361 // TA grid 1 has 20, 25, 30, 35.
362 // TA grid 3 has 30, 35, 40, 45.
363 // >>> [nf(0.75,iat) for iat in [1,10,18,36]]
364 // [15, 24, 28, 35]
365 // >>> [nf(0.75+2*.25,iat) for iat in [1,10,18,36]]
366 // [25, 40, 47, 58]
367 // >>> [nf(0.75+2*.2,iat) for iat in [1,10,18,36]]
368 // [23, 36, 43, 54]
369
370 nRadialPt = (int)((0.75 + ((signed)Params.nLevel-1)*0.2) * 14. * std::pow((double)Atom.AtomicNumber+2., 1./3.));
371 // ^- should produce more or less comparable radial sizes to TA.
372 }
373 // double
374 // EpsTol = std::pow(0.1, 3 + Params.nLevel); // target energy tolerance.
375 // n = std::max(20u, (uint)(-5. * (3. * std::log(EpsTol) - (double)iRow + 6))); // JCP 121 681 eq. 20
376 // _xout0("EpsTol = " << EpsTol << " n = " << n);
377 }
378
379
380 /* this one is Delley's grid [5]/[7]. */
GetAtomRadialGrid(double * r,double * w,uint n,double AtomicScale)381 void FDftGridGenerator::GetAtomRadialGrid( double *r, double *w, uint n, double AtomicScale )
382 {
383 FScalar
384 r_outer = 12.0 * AtomicScale, // given in bohr? sounds a bit small.
385 den = 1/FScalar(1+n),
386 RFac = r_outer/std::log(1 - sqr(n*den));
387
388 for ( uint i = 1; i <= n; ++ i ){ // <- beware of indices starting at 1!
389 FScalar
390 xi,
391 // derivative d[ri]/d[xi] (for weight)
392 dri;
393 xi = RFac * std::log(1-sqr(i*den)); // note: is 0 for i = 0 and +inf for n+1.
394 dri = RFac/((1-sqr(i*den))) * (-2.0*i*sqr(den)); // d/d[s] xi
395
396 r[i-1] = xi;
397 w[i-1] = dri * 4.*M_PI*r[i-1]*r[i-1]; // and that is the radial volume element.
398 // ^- this seems to work fine, but should there not be a scale for the number of
399 // integration points in dri? Or is this implicit in transforming directly from
400 // int[i=0...n] to xi?
401 }
402 }
403
AddAtomGrid(FPointList & Points,uint iAtom,double const * pInvDistAt,FMemoryStack & Mem)404 void FDftGridGenerator::AddAtomGrid( FPointList &Points, uint iAtom, double const *pInvDistAt, FMemoryStack &Mem )
405 {
406 FAtom const
407 &Atom = Atoms.Atoms[iAtom];
408 std::vector<double>
409 ri, wi;
410 std::vector<uint>
411 nAngPts;
412 double
413 AtomicScale;
414 uint
415 nr,
416 iBaseAngGrid;
417
418 GetAtomGridParams(nr, AtomicScale, iBaseAngGrid, iAtom);
419 // _xout0("iAtom " << iAtom << " nr: " << nr << " iBaseL: " << iBaseL );
420 ri.resize(nr);
421 wi.resize(nr);
422 GetAtomRadialGrid( &ri[0], &wi[0], nr, AtomicScale );
423 // _xout0("radial grid:");
424 // for ( uint i = 0; i < nr; ++ i )
425 // _xout0( boost::format("r[%2i] = %10.5f (w[i]=%10.5f)") % i % ri[i] % wi[i] );
426
427 // assign radial grid for each radius and count total number of points.
428 uint
429 nPts = 0,
430 nMaxLMax = 0,
431 nAngPtsMax = 0;
432 double fAtomTotalWeight = 0;
433 for ( uint iShell = 1; iShell <= nr; ++ iShell ){
434 // BEWARE of iShell starting at 1!
435 uint iAngGrid = iBaseAngGrid;
436 // reduce grid size for inner shells ([1], eq. 37 and 38)
437 // (a proper adaptive integration would surely be a better way of dealing with this...)
438 if ( iShell <= nr/3 ) {
439 if (iAngGrid) --iAngGrid;
440 if (iAngGrid) --iAngGrid;
441 if (iAngGrid) --iAngGrid;
442 if (iAngGrid) --iAngGrid;
443 } else if ( iShell <= nr/2 ) {
444 if (iAngGrid) --iAngGrid;
445 if (iAngGrid) --iAngGrid;
446 }
447
448 nMaxLMax = std::max(nMaxLMax, uint(AngularGridInfo[iAngGrid].MaxL));
449 nAngPts.push_back(AngularGridInfo[iAngGrid].nPoints);
450 nPts += nAngPts.back();
451 // _xout0( iShell << " of " << nr << " iL = " << iL << " grd-sz = " << nAngPts.back() );
452 nAngPtsMax = std::max(nAngPtsMax, nAngPts.back());
453 };
454 // std::cout << format(" Atom %3i (%2s): nr = %3i MaxL = %2i (nAng = %s)") % (1+iAtom) % Atom.ElementName() % nr % nMaxLMax % nAngPtsMax << std::endl;
455
456 double
457 (*pAngPts)[4] = reinterpret_cast<double (*)[4]>(::malloc( 4 * sizeof(FScalar) * nAngPtsMax ));
458 uint
459 nLastAng = 0xffffffff;
460
461 for ( uint iShell = 0; iShell < nr; ++ iShell ){
462 double
463 fRad = ri[iShell], // current radius
464 fRadWeight = wi[iShell]; // weight for the radial integration.
465 uint
466 nAng = nAngPts[iShell]; // current number of angular points
467
468 // make lebedev grid if different from the one before.
469 if ( nAng != nLastAng )
470 nAng = MakeAngularGrid(pAngPts, nAng);
471 nLastAng = nAng;
472
473 double
474 fShellTotalWeight = 0;
475
476 // generate output points.
477 for ( uint i = 0; i < nAng; ++ i ){
478 double (&p)[4] = pAngPts[i];
479 fShellTotalWeight += p[3];
480 FPoint out;
481 out.vPos[0] = fRad * p[0] + Atom.vPos[0];
482 out.vPos[1] = fRad * p[1] + Atom.vPos[1];
483 out.vPos[2] = fRad * p[2] + Atom.vPos[2];
484 out.fWeight = p[3] * fRadWeight * GetAtomWeight( out.vPos, iAtom, pInvDistAt, Mem );
485 // _xout0( boost::format("vPos=(%10.4f %10.4f %10.4f) w(1) = %12.6f w(2) = %12.6f sum = %12.6f")
486 // % out.vPos[0] % out.vPos[1] % out.vPos[2] % GetAtomWeight(out.vPos, 0) %
487 // GetAtomWeight(out.vPos, 1) % (GetAtomWeight(out.vPos, 0)+GetAtomWeight(out.vPos, 1))
488 // );
489 fAtomTotalWeight += out.fWeight;
490 if ( std::abs(out.fWeight) > 1e-12 )
491 Points.push_back(out);
492 };
493 assert_rt(std::abs(fShellTotalWeight-1.0)<1e-13);
494 }
495 ::free(pAngPts);
496 // _xout0("Atom " << iAtom << " total weight: " << fAtomTotalWeight);
497 }
498
499
Create()500 void FDftGridGenerator::Create()
501 {
502 size_t
503 nAt = Atoms.size();
504 FMemoryStack2
505 Mem(2000000 + sizeof(double) * nAt*nAt);
506
507 // compute inverse distance between all atoms.
508 double
509 *pInvDistAt;
510 Mem.Alloc(pInvDistAt, nAt*nAt);
511 for (size_t iAt = 0; iAt < nAt; ++ iAt)
512 for (size_t jAt = 0; jAt < nAt; ++ jAt) {
513 double rij = Dist(Atoms[iAt].vPos, Atoms[jAt].vPos);
514 if ( rij != 0 )
515 rij = 1./rij;
516 else
517 rij = 0.;
518 pInvDistAt[nAt*iAt + jAt] = rij;
519 pInvDistAt[nAt*jAt + iAt] = rij;
520 }
521
522
523 double
524 fTotalWeight = 0;
525 for ( uint iAtom = 0; iAtom < Atoms.size(); ++ iAtom ) {
526 uint
527 iOff = Points.size();
528 AddAtomGrid(Points, iAtom, pInvDistAt, Mem);
529 for ( uint i = iOff; i < Points.size(); ++ i )
530 fTotalWeight += Points[i].fWeight;
531 // _xout0(boost::format("Atom #%i: Grid points %i--%i") % iAtom % iOff % Points.size());
532 }
533 Mem.Free(pInvDistAt);
534 BlockifyGridR( GridBlocks, 0, 0, 0 );
535 // _xout0(format("Generated integration grid with %i points") % Points.size());
536 // _xout0("Grid integrated volume^(1/3): " << fmt::ff(std::pow(fTotalWeight,1./3.),12,4));
537 }
538
BlockifyGridR(FGridBlockList & Blocks,uint,FPoint *,FPoint *)539 void FDftGridGenerator::BlockifyGridR( FGridBlockList &Blocks, uint /*iFirst*/, FPoint */*pFirst*/, FPoint */*pLast*/ )
540 {
541 // recursively sub-divide points in [pFirst,pLast) at some axis (in a kd-tree fashion)
542 // such that
543 // (a) the spatial extend of the grid points in the two sub-ranges is minimized
544 // (b) a more-or-less uniform number of total points per block is formed.
545
546 // atm: don't do anything, just make fixed size blocks without any
547 // regard to reasonableness.
548 uint
549 iPt = 0,
550 nTargetPtsPerBlock = 128;
551 while ( iPt < Points.size() ){
552 uint
553 iPtEnd = std::min(iPt + nTargetPtsPerBlock, (uint)Points.size());
554 Blocks.push_back( FGridBlock() );
555 FGridBlock
556 &Block = Blocks.back();
557 Block.iFirst = iPt;
558 Block.iLast = iPtEnd;
559 Block.vCenter = FVector3(0,0,0);
560 Block.fLargestWeight = 0;
561 for ( uint i = iPt; i < iPtEnd; ++ i )
562 Block.vCenter += (1.0/(iPtEnd-iPt)) * Points[i].vPos;
563 Block.fRadius = 0;
564 for ( uint i = iPt; i < iPtEnd; ++ i ) {
565 double fDist1 = (Points[i].vPos - Block.vCenter).LengthSq();
566 Block.fRadius = std::max( Block.fRadius, fDist1 );
567 Block.fLargestWeight = std::max(Block.fLargestWeight, Points[i].fWeight);
568 }
569 Block.fRadius = std::sqrt(Block.fRadius);
570
571 iPt = iPtEnd;
572 }
573 }
574
575
576
FDftGrid(FAtomSet const & Atoms,FDftGridParams const & Params,ct::FLog * pLog)577 FDftGrid::FDftGrid( FAtomSet const &Atoms, FDftGridParams const &Params, ct::FLog *pLog )
578 {
579 FLogStdStream
580 xLog(xout);
581 if (pLog == 0)
582 pLog = &xLog;
583
584 FTimer tDftGrid;
585 FDftGridGenerator(*this, Atoms, Params).Create();
586 MakeAdditionalRepresentations();
587 pLog->Write(" Generated DFT grid with {} points for {} atoms in {:.2} sec.\n", Points.size(), Atoms.size(), (double)tDftGrid);
588 }
589
590
~FDftGrid()591 FDftGrid::~FDftGrid()
592 {
593 }
594
595
MakeAdditionalRepresentations()596 void FDftGrid::MakeAdditionalRepresentations()
597 {
598 Positions.resize(Points.size());
599 Weights.resize(Points.size());
600 for ( uint i = 0; i < Points.size(); ++ i ){
601 Positions[i][0] = Points[i].vPos[0];
602 Positions[i][1] = Points[i].vPos[1];
603 Positions[i][2] = Points[i].vPos[2];
604 Weights[i] = Points[i].fWeight;
605 }
606 }
607
608
609
610 } // namespace ct
611