1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // Species.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include "Species.h"
20 #include "spline.h"
21 #include "sinft.h"
22 #include <cmath>
23 #include <cassert>
24 #include <string>
25 #include <iostream>
26 #include <iomanip>
27 #include <fstream>
28 #include <vector>
29 #include <algorithm>
30
31 using namespace std;
32
33 ////////////////////////////////////////////////////////////////////////////////
simpsn(int n,double * t)34 static double simpsn ( int n, double *t )
35 {
36 const double c0 = 17.0/48.0, c1 = 59.0/48.0,
37 c2 = 43.0/48.0, c3 = 49.0/48.0;
38 double sum = c0 * ( t[0] + t[n-1] ) +
39 c1 * ( t[1] + t[n-2] ) +
40 c2 * ( t[2] + t[n-3] ) +
41 c3 * ( t[3] + t[n-4] );
42
43 for ( int i = 4; i < n-4; i++ )
44 {
45 sum += t[i];
46 }
47 return sum;
48 }
49
50 ////////////////////////////////////////////////////////////////////////////////
Species(string name)51 Species::Species(string name) : name_(name),
52 zval_(-1), mass_(0.0), lmax_(-1), deltar_(0.0), atomic_number_(0),
53 llocal_(-1), nquad_(-1), rquad_(0.0),
54 rcps_(0.0), uri_(""), description_("undefined"), symbol_("Undef")
55 {}
56
57 ////////////////////////////////////////////////////////////////////////////////
initialize(double rcpsval)58 bool Species::initialize(double rcpsval)
59 {
60 rcps_ = rcpsval;
61 // select appropriate initialization method
62 if ( type_ == NCPP )
63 return initialize_ncpp();
64 else if ( type_ == SLPP )
65 return initialize_slpp();
66 else
67 throw SpeciesInitException("potential type unknown");
68 }
69
70 ////////////////////////////////////////////////////////////////////////////////
initialize_ncpp()71 bool Species::initialize_ncpp()
72 {
73 assert(description_ != "undefined");
74
75 const double fpi = 4.0 * M_PI;
76
77 const int np = vps_[0].size();
78 if (zval_ < 0) throw SpeciesInitException("zval_ < 0");
79 if (rcps_ < 0.0) throw SpeciesInitException("rcps_ < 0");
80 if (mass_ < 0.0) throw SpeciesInitException("mass_ < 0");
81 if (lmax_ < 0 || lmax_ > 3) throw SpeciesInitException("lmax_<0 or lmax_>3");
82
83 if (vps_.size() < lmax_+1) throw SpeciesInitException("vps_.size < lmax_+1");
84
85 if (llocal_ < 0 || llocal_ > lmax_)
86 throw SpeciesInitException("llocal_ < 0 || llocal_ > lmax_");
87
88 if ( nquad_ == 0 ) // KB potential
89 {
90 for ( int l = 0; l <= lmax_; l++ )
91 if ( l != llocal_ && phi_[l].size() == 0 )
92 throw SpeciesInitException("phi_[l] undefined for non-local projector");
93 }
94
95 if ( nquad_ < 0 ) throw SpeciesInitException("nquad < 0");
96 if ( nquad_ > 0 && rquad_ <= 0.0 )
97 throw SpeciesInitException("semilocal with rquad_ <= 0");
98
99 // compute number of non-local projectors nlm_
100 nlm_ = 0;
101 nop_ = lmax_ + 1;
102 for ( int l = 0; l <= lmax_; l++ )
103 {
104 if ( l != llocal_ )
105 {
106 nlm_ += 2 * l + 1;
107 }
108 lmap_.push_back(l);
109 }
110
111 // compute ndft_: size of radial FFT array
112 // ndft_ is a power of 2 larger than ( rdftmin / deltar_ )
113 // minimum outer bound in (a.u.)
114 const double rdftmin = 40.0;
115 assert(deltar_ > 0.0);
116 ndft_ = 1;
117 while ( ndft_ * deltar_ < rdftmin )
118 ndft_ *= 2;
119
120 rps_spl_.resize(ndft_);
121 for ( int i = 0; i < ndft_; i++ )
122 rps_spl_[i] = i * deltar_;
123
124 vps_spl_.resize(lmax_+1);
125 vps_spl2_.resize(lmax_+1);
126 phi_spl_.resize(lmax_+1);
127 phi_spl2_.resize(lmax_+1);
128
129 for ( int l = 0; l <= lmax_; l++ )
130 {
131 vps_spl_[l].resize(ndft_);
132 vps_spl2_[l].resize(ndft_);
133 if ( phi_[l].size() > 0 )
134 {
135 phi_spl_[l].resize(ndft_);
136 phi_spl2_[l].resize(ndft_);
137 }
138 }
139
140 // extend rps and vps_ to full mesh (up to i==ndft_-1)
141
142 vector<double> fint(ndft_);
143
144 wsg_.resize(lmax_+1);
145 gspl_.resize(ndft_);
146 vlocg_spl_.resize(ndft_);
147 vlocg_spl2_.resize(ndft_);
148
149 vnlg_spl_.resize(lmax_+1);
150 vnlg_spl2_.resize(lmax_+1);
151
152 vector<double> vlocr(ndft_);
153 vector<vector<double> > vnlr(lmax_+1);
154
155 for ( int l = 0; l <= lmax_; l++ )
156 {
157 vnlr[l].resize(ndft_);
158 vnlg_spl_[l].resize(ndft_+1);
159 vnlg_spl2_[l].resize(ndft_+1);
160 }
161
162 // Extend vps_[l][i] up to ndft_ using -zv/r
163
164 for ( int l = 0; l <= lmax_; l++ )
165 {
166 for ( int i = 0; i < np; i++ )
167 vps_spl_[l][i] = vps_[l][i];
168 for ( int i = np; i < ndft_; i++ )
169 vps_spl_[l][i] = - zval_ / rps_spl_[i];
170 if ( phi_[l].size() > 0 )
171 {
172 for ( int i = 0; i < np; i++ )
173 phi_spl_[l][i] = phi_[l][i];
174 for ( int i = np; i < ndft_; i++ )
175 phi_spl_[l][i] = 0.0;
176 }
177 }
178
179 // compute spline coefficients of vps_spl_ and phi_spl_
180 for ( int l = 0; l <= lmax_; l++ )
181 {
182 const double dvdr = zval_ / (rps_spl_[ndft_-1]*rps_spl_[ndft_-1]);
183 spline(ndft_,&rps_spl_[0],&vps_spl_[l][0],0.0,dvdr,0,0,&vps_spl2_[l][0]);
184 }
185 for ( int l = 0; l <= lmax_; l++ )
186 {
187 if ( l != llocal_ && phi_[l].size() > 0 )
188 {
189 spline(ndft_,&rps_spl_[0],&phi_spl_[l][0],0.0,0.0,0,1,&phi_spl2_[l][0]);
190 }
191 }
192
193 // nonlinear core correction
194 if ( has_nlcc() ) initialize_nlcc();
195 ztot_ = zval_;
196
197 // local potential: subtract the long range part due to the smeared charge
198 // Next line: constant is 2/sqrt(pi)
199 // math.h: # define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */
200 vlocr[0] = vps_spl_[llocal_][0] + (ztot_/rcps_) * M_2_SQRTPI;
201 for ( int i = 1; i < ndft_; i++ )
202 {
203 vlocr[i] = vps_spl_[llocal_][i] + (ztot_/rps_spl_[i]) *
204 erf( rps_spl_[i]/rcps_ );
205 }
206
207 // Prepare the function vlocr to be used later in the Bessel transforms:
208 //
209 // local potential: v(G) = 4 pi \int r^2 vloc(r) sin(Gr)/(Gr) dr
210 // -> store 4 pi r dr vps_(lmax_) in vlocr(i,is)
211 //
212 // the Bessel transform is then:
213 // v(G) = 1/G \sum_r sin(Gr) vlocr
214
215 for ( int i = 0; i < ndft_; i++ )
216 {
217 vlocr[i] *= fpi * rps_spl_[i] * deltar_;
218 }
219 // Local potential
220 // Compute Fourier coefficients of the local potential
221 // vlocr[i] contains 4 pi r dr vpsr(lmax_)
222 // v(G) = 4 pi \int r^2 vpsr(r) sin(Gr)/(Gr) dr
223 // = 1/G \sum_r sin(Gr) vlocr
224 //
225 // v(G=0) by simpson integration
226 // v(G) = 4 pi \int r^2 vpsr(r) dr
227 // = \sum_r r vlocr
228 //
229 // N.B. vlocr[i] contains 4 pi r dr (vps_(lmax_)-v_pseudocharge(r))
230 // Simpson integration up to ndft_ (V is zero beyond that point)
231 // Use fint as temporary array for integration
232
233 for ( int i = 0; i < ndft_; i++ )
234 {
235 fint[i] = vlocr[i] * rps_spl_[i];
236 }
237
238 const double v0 = simpsn(ndft_,&fint[0]);
239 // next line: include correction for v(g=0) to be independent of rcps
240 // const double v0 = simpsn(ndft_,&fint[0]) - M_PI*zval_*rcps_*rcps_;
241
242 for ( int i = 0; i < ndft_; i++ )
243 {
244 vlocg_spl_[i] = vlocr[i];
245 }
246
247 sinft(ndft_,&vlocg_spl_[0]);
248 if ( has_nlcc() ) sinft(ndft_,&nlccg_spl_[0]);
249
250 // Divide by g's
251 gspl_[0] = 0.0;
252 vlocg_spl_[0] = v0;
253 if ( has_nlcc() ) nlccg_spl_[0] = zcore_;
254 double fac = M_PI/(ndft_*deltar_);
255 for ( int i = 1; i < ndft_; i++ )
256 {
257 gspl_[i] = i * fac;
258 vlocg_spl_[i] /= gspl_[i];
259 if ( has_nlcc() ) nlccg_spl_[i] /= gspl_[i];
260 }
261
262 // Initialize cubic spline interpolation for local potential Vloc(G)
263 // Use zero first derivative at G=0 and natural (y"=0) at Gmax
264
265 spline(ndft_,&gspl_[0],&vlocg_spl_[0],0.0,0.0,0,1,&vlocg_spl2_[0]);
266 if ( has_nlcc() )
267 spline(ndft_,&gspl_[0],&nlccg_spl_[0],0.0,0.0,0,1,&nlccg_spl2_[0]);
268
269
270 // Non-local KB projectors
271 if ( nquad_ == 0 && lmax_ > 0)
272 {
273 for ( int l = 0; l <= lmax_; l++ )
274 {
275 wsg_[l] = 0.0;
276
277 if ( l != llocal_ )
278 {
279 // for KB potentials, compute weights wsg[l]
280 // Compute weight wsg_[l] by integration on the linear mesh
281 for ( int i = 0; i < ndft_; i++ )
282 {
283 double tmp = phi_spl_[l][i] * rps_spl_[i];
284 fint[i] = ( vps_spl_[l][i] - vps_spl_[llocal_][i] ) * tmp * tmp;
285 }
286 double tmp = simpsn(ndft_,&fint[0]);
287 assert(tmp != 0.0);
288 // Next lines: store 1/<phi|delta_v| phi > in wsg[is][l]
289 wsg_[l] = 1.0 / ( deltar_ * tmp );
290
291 // compute non-local projectors:
292 // w(G) = Ylm(G) i^l 4 pi \int r^2 phi_l(r) j_l(Gr) v_l(r) dr
293 // -> store 4 pi v_l(r) phi_l(r) dr in vnlr[l][i]
294 // the Bessel transform is then done by
295 // l=0: j_0(Gr) = sin(Gr)/(Gr)
296 // w_0(G) = Ylm(G) 1/G \sum_r sin(Gr) r vnlr
297 // l=1: j_1(Gr) = sin(Gr)/(Gr)^2 - cos(Gr)/Gr
298 // w_1(G) = Ylm(G) i^-1 1/G^2 \sum_r sin(Gr) vnlr -
299 // Ylm(G) i^-1 1/G \sum_r cos(Gr) r vnlr
300 // l=2: j_2(Gr) = sin(Gr)*(3/(Gr)^3 -1/(Gr)) - 3*cos(Gr)/(Gr)^2
301 // w_2(G) = Ylm(G) i^-2 ( 3/G^3 \sum_r sin(Gr)/r vnlr -
302 // 1/G \sum_r sin(Gr)*r vnlr -
303 // 3/G^2 \sum_r cos(Gr) vnlr )
304
305 for ( int i = 0; i < ndft_; i++ )
306 {
307 vnlr[l][i] = fpi * deltar_ *
308 ( vps_spl_[l][i] - vps_spl_[llocal_][i] ) * phi_spl_[l][i];
309 }
310
311 }
312 }
313
314 // compute radial Fourier transforms of vnlr
315
316 // Next line: vnlg_spl_ dimensioned ndft_+1 since it is passed to cosft1
317
318 // Non-local potentials
319 //
320 // w(G) = Ylm(G) i^l 4 pi \int r^2 phi_l(r) j_l(Gr) v_l(r) dr
321 // -> store 4 pi v_l(r) phi_l(r) dr in vnlr[l][i]
322 // the Bessel transform is then done by
323 // s: j_0(Gr) = sin(Gr)/(Gr)
324 // w_0(G) = Ylm(G) 1/G \sum_r sin(Gr) r vnlr
325 // p: j_1(Gr) = sin(Gr)/(Gr)^2 - cos(Gr)/Gr
326 // w_1(G) = Ylm(G) i^-1 1/G^2 \sum_r sin(Gr) vnlr -
327 // Ylm(G) i^-1 1/G \sum_r cos(Gr) r vnlr
328 // d: j_2(Gr) = sin(Gr)*(3/(Gr)^3 -1/(Gr)) - 3*cos(Gr)/(Gr)^2
329 // w_2(G) = Ylm(G) i^-2 ( 3/G^3 \sum_r sin(Gr)/r vnlr -
330 // 1/G \sum_r sin(Gr)*r vnlr -
331 // 3/G^2 \sum_r cos(Gr) vnlr )
332 // vnlr[l][i] contains 4 pi dr phi_l(r) v_l(r)
333
334 for ( int l = 0; l <= lmax_; l++ )
335 {
336 if ( l != llocal_ )
337 {
338 if ( l == 0 )
339 {
340
341 // s projector
342 // w_0(G) = Ylm(G) 1/G \sum_r sin(Gr) r vnlr
343 //
344 // G=0: Simpson integration up to ndft_
345 // Use fint as temporary array for integration
346
347 for ( int i = 0; i < ndft_; i++ )
348 {
349 fint[i] = vnlr[l][i] * rps_spl_[i] * rps_spl_[i];
350 }
351 const double v0 = simpsn(ndft_, &fint[0]);
352
353 for ( int i = 0; i < ndft_; i++ )
354 {
355 vnlg_spl_[l][i] = vnlr[l][i] * rps_spl_[i];
356 }
357
358 sinft(ndft_,&vnlg_spl_[l][0]);
359
360 vnlg_spl_[l][0] = v0;
361 // Divide by g
362 for ( int i = 1; i < ndft_; i++ )
363 {
364 vnlg_spl_[l][i] /= gspl_[i];
365 }
366
367 // Initialize cubic spline interpolation
368 // Use zero first derivative at G=0 and natural (y"=0) at Gmax
369
370 spline(ndft_,&gspl_[0],&vnlg_spl_[l][0],0.0,0.0,0,1,
371 &vnlg_spl2_[l][0]);
372
373 }
374 else if ( l == 1 )
375 {
376 // p projectors
377 // w_1(G) = Ylm(G) i 1/G^2 \sum_r sin(Gr) vnlr -
378 // Ylm(G) i 1/G \sum_r cos(Gr) r vnlr
379 // vnlr(i,is,l) contains 4 pi dr phi_l(r) v_l(r)
380
381 // Next line: v(G=0) is zero (j_1(Gr) -> 0 as G->0 )
382 const double v0 = 0.0;
383
384 // First part: 1/G^2 \sum_r sin(Gr) vnlr
385 for ( int i = 0; i < ndft_; i++ )
386 {
387 vnlg_spl_[l][i] = vnlr[l][i];
388 }
389
390 sinft(ndft_,&vnlg_spl_[l][0]);
391
392 // Divide by g**2 and store in fint */
393 fint[0] = v0;
394 for ( int i = 1; i < ndft_; i++ )
395 {
396 fint[i] = vnlg_spl_[l][i] / ( gspl_[i] * gspl_[i] );
397 }
398
399 // Second part: cosine transform: 1/G \sum_r cos(Gr) r vnlr
400
401 for ( int i = 0; i < ndft_; i++ )
402 {
403 vnlg_spl_[l][i] = vnlr[l][i] * rps_spl_[i];
404 }
405
406 // N.B. Next line: Initialize also vnlg_[l][ndft_] to zero
407 // since it is used and modified by cosft1
408 // vnlg_ was dimensioned ndft_[is]+1
409
410 vnlg_spl_[l][ndft_] = 0.0;
411 cosft1(ndft_,&vnlg_spl_[l][0]);
412
413 // Divide by g and add to fint to get vnlg_
414 vnlg_spl_[l][0] = v0;
415 for ( int i = 1; i < ndft_; i++ )
416 {
417 vnlg_spl_[l][i] = fint[i] - vnlg_spl_[l][i]/gspl_[i];
418 }
419
420 // Initialize spline interpolation
421 // Use natural bc (y"=0) at G=0 and natural bc (y"=0) at Gmax
422
423 spline(ndft_,&gspl_[0],&vnlg_spl_[l][0],0.0,0.0,1,1,
424 &vnlg_spl2_[l][0]);
425 }
426 else if ( l == 2 )
427 {
428 // d projectors
429 // d: j_2(Gr) = sin(Gr)*(3/(Gr)^3 -1/(Gr)) - 3*cos(Gr)/(Gr)^2
430 // w_2(G) = Ylm(G) i^-2 ( 3/G^3 \sum_r sin(Gr)/r vnlr -
431 // 1/G \sum_r sin(Gr)*r vnlr -
432 // 3/G^2 \sum_r cos(Gr) vnlr )
433 // vnlr[l][i] contains 4 pi dr phi_l(r) v_l(r)
434
435 // Next line: v(G=0) is zero (j_2(Gr) -> 0 as G->0 )
436 const double v0 = 0.0;
437
438 // First part: sine transform 3/G^3 \sum_r sin(Gr)/r vnlr
439 // Note: the integrand is linear near r=0 since vnlr(r) ~ r^2
440 vnlg_spl_[l][0] = 0.0;
441 for ( int i = 1; i < ndft_; i++ )
442 {
443 vnlg_spl_[l][i] = vnlr[l][i] / rps_spl_[i];
444 }
445
446 sinft(ndft_,&vnlg_spl_[l][0]);
447
448 // multiply by 3/G^3 and store in fint */
449 fint[0] = v0;
450 for ( int i = 1; i < ndft_; i++ )
451 {
452 fint[i] = 3.0 * vnlg_spl_[l][i] / ( gspl_[i]*gspl_[i]*gspl_[i] );
453 }
454
455 // Second part: sine transform -1/G \sum_r sin(Gr)*r vnlr
456 for ( int i = 0; i < ndft_; i++ )
457 {
458 vnlg_spl_[l][i] = vnlr[l][i] * rps_spl_[i];
459 }
460
461 sinft(ndft_,&vnlg_spl_[l][0]);
462
463 // multiply by -1/G and accumulate in fint */
464 fint[0] += v0;
465 for ( int i = 1; i < ndft_; i++ )
466 {
467 fint[i] += - vnlg_spl_[l][i] / gspl_[i];
468 }
469
470 // Third part: cosine transform: -3/G^2 \sum_r cos(Gr) vnlr
471
472 for ( int i = 0; i < ndft_; i++ )
473 {
474 vnlg_spl_[l][i] = vnlr[l][i];
475 }
476
477 // N.B. Next line: Initialize also vnlg_[l][ndft_] to zero
478 // since it is used and modified by cosft1
479 // vnlg_ was dimensioned ndft_[is]+1
480
481 vnlg_spl_[l][ndft_] = 0.0;
482 cosft1(ndft_,&vnlg_spl_[l][0]);
483
484 // Multiply by -3/G^2 and add to fint
485 fint[0] += v0;
486 for ( int i = 1; i < ndft_; i++ )
487 {
488 fint[i] += - 3.0 * vnlg_spl_[l][i] / (gspl_[i] * gspl_[i]);
489 }
490
491 vnlg_spl_[l][0] = v0;
492 for ( int i = 1; i < ndft_; i++ )
493 {
494 vnlg_spl_[l][i] = fint[i];
495 }
496
497 // Initialize spline interpolation
498 // Use zero first derivative at G=0 and natural (y"=0) at Gmax
499
500 spline(ndft_,&gspl_[0],&vnlg_spl_[l][0],0.0,0.0,0,1,
501 &vnlg_spl2_[l][0]);
502 }
503 } // l != llocal_
504 } // l
505 } // nquad_ == 0 && lmax_ > 0 (KB projectors)
506 return true;
507 }
508
509 ////////////////////////////////////////////////////////////////////////////////
initialize_slpp()510 bool Species::initialize_slpp()
511 {
512 assert(description_ != "undefined");
513
514 const double fpi = 4.0 * M_PI;
515
516 const int np = vlocal_.size();
517 nquad_ = 0;
518
519 // number of projectors
520 nop_ = 0;
521 for ( int l = 0; l < proj_.size(); l++ )
522 for ( int n = 0; n < proj_[l].size(); n++ )
523 {
524 lmap_.push_back(l);
525 nop_++;
526 }
527
528 // check sanity of input
529 if (zval_ < 0) throw SpeciesInitException("zval_ < 0");
530 if (rcps_ < 0.0) throw SpeciesInitException("rcps_ < 0");
531 if (mass_ < 0.0) throw SpeciesInitException("mass_ < 0");
532
533 // check consistency of input
534 // a) for each angular momentum l there has to be a D matrix of
535 // the size n x n, where n are the number of projectors with
536 // this angular momentum l
537 if ( proj_.size() != d_.size() )
538 throw SpeciesInitException("projector and dmatrix inconsistent");
539 for ( int l = 0; l < proj_.size(); l++ )
540 {
541 // trivial case only one projector (n = 1)
542 // if D matrix is not present, create a trivial one D = 1
543 // otherwise leave D matrix unchanged
544 if ( proj_[l].size() == 1 )
545 {
546 if ( d_[l].size() > 1 )
547 throw SpeciesInitException("dmatrix has wrong dimension");
548 d_[l].resize(1);
549 if ( d_[l][0].size() > 1)
550 throw SpeciesInitException("dmatrix has wrong dimension");
551 d_[l][0].resize(1,1);
552 }
553 // general case dim( d[l] ) = n
554 else
555 {
556 const int nmax = proj_[l].size();
557 if ( d_[l].size() != nmax )
558 throw SpeciesInitException("dmatrix has wrong dimension");
559 for ( int n = 0; n < nmax; n++ )
560 {
561 if ( d_[l][n].size() != nmax )
562 throw SpeciesInitException("dmatrix has wrong dimension");
563 }
564 }
565 }
566
567 // compute number of non-local projectors nlm_
568 nlm_ = 0;
569 for ( int l = 0; l < proj_.size(); l++ )
570 {
571 nlm_ += (2 * l + 1) * proj_[l].size();
572 }
573
574 // compute ndft_: size of radial FFT array
575 // ndft_ is a power of 2 larger than ( rdftmin / deltar_ )
576 // minimum outer bound in (a.u.)
577 const double rdftmin = 40.0;
578 // next line: limit small mesh sizes
579 assert(deltar_ > 0.0001);
580 ndft_ = 1;
581 while ( ndft_ * deltar_ < rdftmin )
582 ndft_ *= 2;
583
584 rps_spl_.resize(ndft_);
585 for ( int i = 0; i < ndft_; i++ )
586 rps_spl_[i] = i * deltar_;
587
588 phi_spl_.resize(nop_);
589 phi_spl2_.resize(nop_);
590
591 for ( int iop = 0; iop < nop_; iop++ )
592 {
593 phi_spl_[iop].resize(ndft_);
594 phi_spl2_[iop].resize(ndft_);
595 }
596
597
598 // extend rps and vps_ to full mesh (up to i==ndft_-1)
599
600 vector<double> fint(ndft_);
601
602 gspl_.resize(ndft_);
603 vlocg_spl_.resize(ndft_);
604 vlocg_spl2_.resize(ndft_);
605
606 vnlg_spl_.resize(nop_);
607 vnlg_spl2_.resize(nop_);
608
609 vector<double> vlocr(ndft_);
610 vector<vector<double> > vnlr(nop_);
611
612 for ( int iop = 0; iop < nop_; iop++ )
613 {
614 vnlr[iop].resize(ndft_);
615 vnlg_spl_[iop].resize(ndft_ + 1);
616 vnlg_spl2_[iop].resize(ndft_ + 1);
617 }
618
619 // Extend vps_[l][i] up to ndft_ using -zv/r
620
621 // local potential
622 for ( int i = 0; i < np; i++ )
623 vlocr[i] = vlocal_[i];
624 for ( int i = np; i < ndft_; i++ )
625 vlocr[i] = -zval_ / rps_spl_[i];
626
627 for ( int iop = 0; iop < nop_; iop++ )
628 {
629 // counter for projectors of same l
630 int n = 0;
631 if ( iop == 0 || lmap_[iop] != lmap_[iop - 1] ) n = 0;
632 else n++;
633 for ( int i = 0; i < np; i++ )
634 phi_spl_[iop][i] = proj_[lmap_[iop]][n][i];
635 for ( int i = np; i < ndft_; i++ )
636 phi_spl_[iop][i] = 0.0;
637 }
638
639 for ( int iop = 0; iop < nop_; iop++ )
640 {
641 spline(ndft_,&rps_spl_[0],&phi_spl_[iop][0],0.0,0.0,0,1,&phi_spl2_[iop][0]);
642 }
643
644 // nonlinear core correction
645 if ( has_nlcc() ) initialize_nlcc();
646 ztot_ = zval_;
647
648 // local potential: subtract the long range part due to the smeared charge
649 // Next line: constant is 2/sqrt(pi)
650 // math.h: # define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */
651 vlocr[0] += ( ztot_ / rcps_ ) * M_2_SQRTPI;
652 for ( int i = 1; i < ndft_; i++ )
653 {
654 vlocr[i] += ( ztot_ / rps_spl_[i] ) * erf(rps_spl_[i] / rcps_);
655 }
656
657 // Prepare the function vlocr to be used later in the Bessel transforms:
658 //
659 // local potential: v(G) = 4 pi \int r^2 vloc(r) sin(Gr)/(Gr) dr
660 // -> store 4 pi r dr vps_(lmax_) in vlocr(i,is)
661 //
662 // the Bessel transform is then:
663 // v(G) = 1/G \sum_r sin(Gr) vlocr
664
665 for ( int i = 0; i < ndft_; i++ )
666 {
667 vlocr[i] *= fpi * rps_spl_[i] * deltar_;
668 }
669 // Local potential
670 // Compute Fourier coefficients of the local potential
671 // vlocr[i] contains 4 pi r dr vpsr(lmax_)
672 // v(G) = 4 pi \int r^2 vpsr(r) sin(Gr)/(Gr) dr
673 // = 1/G \sum_r sin(Gr) vlocr
674 //
675 // v(G=0) by simpson integration
676 // v(G) = 4 pi \int r^2 vpsr(r) dr
677 // = \sum_r r vlocr
678 //
679 // N.B. vlocr[i] contains 4 pi r dr (vps_(lmax_)-v_pseudocharge(r))
680 // Simpson integration up to ndft_ (V is zero beyond that point)
681 // Use fint as temporary array for integration
682
683 for ( int i = 0; i < ndft_; i++ )
684 {
685 fint[i] = vlocr[i] * rps_spl_[i];
686 }
687
688 const double v0 = simpsn(ndft_,&fint[0]);
689 // next line: include correction for v(g=0) to be independent of rcps
690 // const double v0 = simpsn(ndft_,&fint[0]) - M_PI*zval_*rcps_*rcps_;
691
692 for ( int i = 0; i < ndft_; i++ )
693 {
694 vlocg_spl_[i] = vlocr[i];
695 }
696
697 sinft(ndft_,&vlocg_spl_[0]);
698 if ( has_nlcc() ) sinft(ndft_,&nlccg_spl_[0]);
699
700 // Divide by g's
701 gspl_[0] = 0.0;
702 vlocg_spl_[0] = v0;
703 if ( has_nlcc() ) nlccg_spl_[0] = zcore_;
704 double fac = M_PI/(ndft_*deltar_);
705 for ( int i = 1; i < ndft_; i++ )
706 {
707 gspl_[i] = i * fac;
708 vlocg_spl_[i] /= gspl_[i];
709 if ( has_nlcc() ) nlccg_spl_[i] /= gspl_[i];
710 }
711
712 // Initialize cubic spline interpolation for local potential Vloc(G)
713 // Use zero first derivative at G=0 and natural (y"=0) at Gmax
714
715 spline(ndft_,&gspl_[0],&vlocg_spl_[0],0.0,0.0,0,1,&vlocg_spl2_[0]);
716 if ( has_nlcc() )
717 spline(ndft_,&gspl_[0],&nlccg_spl_[0],0.0,0.0,0,1,&nlccg_spl2_[0]);
718
719
720 // Non-local projectors
721 // inf
722 // / /
723 // | 3 i G r l | 2
724 // | d r e f (r) Y (r) = 4 pi i Y (G) | dr r j (Gr) f (r)
725 // | l lm lm | l l
726 // / /
727 // 0
728 // where f_l(r) = phi_l(r)
729
730 // multiply projectors with weight
731 // w_l = 1 / < phi_l | V_l | phi_l >
732 // as there is no V_l all weights are 1
733 wsg_.resize(nop_,1.0);
734 for ( int iop = 0; iop < nop_; iop++ )
735 {
736 // compute non-local projectors:
737 // w(G) = Ylm(G) i^l 4 pi \int r^2 j_l(Gr) f_l(r) dr
738 // -> store 4 pi f_l(r) dr in vnlr[l][i]
739 // the Bessel transform is then done by
740 // l=0: j_0(Gr) = sin(Gr)/(Gr)
741 // w_0(G) = Ylm(G) 1/G \sum_r sin(Gr) r vnlr
742 // l=1: j_1(Gr) = sin(Gr)/(Gr)^2 - cos(Gr)/Gr
743 // w_1(G) = Ylm(G) i^-1 1/G^2 \sum_r sin(Gr) vnlr -
744 // Ylm(G) i^-1 1/G \sum_r cos(Gr) r vnlr
745 // l=2: j_2(Gr) = sin(Gr)*(3/(Gr)^3 -1/(Gr)) - 3*cos(Gr)/(Gr)^2
746 // w_2(G) = Ylm(G) i^-2 ( 3/G^3 \sum_r sin(Gr)/r vnlr -
747 // 1/G \sum_r sin(Gr)*r vnlr -
748 // 3/G^2 \sum_r cos(Gr) vnlr )
749
750 for ( int i = 0; i < ndft_; i++ )
751 {
752 vnlr[iop][i] = fpi * deltar_ * phi_spl_[iop][i];
753 }
754
755 } // loop over projector
756
757 // compute radial Fourier transforms of vnlr
758
759 // Next line: vnlg_spl_ dimensioned ndft_+1 since it is passed to cosft1
760
761 // Non-local potentials
762 //
763 // w(G) = Ylm(G) i^l 4 pi \int r^2 j_l(Gr) f_l(r) dr
764 // -> store 4 pi f_l(r) dr in vnlr[l][i]
765 // the Bessel transform is then done by
766 // s: j_0(Gr) = sin(Gr)/(Gr)
767 // w_0(G) = Ylm(G) 1/G \sum_r sin(Gr) r vnlr
768 // p: j_1(Gr) = sin(Gr)/(Gr)^2 - cos(Gr)/Gr
769 // w_1(G) = Ylm(G) i^-1 1/G^2 \sum_r sin(Gr) vnlr -
770 // Ylm(G) i^-1 1/G \sum_r cos(Gr) r vnlr
771 // d: j_2(Gr) = sin(Gr)*(3/(Gr)^3 -1/(Gr)) - 3*cos(Gr)/(Gr)^2
772 // w_2(G) = Ylm(G) i^-2 ( 3/G^3 \sum_r sin(Gr)/r vnlr -
773 // 1/G \sum_r sin(Gr)*r vnlr -
774 // 3/G^2 \sum_r cos(Gr) vnlr )
775 // f: j_3(Gr) = sin(Gr)*(15/(Gr)^4-6/(Gr)^2)) +
776 // cos(Gr)*(1/(Gr)-15/(Gr)^3)
777 // w_3(G) = Ylm(G) i^-3 ( 15/G^4 \sum_r sin(Gr)/r^2 vnlr
778 // - 6/G^2 \sum_r sin(Gr) vnlr
779 // + 1/G \sum_r cos(Gr)*r vnlr
780 // - 15/G^3 \sum_r cos(Gr)/r vnlr )
781
782 for ( int iop = 0; iop < nop_; iop++ )
783 {
784
785 if ( lmap_[iop] == 0 )
786 {
787 // s projector
788 // w_0(G) = Ylm(G) 1/G \sum_r sin(Gr) r vnlr
789 //
790 // G=0: Simpson integration up to ndft_
791 // Use fint as temporary array for integration
792
793 for ( int i = 0; i < ndft_; i++ )
794 {
795 const double r = rps_spl_[i];
796 fint[i] = vnlr[iop][i] * r * r;
797 }
798 const double v0 = simpsn(ndft_,&fint[0]);
799
800 for ( int i = 0; i < ndft_; i++ )
801 {
802 vnlg_spl_[iop][i] = vnlr[iop][i] * rps_spl_[i];
803 }
804
805 sinft(ndft_,&vnlg_spl_[iop][0]);
806
807 vnlg_spl_[iop][0] = v0;
808 // Divide by g
809 for ( int i = 1; i < ndft_; i++ )
810 {
811 vnlg_spl_[iop][i] /= gspl_[i];
812 }
813
814 // Initialize cubic spline interpolation
815 // Use zero first derivative at G=0 and natural (y"=0) at Gmax
816
817 spline(ndft_,&gspl_[0],&vnlg_spl_[iop][0],
818 0.0,0.0,0,1,&vnlg_spl2_[iop][0]);
819
820 }
821 else if ( lmap_[iop] == 1 )
822 {
823 // p projectors
824 // w_1(G) = Ylm(G) i 1/G^2 \sum_r sin(Gr) vnlr -
825 // Ylm(G) i 1/G \sum_r cos(Gr) r vnlr
826
827 // Next line: v(G=0) is zero (j_1(Gr) -> 0 as G->0 )
828 const double v0 = 0.0;
829
830 // First part: 1/G^2 \sum_r sin(Gr) vnlr
831 for ( int i = 0; i < ndft_; i++ )
832 {
833 vnlg_spl_[iop][i] = vnlr[iop][i];
834 }
835
836 sinft(ndft_,&vnlg_spl_[iop][0]);
837
838 // Divide by g**2 and store in fint */
839 fint[0] = v0;
840 for ( int i = 1; i < ndft_; i++ )
841 {
842 const double g = gspl_[i];
843 fint[i] = vnlg_spl_[iop][i] / ( g * g );
844 }
845
846 // Second part: cosine transform: 1/G \sum_r cos(Gr) r vnlr
847
848 for ( int i = 0; i < ndft_; i++ )
849 {
850 vnlg_spl_[iop][i] = vnlr[iop][i] * rps_spl_[i];
851 }
852
853 // N.B. Next line: Initialize also vnlg_[l][ndft_] to zero
854 // since it is used and modified by cosft1
855 // vnlg_ was dimensioned ndft_[is]+1
856
857 vnlg_spl_[iop][ndft_] = 0.0;
858 cosft1(ndft_,&vnlg_spl_[iop][0]);
859
860 // Divide by g and add to fint to get vnlg_
861 vnlg_spl_[iop][0] = v0;
862 for ( int i = 1; i < ndft_; i++ )
863 {
864 vnlg_spl_[iop][i] = fint[i] - vnlg_spl_[iop][i] / gspl_[i];
865 }
866
867 // Initialize spline interpolation
868 // Use natural bc (y"=0) at G=0 and natural bc (y"=0) at Gmax
869
870 spline(ndft_,&gspl_[0],&vnlg_spl_[iop][0],
871 0.0,0.0,1,1,&vnlg_spl2_[iop][0]);
872 }
873 else if ( lmap_[iop] == 2 )
874 {
875 // d projectors
876 // d: j_2(Gr) = sin(Gr)*(3/(Gr)^3 -1/(Gr)) - 3*cos(Gr)/(Gr)^2
877 // w_2(G) = Ylm(G) i^-2 ( 3/G^3 \sum_r sin(Gr)/r vnlr -
878 // 1/G \sum_r sin(Gr)*r vnlr -
879 // 3/G^2 \sum_r cos(Gr) vnlr )
880
881 // Next line: v(G=0) is zero (j_2(Gr) -> 0 as G->0 )
882 const double v0 = 0.0;
883
884 // First part: sine transform 3/G^3 \sum_r sin(Gr)/r vnlr
885 // Note: the integrand is linear near r=0 since vnlr(r) ~ r^2
886 vnlg_spl_[iop][0] = 0.0;
887 for ( int i = 1; i < ndft_; i++ )
888 {
889 vnlg_spl_[iop][i] = vnlr[iop][i] / rps_spl_[i];
890 }
891
892 sinft(ndft_,&vnlg_spl_[iop][0]);
893
894 // multiply by 3/G^3 and store in fint */
895 fint[0] = v0;
896 for ( int i = 1; i < ndft_; i++ )
897 {
898 const double g = gspl_[i];
899 fint[i] = 3.0 * vnlg_spl_[iop][i] / ( g * g * g );
900 }
901
902 // Second part: sine transform -1/G \sum_r sin(Gr)*r vnlr
903 for ( int i = 0; i < ndft_; i++ )
904 {
905 vnlg_spl_[iop][i] = vnlr[iop][i] * rps_spl_[i];
906 }
907
908 sinft(ndft_,&vnlg_spl_[iop][0]);
909
910 // multiply by -1/G and accumulate in fint */
911 for ( int i = 1; i < ndft_; i++ )
912 {
913 fint[i] += -vnlg_spl_[iop][i] / gspl_[i];
914 }
915
916 // Third part: cosine transform: -3/G^2 \sum_r cos(Gr) vnlr
917 for ( int i = 0; i < ndft_; i++ )
918 {
919 vnlg_spl_[iop][i] = vnlr[iop][i];
920 }
921
922 // N.B. Next line: Initialize also vnlg_[l][ndft_] to zero
923 // since it is used and modified by cosft1
924 // vnlg_ was dimensioned ndft_[is]+1
925
926 vnlg_spl_[iop][ndft_] = 0.0;
927 cosft1(ndft_,&vnlg_spl_[iop][0]);
928
929 // Multiply by -3/G^2 and add to fint
930 for ( int i = 1; i < ndft_; i++ )
931 {
932 const double g = gspl_[i];
933 fint[i] += -3.0 * vnlg_spl_[iop][i] / ( g * g );
934 }
935
936 vnlg_spl_[iop][0] = v0;
937 for ( int i = 1; i < ndft_; i++ )
938 {
939 vnlg_spl_[iop][i] = fint[i];
940 }
941
942 // Initialize spline interpolation
943 // Use zero first derivative at G=0 and natural (y"=0) at Gmax
944
945 spline(ndft_,&gspl_[0],&vnlg_spl_[iop][0],
946 0.0,0.0,0,1,&vnlg_spl2_[iop][0]);
947 }
948 else if ( lmap_[iop] == 3 )
949 {
950 // f projectors
951 // f: j_3(Gr) = sin(Gr)*(15/(Gr)^4-6/(Gr)^2)) +
952 // cos(Gr)*(1/(Gr)-15/(Gr)^3)
953 // w_3(G) = Ylm(G) i^-3 ( 15/G^4 \sum_r sin(Gr)/r^2 vnlr
954 // - 6/G^2 \sum_r sin(Gr) vnlr
955 // + 1/G \sum_r cos(Gr)*r vnlr
956 // - 15/G^3 \sum_r cos(Gr)/r vnlr )
957
958 // Next line: v(G=0) is zero (j_3(Gr) -> 0 as G->0 )
959 const double v0 = 0.0;
960
961 // First part: sine transform 15/G^4 \sum_r sin(Gr)/r^2 vnlr
962 // Note: the integrand is linear near r=0 since vnlr(r) ~ r^3
963 vnlg_spl_[iop][0] = 0.0;
964 for ( int i = 1; i < ndft_; i++ )
965 {
966 const double r = rps_spl_[i];
967 vnlg_spl_[iop][i] = vnlr[iop][i] / (r*r);
968 }
969
970 sinft(ndft_,&vnlg_spl_[iop][0]);
971
972 // multiply by 15/G^4 and store in fint */
973 fint[0] = v0;
974 for ( int i = 1; i < ndft_; i++ )
975 {
976 const double g = gspl_[i];
977 fint[i] = 15.0 * vnlg_spl_[iop][i] / (g*g*g*g);
978 }
979
980 // Second part: sine transform -6/G^2 \sum_r sin(Gr) vnlr
981 for ( int i = 0; i < ndft_; i++ )
982 {
983 vnlg_spl_[iop][i] = vnlr[iop][i];
984 }
985
986 sinft(ndft_,&vnlg_spl_[iop][0]);
987
988 // multiply by -6/G^2 and accumulate in fint */
989 for ( int i = 1; i < ndft_; i++ )
990 {
991 const double g = gspl_[i];
992 fint[i] += -6.0 * vnlg_spl_[iop][i] / (g*g);
993 }
994
995 // Third part: cosine transform: 1/G \sum_r cos(Gr)*r vnlr
996 for ( int i = 0; i < ndft_; i++ )
997 {
998 vnlg_spl_[iop][i] = vnlr[iop][i] * rps_spl_[i];
999 }
1000
1001 // N.B. Next line: Initialize also vnlg_[l][ndft_] to zero
1002 // since it is used and modified by cosft1
1003 // vnlg_ was dimensioned ndft_[is]+1
1004
1005 vnlg_spl_[iop][ndft_] = 0.0;
1006 cosft1(ndft_,&vnlg_spl_[iop][0]);
1007
1008 // Multiply by 1/G and add to fint
1009 for ( int i = 1; i < ndft_; i++ )
1010 {
1011 fint[i] += vnlg_spl_[iop][i] / gspl_[i];
1012 }
1013
1014 // Fourth part: cosine transform: -15/G^3 \sum_r cos(Gr)/r vnlr
1015 vnlg_spl_[iop][0] = 0.0;
1016 for ( int i = 1; i < ndft_; i++ )
1017 {
1018 vnlg_spl_[iop][i] = vnlr[iop][i] / rps_spl_[i];
1019 }
1020
1021 // N.B. Next line: Initialize also vnlg_[l][ndft_] to zero
1022 // since it is used and modified by cosft1
1023 // vnlg_ was dimensioned ndft_[is]+1
1024
1025 vnlg_spl_[iop][ndft_] = 0.0;
1026 cosft1(ndft_,&vnlg_spl_[iop][0]);
1027
1028 // Multiply by -15/G^3 and add to fint
1029 for ( int i = 1; i < ndft_; i++ )
1030 {
1031 const double g = gspl_[i];
1032 fint[i] += -15.0 * vnlg_spl_[iop][i] / (g*g*g);
1033 }
1034
1035 vnlg_spl_[iop][0] = v0;
1036 for ( int i = 1; i < ndft_; i++ )
1037 {
1038 vnlg_spl_[iop][i] = fint[i];
1039 }
1040
1041 // Initialize spline interpolation
1042 // Use zero first derivative at G=0 and natural (y"=0) at Gmax
1043
1044 spline(ndft_,&gspl_[0],&vnlg_spl_[iop][0],
1045 0.0,0.0,0,1,&vnlg_spl2_[iop][0]);
1046 }
1047 else
1048 {
1049 throw SpeciesInitException("Species::initialize_slpp: l < 0 or l > 3");
1050 }
1051 } // loop over projector
1052 return true;
1053 }
1054
1055 ////////////////////////////////////////////////////////////////////////////////
initialize_nlcc()1056 void Species::initialize_nlcc()
1057 {
1058 const double fpi = 4.0 * M_PI;
1059
1060 // allocate array
1061 nlcc_spl_.resize(ndft_);
1062 nlcc_spl2_.resize(ndft_);
1063 nlccg_spl_.resize(ndft_);
1064 nlccg_spl2_.resize(ndft_);
1065
1066 fill( nlcc_spl_.begin(), nlcc_spl_.end(), 0.0 );
1067 copy( nlcc_.begin(), nlcc_.end(), nlcc_spl_.begin() );
1068
1069 // compute spline coefficients
1070 spline(ndft_,&rps_spl_[0],&nlcc_spl_[0],0.0,0.0,0,1,&nlcc_spl2_[0]);
1071
1072 // prepare the function nlccr to be used later in the Bessel transforms
1073 //
1074 // local density: nlcc(G) = 4 pi \int r^2 nlcc(r) sin(Gr)/(Gr) dr
1075 // -> store 4 pi r dr nlcc_spl_ in nlccr(i)
1076 //
1077 // the Bessel transform is then:
1078 // nlcc(G) = 1/G \sum_r sin(Gr) nlccr
1079
1080 for ( int i = 0; i < ndft_; i++ )
1081 {
1082 nlccg_spl_[i] = fpi * rps_spl_[i] * deltar_ * nlcc_spl_[i];
1083 }
1084
1085 // integrate core correction density
1086 // / 2
1087 // | dr 4 pi r n(r)
1088 // /
1089 // nlccr already contains factor 4 pi r
1090 // use fint for integration
1091 vector<double> fint(ndft_);
1092 for ( int i = 0; i < ndft_; i++ )
1093 {
1094 fint[i] = nlccg_spl_[i] * rps_spl_[i];
1095 }
1096 zcore_ = simpsn(ndft_,&fint[0]);
1097
1098 }
1099
phi(int l,double r,double & val)1100 void Species::phi(int l, double r, double &val)
1101 {
1102 val = 0.0;
1103 if ( phi_[l].size() == 0 || l > lmax_ || r > rps_spl_[ndft_-1] )
1104 return;
1105 splint(ndft_,&rps_spl_[0],&phi_spl_[l][0],&phi_spl2_[l][0],r,&val);
1106 }
1107
vpsr(int l,double r,double & v)1108 void Species::vpsr(int l, double r, double &v)
1109 {
1110 if ( l > lmax_ || r > rps_spl_[ndft_-1] )
1111 {
1112 v = 0.0;
1113 }
1114 else
1115 {
1116 splint(ndft_,&rps_spl_[0],&vps_spl_[l][0],&vps_spl2_[l][0],r,&v);
1117 }
1118 }
1119
dvpsr(int l,double r,double & v,double & dv)1120 void Species::dvpsr(int l, double r, double &v, double &dv)
1121 {
1122 if ( l > lmax_ || r > rps_spl_[ndft_-1] )
1123 {
1124 v = 0.0;
1125 dv = 0.0;
1126 }
1127 else
1128 {
1129 splintd(ndft_,&rps_spl_[0],&vps_spl_[l][0],&vps_spl2_[l][0],r,&v,&dv);
1130 }
1131 }
1132
vlocg(double g,double & v)1133 void Species::vlocg(double g, double &v)
1134 {
1135 if ( g > gspl_[ndft_-1] )
1136 {
1137 v = 0.0;
1138 }
1139 else
1140 {
1141 splint(ndft_,&gspl_[0],&vlocg_spl_[0],&vlocg_spl2_[0],g,&v);
1142 }
1143 }
1144
dvlocg(double g,double & v,double & dv)1145 void Species::dvlocg(double g, double &v, double &dv)
1146 {
1147 if ( g > gspl_[ndft_-1] )
1148 {
1149 v = 0.0;
1150 dv = 0.0;
1151 }
1152 else
1153 {
1154 splintd(ndft_,&gspl_[0],&vlocg_spl_[0],&vlocg_spl2_[0],g,&v,&dv);
1155 }
1156 }
1157
vnlg(int iop,double g,double & v)1158 void Species::vnlg(int iop, double g, double &v)
1159 {
1160 assert ( iop >= 0 && iop < nop_ );
1161 int l = lmap_[iop];
1162 if ( l == llocal_ || g > gspl_[ndft_-1] )
1163 {
1164 v = 0.0;
1165 }
1166 else
1167 {
1168 splint(ndft_,&gspl_[0],&vnlg_spl_[iop][0],&vnlg_spl2_[iop][0],g,&v);
1169 }
1170 }
1171
dvnlg(int iop,double g,double & v,double & dv)1172 void Species::dvnlg(int iop, double g, double &v, double &dv)
1173 {
1174 assert ( iop >= 0 && iop < nop_ );
1175 int l = lmap_[iop];
1176 if ( l == llocal_ || g > gspl_[ndft_-1] )
1177 {
1178 v = 0.0;
1179 dv = 0.0;
1180 }
1181 else
1182 {
1183 splintd(ndft_,&gspl_[0],&vnlg_spl_[iop][0],&vnlg_spl2_[iop][0],g,&v,&dv);
1184 }
1185 }
1186
rhopsg(double g)1187 double Species::rhopsg( double g )
1188 {
1189 double arg = 0.25 * rcps_ * rcps_ * g * g;
1190 return -ztot_ * exp( -arg );
1191 }
1192
1193 // core density for nonlinear core correction in reciprocal space
rhocoreg(double g,double & rho)1194 void Species::rhocoreg(double g, double &rho)
1195 {
1196 if ( has_nlcc() && g <= gspl_[ndft_-1] )
1197 {
1198 // spline interpolation
1199 splint(ndft_,&gspl_[0],&nlccg_spl_[0],&nlccg_spl2_[0],g,&rho);
1200 }
1201 else
1202 {
1203 rho = 0.0;
1204 }
1205 }
1206
1207 // core density for nonlinear core correction in reciprocal space
drhocoreg(double g,double & rho,double & drho)1208 void Species::drhocoreg(double g, double &rho, double &drho)
1209 {
1210 if ( has_nlcc() && g <= gspl_[ndft_ - 1] )
1211 {
1212 // spline interpolation
1213 splintd(ndft_,&gspl_[0],&nlccg_spl_[0],&nlccg_spl2_[0],g,&rho,&drho);
1214 }
1215 else
1216 {
1217 rho = 0.0;
1218 drho = 0.0;
1219 }
1220 }
1221
operator <<(ostream & os,Species & s)1222 ostream& operator << ( ostream &os, Species &s )
1223 {
1224 s.print(os, false);
1225 return os;
1226 }
1227
print(ostream & os,bool expanded_form)1228 void Species::print(ostream &os, bool expanded_form)
1229 {
1230 // XML output of species
1231 // print in compact form if the uri is known
1232 // and if expanded_form==false
1233 if ( !expanded_form && uri() != "" )
1234 {
1235 os <<"<species name=\"" << name()
1236 << "\" href=\"" << uri() << "\"/>" << endl;
1237 }
1238 else
1239 {
1240 os.setf(ios::scientific,ios::floatfield);
1241 os << setprecision(12);
1242 os <<"<species name=\"" << name() << "\">" << endl;
1243 os << "<description>" << description() << "</description>" << endl;
1244 os << "<symbol>" << symbol() << "</symbol>" << endl;
1245 os << "<atomic_number>" << atomic_number() << "</atomic_number>" << endl;
1246 os << "<mass>" << mass() << "</mass>" << endl;
1247 if ( type_ == NCPP )
1248 {
1249 os << "<norm_conserving_pseudopotential>" << endl;
1250 os << "<valence_charge>" << zval() << "</valence_charge>" << endl;
1251 os << "<lmax>" << lmax() << "</lmax>" << endl;
1252 os << "<llocal>" << llocal() << "</llocal>" << endl;
1253 os << "<nquad>" << nquad() << "</nquad>" << endl;
1254 os << "<rquad>" << rquad() << "</rquad>" << endl;
1255 os << "<mesh_spacing>" << deltar() << "</mesh_spacing>" << endl;
1256 if ( nlcc_.size() > 0 ) print_nlcc(os);
1257 for ( int l = 0; l <= lmax(); l++ )
1258 {
1259 const int size = vps_[l].size();
1260 os << "<projector l=\"" << l << "\" size=\"" << size
1261 << "\">" << endl;
1262 os << "<radial_potential>\n";
1263 for ( int i = 0; i < size; i++ )
1264 os << vps_[l][i] << "\n";
1265 os << "</radial_potential>\n";
1266 os << "<radial_function>\n";
1267 if ( l < phi_.size() && phi_[l].size() == size )
1268 print_radial_function(os,phi_[l]);
1269 os << "</radial_function>\n";
1270 os << "</projector>" << endl;
1271 }
1272 os << "</norm_conserving_pseudopotential>" << endl;
1273 }
1274 else if ( type_ == SLPP )
1275 {
1276 os << "<norm_conserving_semilocal_pseudopotential>" << endl;
1277 os << "<valence_charge>" << zval() << "</valence_charge>" << endl;
1278 os << "<mesh_spacing>" << deltar() << "</mesh_spacing>" << endl;
1279 if ( nlcc_.size() > 0 ) print_nlcc(os);
1280 os << "<local_potential size=\"" << vlocal_.size() << "\">" << endl;
1281 print_radial_function(os,vlocal_);
1282 os << "</local_potential>" << endl;
1283 // print projector
1284 for ( int l = 0; l < proj_.size(); l++ )
1285 {
1286 for ( int n = 0; n < proj_[l].size(); n++ )
1287 {
1288 os << "<projector l=\"" << l << "\" i=\"" << n + 1 << "\" size=\""
1289 << proj_[l][n].size() << "\">" << endl;
1290 // print radial function
1291 print_radial_function(os,proj_[l][n]);
1292 os << "</projector>" << endl;
1293 }
1294 }
1295 // print D matrix
1296 for ( int l = 0; l < d_.size(); l++ )
1297 {
1298 for ( int n = 0; n < d_[l].size(); n++ )
1299 for ( int m = 0; m < d_[l][n].size(); m++ )
1300 os << "<d_ij l=\"" << l << "\" i=\"" << n + 1 << "\" j=\"" << m + 1
1301 << "\">" << d_[l][n][m] << "</d_ij>" << endl;
1302 }
1303 os << "</norm_conserving_semilocal_pseudopotential>" << endl;
1304 }
1305 os << "</species>" << endl;
1306 }
1307 }
1308
1309 // print core density of nonlinear core correction
print_nlcc(std::ostream & os)1310 void Species::print_nlcc(std::ostream& os)
1311 {
1312 os << "<core_density size=\"" << nlcc_.size() << "\">" << endl;
1313 print_radial_function(os,nlcc_);
1314 os << "</core_density>" << endl;
1315 }
1316
1317 // print any radial function
print_radial_function(std::ostream & os,const std::vector<double> & rad_func)1318 void Species::print_radial_function(std::ostream& os,
1319 const std::vector<double>& rad_func)
1320 {
1321 for ( int i = 0; i < rad_func.size(); i++ )
1322 os << rad_func[i] << endl;
1323 }
1324
info(ostream & os)1325 void Species::info(ostream &os)
1326 {
1327 os.setf(ios::left,ios::adjustfield);
1328
1329 if ( uri() != "" )
1330 {
1331 os <<"<species name=\"" << name()
1332 << "\" href=\"" << uri() << "\">" << endl;
1333 }
1334 else
1335 {
1336 os <<"<species name=\"" << name() << "\">" << endl;
1337 }
1338 os << " <description>" << description() << " </description>" << endl;
1339 os << " <symbol>" << symbol() << "</symbol>" << endl;
1340 os << " <atomic_number>" << atomic_number() << "</atomic_number>" << endl;
1341 os << " <mass>" << mass() << "</mass>" << endl;
1342 if ( type_ == NCPP )
1343 {
1344 os << " <norm_conserving_pseudopotential>" << endl;
1345 os << " <valence_charge>" << zval() << "</valence_charge>" << endl;
1346 os << " <lmax>" << lmax() << "</lmax>" << endl;
1347 os << " <llocal>" << llocal() << "</llocal>" << endl;
1348 os << " <nquad>" << nquad() << "</nquad>" << endl;
1349 os << " <rquad>" << rquad() << "</rquad>" << endl;
1350 os << " <mesh_spacing>" << deltar() << "</mesh_spacing>" << endl;
1351 os << " </norm_conserving_pseudopotential>" << endl;
1352 }
1353 else
1354 {
1355 os << " <norm_conserving_semilocal_pseudopotential>" << endl;
1356 os << " <valence_charge>" << zval() << "</valence_charge>" << endl;
1357 os << " <mesh_spacing>" << deltar() << "</mesh_spacing>" << endl;
1358 os << " </norm_conserving_semilocal_pseudopotential>" << endl;
1359 }
1360 os << "</species>" << endl;
1361
1362 // describe type of potential
1363 if ( type_ == NCPP )
1364 {
1365 if ( nquad() == 0 )
1366 {
1367 if ( lmax() == 0 )
1368 os << " local potential" << endl;
1369 else
1370 os << " Kleinman-Bylander potential" << endl;
1371 }
1372 else if ( nquad() > 0 )
1373 {
1374 os << " NCPP potential with " << nquad()
1375 << " quadrature points in [0.0, " << rquad() << "]" << endl;
1376 os << " local (within 1.e-6) beyond r = " << rcut_loc(1.e-6) << endl;
1377 os << " local (within 1.e-5) beyond r = " << rcut_loc(1.e-5) << endl;
1378 os << " local (within 1.e-4) beyond r = " << rcut_loc(1.e-4) << endl;
1379 os << " local (within 1.e-3) beyond r = " << rcut_loc(1.e-3) << endl;
1380 }
1381 }
1382 else if ( type_ == SLPP )
1383 {
1384 os << " SLPP semilocal potential" << endl;
1385 }
1386 os << " rcps_ = " << rcps() << endl;
1387 os.setf(ios::right,ios::adjustfield);
1388 }
1389
1390 ////////////////////////////////////////////////////////////////////////////////
rcut_loc(double epsilon)1391 double Species::rcut_loc(double epsilon)
1392 {
1393 // find radius at which the largest deviation delta_vnl(r) is < epsilon
1394 double delta = 0.0;
1395 int i = ndft_-1;
1396 while ( ( delta < epsilon ) && i > 0 )
1397 {
1398 i--;
1399 for ( int l = 0; l <= lmax_; l++ )
1400 {
1401 if ( l != llocal_ )
1402 {
1403 // compute deviation vps_[l][i]
1404 double dv = fabs( vps_spl_[l][i] - vps_spl_[llocal_][i] );
1405 delta = dv > delta ? dv : delta;
1406 }
1407 }
1408 }
1409 // adjust i so that delta_v[i] < epsilon
1410 if ( i < ndft_-1 ) i++;
1411
1412 return rps_spl_[i];
1413 }
1414
1415 ////////////////////////////////////////////////////////////////////////////////
1416 // helper function that extract l and n from projector index
get_proj_data(int ipr)1417 Species::ProjectorData Species::get_proj_data(int ipr)
1418 {
1419 ProjectorData res;
1420 res.l = -1;
1421 res.n = 0;
1422 int jpr = 0;
1423 for ( int iop = 0; iop < nop_; iop++ )
1424 {
1425 int l = lmap_[iop];
1426 // n labels differentiates between projectors of same l
1427 if ( res.l == l ) res.n++;
1428 // reset n if new angular momentum is found
1429 else
1430 {
1431 res.l = l;
1432 res.n = 0;
1433 }
1434 // there is one projector Y_lm
1435 for ( res.m = -l; res.m <= l; res.m++ )
1436 {
1437 if ( ipr == jpr ) return res;
1438 jpr++;
1439 }
1440 }
1441 // projector out of bounds
1442 assert(false);
1443 res.l = -1;
1444 res.n = 0;
1445 return res;
1446 }
1447
1448 ////////////////////////////////////////////////////////////////////////////////
1449 // extract D matrix
dmatrix(int ipr,int jpr)1450 double Species::dmatrix(int ipr, int jpr)
1451 {
1452 // extract angular momentum of projectors
1453 const ProjectorData iem = get_proj_data(ipr);
1454 const ProjectorData jem = get_proj_data(jpr);
1455
1456 // projectors have only nonzero overlap if diagonal in l and m
1457 if ( iem.l != jem.l || iem.m != jem.m ) return 0;
1458
1459 // if dmatrix is not present identity matrix is assumed
1460 if ( d_.size() <= iem.l || d_[iem.l].size() == 0 ) return 1;
1461
1462 // return dmatrix
1463 return d_[iem.l][iem.n][jem.n];
1464 }
1465