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  */
35 template<typename Matrix>
36 struct KsData {
37  Matrix *excmat;
38  real* dR;
39  real* dZ;
41  KsData(Matrix *m_, size_t s)
42  : excmat(m_), dR(new real[s]), dZ(new real[s]), energy(0)
43  {}
44  ~KsData() {
45  delete []dR;
46  delete []dZ;
47  }
48 };
49 
50 
55 template<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 
98 template<typename Matrix, typename LDADistributor>
99 void
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 
125 template<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 
180 template<typename Matrix, typename GGADistributor>
181 void
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 
219 template<typename Matrix>
220 struct UksData {
221  Matrix *exca, *excb;
223  real* dZa, *dZb, *dZab;
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 
237 template<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 
247 template<typename Matrix>
248 void
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 
320 template<typename Matrix, typename LDADistributor>
321 void
323  int bllen, int blstart, int blend,
324  UksData<Matrix>* d)
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 
351 template<typename Matrix, typename GGADistributor>
352 static void
354  int bllen, int blstart, int blend,
355  UksData<Matrix>* d)
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 }
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
DftIntegratorBl_::bas_bl_cnt
int bas_bl_cnt[8]
Definition: integrator.h:58
BASBLOCK
#define BASBLOCK(grid, isym)
Definition: integrator.h:61
FunDensProp_
Definition: functionals.h:375
FirstDrv::fZ
real fZ
Definition: dft_common.h:62
UksData::exca
Matrix * exca
Definition: xc_evaluators.h:221
XCDistributorLda
distributes a LDA-type xc potential over the XC-matrix elements, with optimization for a closed shell...
Definition: xc_evaluators.h:56
xcCallbackLdaU
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
DftIntegratorBl_::grad
real(* grad)[3]
Definition: integrator.h:82
restrict
#define restrict
Definition: config.h:283
xcCallbackGgaU
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
Functional_::first
FirstOrderFun first
Definition: functionals.h:410
KsData
structure describing the data needed by distributors.
Definition: xc_evaluators.h:36
DftIntegratorBl_::ho
struct DftIntegratorBl_::@0::@2 ho
FunDensProp_::rhoa
real rhoa
Definition: functionals.h:376
dftpot0_
EXTERN_C void dftpot0_(FirstDrv *ds, const real *weight, const FunDensProp *dp)
Definition: dft_common.cc:766
KsData::dR
real * dR
Definition: xc_evaluators.h:38
UksData
Definition: xc_evaluators.h:220
UksData::energy
real energy
Definition: xc_evaluators.h:224
FunFirstFuncDrv::df00001
real df00001
Definition: functionals.h:124
Functional_::func
EnergyFunc func
Definition: functionals.h:409
drv1_clear
void drv1_clear(FunFirstFuncDrv *gga)
Definition: functionals.c:131
real
ergo_real real
Definition: test.cc:46
UksData::excb
Matrix * excb
Definition: xc_evaluators.h:221
XCDistributorGga
distributes a GGA-type xc potential over the XC-matrix elements.
Definition: xc_evaluators.h:126
FirstDrv
A vector of first order derivatives with respect to two parameters: density rho and SQUARE of the gra...
Definition: dft_common.h:60
DftIntegratorBl_
Definition: integrator.h:49
DftIntegratorBl_::curr_point
int curr_point
Definition: integrator.h:89
DftIntegratorBl_::weight
real * weight
Definition: integrator.h:52
UksData::dRb
real * dRb
Definition: xc_evaluators.h:222
KsData::energy
real energy
Definition: xc_evaluators.h:40
DftIntegratorBl_::g
union DftIntegratorBl_::@1 g
UksData::UksData
UksData(Matrix *a_, Matrix *b_, size_t s)
Definition: xc_evaluators.h:225
DftIntegratorBl_::ntypso
int ntypso
Definition: integrator.h:63
UksData::dRa
real * dRa
Definition: xc_evaluators.h:222
FunDensProp_::rhob
real rhob
Definition: functionals.h:376
DftIntegratorBl_::r
union DftIntegratorBl_::@0 r
XCDistributorGga::distribute
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
UksData::dZb
real * dZb
Definition: xc_evaluators.h:223
UksData::~UksData
~UksData()
Definition: xc_evaluators.h:231
DftIntegratorBl_::rad
struct DftIntegratorBl_::@1::@3 rad
UksData::dZab
real * dZab
Definition: xc_evaluators.h:223
UksData::dZa
real * dZa
Definition: xc_evaluators.h:223
FunFirstFuncDrv
Definition: functionals.h:119
KsData::excmat
Matrix * excmat
Definition: xc_evaluators.h:37
KsData::~KsData
~KsData()
Definition: xc_evaluators.h:44
long_real
ergo_long_real long_real
Definition: grid_atomic.h:43
XCDistributorGgaU::distribute
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
DftIntegratorBl_::nbast
int nbast
Definition: integrator.h:72
selected_func
Functional * selected_func
Definition: functionals.c:106
DftIntegratorBl_::rho
real * rho
Definition: integrator.h:76
DftIntegratorBl_::atv
real * atv
Definition: integrator.h:53
FunDensProp_::gradab
real gradab
Definition: functionals.h:378
XCDistributorGgaU
Definition: xc_evaluators.h:238
FirstDrv::fR
real fR
Definition: dft_common.h:61
XCDistributorLda::distribute
static void distribute(DftIntegratorBl *grid, int bllen, int blstart, int blend, real *restrict tmp, real *restrict dR, Matrix &excmat)
Definition: xc_evaluators.h:57
FunDensProp_::grada
real grada
Definition: functionals.h:377
FunFirstFuncDrv::df1000
real df1000
Definition: functionals.h:120
THR
#define THR
Definition: fun-cam.c:75
FunDensProp_::gradb
real gradb
Definition: functionals.h:377
xcCallbackLdaR
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
KsData::KsData
KsData(Matrix *m_, size_t s)
Definition: xc_evaluators.h:41
DftIntegratorBl_::nsym
int nsym
Definition: integrator.h:58
FunFirstFuncDrv::df0100
real df0100
Definition: functionals.h:121
KsData::dZ
real * dZ
Definition: xc_evaluators.h:39
FunFirstFuncDrv::df0001
real df0001
Definition: functionals.h:123
xcCallbackGgaR
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
FunFirstFuncDrv::df0010
real df0010
Definition: functionals.h:122