5#ifndef MADNESS_CHEM_GTH_PSEUDOPOTENTIAL_H__INCLUDED
6#define MADNESS_CHEM_GTH_PSEUDOPOTENTIAL_H__INCLUDED
9#include <madness/external/tinyxml/tinyxml.h>
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);
101 double rsq = x*x + y*y + z*z;
103 if (
rsq > 40.0)
return 0.0;
105 double rr = std::sqrt(
rsq);
110 rval *= std::sqrt(2);
112 else if (
itmp2 == 1) {
115 else if (
itmp2 == 2) {
118 else if (
itmp2 == 3) {
121 else if (
itmp2 == 4) {
124 else if (
itmp2 == 5) {
127 else if (
itmp2 == 6) {
130 else if (
itmp2 == 7) {
135 rval *= (1./2.)*std::sqrt(1./
PI);
138 rval *= std::sqrt(3./4./
PI)*x;
141 rval *= std::sqrt(3./4./
PI)*y;
144 rval *= std::sqrt(3./4./
PI)*z;
151 rval *= (1./4.)*std::sqrt(5./
PI)*(-x*x - y*y + 2*z*z);
154 rval *= (1./2.)*std::sqrt(15./
PI)*(y*z);
157 rval *= (1./2.)*std::sqrt(15./
PI)*(x*z);
160 rval *= (1./2.)*std::sqrt(15./
PI)*(x*y);
163 rval *= (1./4.)*std::sqrt(15./
PI)*(x*x - y*y);
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];
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];
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];
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];
229 for (
int i = 0; i <
npts; i++) {
233 rsq[i] = x[i]*x[i] + y[i]*y[i] + z[i]*z[i];
234 rr[i] = std::sqrt(
rsq[i]);
242 for (
int i = 0; i <
npts; i++) {
243 fvals[i] *= std::sqrt(2);
246 else if (
itmp2 == 1) {
247 for (
int i = 0; i <
npts; i++) {
248 fvals[i] *=
rr[i]*std::sqrt(2);
251 else if (
itmp2 == 2) {
252 for (
int i = 0; i <
npts; i++) {
256 else if (
itmp2 == 3) {
257 for (
int i = 0; i <
npts; i++) {
261 else if (
itmp2 == 4) {
262 for (
int i = 0; i <
npts; i++) {
266 else if (
itmp2 == 5) {
267 for (
int i = 0; i <
npts; i++) {
271 else if (
itmp2 == 6) {
272 for (
int i = 0; i <
npts; i++) {
276 else if (
itmp2 == 7) {
277 for (
int i = 0; i <
npts; i++) {
283 for (
int i = 0; i <
npts; i++) {
288 for (
int i = 0; i <
npts; i++) {
293 for (
int i = 0; i <
npts; i++) {
298 for (
int i = 0; i <
npts; i++) {
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));
312 for (
int i = 0; i <
npts; i++) {
317 for (
int i = 0; i <
npts; i++) {
322 for (
int i = 0; i <
npts; i++) {
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));
371 f1 =
real_factory_3d(world).functor(functor).truncate_on_project().nofence().truncate_mode(0);
392 std::array<real_tensor,118>
hlij;
393 std::array<real_tensor,118>
klij;
453 if (!
doc.LoadFile()) {
466 if (
debug && world.
rank() == 0) std::cout <<
"Loading pseudopotential file " <<
name << std::endl;
468 else if (
strcmp(
node->Value(),
"atom") == 0) {
469 const char* symbol =
node->Attribute(
"symbol");
473 if (
debug && world.
rank() == 0) std::cout <<
" found atomic pseudopotential " << symbol << std::endl;
475 node->Attribute(
"lmax", &lmax);
476 if (
debug && world.
rank() == 0) std::cout <<
" maximum L is " << lmax << std::endl;
503 t_hlij(0, 0, 1) = -1./2.*std::sqrt(3./5.)*
t_hlij(0, 1, 1);
505 t_hlij(0, 0, 2) = 1./2.*std::sqrt(5./21.)*
t_hlij(0, 2, 2);
507 t_hlij(0, 1, 2) = -1./2.*std::sqrt(100./63.)*
t_hlij(0, 2, 2);
510 t_hlij(1, 0, 1) = -1./2.*std::sqrt(5./7.)*
t_hlij(1, 1, 1);
512 t_hlij(1, 0, 2) = 1./6.*std::sqrt(35./11.)*
t_hlij(1, 2, 2);
514 t_hlij(1, 1, 2) = -1./6.*14./std::sqrt(11.)*
t_hlij(1, 2, 2);
517 t_hlij(2, 0, 1) = -1./2.*std::sqrt(7./9.)*
t_hlij(2, 1, 1);
519 t_hlij(2, 0, 2) = 1./2.*std::sqrt(63./143.)*
t_hlij(2, 2, 2);
521 t_hlij(2, 1, 2) = -1./2.*18./std::sqrt(143.)*
t_hlij(2, 2, 2);
548 unsigned int maxLL = 0;
553 for (
unsigned int iatom = 0; iatom <
natoms; iatom++) {
564 for (
unsigned int iatom = 0; iatom <
natoms; iatom++) {
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++) {
598 for (
unsigned int iatom = 0; iatom <
natoms; iatom++) {
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++) {
609 for (
unsigned int j = 1; j <= 3; j++) {
621 double trantol =
vtol2 / std::min(30.0,
double(
localproj.size()));
626 int nocc = occ.
size();
628 for(
int i = 0;i < nocc;++i){
657 for (
unsigned int iatom = 0; iatom <
natoms; iatom++) {
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++) {
669 for (
unsigned int j = 1; j <= 3; j++) {
double potential(const coord_3d &r)
Definition 3dharmonic.cc:132
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:100
static const double & get_thresh()
Returns the default threshold.
Definition funcdefaults.h:176
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:1006
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition mra.h:746
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 final
Override this to return list of special points to be refined more deeply.
Definition gth_pseudopotential.h:346
double t1
Definition gth_pseudopotential.h:82
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 ¢er)
Definition gth_pseudopotential.h:90
double operator()(const coord_3d &r) const final
Definition gth_pseudopotential.h:99
virtual void operator()(const Vector< double *, 3 > &xvals, double *MADNESS_RESTRICT fvals, int npts) const final
Definition gth_pseudopotential.h:220
virtual bool screened(const coord_3d &c1, const coord_3d &c2) const final
Definition gth_pseudopotential.h:173
Level special_level() const final
Override this to change the minimum level of refinement at special points (default is 6)
Definition gth_pseudopotential.h:348
int l
Definition gth_pseudopotential.h:74
double sqrtPI
Definition gth_pseudopotential.h:80
virtual bool supports_vectorized() const final
Does the interface support a vectorized operator()?
Definition gth_pseudopotential.h:86
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 ¢er)
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
Definition gth_pseudopotential.h:40
double C3
Definition gth_pseudopotential.h:42
double C1
Definition gth_pseudopotential.h:42
double C4
Definition gth_pseudopotential.h:42
double operator()(const coord_3d &r) const final
Definition gth_pseudopotential.h:54
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
Level special_level() const final
Override this to change the minimum level of refinement at special points (default is 6)
Definition gth_pseudopotential.h:66
std::vector< coord_3d > special_points() const final
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 ¢er)
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:320
WorldGopInterface & gop
Global operations.
Definition world.h:207
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 final(a, b, c)
Definition lookup3.c:153
#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.
constexpr 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:122
std::vector< real_function_3d > vector_real_function_3d
Definition functypedefs.h:94
unsigned int symbol_to_atomic_number(const std::string &symbol)
Definition atomutil.cc:173
int Level
Definition key.h:57
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:1774
FunctionFactory< double, 3 > real_factory_3d
Definition functypedefs.h:108
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:2431
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
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:246
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
const auto npts
Definition testgconv.cc:52