MADNESS 0.10.1
gth_pseudopotential.h
Go to the documentation of this file.
1/// \file gth_pseudopotential.cc
2/// \brief GTH pseudopotential functionality
3/// \defgroup moldft The molecular density functional and Hartree-Fock code
4
5#ifndef MADNESS_CHEM_GTH_PSEUDOPOTENTIAL_H__INCLUDED
6#define MADNESS_CHEM_GTH_PSEUDOPOTENTIAL_H__INCLUDED
7
8#include <madness/mra/mra.h>
9#include <madness/external/tinyxml/tinyxml.h>
10
14
15namespace madness {
16
17typedef Tensor<double> tensorT;
18
19double get_charge_from_file(const std::string filename, unsigned int atype);
20
21/*template <typename Q, int NDIM>
22struct function_real2complex_op
23{
24 typedef std::complex<Q> resultT;
25 Tensor<resultT> operator()(const Key<NDIM>& key, const Tensor<Q>& t) const
26 {
27 Tensor<resultT> result(t.ndim(), t.dims());
28 BINARY_OPTIMIZED_ITERATOR(const Q, t, resultT, result, *_p1 = resultT(*_p0,0.0););
29 return result;
30 }
31 template <typename Archive>
32 void serialize(Archive& ar) {}
33};
34
35Function<std::complex<double>,3> function_real2complex(const Function<double,3>& r)
36{
37 return unary_op_coeffs(r, function_real2complex_op<double,3>());
38}*/
39
40class VLocalFunctor : public FunctionFunctorInterface<double,3> {
41private:
42 double Zeff, zi, C1, C2, C3, C4;
44 std::vector<coord_3d> specialpts;
45
46public:
47 VLocalFunctor(double Zeff, double zi,
48 double C1, double C2, double C3, double C4, const coord_3d& center)
49 : Zeff(Zeff), zi(zi), C1(C1), C2(C2),
50 C3(C3), C4(C4), center(center) {
51 specialpts.push_back(center);
52 }
53
54 double operator()(const coord_3d& r) const {
55 const double x = r[0]-center[0]; const double y = r[1]-center[1]; const double z = r[2]-center[2];
56 double rr = std::sqrt(x*x + y*y + z*z);
57 double rs = rr/zi;
58 double rs2 = rs*rs; double rs4 = rs2*rs2; double rs6 = rs2*rs4;
59 return -(Zeff/rr)*erf(rs/std::sqrt(2.0))
60 + std::exp(-0.5*rs2)*
61 (C1 + C2*rs2 + C3*rs4 + C4*rs6);
62 }
63
64 std::vector<coord_3d> special_points() const {return specialpts;}
65
67 return 6;
68 }
69};
70
71class ProjRLMFunctor : public FunctionFunctorInterface<double,3> {
72private:
73 double alpha; // radius
74 int l, m; // i = 1,2,3 and m = 0,1,2,3 and l = 0,1,2,3 (angular momentum) (i just used in constructor)
76 std::vector<coord_3d> specialpts;
77 // gamma of half-integers starting and 1 (which means 1/2)
78 /*const double gamma_data[17] = {1.0, 0.0, 1.0/2.0, 0.0, 3.0/4.0, 0.0, 15.0/8.0, 0.0, 105.0/16.0, 0.0,
79 945.0/32.0, 0.0, 10395.0/64.0, 0.0, 135135.0/128.0, 0.0, 2027025.0/256.0};*/
80 double sqrtPI;
81 int itmp, itmp2;
82 double t1;
83
84public:
85
86 virtual bool supports_vectorized() const {return false;}
87
88 static const double gamma_data[17];
89
90 ProjRLMFunctor(double alpha, int l, int m, int i, const coord_3d& center)
91 : alpha(alpha), l(l), m(m), center(center) {
92 specialpts.push_back(coord_3d(0.0));
93 sqrtPI = std::sqrt(constants::pi);
94 itmp = 2*l + (4*i-1);
95 itmp2 = 2*(i-1);
96 t1 = 1./std::pow(alpha, 0.5*(double)itmp)/std::sqrt(gamma_data[itmp-1]*sqrtPI);
97 }
98
99 double operator()(const coord_3d& r) const {
100 double x = r[0]-center[0]; double y = r[1]-center[1]; double z = r[2]-center[2];
101 double rsq = x*x + y*y + z*z;
102
103 if (rsq > 40.0) return 0.0;
104
105 double rr = std::sqrt(rsq);
106 double rval = t1;
107 const double PI = constants::pi;
108 // Radial part
109 if (itmp2 == 0) {
110 rval *= std::sqrt(2);
111 }
112 else if (itmp2 == 1) {
113 rval *= rr*std::sqrt(2);
114 }
115 else if (itmp2 == 2) {
116 rval *= rsq*std::sqrt(2);
117 }
118 else if (itmp2 == 3) {
119 rval *= rr*rsq*std::sqrt(2);
120 }
121 else if (itmp2 == 4) {
122 rval *= rsq*rsq*std::sqrt(2);
123 }
124 else if (itmp2 == 5) {
125 rval *= rr*rsq*rsq*std::sqrt(2);
126 }
127 else if (itmp2 == 6) {
128 rval *= rsq*rsq*rsq*std::sqrt(2);
129 }
130 else if (itmp2 == 7) {
131 rval *= rr*rsq*rsq*rsq*std::sqrt(2);
132 }
133 // Angular part
134 if (l == 0) {
135 rval *= (1./2.)*std::sqrt(1./PI);
136 } else if (l == 1) {
137 if (m == 0) {
138 rval *= std::sqrt(3./4./PI)*x;
139 }
140 else if (m == 1) {
141 rval *= std::sqrt(3./4./PI)*y;
142 }
143 else if (m == 2) {
144 rval *= std::sqrt(3./4./PI)*z;
145 }
146 else {
147 MADNESS_EXCEPTION("m out of range for l = 1", 0);
148 }
149 } else if (l == 2) {
150 if (m == 0) {
151 rval *= (1./4.)*std::sqrt(5./PI)*(-x*x - y*y + 2*z*z);
152 }
153 else if (m == 1) {
154 rval *= (1./2.)*std::sqrt(15./PI)*(y*z);
155 }
156 else if (m == 2) {
157 rval *= (1./2.)*std::sqrt(15./PI)*(x*z);
158 }
159 else if (m == 3) {
160 rval *= (1./2.)*std::sqrt(15./PI)*(x*y);
161 }
162 else if (m == 4) {
163 rval *= (1./4.)*std::sqrt(15./PI)*(x*x - y*y);
164 }
165 else {
166 MADNESS_EXCEPTION("m out of range for l = 2", 0);
167 }
168 }
169 rval *= std::exp(-0.5*(rsq/alpha/alpha));
170 return rval;
171 }
172
173 virtual bool screened(const coord_3d& c1, const coord_3d& c2) const {
174 double ftol = 1e-12;
175
176 double x1 = c1[0]; double y1 = c1[1]; double z1 = c1[2];
177 double x2 = c2[0]; double y2 = c2[1]; double z2 = c2[2];
178
179 // if center is inside box, then return false
180 // otherwise, look for the closest point and check
181 bool inside = (center[0] >= x1) && (center[0] <= x2) &&
182 (center[1] >= y1) && (center[1] <= y2) &&
183 (center[2] >= z1) && (center[2] <= z2);
184 if (inside) {
185 return false;
186 }
187 else {
188 //printf("GTH_pseudopotential: (point not inside)\n");
189 //print(" c1: ", c1, " c2: ", c2);
190 double minr = 1e10;
191 int ii = -1; int jj = -1; int kk = -1;
192 for (int i = 0; i <= 1; i++) {
193 for (int j = 0; j <= 1; j++) {
194 for (int k = 0; k <= 1; k++) {
195 double x = (i == 0) ? c1[0] : c2[0];
196 double y = (j == 0) ? c1[1] : c2[1];
197 double z = (k == 0) ? c1[2] : c2[2];
198 coord_3d rr = coord_3d{x, y, z} - center;
199 double rsq = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
200 // print("current minr: ", minr, " point: ", {x,y,z}, " center: ", center, " p: ", {x,y,z}, " rsq: ", rsq);
201 if (minr > rsq) {
202 minr = rsq;
203 ii = i; jj = j; kk = k;
204 }
205 }
206 }
207 }
208 //print("ii: ", ii, "jj: ", jj, "kk: ", kk);
209 //printf("\n");
210 if ((ii < 0) || (jj < 0) || (kk < 0)) MADNESS_EXCEPTION("GTH_Pseudopotential: failed to find suitable minimum point\n", 0);
211 double x = (ii == 0) ? c1[0] : c2[0];
212 double y = (jj == 0) ? c1[1] : c2[1];
213 double z = (kk == 0) ? c1[2] : c2[2];
214 double fval = this->operator()({x, y, z});
215 if (fabs(fval) < ftol) return true;
216 else return false;
217 }
218 }
219
220 virtual void operator()(const Vector<double*,3>& xvals, double* MADNESS_RESTRICT fvals, int npts) const {
221
222 double* x = new double[npts];
223 double* y = new double[npts];
224 double* z = new double[npts];
225 double* rsq = new double[npts];
226 double* rr = new double[npts];
227
228 double* x1 = xvals[0]; double* x2 = xvals[1]; double* x3 = xvals[2];
229 for (int i = 0; i < npts; i++) {
230 x[i] = x1[i]-center[0];
231 y[i] = x2[i]-center[1];
232 z[i] = x3[i]-center[2];
233 rsq[i] = x[i]*x[i] + y[i]*y[i] + z[i]*z[i];
234 rr[i] = std::sqrt(rsq[i]);
235 fvals[i] = t1;
236 }
237
238 const double PI = constants::pi;
239
240 // Radial part
241 if (itmp2 == 0) {
242 for (int i = 0; i < npts; i++) {
243 fvals[i] *= std::sqrt(2);
244 }
245 }
246 else if (itmp2 == 1) {
247 for (int i = 0; i < npts; i++) {
248 fvals[i] *= rr[i]*std::sqrt(2);
249 }
250 }
251 else if (itmp2 == 2) {
252 for (int i = 0; i < npts; i++) {
253 fvals[i] *= rsq[i]*std::sqrt(2);
254 }
255 }
256 else if (itmp2 == 3) {
257 for (int i = 0; i < npts; i++) {
258 fvals[i] *= rr[i]*rsq[i]*std::sqrt(2);
259 }
260 }
261 else if (itmp2 == 4) {
262 for (int i = 0; i < npts; i++) {
263 fvals[i] *= rsq[i]*rsq[i]*std::sqrt(2);
264 }
265 }
266 else if (itmp2 == 5) {
267 for (int i = 0; i < npts; i++) {
268 fvals[i] *= rr[i]*rsq[i]*rsq[i]*std::sqrt(2);
269 }
270 }
271 else if (itmp2 == 6) {
272 for (int i = 0; i < npts; i++) {
273 fvals[i] *= rsq[i]*rsq[i]*rsq[i]*std::sqrt(2);
274 }
275 }
276 else if (itmp2 == 7) {
277 for (int i = 0; i < npts; i++) {
278 fvals[i] *= rr[i]*rsq[i]*rsq[i]*rsq[i];
279 }
280 }
281 // Angular part
282 if (l == 0) {
283 for (int i = 0; i < npts; i++) {
284 fvals[i] *= (1./2.)*std::sqrt(1./PI)*std::exp(-0.5*(rsq[i]/alpha/alpha));
285 }
286 } else if (l == 1) {
287 if (m == 0) {
288 for (int i = 0; i < npts; i++) {
289 fvals[i] *= std::sqrt(3./4./PI)*x[i]*std::exp(-0.5*(rsq[i]/alpha/alpha));
290 }
291 }
292 else if (m == 1) {
293 for (int i = 0; i < npts; i++) {
294 fvals[i] *= std::sqrt(3./4./PI)*y[i]*std::exp(-0.5*(rsq[i]/alpha/alpha));
295 }
296 }
297 else if (m == 2) {
298 for (int i = 0; i < npts; i++) {
299 fvals[i] *= std::sqrt(3./4./PI)*z[i]*std::exp(-0.5*(rsq[i]/alpha/alpha));
300 }
301 }
302 else {
303 MADNESS_EXCEPTION("m out of range for l = 1", 0);
304 }
305 } else if (l == 2) {
306 if (m == 0) {
307 for (int i = 0; i < npts; i++) {
308 fvals[i] *= (1./4.)*std::sqrt(5./PI)*(-x[i]*x[i] - y[i]*y[i] + 2*z[i]*z[i])*std::exp(-0.5*(rsq[i]/alpha/alpha));
309 }
310 }
311 else if (m == 1) {
312 for (int i = 0; i < npts; i++) {
313 fvals[i] *= (1./2.)*std::sqrt(15./PI)*(y[i]*z[i])*std::exp(-0.5*(rsq[i]/alpha/alpha));
314 }
315 }
316 else if (m == 2) {
317 for (int i = 0; i < npts; i++) {
318 fvals[i] *= (1./2.)*std::sqrt(15./PI)*(x[i]*z[i])*std::exp(-0.5*(rsq[i]/alpha/alpha));
319 }
320 }
321 else if (m == 3) {
322 for (int i = 0; i < npts; i++) {
323 fvals[i] *= (1./2.)*std::sqrt(15./PI)*(x[i]*y[i])*std::exp(-0.5*(rsq[i]/alpha/alpha));
324 }
325 }
326 else if (m == 4) {
327 for (int i = 0; i < npts; i++) {
328 fvals[i] *= (1./4.)*std::sqrt(15./PI)*(x[i]*x[i] - y[i]*y[i])*std::exp(-0.5*(rsq[i]/alpha/alpha));
329 }
330 }
331 else {
332 MADNESS_EXCEPTION("m out of range for l = 2", 0);
333 }
334 }
335 //for (int i = 0; i < npts; i++) {
336 // fvals[i] *= std::exp(-0.5*(rsq[i]/alpha/alpha));
337 //}
338
339 delete [] x;
340 delete [] y;
341 delete [] z;
342 delete [] rsq;
343 delete [] rr;
344 }
345
346 std::vector<coord_3d> special_points() const {return specialpts;}
347
349 return 6;
350 }
351};
352
354private:
355 int maxL;
358
359public:
362
363 real_function_3d nlmproj(World& world, int l, int m, int i) {
364 // real_function_3d f1 = (m < 2*l+1) ?
365 // real_factory_3d(world).functor(real_functor_3d(new ProjRLMFunctor(radii(l), l, m, i, center))).
366 // truncate_on_project().nofence().truncate_mode(0) : real_factory_3d(world);
367
369 if (m < 2*l+1) {
370 auto functor = real_functor_3d(new ProjRLMFunctor(radii(l), l, m, i, center));
371 f1 = real_factory_3d(world).functor(functor).truncate_on_project().nofence().truncate_mode(0);
372 }
373 else {
374 f1 = real_factory_3d(world);
375 }
376
377 return f1;
378 }
379
380 ProjRLMFunctor nlmproj_functor(World& world, int l, int m, int i) {
381 return ProjRLMFunctor(radii(l), l, m, i, center);
382 }
383};
384
385template <typename Q>
387private:
388public:
390 std::array<real_tensor,118> localp;
391 std::array<real_tensor,118> radii;
392 std::array<real_tensor,118> hlij;
393 std::array<real_tensor,118> klij;
395 std::vector<unsigned int> atoms_with_projectors;
396
397
398public:
400
402 // Load info from file
403 load_pseudo_from_file(world, "gth.xml");
404 atoms_with_projectors.clear();
405
406 // fill list with atoms-with-projectors (i.e. not H or He)
407 for (size_t iatom = 0; iatom < molecule.natom(); iatom++) {
408 Atom atom = molecule.get_atom(iatom);
409
410 //make sure this is actually a pseudo-atom
411 if (!atom.pseudo_atom) continue;
412
413 unsigned int atype = atom.atomic_number;
414 if (radii[atype-1].dim(0) > 0)
415 atoms_with_projectors.push_back(iatom);
416 }
417
418 vlocalp = real_factory_3d(world);
420 for (size_t iatom = 0; iatom < molecule.natom(); iatom++) {
421 // Get atom and its associated GTH tensors
422 Atom atom = molecule.get_atom(iatom);
423
424 //make sure this is actually a pseudo-atom
425 if (!atom.pseudo_atom) continue;
426
427 coord_3d center = atom.get_coords();
428 unsigned int atype = atom.atomic_number;
429 // do local part
430 real_tensor atom_localp = localp[atype-1];
431 real_function_3d temp = real_factory_3d(world).functor(
432 real_functor_3d(new
433 VLocalFunctor(atom_localp[0], atom_localp[1], atom_localp[2], atom_localp[3], atom_localp[4], atom_localp[5], center))).
434 truncate_mode(0).truncate_on_project();
435 temp.compress();
436 //vlocalp += temp;
437 vlocalp.gaxpy(1.0, temp, 1.0, true);
438 }
439 }
440
442 return vlocalp;
443 }
444
445 void reproject(int k, double thresh) {
447 }
448
449 void load_pseudo_from_file(World& world, const std::string filename) {
450 bool debug = true;
451
452 TiXmlDocument doc(filename);
453 if (!doc.LoadFile()) {
454 MADNESS_EXCEPTION("Failed to load GTH pseudopotential file", 0);
455 }
456
457 for (size_t iatom = 0; iatom < molecule.natom(); iatom++) {
458 Atom atom = molecule.get_atom(iatom);
459 unsigned int atype = atom.atomic_number;
460 if (debug && world.rank() == 0) {printf("atom atomic_number = %d\n", atype);}
461
462 bool success = false;
463 for (TiXmlElement* node=doc.FirstChildElement(); node && !success; node=node->NextSiblingElement()) {
464 if (strcmp(node->Value(),"name") == 0) {
465 std::string name = node->GetText();
466 if (debug && world.rank() == 0) std::cout << "Loading pseudopotential file " << name << std::endl;
467 }
468 else if (strcmp(node->Value(), "atom") == 0) {
469 const char* symbol = node->Attribute("symbol");
470 unsigned int atn = symbol_to_atomic_number(symbol);
471 if (atype == atn) {
472 success = true;
473 if (debug && world.rank() == 0) std::cout << " found atomic pseudopotential " << symbol << std::endl;
474 int lmax = -1;
475 node->Attribute("lmax", &lmax);
476 if (debug && world.rank() == 0) std::cout << " maximum L is " << lmax << std::endl;
477 real_tensor t_radii((long)lmax+1);
478 real_tensor t_hlij((long)lmax+1, (long)3, (long)3);
479 real_tensor t_klij((long)lmax+1, (long)3, (long)3);
480 // local part
481 TiXmlElement* xmlVLocal = node->FirstChildElement();
482 real_tensor t_localp((long)6);
483 double zeff = 0.0; xmlVLocal->Attribute("Zeff", &zeff); t_localp[0] = zeff;
484 double lradius = 0.0; xmlVLocal->Attribute("radius", &lradius); t_localp(1) = lradius;
485 double C1 = 0.0; xmlVLocal->Attribute("C1", &C1); t_localp[2] = C1;
486 double C2 = 0.0; xmlVLocal->Attribute("C2", &C2); t_localp[3] = C2;
487 double C3 = 0.0; xmlVLocal->Attribute("C3", &C3); t_localp[4] = C3;
488 double C4 = 0.0; xmlVLocal->Attribute("C4", &C4); t_localp[5] = C4;
489 // loop through nonlocal part
490 for (TiXmlElement* xmlLnlproj = xmlVLocal->NextSiblingElement();
491 xmlLnlproj; xmlLnlproj=xmlLnlproj->NextSiblingElement()) {
492 int lvalue = -1; xmlLnlproj->Attribute("l", &lvalue);
493 double radius = 0.0; xmlLnlproj->Attribute("radius", &radius); t_radii[lvalue] = radius;
494 double h00 = 0.0; xmlLnlproj->Attribute("h00", &h00); t_hlij(lvalue, 0, 0) = h00;
495 double h11 = 0.0; xmlLnlproj->Attribute("h11", &h11); t_hlij(lvalue, 1, 1) = h11;
496 double h22 = 0.0; xmlLnlproj->Attribute("h22", &h22); t_hlij(lvalue, 2, 2) = h22;
497 double k00 = 0.0; xmlLnlproj->Attribute("k00", &k00); t_klij(lvalue, 0, 0) = k00;
498 double k11 = 0.0; xmlLnlproj->Attribute("k11", &k11); t_klij(lvalue, 1, 1) = k11;
499 double k22 = 0.0; xmlLnlproj->Attribute("k22", &k22); t_klij(lvalue, 2, 2) = k22;
500 }
501 // off-diagonal elements
502 if (lmax >= 0) {
503 t_hlij(0, 0, 1) = -1./2.*std::sqrt(3./5.)*t_hlij(0, 1, 1);
504 t_hlij(0, 1, 0) = t_hlij(0, 0, 1);
505 t_hlij(0, 0, 2) = 1./2.*std::sqrt(5./21.)*t_hlij(0, 2, 2);
506 t_hlij(0, 2, 0) = t_hlij(0, 0, 2);
507 t_hlij(0, 1, 2) = -1./2.*std::sqrt(100./63.)*t_hlij(0, 2, 2);
508 t_hlij(0, 2, 1) = t_hlij(0, 1, 2);
509 } if (lmax >= 1) {
510 t_hlij(1, 0, 1) = -1./2.*std::sqrt(5./7.)*t_hlij(1, 1, 1);
511 t_hlij(1, 1, 0) = t_hlij(1, 0, 1);
512 t_hlij(1, 0, 2) = 1./6.*std::sqrt(35./11.)*t_hlij(1, 2, 2);
513 t_hlij(1, 2, 0) = t_hlij(1, 0, 2);
514 t_hlij(1, 1, 2) = -1./6.*14./std::sqrt(11.)*t_hlij(1, 2, 2);
515 t_hlij(1, 2, 1) = t_hlij(1, 1, 2);
516 } if (lmax >= 2) {
517 t_hlij(2, 0, 1) = -1./2.*std::sqrt(7./9.)*t_hlij(2, 1, 1);
518 t_hlij(2, 1, 0) = t_hlij(2, 0, 1);
519 t_hlij(2, 0, 2) = 1./2.*std::sqrt(63./143.)*t_hlij(2, 2, 2);
520 t_hlij(2, 2, 0) = t_hlij(2, 0, 2);
521 t_hlij(2, 1, 2) = -1./2.*18./std::sqrt(143.)*t_hlij(2, 2, 2);
522 t_hlij(2, 2, 1) = t_hlij(2, 1, 2);
523 }
524
525 // Copy to main array
526 localp[atype-1] = t_localp;
527 radii[atype-1] = t_radii;
528 hlij[atype-1] = t_hlij;
529 klij[atype-1] = t_klij;
530 }
531 }
532 }
533 }
534 }
535
536
537 std::vector<Function<Q,3> > apply_potential(World& world, const real_function_3d& potential, const std::vector<Function<Q,3> >& psi, const tensorT & occ, Q & enl) {
539 double vtol = 1e-2*thresh;
540 std::vector<Function<Q,3> > vpsi = mul_sparse(world,(potential), psi, vtol);
541
542 unsigned int norbs = psi.size();
543 unsigned int natoms = atoms_with_projectors.size();
544
545 // NEW (VECTORIZED) ... hopefully
546 vector_real_function_3d localproj;
547 int lidx = 0;
548 unsigned int maxLL = 0;
549 // loop through all of the atom types in the molecule and get the maximum L value
550 // we need this because we are going to create a fixed sized tensor to store
551 // mapping between a linear index and (ias=which atom, i=which projector, l=angular momentum,
552 // m=angular momentum projection)
553 for (unsigned int iatom = 0; iatom < natoms; iatom++) {
554 // Get atom and its associated GTH tensors
556 unsigned int atype = atom.atomic_number;
557 real_tensor& atom_radii = radii[atype-1];
558 if (atom_radii.dim(0) > 0)
559 maxLL = std::max(maxLL,(unsigned int)atom_radii.dim(0)-1);
560 }
561
562 // Pilm_lookup is a mapping between a linear index and (ias,i,l,m)
563 Tensor<int> Pilm_lookup((unsigned int) natoms, (unsigned long) 3, (unsigned long) maxLL+1, (unsigned long) 2*maxLL+1);
564 for (unsigned int iatom = 0; iatom < natoms; iatom++) {
565 // Get atom and its associated GTH tensors
567 coord_3d center = atom.get_coords();
568 unsigned int atype = atom.atomic_number;
569 real_tensor& atom_radii = radii[atype-1];
570 //real_tensor& atom_hlij = hlij[atype-1];
571
572 // Create function stores for projectors
573 ProjRLMStore prlmstore(atom_radii, center);
574 for (unsigned int j = 1; j <= 3; j++) {
575 for (unsigned int l = 0; l <= maxLL; l++) {
576 for (unsigned int m = 0; m < 2*maxLL+1; m++) {
577 Pilm_lookup(iatom, j-1, l, m) = lidx++;
578 if (m < 2*l+1) localproj.push_back(prlmstore.nlmproj(world,l,m,j));
579 else localproj.push_back(real_factory_3d(world));
580 }
581 }
582 }
583 // Somehow this scares me ... thread-safety of the container localproj (???)
584 world.gop.fence();
585 }
586 truncate(world, localproj, FunctionDefaults<3>::get_thresh());
587 compress(world, localproj);
588 //truncate(world, psi, FunctionDefaults<3>::get_thresh());
589 compress(world, psi);
590 //truncate(world, vpsi, FunctionDefaults<3>::get_thresh());
591 compress(world, vpsi);
592
593 Tensor<Q> Pilm = matrix_inner(world, localproj, psi);
594 Pilm = Pilm.reshape(natoms, 3, maxLL+1, 2*maxLL+1, norbs);
595
596 Tensor<Q> Qilm((unsigned int) natoms, (unsigned long) 3, (unsigned long) maxLL+1, (unsigned long) 2*maxLL+1, (unsigned int) norbs);
597 for (unsigned int iorb=0; iorb<psi.size(); iorb++) {
598 for (unsigned int iatom = 0; iatom < natoms; iatom++) {
599 // Get atom and its associated GTH tensors
601 unsigned int atype = atom.atomic_number;
602 real_tensor& atom_radii = radii[atype-1];
603 real_tensor& atom_hlij = hlij[atype-1];
604 int maxL = atom_radii.dim(0)-1;
605 for (unsigned int i = 1; i <= 3; i++) {
606 for (int l = 0; l <= maxL; l++) {
607 for (int m = 0; m < 2*l+1; m++) {
608 Q s = 0.0;
609 for (unsigned int j = 1; j <= 3; j++) {
610 s += atom_hlij(l,i-1,j-1)*Pilm(iatom, j-1,l,m,iorb);
611 }
612 Qilm(iatom, i-1,l,m,iorb) = s;
613 }
614 }
615 }
616 }
617 }
618 Qilm = Qilm.reshape(natoms*3*(maxLL+1)*(2*maxLL+1),norbs);
619
620 double vtol2 = 1e-4*thresh;
621 double trantol = vtol2 / std::min(30.0, double(localproj.size()));
622 vector_real_function_3d dpsi = transform(world, localproj, Qilm, trantol, true);
623
624 // calculate non-local energy
625 tensorT nlmat = matrix_inner(world, dpsi, psi, true);
626 int nocc = occ.size();
627 enl = 0.0;
628 for(int i = 0;i < nocc;++i){
629 enl += occ[i] * nlmat(i, i);
630 }
631
632 //debug printing
633 /*tensorT lmat = matrix_inner(world, vpsi, psi, true);
634 Q el = 0.0;
635 for(int i = 0;i < nocc;++i){
636 el += occ[i] * lmat(i, i);
637 std::cout << "nloc/loc " << i << " " << occ[i] << " " << nlmat(i,i) << " "<< lmat(i,i) << std::endl;
638 }
639
640 if(world.rank() == 0){
641 printf("\n enl, el, epot %16.8f %16.8f %16.8f\n", enl, el, enl+el);
642 }*/
643
644 gaxpy(world, 1.0, vpsi, 1.0, dpsi);
645
646 return vpsi;
647 }
648
649 std::vector<Function<Q,3> > apply_potential_simple(World& world, const real_function_3d& potential, const std::vector<Function<Q,3> >& psi, const tensorT & occ, Q & enl) {
651 double vtol = 1e-2*thresh;
652 std::vector<Function<Q,3> > vpsi = mul_sparse(world,(potential), psi, vtol);
653
654 //unsigned int norbs = psi.size();
655 unsigned int natoms = atoms_with_projectors.size();
656
657 for (unsigned int iatom = 0; iatom < natoms; iatom++) {
659 unsigned int atype = atom.atomic_number;
660 real_tensor& atom_radii = radii[atype-1];
661 real_tensor& atom_hlij = hlij[atype-1];
662 coord_3d center = atom.get_coords();
663 ProjRLMStore prlmstore(atom_radii, center);
664 unsigned int maxLL = atom_radii.dim(0)-1;
665 for (unsigned int l = 0; l <= maxLL; l++) {
666 for (unsigned int m = 0; m < 2*l+1; m++) {
667 for (unsigned int i = 1; i <= 3; i++) {
668 real_function_3d fproji = prlmstore.nlmproj(world,l,m,i);
669 for (unsigned int j = 1; j <= 3; j++) {
670 real_function_3d fprojj = prlmstore.nlmproj(world,l,m,j);
671 for (unsigned int iorb = 0; iorb < psi.size(); iorb++) {
672 double val = atom_hlij(l,i-1,j-1)*(fprojj.inner(psi[iorb]));
673 if (std::abs(val) > vtol*1e-2) {
674 vpsi[iorb] += val*fproji;
675 }
676 }
677 }
678 }
679 }
680 }
681 }
682 return vpsi;
683 }
684};
685
686
687
688
689}
690
691#endif // MADNESS_CHEM_GTH_PSEUDOPOTENTIAL_H__INCLUDED
double potential(const coord_3d &r)
Definition 3dharmonic.cc:132
Definition molecule.h:58
unsigned int atomic_number
Atomic number.
Definition molecule.h:61
madness::Vector< double, 3 > get_coords() const
Definition molecule.h:99
bool pseudo_atom
Indicates if this atom uses a pseudopotential.
Definition molecule.h:63
long dim(int i) const
Returns the size of dimension i.
Definition basetensor.h:147
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:204
static const double & get_thresh()
Returns the default threshold.
Definition funcdefaults.h:279
Abstract base class interface required for functors used as input to Functions.
Definition function_interface.h:68
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition mra.h:981
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition mra.h:721
Definition gth_pseudopotential.h:386
std::vector< Function< Q, 3 > > apply_potential(World &world, const real_function_3d &potential, const std::vector< Function< Q, 3 > > &psi, const tensorT &occ, Q &enl)
Definition gth_pseudopotential.h:537
void make_pseudo_potential(World &world)
Definition gth_pseudopotential.h:401
std::array< real_tensor, 118 > klij
Definition gth_pseudopotential.h:393
void reproject(int k, double thresh)
Definition gth_pseudopotential.h:445
real_function_3d vlocalp
Definition gth_pseudopotential.h:394
std::vector< unsigned int > atoms_with_projectors
Definition gth_pseudopotential.h:395
GTHPseudopotential(World &world, Molecule molecule)
Definition gth_pseudopotential.h:399
std::array< real_tensor, 118 > hlij
Definition gth_pseudopotential.h:392
std::array< real_tensor, 118 > localp
Definition gth_pseudopotential.h:390
std::vector< Function< Q, 3 > > apply_potential_simple(World &world, const real_function_3d &potential, const std::vector< Function< Q, 3 > > &psi, const tensorT &occ, Q &enl)
Definition gth_pseudopotential.h:649
void load_pseudo_from_file(World &world, const std::string filename)
Definition gth_pseudopotential.h:449
std::array< real_tensor, 118 > radii
Definition gth_pseudopotential.h:391
real_function_3d vlocalpot()
Definition gth_pseudopotential.h:441
Molecule molecule
Definition gth_pseudopotential.h:389
Definition molecule.h:124
const Atom & get_atom(unsigned int i) const
Definition molecule.cc:447
size_t natom() const
Definition molecule.h:387
Definition gth_pseudopotential.h:71
int itmp
Definition gth_pseudopotential.h:81
static const double gamma_data[17]
Definition gth_pseudopotential.h:88
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition gth_pseudopotential.h:346
virtual bool supports_vectorized() const
Does the interface support a vectorized operator()?
Definition gth_pseudopotential.h:86
double t1
Definition gth_pseudopotential.h:82
virtual bool screened(const coord_3d &c1, const coord_3d &c2) const
Definition gth_pseudopotential.h:173
double alpha
Definition gth_pseudopotential.h:73
std::vector< coord_3d > specialpts
Definition gth_pseudopotential.h:76
ProjRLMFunctor(double alpha, int l, int m, int i, const coord_3d &center)
Definition gth_pseudopotential.h:90
int l
Definition gth_pseudopotential.h:74
double operator()(const coord_3d &r) const
Definition gth_pseudopotential.h:99
double sqrtPI
Definition gth_pseudopotential.h:80
virtual void operator()(const Vector< double *, 3 > &xvals, double *MADNESS_RESTRICT fvals, int npts) const
Definition gth_pseudopotential.h:220
Level special_level()
Override this change level refinement for special points (default is 6)
Definition gth_pseudopotential.h:348
int m
Definition gth_pseudopotential.h:74
int itmp2
Definition gth_pseudopotential.h:81
coord_3d center
Definition gth_pseudopotential.h:75
Definition gth_pseudopotential.h:353
coord_3d center
Definition gth_pseudopotential.h:357
real_function_3d nlmproj(World &world, int l, int m, int i)
Definition gth_pseudopotential.h:363
ProjRLMStore(const real_tensor &radii, const coord_3d &center)
Definition gth_pseudopotential.h:360
real_tensor radii
Definition gth_pseudopotential.h:356
ProjRLMFunctor nlmproj_functor(World &world, int l, int m, int i)
Definition gth_pseudopotential.h:380
int maxL
Definition gth_pseudopotential.h:355
Tensor< T > reshape(int ndimnew, const long *d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition tensor.h:1384
Definition gth_pseudopotential.h:40
double operator()(const coord_3d &r) const
Definition gth_pseudopotential.h:54
double C3
Definition gth_pseudopotential.h:42
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition gth_pseudopotential.h:64
double C1
Definition gth_pseudopotential.h:42
Level special_level()
Override this change level refinement for special points (default is 6)
Definition gth_pseudopotential.h:66
double C4
Definition gth_pseudopotential.h:42
double Zeff
Definition gth_pseudopotential.h:42
double zi
Definition gth_pseudopotential.h:42
double C2
Definition gth_pseudopotential.h:42
coord_3d center
Definition gth_pseudopotential.h:43
std::vector< coord_3d > specialpts
Definition gth_pseudopotential.h:44
VLocalFunctor(double Zeff, double zi, double C1, double C2, double C3, double C4, const coord_3d &center)
Definition gth_pseudopotential.h:47
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition worldgop.cc:161
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:318
WorldGopInterface & gop
Global operations.
Definition world.h:205
double(* f1)(const coord_3d &)
Definition derivatives.cc:55
static bool debug
Definition dirac-hatom.cc:16
double psi(const Vector< double, 3 > &r)
Definition hatom_energy.cc:78
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition madness_exception.h:119
Main include file for MADNESS and defines Function interface.
const double pi
Mathematical constant .
Definition constants.h:48
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
static const char * filename
Definition legendre.cc:96
void truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition vmra.h:689
Tensor< double > tensorT
Definition distpm.cc:21
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition vmra.h:133
std::shared_ptr< FunctionFunctorInterface< double, 3 > > real_functor_3d
Definition functypedefs.h:107
std::vector< real_function_3d > vector_real_function_3d
Definition functypedefs.h:79
unsigned int symbol_to_atomic_number(const std::string &symbol)
Definition atomutil.cc:173
int Level
Definition key.h:55
Function< TENSOR_RESULT_TYPE(L, R), NDIM > mul_sparse(const Function< L, NDIM > &left, const Function< R, NDIM > &right, double tol, bool fence=true)
Sparse multiplication — left and right must be reconstructed and if tol!=0 have tree of norms already...
Definition mra.h:1742
FunctionFactory< double, 3 > real_factory_3d
Definition functypedefs.h:93
double get_charge_from_file(const std::string filename, unsigned int atype)
Definition gth_pseudopotential.cc:8
Vector< double, 3 > coord_3d
Definition funcplot.h:1042
Function< T, NDIM > project(const Function< T, NDIM > &other, int k=FunctionDefaults< NDIM >::get_k(), double thresh=FunctionDefaults< NDIM >::get_thresh(), bool fence=true)
Definition mra.h:2399
std::string name(const FuncType &type, const int ex=-1)
Definition ccpairfunction.h:28
void matrix_inner(DistributedMatrix< T > &A, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition distpm.cc:46
void gaxpy(const double a, ScalarResult< T > &left, const double b, const T &right, const bool fence=true)
the result type of a macrotask must implement gaxpy
Definition macrotaskq.h:140
static long abs(long a)
Definition tensor.h:218
double Q(double a)
Definition relops.cc:20
static const double PI
Definition relops.cc:12
static const double m
Definition relops.cc:9
static const double thresh
Definition rk.cc:45
static const long k
Definition rk.cc:44
void e()
Definition test_sig.cc:75
static const int truncate_mode
Definition testcosine.cc:14