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 
15 namespace madness {
16 
17 typedef Tensor<double> tensorT;
18 
19 double get_charge_from_file(const std::string filename, unsigned int atype);
20 
21 /*template <typename Q, int NDIM>
22 struct 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 
35 Function<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 
40 class VLocalFunctor : public FunctionFunctorInterface<double,3> {
41 private:
42  double Zeff, zi, C1, C2, C3, C4;
44  std::vector<coord_3d> specialpts;
45 
46 public:
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 
71 class ProjRLMFunctor : public FunctionFunctorInterface<double,3> {
72 private:
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 
84 public:
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 
354 private:
355  int maxL;
358 
359 public:
361  : maxL(radii.dim(0)), radii(radii), center(center) {}
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 
385 template <typename Q>
387 private:
388 public:
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 
398 public:
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);
419  vlocalp.compress();
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
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition: mra.h:721
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
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
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
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
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
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
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
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
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 > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: gth_pseudopotential.h:64
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
const double m
Definition: gfit.cc:199
double psi(const Vector< double, 3 > &r)
Definition: hatom_energy.cc:78
static double pow(const double *a, const double *b)
Definition: lda.h:74
#define max(a, b)
Definition: lda.h:51
#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
File holds all helper structures necessary for the CC_Operator and CC2 class.
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
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
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
FunctionFactory< double, 3 > real_factory_3d
Definition: functypedefs.h:93
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
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
std::string name(const FuncType &type, const int ex=-1)
Definition: ccpairfunction.h:28
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
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 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