Skip to content

Commit 28d678c

Browse files
lratclifflratcliff
committed
fixing PSP bug and enabling PBE PSPs (previously was LDA style only)
1 parent 5db53b6 commit 28d678c

1 file changed

Lines changed: 39 additions & 24 deletions

File tree

‎src/madness/chem/gth_pseudopotential.h‎

Lines changed: 39 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,7 @@ class ProjRLMFunctor : public FunctionFunctorInterface<double,3> {
9191

9292
ProjRLMFunctor(double alpha, int l, int m, int i, const coord_3d& center)
9393
: alpha(alpha), l(l), m(m), i(i), center(center) {
94-
specialpts.push_back(coord_3d(0.0));
94+
specialpts.push_back(center);
9595
sqrtPI = std::sqrt(constants::pi);
9696
itmp = 2*l + (4*i-1);
9797
itmp2 = 2*(i-1);
@@ -486,33 +486,44 @@ class GTHPseudopotential {
486486
double h00 = 0.0; xmlLnlproj->Attribute("h00", &h00); t_hlij(lvalue, 0, 0) = h00;
487487
double h11 = 0.0; xmlLnlproj->Attribute("h11", &h11); t_hlij(lvalue, 1, 1) = h11;
488488
double h22 = 0.0; xmlLnlproj->Attribute("h22", &h22); t_hlij(lvalue, 2, 2) = h22;
489+
490+
// read off-diag terms in directly, since PBE PSPs don't follow the LDA relations
491+
double h01 = 0.0; xmlLnlproj->Attribute("h01", &h01); t_hlij(lvalue, 0, 1) = h01;
492+
t_hlij(lvalue, 1, 0) = t_hlij(lvalue, 0, 1);
493+
double h02 = 0.0; xmlLnlproj->Attribute("h02", &h02); t_hlij(lvalue, 0, 2) = h02;
494+
t_hlij(lvalue, 2, 0) = t_hlij(lvalue, 0, 2);
495+
double h12 = 0.0; xmlLnlproj->Attribute("h12", &h12); t_hlij(lvalue, 1, 2) = h12;
496+
t_hlij(lvalue, 2, 1) = t_hlij(lvalue, 1, 2);
497+
498+
// we don't use k terms, so for now these are just set to zero and not included in the file
489499
double k00 = 0.0; xmlLnlproj->Attribute("k00", &k00); t_klij(lvalue, 0, 0) = k00;
490500
double k11 = 0.0; xmlLnlproj->Attribute("k11", &k11); t_klij(lvalue, 1, 1) = k11;
491501
double k22 = 0.0; xmlLnlproj->Attribute("k22", &k22); t_klij(lvalue, 2, 2) = k22;
492502
}
493503
// off-diagonal elements
494-
if (lmax >= 0) {
495-
t_hlij(0, 0, 1) = -1./2.*std::sqrt(3./5.)*t_hlij(0, 1, 1);
496-
t_hlij(0, 1, 0) = t_hlij(0, 0, 1);
497-
t_hlij(0, 0, 2) = 1./2.*std::sqrt(5./21.)*t_hlij(0, 2, 2);
498-
t_hlij(0, 2, 0) = t_hlij(0, 0, 2);
499-
t_hlij(0, 1, 2) = -1./2.*std::sqrt(100./63.)*t_hlij(0, 2, 2);
500-
t_hlij(0, 2, 1) = t_hlij(0, 1, 2);
501-
} if (lmax >= 1) {
502-
t_hlij(1, 0, 1) = -1./2.*std::sqrt(5./7.)*t_hlij(1, 1, 1);
503-
t_hlij(1, 1, 0) = t_hlij(1, 0, 1);
504-
t_hlij(1, 0, 2) = 1./6.*std::sqrt(35./11.)*t_hlij(1, 2, 2);
505-
t_hlij(1, 2, 0) = t_hlij(1, 0, 2);
506-
t_hlij(1, 1, 2) = -1./6.*14./std::sqrt(11.)*t_hlij(1, 2, 2);
507-
t_hlij(1, 2, 1) = t_hlij(1, 1, 2);
508-
} if (lmax >= 2) {
509-
t_hlij(2, 0, 1) = -1./2.*std::sqrt(7./9.)*t_hlij(2, 1, 1);
510-
t_hlij(2, 1, 0) = t_hlij(2, 0, 1);
511-
t_hlij(2, 0, 2) = 1./2.*std::sqrt(63./143.)*t_hlij(2, 2, 2);
512-
t_hlij(2, 2, 0) = t_hlij(2, 0, 2);
513-
t_hlij(2, 1, 2) = -1./2.*18./std::sqrt(143.)*t_hlij(2, 2, 2);
514-
t_hlij(2, 2, 1) = t_hlij(2, 1, 2);
515-
}
504+
// this is only correct for LDA PSPs, therefore read directly from file
505+
//if (lmax >= 0) {
506+
// t_hlij(0, 0, 1) = -1./2.*std::sqrt(3./5.)*t_hlij(0, 1, 1);
507+
// t_hlij(0, 1, 0) = t_hlij(0, 0, 1);
508+
// t_hlij(0, 0, 2) = 1./2.*std::sqrt(5./21.)*t_hlij(0, 2, 2);
509+
// t_hlij(0, 2, 0) = t_hlij(0, 0, 2);
510+
// t_hlij(0, 1, 2) = -1./2.*std::sqrt(100./63.)*t_hlij(0, 2, 2);
511+
// t_hlij(0, 2, 1) = t_hlij(0, 1, 2);
512+
//} if (lmax >= 1) {
513+
// t_hlij(1, 0, 1) = -1./2.*std::sqrt(5./7.)*t_hlij(1, 1, 1);
514+
// t_hlij(1, 1, 0) = t_hlij(1, 0, 1);
515+
// t_hlij(1, 0, 2) = 1./6.*std::sqrt(35./11.)*t_hlij(1, 2, 2);
516+
// t_hlij(1, 2, 0) = t_hlij(1, 0, 2);
517+
// t_hlij(1, 1, 2) = -1./6.*14./std::sqrt(11.)*t_hlij(1, 2, 2);
518+
// t_hlij(1, 2, 1) = t_hlij(1, 1, 2);
519+
//} if (lmax >= 2) {
520+
// t_hlij(2, 0, 1) = -1./2.*std::sqrt(7./9.)*t_hlij(2, 1, 1);
521+
// t_hlij(2, 1, 0) = t_hlij(2, 0, 1);
522+
// t_hlij(2, 0, 2) = 1./2.*std::sqrt(63./143.)*t_hlij(2, 2, 2);
523+
// t_hlij(2, 2, 0) = t_hlij(2, 0, 2);
524+
// t_hlij(2, 1, 2) = -1./2.*18./std::sqrt(143.)*t_hlij(2, 2, 2);
525+
// t_hlij(2, 2, 1) = t_hlij(2, 1, 2);
526+
//}
516527

517528
// Copy to main array
518529
localp[atype-1] = t_localp;
@@ -624,13 +635,17 @@ class GTHPseudopotential {
624635
//debug printing
625636
/*tensorT lmat = matrix_inner(world, vpsi, psi, true);
626637
Q el = 0.0;
638+
Q elt = 0.0;
639+
Q enlt = 0.0;
627640
for(int i = 0;i < nocc;++i){
628641
el += occ[i] * lmat(i, i);
642+
elt += lmat(i, i);
643+
enlt += nlmat(i, i);
629644
std::cout << "nloc/loc " << i << " " << occ[i] << " " << nlmat(i,i) << " "<< lmat(i,i) << std::endl;
630645
}
631646
632647
if(world.rank() == 0){
633-
printf("\n enl, el, epot %16.8f %16.8f %16.8f\n", enl, el, enl+el);
648+
printf("\n enl, el, epot, enltot, eltot, %16.8f %16.8f %16.8f %16.8f %16.8f\n", enl, el, enl+el, enlt, elt);
634649
}*/
635650

636651
gaxpy(world, 1.0, vpsi, 1.0, dpsi);

0 commit comments

Comments
 (0)
⚡