00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00035 template<typename Matrix>
00036 struct KsData {
00037 Matrix *excmat;
00038 real* dR;
00039 real* dZ;
00040 real energy;
00041 KsData(Matrix *m_, size_t s)
00042 : excmat(m_), dR(new real[s]), dZ(new real[s]), energy(0)
00043 {}
00044 ~KsData() {
00045 delete []dR;
00046 delete []dZ;
00047 }
00048 };
00049
00050
00055 template<typename Matrix>
00056 struct XCDistributorLda {
00057 static void distribute(DftIntegratorBl *grid,
00058 int bllen, int blstart, int blend,
00059 real * restrict tmp, real *restrict dR,
00060 Matrix& excmat)
00061 {
00062 static const int isym = 0;
00063 int jbl, j, ibl, i, k;
00064 real * restrict aos = grid->atv;
00065
00066 int (*restrict blocks)[2] = BASBLOCK(grid,isym);
00067 int bl_cnt = grid->bas_bl_cnt[isym];
00068
00069 for(jbl=0; jbl<bl_cnt; jbl++) {
00070 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
00071 int joff = j*bllen;
00072 for(k=blstart; k<blend; k++)
00073 tmp[k] = aos[k+joff]*dR[k]*0.5;
00074
00075 for(ibl=0; ibl<bl_cnt; ibl++) {
00076 int top = blocks[ibl][1] < j
00077 ? blocks[ibl][1] : j;
00078 for(i=blocks[ibl][0]; i<top; i++) {
00079 real *restrict aosi = aos + i*bllen;
00080 long_real s = 0;
00081 for(k=blstart; k<blend; k++)
00082 s += aosi[k]*tmp[k];
00083 excmat.add(j,i, s);
00084 }
00085 }
00086 long_real s = 0;
00087 for(k=blstart; k<blend; k++)
00088 s += aos[k+j*bllen]*tmp[k];
00089 excmat.add(j,j, s);
00090 }
00091 }
00092 }
00093 };
00094
00098 template<typename Matrix, typename LDADistributor>
00099 void
00100 xcCallbackLdaR(DftIntegratorBl *grid, real * restrict tmp,
00101 int bllen, int blstart, int blend,
00102 KsData<Matrix>* data)
00103 {
00104 FirstDrv drvs;
00105 int k;
00106 real * restrict dR = data->dR;
00107 FunDensProp dp = { 0 };
00108
00109 assert(grid->ntypso >0);
00110 for(k=blstart; k<blend; k++) {
00111 real weight = grid->weight[grid->curr_point+k];
00112 dp.rhoa = dp. rhob = 0.5*grid->r.rho[k];
00113 data->energy += selected_func->func(&dp)*weight;
00114 dftpot0_(&drvs, &weight, &dp);
00115 dR[k] = 2*drvs.fR;
00116 }
00117 LDADistributor::distribute(grid, bllen, blstart, blend, tmp, dR,
00118 *data->excmat);
00119 }
00120
00125 template<typename Matrix>
00126 struct XCDistributorGga {
00127 static void distribute(DftIntegratorBl *grid,
00128 int bllen, int blstart, int blend,
00129 real * restrict tmp,
00130 const real *dR, const real *dZ,
00131 Matrix& excmat)
00132 {
00133 static const int isym = 0;
00134 int jbl, j, ibl, i, k;
00135 const real * aox = grid->atv+bllen*grid->nbast;
00136 const real * aoy = grid->atv+bllen*grid->nbast*2;
00137 const real * aoz = grid->atv+bllen*grid->nbast*3;
00138 const real * aos = grid->atv;
00139
00140 const int (*blocks)[2] = BASBLOCK(grid,isym);
00141 int nblocks = grid->bas_bl_cnt[isym];
00142 for(jbl=0; jbl<nblocks; jbl++)
00143 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
00144 int joff = j*bllen;
00145 for(k=blstart; k<blend; k++)
00146 tmp[k+joff] =
00147 (dR[k]* aos[k+j*bllen] +
00148 dZ[k]*(aox[k+j*bllen]*grid->g.rad.a[k][0]+
00149 aoy[k+j*bllen]*grid->g.rad.a[k][1]+
00150 aoz[k+j*bllen]*grid->g.rad.a[k][2]));
00151 }
00152
00153 for(jbl=0; jbl<nblocks; jbl++) {
00154 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
00155 const real * tmpj = tmp + j*bllen;
00156 for(ibl=0; ibl<nblocks; ibl++) {
00157 int top = blocks[ibl][1] < j
00158 ? blocks[ibl][1] : j;
00159 for(i=blocks[ibl][0]; i<top; i++) {
00160 long_real s = 0;
00161 const real * tmpi = tmp + i*bllen;
00162 for(k=blstart; k<blend; k++)
00163 s += aos[k+i*bllen]*tmpj[k] +
00164 aos[k+j*bllen]*tmpi[k];
00165 excmat.add(i, j, s);
00166 }
00167 }
00168 long_real s = 0;
00169 for(k=blstart; k<blend; k++)
00170 s += 2*aos[k+j*bllen]*tmpj[k];
00171 excmat.add(j,j, s);
00172 }
00173 }
00174 }
00175 };
00176
00180 template<typename Matrix, typename GGADistributor>
00181 void
00182 xcCallbackGgaR(DftIntegratorBl* grid, real * restrict tmp,
00183 int bllen, int blstart, int blend,
00184 KsData<Matrix>* data)
00185 {
00186 FirstDrv drvs;
00187 int k;
00188 real * restrict dR = data->dR;
00189 real * restrict dZ = data->dZ;
00190 FunDensProp dp = { 0 };
00191
00192 assert(grid->ntypso >0);
00193 for(k=blstart; k<blend; k++) {
00194 real weight = grid->weight[grid->curr_point+k];
00195 dp.grada = 0.5*template_blas_sqrt(grid->g.grad[k][0]*grid->g.grad[k][0]+
00196 grid->g.grad[k][1]*grid->g.grad[k][1]+
00197 grid->g.grad[k][2]*grid->g.grad[k][2]);
00198 dp. rhoa = dp.rhob = 0.5*grid->r.rho[k];
00199 dp.gradb = dp.grada;
00200 dp.gradab = dp.grada*dp.gradb;
00201 if(dp.rhoa>1e-14) {
00202 if(dp.grada<1e-35) dp.grada = 1e-35;
00203 data->energy += selected_func->func(&dp)*weight;
00204 dftpot0_(&drvs, &weight, &dp);
00205 dR[k] = 0.5*drvs.fR;
00206 dZ[k] = 0.5*drvs.fZ/dp.grada;
00207 }else {
00208 dR[k] = dZ[k] = 0;
00209 }
00210 }
00211
00212 GGADistributor::distribute(grid, bllen, blstart, blend, tmp, dR, dZ,
00213 *data->excmat);
00214 }
00215
00216
00217
00218
00219 template<typename Matrix>
00220 struct UksData {
00221 Matrix *exca, *excb;
00222 real* dRa, *dRb;
00223 real* dZa, *dZb, *dZab;
00224 real energy;
00225 UksData(Matrix * a_, Matrix * b_, size_t s)
00226 : exca(a_), excb(b_),
00227 dRa(new real[s]), dRb(new real[s]),
00228 dZa(new real[s]), dZb(new real[s]), dZab(new real[s]),
00229 energy(0.0)
00230 {}
00231 ~UksData() {
00232 delete []dRa; delete []dRb;
00233 delete []dZa; delete []dZb; delete []dZab;
00234 }
00235 };
00236
00237 template<typename Matrix>
00238 struct XCDistributorGgaU {
00239 static void
00240 distribute(DftIntegratorBl *grid, int bllen, int blstart, int blend,
00241 real * restrict tmp, const real * dR,
00242 const real *dZ1, const real (*grad1)[3],
00243 const real *dZ2, const real (*grad2)[3],
00244 Matrix& excmat);
00245 };
00246
00247 template<typename Matrix>
00248 void
00249 XCDistributorGgaU<Matrix>::distribute(DftIntegratorBl *grid,
00250 int bllen, int blstart, int blend,
00251 real * restrict tmp, const real * dR,
00252 const real *dZ1,
00253 const real (*grad1)[3],
00254 const real *dZ2,
00255 const real (*grad2)[3],
00256 Matrix& excmat)
00257 {
00258 int isym, jbl, j, ibl, i, k;
00259 real * restrict aox = grid->atv+bllen*grid->nbast;
00260 real * restrict aoy = grid->atv+bllen*grid->nbast*2;
00261 real * restrict aoz = grid->atv+bllen*grid->nbast*3;
00262 real * restrict aos = grid->atv;
00263
00264 for(isym=0; isym<grid->nsym; isym++) {
00265 int (*restrict blocks)[2] = BASBLOCK(grid,isym);
00266 int nblocks = grid->bas_bl_cnt[isym];
00267 for(jbl=0; jbl<nblocks; jbl++)
00268 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
00269 int joff = j*bllen;
00270 for(k=blstart; k<blend; k++)
00271 tmp[k+joff] =
00272 dR[k]* aos[k+j*bllen] +
00273 dZ1[k]*(aox[k+j*bllen]*grad1[k][0]+
00274 aoy[k+j*bllen]*grad1[k][1]+
00275 aoz[k+j*bllen]*grad1[k][2])+
00276 dZ2[k]*(aox[k+j*bllen]*grad2[k][0]+
00277 aoy[k+j*bllen]*grad2[k][1]+
00278 aoz[k+j*bllen]*grad2[k][2]);
00279 }
00280
00281 for(jbl=0; jbl<nblocks; jbl++) {
00282 for(j=blocks[jbl][0]; j<blocks[jbl][1]; j++) {
00283 real *restrict tmpj = tmp + j*bllen;
00284 for(ibl=0; ibl<nblocks; ibl++) {
00285
00286 #if defined(FULL_ONLY)
00287 for(i=blocks[ibl][0]; i<blocks[ibl][1]; i++) {
00288 long_real s = 0;
00289 for(k=blstart; k<blend; k++)
00290 s += aos[k+i*bllen]*tmpj[k];
00291 excmat.add(i, j, s);
00292 }
00293 #else
00294 int top = blocks[ibl][1] < j
00295 ? blocks[ibl][1] : j;
00296 for(i=blocks[ibl][0]; i<top; i++) {
00297 long_real s = 0;
00298 const real * tmpi = tmp + i*bllen;
00299 for(k=blstart; k<blend; k++)
00300 s += aos[k+i*bllen]*tmpj[k] +
00301 aos[k+j*bllen]*tmpi[k];
00302 excmat.add(i, j, s);
00303 }
00304 #endif
00305 }
00306 #if !defined(FULL_ONLY)
00307 long_real s = 0;
00308 for(k=blstart; k<blend; k++)
00309 s += 2*aos[k+j*bllen]*tmpj[k];
00310 excmat.add(j,j, s);
00311 #endif
00312 }
00313 }
00314 }
00315 }
00316
00320 template<typename Matrix, typename LDADistributor>
00321 void
00322 xcCallbackLdaU(DftIntegratorBl *grid, real * restrict tmp,
00323 int bllen, int blstart, int blend,
00324 UksData<Matrix>* d)
00325 {
00326 FunFirstFuncDrv drvs;
00327 FunDensProp dp = { 0 };
00328 int k;
00329
00330 assert(grid->ntypso >0);
00331 for(k=blstart; k<blend; k++) {
00332 real weight = grid->weight[grid->curr_point+k];
00333 dp.rhoa = grid->r.ho.a[k];
00334 dp.rhob = grid->r.ho.b[k];
00335 d->energy += selected_func->func(&dp)*weight;
00336 drv1_clear(&drvs);
00337 selected_func->first(&drvs, weight, &dp);
00338 d->dRa[k] = 2*drvs.df1000;
00339 d->dRb[k] = 2*drvs.df0100;
00340 }
00341
00342 LDADistributor::distribute(grid, bllen, blstart, blend,
00343 tmp, d->dRa, *d->exca);
00344 LDADistributor::distribute(grid, bllen, blstart, blend,
00345 tmp, d->dRb, *d->excb);
00346 }
00347
00351 template<typename Matrix, typename GGADistributor>
00352 static void
00353 xcCallbackGgaU(DftIntegratorBl* grid, real * restrict tmp,
00354 int bllen, int blstart, int blend,
00355 UksData<Matrix>* d)
00356 {
00357 FunFirstFuncDrv drvs;
00358 int k;
00359 FunDensProp dp = { 0 };
00360
00361 assert(grid->ntypso >0);
00362 for(k=blstart; k<blend; k++) {
00363 const real THR = 1e-40;
00364 real weight = grid->weight[grid->curr_point+k];
00365 dp.rhoa = grid->r.ho.a[k];
00366 dp.rhob = grid->r.ho.b[k];
00367 dp.grada = template_blas_sqrt(grid->g.rad.a[k][0]*grid->g.rad.a[k][0]+
00368 grid->g.rad.a[k][1]*grid->g.rad.a[k][1]+
00369 grid->g.rad.a[k][2]*grid->g.rad.a[k][2]);
00370 dp.gradb = template_blas_sqrt(grid->g.rad.b[k][0]*grid->g.rad.b[k][0]+
00371 grid->g.rad.b[k][1]*grid->g.rad.b[k][1]+
00372 grid->g.rad.b[k][2]*grid->g.rad.b[k][2]);
00373 dp.gradab = grid->g.rad.a[k][0]*grid->g.rad.b[k][0]+
00374 grid->g.rad.a[k][1]*grid->g.rad.b[k][1]+
00375 grid->g.rad.a[k][2]*grid->g.rad.b[k][2];
00376 if(dp.rhoa+dp.rhob>1e-17) {
00377 if(dp.grada<THR) dp.grada = THR;
00378 if(dp.rhob<THR) dp.rhob = THR;
00379 if(dp.gradb<THR) dp.gradb = THR;
00380 if(dp.gradab<THR) dp.gradab = THR;
00381 d->energy += selected_func->func(&dp)*weight;
00382 drv1_clear(&drvs);
00383 selected_func->first(&drvs, weight, &dp);
00384 d->dRa[k] = drvs.df1000*0.5;
00385 d->dRb[k] = drvs.df0100*0.5;
00386 d->dZa[k] = drvs.df0010/dp.grada;
00387 d->dZb[k] = drvs.df0001/dp.gradb;
00388 d->dZab[k] = drvs.df00001;
00389 } else {
00390 d->dRa[k] = d->dRb[k] = d->dZa[k] = d->dZb[k] = d->dZab[k] = 0;
00391 }
00392 }
00393
00394 GGADistributor::distribute(grid, bllen, blstart, blend, tmp, d->dRa,
00395 d->dZa, grid->g.rad.a, d->dZab,
00396 grid->g.rad.b, *d->exca);
00397 GGADistributor::distribute(grid, bllen, blstart, blend, tmp, d->dRb,
00398 d->dZb, grid->g.rad.b, d->dZab,
00399 grid->g.rad.a, *d->excb);
00400 }