ergo
xc_evaluators.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
35template<typename Matrix>
36struct KsData {
37 Matrix *excmat;
41 KsData(Matrix *m_, size_t s)
42 : excmat(m_), dR(new real[s]), dZ(new real[s]), energy(0)
43 {}
45 delete []dR;
46 delete []dZ;
47 }
48};
49
50
55template<typename Matrix>
57 static void distribute(DftIntegratorBl *grid,
58 int bllen, int blstart, int blend,
59 real * restrict tmp, real *restrict dR,
60 Matrix& excmat)
61 {
62 static const int isym = 0;
63 int jbl, j, ibl, i, k;
64 real * restrict aos = grid->atv;
65
66 int (*restrict blocks)[2] = BASBLOCK(grid,isym);
67 int bl_cnt = grid->bas_bl_cnt[isym];
68
69 for(jbl=0; jbl<bl_cnt; jbl++) {
70 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
71 int joff = j*bllen;
72 for(k=blstart; k<blend; k++)
73 tmp[k] = aos[k+joff]*dR[k]*0.5;
74
75 for(ibl=0; ibl<bl_cnt; ibl++) {
76 int top = blocks[ibl][1] < j
77 ? blocks[ibl][1] : j;
78 for(i=blocks[ibl][0]; i<top; i++) {
79 real *restrict aosi = aos + i*bllen; /* ith orbital */
80 long_real s = 0;
81 for(k=blstart; k<blend; k++)
82 s += aosi[k]*tmp[k];
83 excmat.add(j,i, s);
84 }
85 }
86 long_real s = 0;
87 for(k=blstart; k<blend; k++)
88 s += aos[k+j*bllen]*tmp[k];
89 excmat.add(j,j, s);
90 }
91 }
92 }
93};
94
98template<typename Matrix, typename LDADistributor>
99void
101 int bllen, int blstart, int blend,
102 KsData<Matrix>* data)
103{
104 FirstDrv drvs;
105 int k;
106 real * restrict dR = data->dR;
107 FunDensProp dp = { 0 };
108
109 assert(grid->ntypso >0);
110 for(k=blstart; k<blend; k++) {
111 real weight = grid->weight[grid->curr_point+k];
112 dp.rhoa = dp. rhob = 0.5*grid->r.rho[k];
113 data->energy += selected_func->func(&dp)*weight;
114 dftpot0_(&drvs, &weight, &dp);
115 dR[k] = 2*drvs.fR;
116 }
117 LDADistributor::distribute(grid, bllen, blstart, blend, tmp, dR,
118 *data->excmat);
119}
120
125template<typename Matrix>
127 static void distribute(DftIntegratorBl *grid,
128 int bllen, int blstart, int blend,
129 real * restrict tmp,
130 const real *dR, const real *dZ,
131 Matrix& excmat)
132 {
133 static const int isym = 0;
134 int jbl, j, ibl, i, k;
135 const real * aox = grid->atv+bllen*grid->nbast;
136 const real * aoy = grid->atv+bllen*grid->nbast*2;
137 const real * aoz = grid->atv+bllen*grid->nbast*3;
138 const real * aos = grid->atv;
139
140 const int (*blocks)[2] = BASBLOCK(grid,isym);
141 int nblocks = grid->bas_bl_cnt[isym];
142 for(jbl=0; jbl<nblocks; jbl++)
143 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
144 int joff = j*bllen;
145 for(k=blstart; k<blend; k++)
146 tmp[k+joff] =
147 (dR[k]* aos[k+j*bllen] +
148 dZ[k]*(aox[k+j*bllen]*grid->g.rad.a[k][0]+
149 aoy[k+j*bllen]*grid->g.rad.a[k][1]+
150 aoz[k+j*bllen]*grid->g.rad.a[k][2]));
151 }
152
153 for(jbl=0; jbl<nblocks; jbl++) {
154 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
155 const real * tmpj = tmp + j*bllen; /* jth orbital */
156 for(ibl=0; ibl<nblocks; ibl++) {
157 int top = blocks[ibl][1] < j
158 ? blocks[ibl][1] : j;
159 for(i=blocks[ibl][0]; i<top; i++) {
160 long_real s = 0;
161 const real * tmpi = tmp + i*bllen; /* ith orbital */
162 for(k=blstart; k<blend; k++)
163 s += aos[k+i*bllen]*tmpj[k] +
164 aos[k+j*bllen]*tmpi[k];
165 excmat.add(i, j, s);
166 }
167 }
168 long_real s = 0;
169 for(k=blstart; k<blend; k++)
170 s += 2*aos[k+j*bllen]*tmpj[k];
171 excmat.add(j,j, s);
172 }
173 }
174 }
175};
176
180template<typename Matrix, typename GGADistributor>
181void
183 int bllen, int blstart, int blend,
184 KsData<Matrix>* data)
185{
186 FirstDrv drvs;
187 int k;
188 real * restrict dR = data->dR;
189 real * restrict dZ = data->dZ;
190 FunDensProp dp = { 0 };
191
192 assert(grid->ntypso >0);
193 for(k=blstart; k<blend; k++) {
194 real weight = grid->weight[grid->curr_point+k];
195 dp.grada = 0.5*template_blas_sqrt(grid->g.grad[k][0]*grid->g.grad[k][0]+
196 grid->g.grad[k][1]*grid->g.grad[k][1]+
197 grid->g.grad[k][2]*grid->g.grad[k][2]);
198 dp. rhoa = dp.rhob = 0.5*grid->r.rho[k];
199 dp.gradb = dp.grada;
200 dp.gradab = dp.grada*dp.gradb;
201 if(dp.rhoa>1e-14) {
202 if(dp.grada<1e-35) dp.grada = 1e-35;
203 data->energy += selected_func->func(&dp)*weight;
204 dftpot0_(&drvs, &weight, &dp);
205 dR[k] = 0.5*drvs.fR;
206 dZ[k] = 0.5*drvs.fZ/dp.grada;
207 }else {
208 dR[k] = dZ[k] = 0;
209 }
210 }
211
212 GGADistributor::distribute(grid, bllen, blstart, blend, tmp, dR, dZ,
213 *data->excmat);
214}
215
216/* =================================================================== */
217/* Unrestricted evaluators. */
218
219template<typename Matrix>
220struct UksData {
221 Matrix *exca, *excb;
225 UksData(Matrix * a_, Matrix * b_, size_t s)
226 : exca(a_), excb(b_),
227 dRa(new real[s]), dRb(new real[s]),
228 dZa(new real[s]), dZb(new real[s]), dZab(new real[s]),
229 energy(0.0)
230 {}
232 delete []dRa; delete []dRb;
233 delete []dZa; delete []dZb; delete []dZab;
234 }
235};
236
237template<typename Matrix>
239 static void
240 distribute(DftIntegratorBl *grid, int bllen, int blstart, int blend,
241 real * restrict tmp, const real * dR,
242 const real *dZ1, const real (*grad1)[3],
243 const real *dZ2, const real (*grad2)[3],
244 Matrix& excmat);
245};
246
247template<typename Matrix>
248void
250 int bllen, int blstart, int blend,
251 real * restrict tmp, const real * dR,
252 const real *dZ1,
253 const real (*grad1)[3],
254 const real *dZ2,
255 const real (*grad2)[3],
256 Matrix& excmat)
257{
258 int isym, jbl, j, ibl, i, k;
259 real * restrict aox = grid->atv+bllen*grid->nbast;
260 real * restrict aoy = grid->atv+bllen*grid->nbast*2;
261 real * restrict aoz = grid->atv+bllen*grid->nbast*3;
262 real * restrict aos = grid->atv;
263
264 for(isym=0; isym<grid->nsym; isym++) {
265 int (*restrict blocks)[2] = BASBLOCK(grid,isym);
266 int nblocks = grid->bas_bl_cnt[isym];
267 for(jbl=0; jbl<nblocks; jbl++)
268 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
269 int joff = j*bllen;
270 for(k=blstart; k<blend; k++)
271 tmp[k+joff] =
272 dR[k]* aos[k+j*bllen] +
273 dZ1[k]*(aox[k+j*bllen]*grad1[k][0]+
274 aoy[k+j*bllen]*grad1[k][1]+
275 aoz[k+j*bllen]*grad1[k][2])+
276 dZ2[k]*(aox[k+j*bllen]*grad2[k][0]+
277 aoy[k+j*bllen]*grad2[k][1]+
278 aoz[k+j*bllen]*grad2[k][2]);
279 }
280
281 for(jbl=0; jbl<nblocks; jbl++) {
282 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
283 real *restrict tmpj = tmp + j*bllen; /* jth orbital */
284 for(ibl=0; ibl<nblocks; ibl++) {
285//#define FULL_ONLY 1
286#if defined(FULL_ONLY)
287 for(i=blocks[ibl][0]; i<blocks[ibl][1]; i++) {
288 long_real s = 0;
289 for(k=blstart; k<blend; k++)
290 s += aos[k+i*bllen]*tmpj[k];
291 excmat.add(i, j, s);
292 }
293#else /* FULL_ONLY */
294 int top = blocks[ibl][1] < j
295 ? blocks[ibl][1] : j;
296 for(i=blocks[ibl][0]; i<top; i++) {
297 long_real s = 0;
298 const real * tmpi = tmp + i*bllen; /* ith orbital */
299 for(k=blstart; k<blend; k++)
300 s += aos[k+i*bllen]*tmpj[k] +
301 aos[k+j*bllen]*tmpi[k];
302 excmat.add(i, j, s);
303 }
304#endif /* FULL_ONLY */
305 }
306#if !defined(FULL_ONLY)
307 long_real s = 0;
308 for(k=blstart; k<blend; k++)
309 s += 2*aos[k+j*bllen]*tmpj[k];
310 excmat.add(j,j, s);
311#endif /* FULL_ONLY */
312 }
313 }
314 }
315}
316
320template<typename Matrix, typename LDADistributor>
321void
323 int bllen, int blstart, int blend,
325{
326 FunFirstFuncDrv drvs;
327 FunDensProp dp = { 0 };
328 int k;
329
330 assert(grid->ntypso >0);
331 for(k=blstart; k<blend; k++) {
332 real weight = grid->weight[grid->curr_point+k];
333 dp.rhoa = grid->r.ho.a[k];
334 dp.rhob = grid->r.ho.b[k];
335 d->energy += selected_func->func(&dp)*weight;
336 drv1_clear(&drvs);
337 selected_func->first(&drvs, weight, &dp);
338 d->dRa[k] = 2*drvs.df1000;
339 d->dRb[k] = 2*drvs.df0100;
340 }
341
342 LDADistributor::distribute(grid, bllen, blstart, blend,
343 tmp, d->dRa, *d->exca);
344 LDADistributor::distribute(grid, bllen, blstart, blend,
345 tmp, d->dRb, *d->excb);
346}
347
351template<typename Matrix, typename GGADistributor>
352static void
354 int bllen, int blstart, int blend,
356{
357 FunFirstFuncDrv drvs;
358 int k;
359 FunDensProp dp = { 0 };
360
361 assert(grid->ntypso >0);
362 for(k=blstart; k<blend; k++) {
363 const real THR = 1e-40;
364 real weight = grid->weight[grid->curr_point+k];
365 dp.rhoa = grid->r.ho.a[k];
366 dp.rhob = grid->r.ho.b[k];
367 dp.grada = template_blas_sqrt(grid->g.rad.a[k][0]*grid->g.rad.a[k][0]+
368 grid->g.rad.a[k][1]*grid->g.rad.a[k][1]+
369 grid->g.rad.a[k][2]*grid->g.rad.a[k][2]);
370 dp.gradb = template_blas_sqrt(grid->g.rad.b[k][0]*grid->g.rad.b[k][0]+
371 grid->g.rad.b[k][1]*grid->g.rad.b[k][1]+
372 grid->g.rad.b[k][2]*grid->g.rad.b[k][2]);
373 dp.gradab = grid->g.rad.a[k][0]*grid->g.rad.b[k][0]+
374 grid->g.rad.a[k][1]*grid->g.rad.b[k][1]+
375 grid->g.rad.a[k][2]*grid->g.rad.b[k][2];
376 if(dp.rhoa+dp.rhob>1e-17) {
377 if(dp.grada<THR) dp.grada = THR;
378 if(dp.rhob<THR) dp.rhob = THR;
379 if(dp.gradb<THR) dp.gradb = THR;
380 if(dp.gradab<THR) dp.gradab = THR;
381 d->energy += selected_func->func(&dp)*weight;
382 drv1_clear(&drvs);
383 selected_func->first(&drvs, weight, &dp);
384 d->dRa[k] = drvs.df1000*0.5;
385 d->dRb[k] = drvs.df0100*0.5;
386 d->dZa[k] = drvs.df0010/dp.grada;
387 d->dZb[k] = drvs.df0001/dp.gradb;
388 d->dZab[k] = drvs.df00001;
389 } else {
390 d->dRa[k] = d->dRb[k] = d->dZa[k] = d->dZb[k] = d->dZab[k] = 0;
391 }
392 }
393
394 GGADistributor::distribute(grid, bllen, blstart, blend, tmp, d->dRa,
395 d->dZa, grid->g.rad.a, d->dZab,
396 grid->g.rad.b, *d->exca);
397 GGADistributor::distribute(grid, bllen, blstart, blend, tmp, d->dRb,
398 d->dZb, grid->g.rad.b, d->dZab,
399 grid->g.rad.a, *d->excb);
400}
#define restrict
Definition: config.h:283
ergo_real real
Definition: test.cc:46
EXTERN_C void dftpot0_(FirstDrv *ds, const real *weight, const FunDensProp *dp)
Definition: dft_common.cc:766
#define THR
Definition: fun-cam.c:75
Functional * selected_func
Definition: functionals.c:106
void drv1_clear(FunFirstFuncDrv *gga)
Definition: functionals.c:131
ergo_long_real long_real
Definition: grid_atomic.h:43
#define BASBLOCK(grid, isym)
Definition: integrator.h:61
Definition: integrator.h:49
int nbast
Definition: integrator.h:72
real * rho
Definition: integrator.h:76
int nsym
Definition: integrator.h:58
real * a
Definition: integrator.h:78
real * weight
Definition: integrator.h:52
struct DftIntegratorBl_::@1::@3 rad
struct DftIntegratorBl_::@0::@2 ho
int ntypso
Definition: integrator.h:63
real * atv
Definition: integrator.h:53
union DftIntegratorBl_::@0 r
union DftIntegratorBl_::@1 g
int curr_point
Definition: integrator.h:89
real(* grad)[3]
Definition: integrator.h:82
int bas_bl_cnt[8]
Definition: integrator.h:58
real * b
Definition: integrator.h:78
A vector of first order derivatives with respect to two parameters: density rho and SQUARE of the gra...
Definition: dft_common.h:60
real fR
Definition: dft_common.h:61
real fZ
Definition: dft_common.h:62
Definition: functionals.h:375
real rhob
Definition: functionals.h:376
real gradab
Definition: functionals.h:378
real rhoa
Definition: functionals.h:376
real gradb
Definition: functionals.h:377
real grada
Definition: functionals.h:377
Definition: functionals.h:119
real df0010
Definition: functionals.h:122
real df0100
Definition: functionals.h:121
real df00001
Definition: functionals.h:124
real df0001
Definition: functionals.h:123
real df1000
Definition: functionals.h:120
FirstOrderFun first
Definition: functionals.h:410
EnergyFunc func
Definition: functionals.h:409
structure describing the data needed by distributors.
Definition: xc_evaluators.h:36
real * dZ
Definition: xc_evaluators.h:39
~KsData()
Definition: xc_evaluators.h:44
KsData(Matrix *m_, size_t s)
Definition: xc_evaluators.h:41
real energy
Definition: xc_evaluators.h:40
real * dR
Definition: xc_evaluators.h:38
Matrix * excmat
Definition: xc_evaluators.h:37
Definition: xc_evaluators.h:220
real energy
Definition: xc_evaluators.h:224
real * dZa
Definition: xc_evaluators.h:223
Matrix * excb
Definition: xc_evaluators.h:221
real * dRb
Definition: xc_evaluators.h:222
~UksData()
Definition: xc_evaluators.h:231
real * dZb
Definition: xc_evaluators.h:223
Matrix * exca
Definition: xc_evaluators.h:221
real * dRa
Definition: xc_evaluators.h:222
real * dZab
Definition: xc_evaluators.h:223
UksData(Matrix *a_, Matrix *b_, size_t s)
Definition: xc_evaluators.h:225
Definition: xc_evaluators.h:238
static void distribute(DftIntegratorBl *grid, int bllen, int blstart, int blend, real *restrict tmp, const real *dR, const real *dZ1, const real(*grad1)[3], const real *dZ2, const real(*grad2)[3], Matrix &excmat)
Definition: xc_evaluators.h:249
distributes a GGA-type xc potential over the XC-matrix elements.
Definition: xc_evaluators.h:126
static void distribute(DftIntegratorBl *grid, int bllen, int blstart, int blend, real *restrict tmp, const real *dR, const real *dZ, Matrix &excmat)
Definition: xc_evaluators.h:127
distributes a LDA-type xc potential over the XC-matrix elements, with optimization for a closed shell...
Definition: xc_evaluators.h:56
static void distribute(DftIntegratorBl *grid, int bllen, int blstart, int blend, real *restrict tmp, real *restrict dR, Matrix &excmat)
Definition: xc_evaluators.h:57
Treal template_blas_sqrt(Treal x)
void xcCallbackLdaR(DftIntegratorBl *grid, real *restrict tmp, int bllen, int blstart, int blend, KsData< Matrix > *data)
modifies data->excmat by adding LDA-type contributions from a given set of bllen grid points as saved...
Definition: xc_evaluators.h:100
static void xcCallbackGgaU(DftIntegratorBl *grid, real *restrict tmp, int bllen, int blstart, int blend, UksData< Matrix > *d)
modifes data->excmat by adding GGA-type xc contributions, for a general unrestricted case.
Definition: xc_evaluators.h:353
void xcCallbackLdaU(DftIntegratorBl *grid, real *restrict tmp, int bllen, int blstart, int blend, UksData< Matrix > *d)
modifies data->excmat by adding LDA-type xc contributions, for a general unrestricted case.
Definition: xc_evaluators.h:322
void xcCallbackGgaR(DftIntegratorBl *grid, real *restrict tmp, int bllen, int blstart, int blend, KsData< Matrix > *data)
modifies data->excmat by adding GGA-type contributions from a given set of bllen grid points as saved...
Definition: xc_evaluators.h:182