1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * HF/DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2013
9 * Copyright (c) 2010-2013, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16
17 #include "stockholder.h"
18 #include "timer.h"
19 #include "dftgrid.h"
20 #include "lebedev.h"
21 #include <cfloat>
22
23 /// Maximum amount of iterations
24 #define MAXITER 10000
25
StockholderAtom()26 StockholderAtom::StockholderAtom() {
27 }
28
~StockholderAtom()29 StockholderAtom::~StockholderAtom() {
30 }
31
compute(const BasisSet & basis,const arma::mat & P,const std::vector<double> & shran,const std::vector<size_t> & compute_shells,double dr,size_t irad,int lmax)32 void StockholderAtom::compute(const BasisSet & basis, const arma::mat & P, const std::vector<double> & shran, const std::vector<size_t> & compute_shells, double dr, size_t irad, int lmax) {
33 // Current radius
34 double rad=irad*dr;
35
36 // Get rule
37 std::vector<lebedev_point_t> leb=lebedev_sphere(lmax);
38
39 // Nuclear coordinate
40 coords_t nuc=basis.get_nucleus(atind).r;
41
42 // Allocate memory
43 rho[irad].resize(leb.size());
44 weights[irad].resize(leb.size());
45 grid[irad].resize(leb.size());
46
47 #ifdef _OPENMP
48 #pragma omp parallel for
49 #endif
50 for(size_t ip=0;ip<grid[irad].size();ip++) {
51 // Store weight
52 weights[irad][ip]=leb[ip].w;
53
54 // and coordinates
55 grid[irad][ip].x=nuc.x+rad*leb[ip].x;
56 grid[irad][ip].y=nuc.y+rad*leb[ip].y;
57 grid[irad][ip].z=nuc.z+rad*leb[ip].z;
58
59 // Compute functions at grid point
60 std::vector<bf_f_t> flist;
61 for(size_t i=0;i<compute_shells.size();i++) {
62 // Shell is
63 size_t ish=compute_shells[i];
64
65 // Center of shell is
66 coords_t shell_center=basis.get_shell_center(ish);
67 // Compute distance of point to center of shell
68 double shell_dist=norm(shell_center-grid[irad][ip]);
69
70 // Add shell to point if it is within the range of the shell
71 if(shell_dist<shran[ish]) {
72 // Index of first function on shell is
73 size_t ind0=basis.get_first_ind(ish);
74
75 // Compute values of basis functions
76 arma::vec fval=basis.eval_func(ish,grid[irad][ip].x,grid[irad][ip].y,grid[irad][ip].z);
77
78 // and add them to the list
79 bf_f_t hlp;
80 for(size_t ifunc=0;ifunc<fval.n_elem;ifunc++) {
81 // Index of function is
82 hlp.ind=ind0+ifunc;
83 // Value is
84 hlp.f=fval(ifunc);
85 // Add to stack
86 flist.push_back(hlp);
87 }
88 }
89 }
90
91 // Compute density
92 double d=0.0;
93 for(size_t ii=0;ii<flist.size();ii++) {
94 // Index of function is
95 size_t i=flist[ii].ind;
96 for(size_t jj=0;jj<ii;jj++) {
97 // Index of function is
98 size_t j=flist[jj].ind;
99 d+=2.0*P(i,j)*flist[ii].f*flist[jj].f;
100 }
101 d+=P(i,i)*flist[ii].f*flist[ii].f;
102 }
103
104 // Store the density
105 rho[irad][ip]=d;
106 }
107 }
108
fill_adaptive(const BasisSet & basis,const arma::mat & P,const Hirshfeld & hirsh,size_t atindv,double dr,int nrad,int lmax,double tol,bool verbose)109 void StockholderAtom::fill_adaptive(const BasisSet & basis, const arma::mat & P, const Hirshfeld & hirsh, size_t atindv, double dr, int nrad, int lmax, double tol, bool verbose) {
110 // Allocate memory for radial shells
111 rho.resize(nrad);
112 weights.resize(nrad);
113 grid.resize(nrad);
114
115 // Store atom index
116 atind=atindv;
117
118 // Nuclear distances
119 std::vector<double> nucdist=basis.get_nuclear_distances(atind);
120 // Shell ranges
121 std::vector<double> shran=basis.get_shell_ranges();
122
123 // Add points
124 for(int irad=0;irad<nrad;irad++) {
125 // Radius is
126 double rad=irad*dr;
127
128 // Indices of shells to compute
129 std::vector<size_t> compute_shells;
130
131 // Determine which shells might contribute to this radial shell
132 for(size_t inuc=0;inuc<basis.get_Nnuc();inuc++) {
133 // Determine closest distance of nucleus
134 double dist=fabs(nucdist[inuc]-rad);
135 // Get indices of shells centered on nucleus
136 std::vector<size_t> shellinds=basis.get_shell_inds(inuc);
137
138 // Loop over shells on nucleus
139 for(size_t ish=0;ish<shellinds.size();ish++) {
140
141 // Shell is relevant if range is larger than minimal distance
142 if(dist<=shran[shellinds[ish]]) {
143 // Add shell to list of shells to compute
144 compute_shells.push_back(shellinds[ish]);
145 }
146 }
147 }
148
149 // Initial order
150 int l=3;
151 // Compute radial shell
152 compute(basis,P,shran,compute_shells,dr,irad,l);
153 // and get the average
154 double avg=average(hirsh,irad);
155
156 // Next rule and new average
157 int lnext;
158 double avgnext;
159
160 std::vector<double> oldrho, oldweights;
161 std::vector<coords_t> oldgrid;
162
163 while(l<lmax) {
164 // Get old rho, weight and grid
165 oldrho=rho[irad];
166 oldweights=weights[irad];
167 oldgrid=grid[irad];
168
169 // Increment rule
170 lnext=next_lebedev(l);
171 compute(basis,P,shran,compute_shells,dr,irad,lnext);
172 avgnext=average(hirsh,irad);
173
174 // Difference is
175 double diff=rad*rad*fabs(avg-avgnext)*dr;
176
177 // Check if convergence is fulfilled.
178 if(diff < tol/nrad)
179 break;
180 else {
181 // Switch variables
182 l=lnext;
183 avg=avgnext;
184 }
185 }
186
187 // Revert to converged grid
188 rho[irad]=oldrho;
189 weights[irad]=oldweights;
190 grid[irad]=oldgrid;
191 }
192
193 // Count size of grid
194 size_t N=0;
195 for(size_t ir=0;ir<grid.size();ir++)
196 N+=grid[ir].size();
197
198 if(verbose) {
199 printf("%4i %7i\n",(int) atind+1,(int) N);
200 fflush(stdout);
201 }
202 }
203
fill_static(const BasisSet & basis,const arma::mat & P,size_t atindv,double dr,int nrad,int l,bool verbose)204 void StockholderAtom::fill_static(const BasisSet & basis, const arma::mat & P, size_t atindv, double dr, int nrad, int l, bool verbose) {
205 // Allocate memory for radial shells
206 rho.resize(nrad);
207 weights.resize(nrad);
208 grid.resize(nrad);
209
210 // Store atom index
211 atind=atindv;
212
213 // Nuclear distances
214 std::vector<double> nucdist=basis.get_nuclear_distances(atind);
215 // Shell ranges
216 std::vector<double> shran=basis.get_shell_ranges();
217
218 // Add points
219 for(int irad=0;irad<nrad;irad++) {
220 // Radius is
221 double rad=irad*dr;
222
223 // Indices of shells to compute
224 std::vector<size_t> compute_shells;
225
226 // Determine which shells might contribute to this radial shell
227 for(size_t inuc=0;inuc<basis.get_Nnuc();inuc++) {
228 // Determine closest distance of nucleus
229 double dist=fabs(nucdist[inuc]-rad);
230 // Get indices of shells centered on nucleus
231 std::vector<size_t> shellinds=basis.get_shell_inds(inuc);
232
233 // Loop over shells on nucleus
234 for(size_t ish=0;ish<shellinds.size();ish++) {
235
236 // Shell is relevant if range is larger than minimal distance
237 if(dist<=shran[shellinds[ish]]) {
238 // Add shell to list of shells to compute
239 compute_shells.push_back(shellinds[ish]);
240 }
241 }
242 }
243
244 // Compute radial shell
245 compute(basis,P,shran,compute_shells,dr,irad,l);
246 }
247
248 // Count size of grid
249 size_t N=0;
250 for(size_t ir=0;ir<grid.size();ir++)
251 N+=grid[ir].size();
252
253 if(verbose) {
254 printf("%4i %7i\n",(int) atind+1,(int) N);
255 fflush(stdout);
256 }
257 }
258
average(const Hirshfeld & hirsh,size_t irad) const259 double StockholderAtom::average(const Hirshfeld & hirsh, size_t irad) const {
260 // Initialize density
261 double d=0.0;
262 // and total weight
263 double w=0.0;
264
265 // Loop over grid
266 for(size_t ip=0;ip<grid[irad].size();ip++) {
267 // Increment total angular weight
268 w+=weights[irad][ip];
269 // and spherically averaged density
270 double c=weights[irad][ip]*rho[irad][ip]*hirsh.get_weight(atind,grid[irad][ip]);
271 d+=c;
272 }
273
274 // Store spherically averaged density
275 return d/w;
276 }
277
update(const Hirshfeld & hirsh,std::vector<double> & dens)278 void StockholderAtom::update(const Hirshfeld & hirsh, std::vector<double> & dens) {
279 // Returned densities
280 dens.resize(rho.size());
281
282 // Loop over radii
283 #ifdef _OPENMP
284 #pragma omp parallel for
285 #endif
286 for(size_t irad=0;irad<grid.size();irad++) {
287 dens[irad]=average(hirsh,irad);
288 }
289 }
290
Stockholder(const BasisSet & basis,const arma::mat & P,double finaltol,double dr,int nrad,int l0,int lmax,bool verbose)291 Stockholder::Stockholder(const BasisSet & basis, const arma::mat & P, double finaltol, double dr, int nrad, int l0, int lmax, bool verbose) {
292 Timer t, ttot;
293
294 // Allocate atomic grids
295 atoms.resize(basis.get_Nnuc());
296
297 // Get centers
298 cen.resize(basis.get_Nnuc());
299 for(size_t i=0;i<basis.get_Nnuc();i++)
300 cen[i]=basis.get_nuclear_coords(i);
301
302 // Initial weight.
303 std::vector<double> w0(nrad,1.0);
304 for(int ir=0;ir<nrad;ir++)
305 w0[ir]=exp(-ir*dr);
306
307 std::vector< std::vector<double> > oldrho(cen.size(),w0);
308 std::vector< std::vector<double> > newrho(cen.size(),w0);
309
310 // Update weights
311 ISA.set(cen,dr,oldrho);
312 (void) l0;
313
314 // Current tolerance
315 double tol=1e-2;
316
317 if(verbose) {
318 printf("Filling initial Stockholder molecular density grid.\n");
319 printf("%4s %7s\n","atom","Npoints");
320 fflush(stdout);
321 }
322 for(size_t i=0;i<basis.get_Nnuc();i++)
323 atoms[i].fill_static(basis,P,i,dr,nrad,l0,verbose);
324 if(verbose) {
325 printf("Initial fill done in %s.\n",t.elapsed().c_str());
326 fflush(stdout);
327 t.set();
328 }
329
330 // Evaluate new spherically averaged densities
331 for(size_t i=0;i<atoms.size();i++)
332 atoms[i].update(ISA,oldrho[i]);
333
334 if(verbose) {
335 printf("Spherically averaged densities updated in %s.\n",t.elapsed().c_str());
336 fflush(stdout);
337 t.set();
338 }
339
340
341 /*
342 {
343 // Print out densities
344 std::ostringstream fname;
345 fname << "atoms_start.dat";
346 FILE *out=fopen(fname.str().c_str(),"w");
347 for(size_t irad=0;irad<oldrho[0].size();irad++) {
348 fprintf(out,"%e",irad*dr);
349 for(size_t iat=0;iat<oldrho.size();iat++)
350 fprintf(out," %e",oldrho[iat][irad]);
351 fprintf(out,"\n");
352 }
353 fclose(out);
354 }
355 */
356
357 // and use these as starting weights for the self-consistent iteration
358 ISA.set(cen,dr,oldrho);
359
360 while(true) {
361 // Compute molecular density
362 if(verbose) {
363 printf("\nFilling Stockholder molecular density grid.\n");
364 printf("%4s %7s\n","atom","Npoints");
365 fflush(stdout);
366 }
367
368 // Adaptive generation of grid
369 for(size_t i=0;i<basis.get_Nnuc();i++)
370 atoms[i].fill_adaptive(basis,P,ISA,i,dr,nrad,lmax,tol,verbose);
371
372 if(verbose) {
373 printf("Grid filled in %s. Grid iteration\n",t.elapsed().c_str());
374 printf("%5s %12s %12s\n","iter","max","mean");
375 fflush(stdout);
376 t.set();
377 }
378
379 size_t iiter;
380 for(iiter=0;iiter<MAXITER;iiter++) {
381 // Evaluate new spherically averaged densities
382 for(size_t i=0;i<atoms.size();i++)
383 atoms[i].update(ISA,newrho[i]);
384
385 // Check for convergence
386 arma::vec diff(atoms.size());
387 diff.zeros();
388
389 for(size_t irad=0;irad<newrho[0].size();irad++) {
390 // Radius
391 double r=irad*dr;
392
393 // Compute density difference
394 for(size_t iat=0;iat<atoms.size();iat++) {
395 double d=fabs(newrho[iat][irad]-oldrho[iat][irad]);
396 diff(iat)+=r*r*d*dr;
397 }
398 }
399
400 // Swap densities
401 std::swap(newrho,oldrho);
402
403 double maxdiff=arma::max(diff);
404 double meandiff=arma::mean(diff);
405
406 if(verbose) {
407 printf("%5i %e %e\n",(int) iiter+1,maxdiff,meandiff);
408 fflush(stdout);
409 }
410
411 if(maxdiff<tol)
412 break;
413
414 // Update weights
415 ISA.set(cen,dr,oldrho);
416 }
417
418 if(iiter==MAXITER)
419 throw std::runtime_error("Stockholder analysis did not converge!\n");
420 else if(verbose) {
421 printf("Iteration converged within %e in %s.\n",tol,t.elapsed().c_str());
422 fflush(stdout);
423
424 /*
425 // Print out densities
426 char fname[80];
427 sprintf(fname,"atoms_%.3e.dat",tol);
428 FILE *out=fopen(fname,"w");
429 for(size_t irad=0;irad<newrho[0].size();irad++) {
430 fprintf(out,"%e",irad*dr);
431 for(size_t iat=0;iat<oldrho.size();iat++)
432 fprintf(out," %e",oldrho[iat][irad]);
433 fprintf(out,"\n");
434 }
435 fclose(out);
436 */
437 }
438
439 // Reduce tolerance and check convergence
440 tol/=sqrt(10.0);
441 if(tol < (1-sqrt(DBL_EPSILON))*finaltol)
442 break;
443 }
444
445 if(verbose) {
446 printf("Stockholder atoms solved in %s.\n",ttot.elapsed().c_str());
447 fflush(stdout);
448 }
449 }
450
~Stockholder()451 Stockholder::~Stockholder() {
452 }
453
get() const454 Hirshfeld Stockholder::get() const {
455 return ISA;
456 }
457