32 #ifndef MOLECULAR_BASIS_H
33 #define MOLECULAR_BASIS_H
38 #include <madness/external/tinyxml/tinyxml.h>
63 int l_lim = 2*
type - 1;
65 for (
int n=l_lim; n>1; n-=2)
f *= n;
67 for (
int n=0; n<
np; n++)
71 for (
int n1=0; n1<
np; n1++) {
72 for (
int n2=0; n2<
np; n2++) {
80 for (
int n=0; n<
np; n++)
coeff[n] *=
f;
88 const std::vector<double>&
coeff,
89 const std::vector<double>&
expnt,
93 double minexpnt =
expnt[0];
94 for (
unsigned int i=1; i<
expnt.size(); i++)
95 minexpnt = std::min(minexpnt,
expnt[i]);
108 if (rsq >
rsqmax)
return 0.0;
110 for (
unsigned int i=0; i<
coeff.size(); i++) {
111 double ersq =
expnt[i]*rsq;
112 if (ersq < 18.4)
sum +=
coeff[i]*exp(-ersq);
119 double*
eval(
double rsq,
double x,
double y,
double z,
double* bf)
const {
121 if (fabs(
R) < 1
e-8) {
122 for (
int i=0; i<
numbf; i++) bf[i] = 0.0;
157 throw "UNKNOWN ANGULAR MOMENTUM";
191 static const char* tags[4][10] = {
192 {
"s" ,
"" ,
"" ,
"" ,
"" ,
"" ,
"" ,
"" ,
"" ,
"" } ,
193 {
"px" ,
"py" ,
"pz" ,
"" ,
"" ,
"" ,
"" ,
"" ,
"" ,
"" } ,
194 {
"dxx" ,
"dxy" ,
"dxz" ,
"dyy" ,
"dyz" ,
"dzz" ,
"" ,
"" ,
"" ,
"" } ,
195 {
"fxxx",
"fxxy",
"fxxz",
"fxyy",
"fxyz",
"fxzz",
"fxzz",
"fyyy",
"fyzz",
"fzzz"}
198 return tags[
type][ibf];
201 template <
typename Archive>
209 static const char* tag[] = {
"s",
"p",
"d",
"f",
"g"};
212 const std::vector<double>& coeff =
c.get_coeff();
213 const std::vector<double>& expnt =
c.get_expnt();
215 p += sprintf(
p,
"%s [",tag[
c.angular_momentum()]);
216 for (
int i=0; i<
c.nprim(); i++) {
217 p += sprintf(
p,
"%.6f(%.6f)",coeff[i],expnt[i]);
218 if (i != (
c.nprim()-1))
p += sprintf(
p,
", ");
220 p += sprintf(
p,
"]");
228 std::vector<ContractedGaussianShell>
g;
240 for (
unsigned int i=0; i<
g.size(); i++) {
247 const Tensor<double>&
avec,
const Tensor<double>&
bvec) {
264 const std::vector<ContractedGaussianShell>&
get_shells()
const {
273 double*
eval(
double x,
double y,
double z,
double* bf)
const {
274 double rsq = x*x + y*y + z*z;
276 for (
int i=0; i<
numbf; i++) bf[i] = 0.0;
280 double* bfstart = bf;
281 for (
unsigned int i=0; i<
g.size(); i++) {
282 bf =
g[i].eval(rsq, x, y, z, bf);
292 double rsq = x*x + y*y + z*z;
293 if (rsq >
rmaxsq)
return 0.0;
301 for (
int j=0; j<
numbf; ++j)
311 for (
unsigned int i=0; i<
g.size(); i++) {
312 int nbf_in_shell =
g[i].nbf();
313 if (ibf>=n && ibf<(n+nbf_in_shell)) {
314 ibf_in_shell = ibf-n;
340 template <
typename Archive>
348 const std::vector<ContractedGaussianShell>& shells =
c.get_shells();
349 for (
int i=0; i<
c.nshell(); i++) {
350 s <<
" " << shells[i] << std::endl;
352 if (
c.has_guess_info()) {
353 s <<
" " <<
"Guess density matrix" << std::endl;
387 double rsq = x*x + y*y + z*z;
393 s <<
"atomic basis function: center " <<
xx <<
" " <<
yy <<
" " <<
zz <<
" : ibf " <<
ibf <<
" nbf " <<
nbf <<
" : shell " <<
shell << std::endl;
426 std::vector<AtomicBasis>
ag;
428 template <
typename T>
430 TiXmlElement* child = node->FirstChildElement(
name);
432 std::istringstream s(child->GetText());
434 for (
int i=0; i<n; i++) {
440 template <
typename T>
442 TiXmlElement* child = node->FirstChildElement(
name);
444 std::istringstream s(child->GetText());
446 for (
int i=0; i<n; i++) {
447 for (
int j=0; j<
m; j++) {
463 static const bool debug =
false;
465 if (!doc.LoadFile()) {
466 std::cout <<
"AtomicBasisSet: Failed loading from file " <<
filename
467 <<
" : ErrorDesc " << doc.ErrorDesc()
468 <<
" : Row " << doc.ErrorRow()
469 <<
" : Col " << doc.ErrorCol() << std::endl;
472 for (TiXmlElement* node=doc.FirstChildElement(); node; node=node->NextSiblingElement()) {
473 if (strcmp(node->Value(),
"name") == 0) {
474 name = node->GetText();
475 if (
debug) std::cout <<
"Loading basis set " <<
name << std::endl;
477 else if (strcmp(node->Value(),
"basis") == 0) {
478 const char* symbol = node->Attribute(
"symbol");
479 if (
debug) std::cout <<
" found basis set for " << symbol << std::endl;
481 std::vector<ContractedGaussianShell>
g;
482 for (TiXmlElement* shell=node->FirstChildElement(); shell; shell=shell->NextSiblingElement()) {
483 const char*
type = shell->Attribute(
"type");
485 shell->Attribute(
"nprim",&nprim);
486 if (
debug) std::cout <<
" found shell " <<
type <<
" " << nprim << std::endl;
487 std::vector<double> expnt = load_tixml_vector<double>(shell, nprim,
"exponents");
488 if (strcmp(
type,
"L") == 0) {
489 std::vector<double> scoeff = load_tixml_vector<double>(shell, nprim,
"scoefficients");
490 std::vector<double> pcoeff = load_tixml_vector<double>(shell, nprim,
"pcoefficients");
495 static const char* tag[] = {
"S",
"P",
"D",
"F",
"G"};
497 for (i=0; i<5; i++) {
498 if (strcmp(
type,tag[i]) == 0)
goto foundit;
502 std::vector<double> coeff = load_tixml_vector<double>(shell, nprim,
"coefficients");
508 else if (strcmp(node->Value(),
"atomicguess") == 0) {
509 const char* symbol = node->Attribute(
"symbol");
510 if (
debug) std::cout <<
" atomic guess info for " << symbol << std::endl;
513 int nbf =
ag[atn].nbf();
514 Tensor<double> dmat = load_tixml_matrix<double>(node,
nbf,
nbf,
"guessdensitymatrix");
515 Tensor<double> avec = load_tixml_matrix<double>(node,
nbf,
nbf,
"alphavectors");
516 Tensor<double> bvec = load_tixml_matrix<double>(node,
nbf,
nbf,
"betavectors");
517 ag[atn].set_guess_info(dmat, avec, bvec);
533 for (
int i=0; i<mentity.
natom(); i++) {
548 for (
int i=0; i<mentity.
natom(); i++) {
553 const int nbf_on_atom =
ag[atn].nbf();
554 if (ibf >= n && (n+nbf_on_atom) > ibf) {
568 for (
int i=0; i<mentity.
natom(); i++) {
573 const int nbf_on_atom =
ag[atn].nbf();
574 if (ibf >= n && (n+nbf_on_atom) > ibf) {
577 ag[atn].get_shell_from_basis_function(ibf-n, index);
591 for (
int i=0; i<mentity.
natom(); i++) {
602 for (
int i=0; i<mentity.
natom(); i++) {
605 bf =
ag[atn].eval(x-atom.
x, y-atom.
y, z-atom.
z, bf);
613 for (
int i=0; i<mentity.
natom(); i++) {
616 sum +=
ag[atn].eval_guess_density(x-atom.
x, y-atom.
y, z-atom.
z);
622 return ag[atomic_number].nbf() > 0;
627 std::cout <<
"\n " <<
name <<
" atomic basis set" << std::endl;
628 for (
int i=0; i<mentity.
natom(); i++) {
631 for (
int j=0; j<i; j++) {
635 std::cout << std::endl;
637 std::cout <<
ag[atn];
643 template <
typename T>
661 template <
typename T>
663 const double thresh = 0.2*
v.normf();
665 printf(
" zero vector\n");
668 long nbf = int(
v.dim[0]);
671 for (
long i=0; i<
nbf; i++) {
679 if (mentity.
natom() < 10) {
680 format =
" %2s(%1d)%4s(%2ld)%6.3f ";
682 else if (mentity.
natom() < 100) {
683 format =
" %2s(%2d)%4s(%3ld)%6.3f ";
685 else if (mentity.
natom() < 1000) {
686 format =
" %2s(%3d)%4s(%4ld)%6.3f ";
689 format =
" %2s(%4d)%4s(%5ld)%6.3f ";
691 for (
long ii=0; ii<ngot; ii++) {
701 printf(format, element, iat, desc, ibf,
v[ibf]);
708 std::cout <<
"\n " <<
name <<
" atomic basis set" << std::endl;
709 for (
unsigned int i=0; i<
ag.size(); i++) {
710 if (
ag[i].
nbf() > 0) {
717 template <
typename Archive>
std::ostream & operator<<(std::ostream &s, const ContractedGaussianShell &c)
Definition: apps/periodic_old/molecularbasis.h:208
double x
Definition: mentity.h:73
unsigned int atomic_number
Atomic number.
Definition: mentity.h:74
double y
Definition: mentity.h:73
double z
Definition: mentity.h:73
Used to represent one basis function from a shell on a specific center.
Definition: apps/periodic_old/molecularbasis.h:361
void print_me(std::ostream &s) const
Definition: apps/periodic_old/molecularbasis.h:392
const int ibf
Definition: apps/periodic_old/molecularbasis.h:365
const double yy
Definition: apps/periodic_old/molecularbasis.h:363
const char * get_desc() const
Definition: apps/periodic_old/molecularbasis.h:404
const int nbf
Definition: apps/periodic_old/molecularbasis.h:366
const double xx
Definition: apps/periodic_old/molecularbasis.h:363
const double zz
Definition: apps/periodic_old/molecularbasis.h:363
AtomicBasisFunction(const AtomicBasisFunction &aofunc)
Definition: apps/periodic_old/molecularbasis.h:374
int get_index() const
Definition: apps/periodic_old/molecularbasis.h:400
const ContractedGaussianShell & shell
Definition: apps/periodic_old/molecularbasis.h:364
double operator()(double x, double y, double z) const
Definition: apps/periodic_old/molecularbasis.h:382
void get_coords(double &x, double &y, double &z) const
Definition: apps/periodic_old/molecularbasis.h:408
const ContractedGaussianShell & get_shell() const
Definition: apps/periodic_old/molecularbasis.h:396
AtomicBasisFunction(double x, double y, double z, const ContractedGaussianShell &shell, int ibf)
Definition: apps/periodic_old/molecularbasis.h:369
double rangesq() const
Definition: apps/periodic_old/molecularbasis.h:413
Definition: apps/periodic_old/molecularbasis.h:644
AnalysisSorter(const Tensor< T > &v)
Definition: apps/periodic_old/molecularbasis.h:647
bool operator()(long i, long j) const
Definition: apps/periodic_old/molecularbasis.h:648
const Tensor< T > v
Definition: apps/periodic_old/molecularbasis.h:645
Contracted Gaussian basis.
Definition: apps/periodic_old/molecularbasis.h:424
double eval_guess_density(const MolecularEntity &mentity, double x, double y, double z) const
Evaluates the guess density.
Definition: apps/periodic_old/molecularbasis.h:611
AtomicBasisSet()
Definition: apps/periodic_old/molecularbasis.h:455
void read_file(std::string filename)
Definition: apps/periodic_old/molecularbasis.h:462
int basisfn_to_atom(const MolecularEntity &mentity, int ibf) const
Returns the number of the atom the ibf'th basis function is on.
Definition: apps/periodic_old/molecularbasis.h:545
Tensor< T > load_tixml_matrix(TiXmlElement *node, int n, int m, const char *name)
Definition: apps/periodic_old/molecularbasis.h:441
int nbf(const MolecularEntity &mentity) const
Given a molecule count the number of basis functions.
Definition: apps/periodic_old/molecularbasis.h:589
void eval(const MolecularEntity &mentity, double x, double y, double z, double *bf) const
Evaluates the basis functions.
Definition: apps/periodic_old/molecularbasis.h:601
void print(const MolecularEntity &mentity) const
Print basis info for atoms in the molecule (once for each unique atom type)
Definition: apps/periodic_old/molecularbasis.h:626
void atoms_to_bfn(const MolecularEntity &mentity, std::vector< int > &at_to_bf, std::vector< int > &at_nbf)
Makes map from atoms to first basis function on atom and number of basis functions on atom.
Definition: apps/periodic_old/molecularbasis.h:528
std::vector< T > load_tixml_vector(TiXmlElement *node, int n, const char *name)
Definition: apps/periodic_old/molecularbasis.h:429
AtomicBasisSet(std::string filename)
Definition: apps/periodic_old/molecularbasis.h:458
std::vector< AtomicBasis > ag
Definition: apps/periodic_old/molecularbasis.h:426
void print_anal(const MolecularEntity &mentity, const Tensor< T > &v)
Given a vector of AO coefficients prints an analysis.
Definition: apps/periodic_old/molecularbasis.h:662
void print_all() const
Print basis info for all supported atoms.
Definition: apps/periodic_old/molecularbasis.h:707
AtomicBasisFunction get_atomic_basis_function(const MolecularEntity &mentity, int ibf) const
Returns the ibf'th atomic basis function.
Definition: apps/periodic_old/molecularbasis.h:565
void serialize(Archive &ar)
Definition: apps/periodic_old/molecularbasis.h:718
bool is_supported(int atomic_number) const
Definition: apps/periodic_old/molecularbasis.h:621
std::string name
Definition: apps/periodic_old/molecularbasis.h:425
Represents multiple shells of contracted gaussians on a single center.
Definition: apps/periodic_old/molecularbasis.h:227
const ContractedGaussianShell & get_shell_from_basis_function(int ibf, int &ibf_in_shell) const
Return shell that contains basis function ibf and also return index of function in the shell.
Definition: apps/periodic_old/molecularbasis.h:309
std::vector< ContractedGaussianShell > g
Definition: apps/periodic_old/molecularbasis.h:228
double eval_guess_density(double x, double y, double z) const
Evaluates the guess atomic density at point x, y, z relative to atomic center.
Definition: apps/periodic_old/molecularbasis.h:290
AtomicBasis()
Definition: apps/periodic_old/molecularbasis.h:234
int numbf
Definition: apps/periodic_old/molecularbasis.h:230
const Tensor< double > & get_dmat() const
Definition: apps/periodic_old/molecularbasis.h:328
const Tensor< double > & get_bvec() const
Definition: apps/periodic_old/molecularbasis.h:336
Tensor< double > bvec
Definition: apps/periodic_old/molecularbasis.h:231
void set_guess_info(const Tensor< double > &dmat, const Tensor< double > &avec, const Tensor< double > &bvec)
Definition: apps/periodic_old/molecularbasis.h:246
AtomicBasis(const std::vector< ContractedGaussianShell > &g)
Definition: apps/periodic_old/molecularbasis.h:236
Tensor< double > dmat
Definition: apps/periodic_old/molecularbasis.h:231
const Tensor< double > & get_avec() const
Definition: apps/periodic_old/molecularbasis.h:332
void serialize(Archive &ar)
Definition: apps/periodic_old/molecularbasis.h:341
double * eval(double x, double y, double z, double *bf) const
Evaluates the basis functions at point x, y, z relative to atomic center.
Definition: apps/periodic_old/molecularbasis.h:273
bool has_guess_info() const
Definition: apps/periodic_old/molecularbasis.h:324
Tensor< double > avec
Definition: apps/periodic_old/molecularbasis.h:231
const std::vector< ContractedGaussianShell > & get_shells() const
Returns a const reference to the shells.
Definition: apps/periodic_old/molecularbasis.h:264
double rmaxsq
Definition: apps/periodic_old/molecularbasis.h:229
int nshell() const
Returns the number of shells on the center.
Definition: apps/periodic_old/molecularbasis.h:259
int nbf() const
Returns the number of basis functions on the center.
Definition: apps/periodic_old/molecularbasis.h:254
Represents a single shell of contracted, Cartesian, Gaussian primitives.
Definition: apps/periodic_old/molecularbasis.h:49
int type
Angular momentum = 0, 1, 2, ...
Definition: apps/periodic_old/molecularbasis.h:50
ContractedGaussianShell()
Definition: apps/periodic_old/molecularbasis.h:84
std::vector< double > expnt
Definition: apps/periodic_old/molecularbasis.h:52
const char * get_desc(int ibf) const
Returns a string description of the basis function type.
Definition: apps/periodic_old/molecularbasis.h:190
double rsqmax
Definition: apps/periodic_old/molecularbasis.h:53
double eval_radial(double rsq) const
Evaluates the radial part of the contracted function.
Definition: apps/periodic_old/molecularbasis.h:107
double rangesq() const
Returns square of the distance beyond which function is less than 1e-8.
Definition: apps/periodic_old/molecularbasis.h:101
ContractedGaussianShell(int type, const std::vector< double > &coeff, const std::vector< double > &expnt, bool donorm=true)
Definition: apps/periodic_old/molecularbasis.h:87
double * eval(double rsq, double x, double y, double z, double *bf) const
Evaluates the entire shell returning the incremented result pointer.
Definition: apps/periodic_old/molecularbasis.h:119
const std::vector< double > & get_coeff() const
Returns a const reference to the coefficients.
Definition: apps/periodic_old/molecularbasis.h:180
void normalize()
Definition: apps/periodic_old/molecularbasis.h:56
int nprim() const
Returns the number of primitives in the contraction.
Definition: apps/periodic_old/molecularbasis.h:175
void serialize(Archive &ar)
Definition: apps/periodic_old/molecularbasis.h:202
std::vector< double > coeff
Definition: apps/periodic_old/molecularbasis.h:51
const std::vector< double > & get_expnt() const
Returns a const reference to the exponents.
Definition: apps/periodic_old/molecularbasis.h:185
int angular_momentum() const
Returns the shell angular momentum.
Definition: apps/periodic_old/molecularbasis.h:165
int nbf() const
Returns the number of basis functions in the shell.
Definition: apps/periodic_old/molecularbasis.h:170
int numbf
Number of basis functions in shell (type+1)*(type+2)/2.
Definition: apps/periodic_old/molecularbasis.h:54
int natom() const
Definition: mentity.h:112
const Atom & get_atom(unsigned int i) const
Definition: mentity.cc:369
long size() const
Returns the number of elements in the tensor.
Definition: basetensor.h:138
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
Defines common mathematical and physical constants.
static const double R
Definition: csqrt.cc:46
double(* f)(const coord_3d &)
Definition: derivatives.cc:54
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition: derivatives.cc:72
static bool debug
Definition: dirac-hatom.cc:16
Fcwf copy(Fcwf psi)
Definition: fcwf.cc:338
const double m
Definition: gfit.cc:199
static const double v
Definition: hatom_sf_dirac.cc:20
static double pow(const double *a, const double *b)
Definition: lda.h:74
#define max(a, b)
Definition: lda.h:51
Macros and tools pertaining to the configuration of MADNESS.
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
unsigned int symbol_to_atomic_number(const std::string &symbol)
Definition: mentity.cc:181
const AtomicData & get_atomic_data(unsigned int atomic_number)
Definition: mentity.cc:175
const double pi
Mathematical constant .
Definition: constants.h:48
static const char * filename
Definition: legendre.cc:96
std::string type(const PairType &n)
Definition: PNOParameters.h:18
static long abs(long a)
Definition: tensor.h:218
static const double a
Definition: nonlinschro.cc:118
static std::vector< int > at_to_bf
Definition: preal.cc:19
static std::vector< int > at_nbf
Definition: preal.cc:19
static const double c
Definition: relops.cc:10
static const double thresh
Definition: rk.cc:45
const char *const symbol
Definition: mentity.h:54
int np
Definition: tdse1d.cc:165
Defines and implements most of Tensor.
AtomicInt sum
Definition: test_atomicint.cc:46
void e()
Definition: test_sig.cc:75
double g(const coord_1d &r)
Definition: testgconv.cc:49