MADNESS 0.10.1
madness/chem/molecularbasis.h
Go to the documentation of this file.
1/*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30*/
31
32#ifndef MADNESS_CHEM_MOLECULAR_BASIS_H__INCLUDED
33#define MADNESS_CHEM_MOLECULAR_BASIS_H__INCLUDED
34
36#include <madness/constants.h>
39#include <madness/external/tinyxml/tinyxml.h>
41
42#include <vector>
43#include <algorithm>
44#include <iostream>
45#include <sstream>
46#include <iomanip>
47#include <cstdio>
48
49namespace madness {
50/// Represents a single shell of contracted, Cartesian, Gaussian primitives
52public:
53 static const int maxtype=6; ///< Maximum angular momentum supported
54 static const int maxbf=(maxtype+1)*(maxtype+2)/2; ///< Maximum number of basis functions in a shell
55private:
56 int type; ///< Angular momentum = 0, 1, 2, ...
57 std::vector<double> coeff;
58 std::vector<double> expnt;
59 double rsqmax;
60 int numbf; ///< Number of basis functions in shell (type+1)*(type+2)/2
61
62 void normalize() {
63 // nwcchem cartesian normalization conventions
64 // translation of nmcoeff.F into python and thence to c++
65 int np = coeff.size();
66 if (np == 1) coeff[0] = 1.0e0;
67
68 double pi32=pow(madness::constants::pi,1.5);
69 int l_lim = 2*type - 1;
70 double f = 1.0e00;
71 for (int n=l_lim; n>1; n-=2) f *= n;
72
73 for (int n=0; n<np; ++n)
74 coeff[n] *= pow(2.e0*expnt[n]/madness::constants::pi,0.75e0)*pow(4.e0*expnt[n],0.5E0*type)/sqrt(f);
75
76 double sum = 0.0;
77 for (int n1=0; n1<np; ++n1) {
78 for (int n2=0; n2<np; ++n2) {
79 double S =pi32/pow(expnt[n1]+expnt[n2],1.5e0+type)/pow(2e0,type);
80 sum = sum + coeff[n1]*coeff[n2]*S;
81 }
82 }
83 sum *= f;
84
85 f = 1e0/sqrt(sum);
86 for (int n=0; n<np; ++n) coeff[n] *= f;
87 }
88
89public:
91 : type(-1), coeff(), expnt(), rsqmax(0.0), numbf(0) {};
92
94 const std::vector<double>& coeff,
95 const std::vector<double>& expnt,
96 bool donorm=true)
97 : type(type), coeff(coeff), expnt(expnt), numbf((type+1)*(type+2)/2) {
98 if (donorm) normalize();
99 double minexpnt = expnt[0];
100 for (unsigned int i=1; i<expnt.size(); ++i)
101 minexpnt = std::min(minexpnt,expnt[i]);
102 rsqmax = 27.6/minexpnt; // 27.6 = log(1e12)
103 }
104
105
106 /// Returns square of the distance beyond which function is less than 1e-8.
107 double rangesq() const {
108 return rsqmax;
109 }
110
111
112 /// Evaluates the radial part of the contracted function
113 double eval_radial(double rsq) const {
114 if (rsq > rsqmax) return 0.0;
115 double sum = 0.0;
116 for (unsigned int i=0; i<coeff.size(); ++i) {
117 double ersq = expnt[i]*rsq;
118 if (ersq < 27.6) sum += coeff[i]*exp(-ersq); // 27.6 = log(1e12)
119 }
120 return sum;
121 }
122
123
124 /// Evaluates the entire shell returning the incremented result pointer
125 double* eval(double rsq, double x, double y, double z, double* bf) const {
126 double R = eval_radial(rsq);
127 if (fabs(R) < 1e-12) {
128 for (int i=0; i<numbf; ++i) bf[i] = 0.0;
129
130 }
131 else {
132 switch (type) {
133 case 0:
134 bf[0] = R;
135 break;
136 case 1:
137 bf[0] = R*x;
138 bf[1] = R*y;
139 bf[2] = R*z;
140 break;
141 case 2:
142 { // braces need by some compilers to limit scope of fac
143 static const double fac = 1.0; //sqrt(3.0);
144 bf[0] = R*x*x;
145 bf[1] = R*x*y*fac;
146 bf[2] = R*x*z*fac;
147 bf[3] = R*y*y;
148 bf[4] = R*y*z*fac;
149 bf[5] = R*z*z;
150 }
151 break;
152 case 3:
153 bf[0] = R*x*x*x;
154 bf[1] = R*x*x*y;
155 bf[2] = R*x*x*z;
156 bf[3] = R*x*y*y;
157 bf[4] = R*x*y*z;
158 bf[5] = R*x*z*z;
159 bf[6] = R*y*y*y;
160 bf[7] = R*y*y*z;
161 bf[8] = R*y*z*z;
162 bf[9] = R*z*z*z;
163 break;
164
165 default:
166 throw "UNKNOWN ANGULAR MOMENTUM";
167 }
168 }
169 return bf+numbf;
170 }
171
172
173 /// Returns the shell angular momentum
174 int angular_momentum() const {
175 return type;
176 }
177
178 /// Returns the number of basis functions in the shell
179 int nbf() const {
180 return numbf;
181 }
182
183 /// Returns the number of primitives in the contraction
184 int nprim() const {
185 return coeff.size();
186 }
187
188 /// Returns a const reference to the coefficients
189 const std::vector<double>& get_coeff() const {
190 return coeff;
191 }
192
193 /// Returns a const reference to the exponents
194 const std::vector<double>& get_expnt() const {
195 return expnt;
196 }
197
198 /// Returns a string description of the basis function type
199 const char* get_desc(int ibf) const {
200 static const char* tags[4][10] = {
201 {"s" ,"" ,"" ,"" ,"" ,"" ,"" ,"" ,"" ,"" } ,
202 {"px" ,"py" ,"pz" ,"" ,"" ,"" ,"" ,"" ,"" ,"" } ,
203 {"dxx" ,"dxy" ,"dxz" ,"dyy" ,"dyz" ,"dzz" ,"" ,"" ,"" ,"" } ,
204 {"fxxx","fxxy","fxxz","fxyy","fxyz","fxzz","fxzz","fyyy","fyzz","fzzz"}
205 };
206 MADNESS_ASSERT(ibf<numbf && ibf >= 0);
207 return tags[type][ibf];
208 }
209
210 template <typename Archive>
211 void serialize(Archive& ar) {
212 ar & type & coeff & expnt & rsqmax & numbf;
213 }
214};
215
216/// Represents multiple shells of contracted gaussians on a single center
218 std::vector<ContractedGaussianShell> g;
219 double rmaxsq;
220 int numbf;
222
223public:
224 AtomicBasis() : g(), rmaxsq(0.0), numbf(0) {};
225
226 AtomicBasis(const std::vector<ContractedGaussianShell>& g)
227 : g(g) {
228 rmaxsq = 0.0;
229 numbf = 0;
230 for (unsigned int i=0; i<g.size(); ++i) {
231 rmaxsq = std::max(rmaxsq, g[i].rangesq());
232 numbf += g[i].nbf();
233 }
234 }
235
237 const Tensor<double>& avec, const Tensor<double>& bvec,
238 const Tensor<double>& aocc, const Tensor<double>& bocc,
239 const Tensor<double>& aeps, const Tensor<double>& beps,
241 this->dmat = copy(dmat);
242 this->dmatpsp = copy(dmatpsp);
243 this->avec = copy(avec);
244 this->bvec = copy(bvec);
245 this->aocc = copy(aocc);
246 this->beps = copy(beps);
247 this->aeps = copy(aeps);
248 this->bocc = copy(bocc);
249 this->aoccpsp = copy(aoccpsp);
250 this->boccpsp = copy(boccpsp);
251 }
252
253 /// Returns the number of basis functions on the center
254 int nbf() const {
255 return numbf;
256 }
257
258 /// Returns the number of shells on the center
259 int nshell() const {
260 return g.size();
261 }
262
263 /// Returns a const reference to the shells
264 const std::vector<ContractedGaussianShell>& get_shells() const {
265 return g;
266 };
267
268 /// Evaluates the basis functions at point x, y, z relative to atomic center
269
270 /// The array bf[] must be large enough to hold nbf() values.
271 ///
272 /// Returned is the incremented pointer.
273 double* eval(double x, double y, double z, double* bf) const {
274 double rsq = x*x + y*y + z*z;
275 if (rsq > rmaxsq) {
276 for (int i=0; i<numbf; ++i) bf[i] = 0.0;
277 return bf+numbf;
278 }
279
280 double* bfstart = bf;
281 for (unsigned int i=0; i<g.size(); ++i) {
282 bf = g[i].eval(rsq, x, y, z, bf);
283 }
284 // paranoia is good
285 MADNESS_CHECK(bf-bfstart == numbf);
286 return bf;
287 }
288
289 /// Evaluates the guess atomic density at point x, y, z relative to atomic center
290 double eval_guess_density(double x, double y, double z, bool pspat) const {
292 double rsq = x*x + y*y + z*z;
293 if (rsq > rmaxsq) return 0.0;
294
296 eval(x, y, z, bf);
297 const double* p;
298 // check if pseudo-atom
299 if (pspat){
300 p = dmatpsp.ptr();}
301 else{
302 p = dmat.ptr();}
303 double sum = 0.0;
304 for (int i=0; i<numbf; ++i, p+=numbf) {
305 double sumj = 0.0;
306 for (int j=0; j<numbf; ++j)
307 sumj += p[j]*bf[j];
308 sum += bf[i]*sumj;
309 }
310 return sum;
311 }
312
313 /// Return shell that contains basis function ibf and also return index of function in the shell
314 const ContractedGaussianShell& get_shell_from_basis_function(int ibf, int& ibf_in_shell) const {
315 int n=0;
316 for (unsigned int i=0; i<g.size(); ++i) {
317 int nbf_in_shell = g[i].nbf();
318 if (ibf>=n && ibf<(n+nbf_in_shell)) {
319 ibf_in_shell = ibf-n;
320 return g[i];
321 }
322 else {
323 n += g[i].nbf();
324 }
325 }
326 MADNESS_EXCEPTION("AtomicBasis: get_shell_from_basis_function", ibf*100000 + nbf());
327 }
328
329 bool has_guess_info() const {
330 return dmat.size()>0;
331 }
332
333 const Tensor<double>& get_dmat() const {
334 return dmat;
335 };
336
338 dmat = mat;
339 };
340
341 bool has_guesspsp_info() const {
342 return dmatpsp.size()>0;
343 }
344
346 return dmatpsp;
347 };
348
350 dmatpsp = mat;
351 };
352
353 const Tensor<double>& get_avec() const {
354 return avec;
355 };
356
357 const Tensor<double>& get_bvec() const {
358 return bvec;
359 };
360
361 const Tensor<double>& get_aeps() const {
362 return aeps;
363 };
364
365 const Tensor<double>& get_beps() const {
366 return beps;
367 };
368
369 const Tensor<double>& get_aocc() const {
370 return aocc;
371 };
372
373 const Tensor<double>& get_bocc() const {
374 return bocc;
375 };
376
378 aocc = occ;
379 };
380
382 bocc = occ;
383 };
384
386 return aoccpsp;
387 };
388
390 return boccpsp;
391 };
392
394 aoccpsp = occ;
395 };
396
398 boccpsp = occ;
399 };
400
401 template <typename Archive>
402 void serialize(Archive& ar) {
403 ar & g & rmaxsq & numbf & dmat & dmatpsp & avec & bvec & aocc & bocc & aeps & beps & aoccpsp & boccpsp;
404 }
405
406};
407
408/// Used to represent one basis function from a shell on a specific center
410private:
411 const double xx, yy, zz; // Coordinates of the center
412 const ContractedGaussianShell& shell; // Reference to the underlying atomic shell
413 const int ibf; // Index of basis function in the shell (0, 1, ...)
414 const int nbf; // Number of functions in the shell
415
416public:
417 AtomicBasisFunction(double x, double y, double z,
419 : xx(x), yy(y), zz(z), shell(shell), ibf(ibf), nbf(shell.maxbf) {}
420
421
423 : xx(aofunc.xx)
424 , yy(aofunc.yy)
425 , zz(aofunc.zz)
426 , shell(aofunc.shell)
427 , ibf(aofunc.ibf)
428 , nbf(aofunc.nbf) {}
429
430 double operator()(double x, double y, double z) const {
432 x-=xx;
433 y-=yy;
434 z-=zz;
435 double rsq = x*x + y*y + z*z;
436 shell.eval(rsq, x, y, z, bf);
437 return bf[ibf];
438 }
439
440 void print_me(std::ostream& s) const;
441
443 return shell;
444 }
445
446 int get_index() const {
447 return ibf;
448 }
449
450 const char* get_desc() const {
451 return shell.get_desc(ibf);
452 }
453
454 void get_coords(double& x, double& y, double& z) const {
455 x=xx; y=yy; z=zz;
456 return;
457 }
458
462
463 double rangesq() const {
464 return shell.rangesq();
465 }
466};
467
468/// Contracted Gaussian basis
470 std::string name;
471 std::vector<AtomicBasis> ag; ///< Basis associated by atomic number = 1, 2, ...; 0=Bq.
472
473 template <typename T>
474 std::vector<T> load_tixml_vector(TiXmlElement* node, int n, const char* name) {
475 TiXmlElement* child = node->FirstChildElement(name);
476 MADNESS_ASSERT(child);
477 std::istringstream s(child->GetText());
478 std::vector<T> r(n);
479 for (int i=0; i<n; ++i) {
480 s >> r[i];
482 }
483 return r;
484 }
485
486 template <typename T>
487 Tensor<T> load_tixml_matrix(TiXmlElement* node, int n, int m, const char* name) {
488 TiXmlElement* child = node->FirstChildElement(name);
489 MADNESS_ASSERT(child);
490 std::istringstream s(child->GetText());
491 Tensor<T> r(n,m);
492 for (int i=0; i<n; ++i) {
493 for (int j=0; j<m; ++j) {
494 s >> r(i,j);
496 }
497 }
498 return r;
499 }
500
501public:
502 AtomicBasisSet() : name("unknown"), ag(110) {}
503
504
505 AtomicBasisSet(std::string filename) : name(""), ag(110) {
507 }
508
509 std::string get_name() const {return name;}
510
511 /// read the atomic basis set from file
512
513 /// use the default location MRA_CHEMDATA_DIR as defined in the Makefile.am
514 /// unless it is overridden by the environment variable MRA_CHEMDATA_DIR
515 /// @param[in] filename the name of the basis set (sto-3g, 6-31g, etc)
516 void read_file(std::string filename);
517
518 /// read the atomic basis set from file
519
520 /// @param[in] filename the base name of nwchem files (.out and .movecs)
521 void read_nw_file(std::string filename);
522
523
524
525 /// Makes map from atoms to first basis function on atom and number of basis functions on atom
526 void atoms_to_bfn(const Molecule& molecule, std::vector<int>& at_to_bf, std::vector<int>& at_nbf) const {
527 at_to_bf = std::vector<int>(molecule.natom());
528 at_nbf = std::vector<int>(molecule.natom());
529
530 int n = 0;
531 for (size_t i=0; i<molecule.natom(); ++i) {
532 const Atom& atom = molecule.get_atom(i);
533 const int atn = atom.atomic_number;
535 at_to_bf[i] = n;
536 at_nbf[i] = ag[atn].nbf();
537 n += at_nbf[i];
538 }
539 }
540
541 /// Makes map from shells to first basis function on she and number of basis functions on sh
542 void shells_to_bfn(const Molecule& molecule, std::vector<int>& sh_to_bf, std::vector<int>& sh_nbf) const {
543 sh_to_bf = std::vector<int>();
544 sh_nbf = std::vector<int>();
545
546 int nbf = 0;
547 for (size_t i=0; i<molecule.natom(); ++i) {
548 const Atom& atom = molecule.get_atom(i);
549 const int atn = atom.atomic_number;
551 const auto& shells = ag[atn].get_shells();
552 for (const auto& sh : shells) {
553 int n = sh.nbf();
554 sh_nbf.push_back(n);
555 sh_to_bf.push_back(nbf);
556 nbf += n;
557 }
558 }
559 }
560
561 /// Returns the atomic alpha eigenvectors for atom iat
562 const Tensor<double>& get_avec(const Molecule& molecule, size_t iat) const {
563 MADNESS_ASSERT(iat>=0 && iat<molecule.natom());
564 const Atom& atom = molecule.get_atom(iat);
565 const int atn = atom.atomic_number;
567 return ag[atn].get_avec();
568 }
569
570 /// Returns the atomic alpha eigenvalues for atom iat
571 const Tensor<double>& get_aeps(const Molecule& molecule, size_t iat) const {
572 MADNESS_ASSERT(iat>=0 && iat<molecule.natom());
573 const Atom& atom = molecule.get_atom(iat);
574 const int atn = atom.atomic_number;
576 return ag[atn].get_aeps();
577 }
578
579 /// Returns the number of the atom the ibf'th basis function is on
580 int basisfn_to_atom(const Molecule& molecule, size_t ibf) const {
581 MADNESS_ASSERT(ibf >= 0);
582 size_t n = 0;
583 for (size_t i=0; i<molecule.natom(); ++i) {
584 // Is the desired function on this atom?
585 const Atom& atom = molecule.get_atom(i);
586 const int atn = atom.atomic_number;
588 const int nbf_on_atom = ag[atn].nbf();
589 if (ibf >= n && (n+nbf_on_atom) > ibf) {
590 return i;
591 }
592 else {
593 n += nbf_on_atom;
594 }
595 }
596 MADNESS_EXCEPTION("AtomicBasisSet: get_atomic_basis_function: confused?", ibf);
597 }
598
599 /// Returns the ibf'th atomic basis function
601 MADNESS_ASSERT(ibf >= 0);
602 size_t n = 0;
603 for (size_t i=0; i<molecule.natom(); ++i) {
604 // Is the desired function on this atom?
605 const Atom& atom = molecule.get_atom(i);
606 const int atn = atom.atomic_number;
608 const int nbf_on_atom = ag[atn].nbf();
609 if (ibf >= n && (n+nbf_on_atom) > ibf) {
610 int index;
611 const ContractedGaussianShell& shell =
612 ag[atn].get_shell_from_basis_function(ibf-n, index);
613 return AtomicBasisFunction(atom.x, atom.y, atom.z, shell, index);
614 }
615 else {
616 n += nbf_on_atom;
617 }
618 }
619 MADNESS_EXCEPTION("AtomicBasisSet: get_atomic_basis_function: confused?", ibf);
620 }
621
622
623 /// Given a molecule count the number of basis functions
624 int nbf(const Molecule& molecule) const {
625 int n = 0;
626 for (size_t i=0; i<molecule.natom(); ++i) {
627 const Atom& atom = molecule.get_atom(i);
628 const int atn = atom.atomic_number;
630 n += ag[atn].nbf();
631 }
632 return n;
633 }
634
635 /// Evaluates the basis functions
636 void eval(const Molecule& molecule, double x, double y, double z, double *bf) const {
637 for (size_t i=0; i<molecule.natom(); ++i) {
638 const Atom& atom = molecule.get_atom(i);
639 const int atn = atom.atomic_number;
640 bf = ag[atn].eval(x-atom.x, y-atom.y, z-atom.z, bf);
641 }
642 }
643
644
645 /// Evaluates the guess density
646 double eval_guess_density(const Molecule& molecule, double x, double y, double z) const {
647 double sum = 0.0;
648 bool pspat;
649 for (size_t i=0; i<molecule.natom(); ++i) {
650 const Atom& atom = molecule.get_atom(i);
651 if (atom.pseudo_atom){
652 pspat=true;}
653 else{
654 pspat=false;}
655 const int atn = atom.atomic_number;
656 sum += ag[atn].eval_guess_density(x-atom.x, y-atom.y, z-atom.z, pspat);
657 }
658 return sum;
659 }
660
661 bool is_supported(int atomic_number) const {
662 return ag[atomic_number].nbf() > 0;
663 }
664
665 /// Print basis info for atoms in the molecule (once for each unique atom type)
666 void print(const Molecule& molecule) const;
667
668 /// Eliminates core orbitals from the density matrix for pseudopotential calculations
669 void modify_dmat_psp(int atn, double zeff);
670
671 template <typename T>
674 public:
676 bool operator()(long i, long j) const {
677 return std::abs(v[i]) > std::abs(v[j]);
678 }
679 };
680
681 /// Given a vector of AO coefficients prints an analysis
682
683 /// For each significant coeff it prints
684 /// - atomic symbol
685 /// - atom number
686 /// - basis function type (e.g., dxy)
687 /// - basis function number
688 /// - MO coeff
689 template <typename T>
690 void print_anal(const Molecule& molecule, const Tensor<T>& v) const {
691 const double thresh = 0.2*v.normf();
692 if (thresh == 0.0) {
693 printf(" zero vector\n");
694 return;
695 }
696 long nbf = int(v.dim(0));
697 std::vector<long> list(nbf);
698 long ngot=0;
699 for (long i=0; i<nbf; ++i) {
700 if (std::abs(v(i)) > thresh) {
701 list[ngot++] = i;
702 }
703 }
704 std::sort(list.begin(),list.begin()+ngot,AnalysisSorter<T>(v));
705
706 const char* format;
707 if (molecule.natom() < 10) {
708 format = " %2s(%1d)%4s(%2ld)%6.3f ";
709 }
710 else if (molecule.natom() < 100) {
711 format = " %2s(%2d)%4s(%3ld)%6.3f ";
712 }
713 else if (molecule.natom() < 1000) {
714 format = " %2s(%3d)%4s(%4ld)%6.3f ";
715 }
716 else {
717 format = " %2s(%4d)%4s(%5ld)%6.3f ";
718 }
719 printf(" ");
720 for (long ii=0; ii<ngot; ++ii) {
721 long ibf = list[ii];
722
723 const int iat = basisfn_to_atom(molecule, ibf);
724 const Atom& atom = molecule.get_atom(iat);
726 const char* desc = ao.get_desc();
727 const char* element = get_atomic_data(atom.atomic_number).symbol;
728
729 // This will need wrapping in a template for a complex MO vector
730 printf(format, element, iat, desc, ibf, v[ibf]);
731 }
732 printf("\n");
733 }
734
735 /// Print basis info for all supported atoms
736 void print_all() const;
737
738 template <typename Archive>
739 void serialize(Archive& ar) {
740 ar & name & ag;
741 }
742};
743}
744
745
746#endif
Declaration of utility class and functions for atom.
Definition molecule.h:60
double y
Definition molecule.h:62
unsigned int atomic_number
Atomic number.
Definition molecule.h:63
double x
Definition molecule.h:62
double z
Definition molecule.h:62
bool pseudo_atom
Indicates if this atom uses a pseudopotential.
Definition molecule.h:65
Used to represent one basis function from a shell on a specific center.
Definition madness/chem/molecularbasis.h:409
void print_me(std::ostream &s) const
Definition molecularbasis.cc:85
int get_index() const
Definition madness/chem/molecularbasis.h:446
const double yy
Definition madness/chem/molecularbasis.h:411
const ContractedGaussianShell & get_shell() const
Definition madness/chem/molecularbasis.h:442
AtomicBasisFunction(double x, double y, double z, const ContractedGaussianShell &shell, int ibf)
Definition madness/chem/molecularbasis.h:417
madness::Vector< double, 3 > get_coords_vec() const
Definition madness/chem/molecularbasis.h:459
const ContractedGaussianShell & shell
Definition madness/chem/molecularbasis.h:412
void get_coords(double &x, double &y, double &z) const
Definition madness/chem/molecularbasis.h:454
const double xx
Definition madness/chem/molecularbasis.h:411
const int nbf
Definition madness/chem/molecularbasis.h:414
double operator()(double x, double y, double z) const
Definition madness/chem/molecularbasis.h:430
double rangesq() const
Definition madness/chem/molecularbasis.h:463
const char * get_desc() const
Definition madness/chem/molecularbasis.h:450
const double zz
Definition madness/chem/molecularbasis.h:411
const int ibf
Definition madness/chem/molecularbasis.h:413
AtomicBasisFunction(const AtomicBasisFunction &aofunc)
Definition madness/chem/molecularbasis.h:422
Definition madness/chem/molecularbasis.h:672
bool operator()(long i, long j) const
Definition madness/chem/molecularbasis.h:676
AnalysisSorter(const Tensor< T > &v)
Definition madness/chem/molecularbasis.h:675
const Tensor< T > v
Definition madness/chem/molecularbasis.h:673
Contracted Gaussian basis.
Definition madness/chem/molecularbasis.h:469
AtomicBasisSet(std::string filename)
Definition madness/chem/molecularbasis.h:505
void read_nw_file(std::string filename)
read the atomic basis set from file
Definition molecularbasis.cc:206
double eval_guess_density(const Molecule &molecule, double x, double y, double z) const
Evaluates the guess density.
Definition madness/chem/molecularbasis.h:646
AtomicBasisFunction get_atomic_basis_function(const Molecule &molecule, size_t ibf) const
Returns the ibf'th atomic basis function.
Definition madness/chem/molecularbasis.h:600
void read_file(std::string filename)
read the atomic basis set from file
Definition molecularbasis.cc:119
std::vector< T > load_tixml_vector(TiXmlElement *node, int n, const char *name)
Definition madness/chem/molecularbasis.h:474
bool is_supported(int atomic_number) const
Definition madness/chem/molecularbasis.h:661
int nbf(const Molecule &molecule) const
Given a molecule count the number of basis functions.
Definition madness/chem/molecularbasis.h:624
AtomicBasisSet()
Definition madness/chem/molecularbasis.h:502
const Tensor< double > & get_avec(const Molecule &molecule, size_t iat) const
Returns the atomic alpha eigenvectors for atom iat.
Definition madness/chem/molecularbasis.h:562
std::vector< AtomicBasis > ag
Basis associated by atomic number = 1, 2, ...; 0=Bq.
Definition madness/chem/molecularbasis.h:471
void print(const Molecule &molecule) const
Print basis info for atoms in the molecule (once for each unique atom type)
Definition molecularbasis.cc:90
Tensor< T > load_tixml_matrix(TiXmlElement *node, int n, int m, const char *name)
Definition madness/chem/molecularbasis.h:487
std::string name
Definition madness/chem/molecularbasis.h:470
void serialize(Archive &ar)
Definition madness/chem/molecularbasis.h:739
std::string get_name() const
Definition madness/chem/molecularbasis.h:509
void modify_dmat_psp(int atn, double zeff)
Eliminates core orbitals from the density matrix for pseudopotential calculations.
Definition molecularbasis.cc:237
void atoms_to_bfn(const Molecule &molecule, std::vector< int > &at_to_bf, std::vector< int > &at_nbf) const
Makes map from atoms to first basis function on atom and number of basis functions on atom.
Definition madness/chem/molecularbasis.h:526
int basisfn_to_atom(const Molecule &molecule, size_t ibf) const
Returns the number of the atom the ibf'th basis function is on.
Definition madness/chem/molecularbasis.h:580
void print_anal(const Molecule &molecule, const Tensor< T > &v) const
Given a vector of AO coefficients prints an analysis.
Definition madness/chem/molecularbasis.h:690
void print_all() const
Print basis info for all supported atoms.
Definition molecularbasis.cc:109
void eval(const Molecule &molecule, double x, double y, double z, double *bf) const
Evaluates the basis functions.
Definition madness/chem/molecularbasis.h:636
const Tensor< double > & get_aeps(const Molecule &molecule, size_t iat) const
Returns the atomic alpha eigenvalues for atom iat.
Definition madness/chem/molecularbasis.h:571
void shells_to_bfn(const Molecule &molecule, std::vector< int > &sh_to_bf, std::vector< int > &sh_nbf) const
Makes map from shells to first basis function on she and number of basis functions on sh.
Definition madness/chem/molecularbasis.h:542
Represents multiple shells of contracted gaussians on a single center.
Definition madness/chem/molecularbasis.h:217
const Tensor< double > & get_dmatpsp() const
Definition madness/chem/molecularbasis.h:345
Tensor< double > beps
Definition madness/chem/molecularbasis.h:221
const std::vector< ContractedGaussianShell > & get_shells() const
Returns a const reference to the shells.
Definition madness/chem/molecularbasis.h:264
void set_aocc(Tensor< double > &occ)
Definition madness/chem/molecularbasis.h:377
Tensor< double > aocc
Definition madness/chem/molecularbasis.h:221
void set_dmatpsp(Tensor< double > &mat)
Definition madness/chem/molecularbasis.h:349
std::vector< ContractedGaussianShell > g
Definition madness/chem/molecularbasis.h:218
double rmaxsq
Definition madness/chem/molecularbasis.h:219
int nshell() const
Returns the number of shells on the center.
Definition madness/chem/molecularbasis.h:259
Tensor< double > bocc
Definition madness/chem/molecularbasis.h:221
void set_aoccpsp(Tensor< double > &occ)
Definition madness/chem/molecularbasis.h:393
int numbf
Definition madness/chem/molecularbasis.h:220
const Tensor< double > & get_aoccpsp() const
Definition madness/chem/molecularbasis.h:385
const Tensor< double > & get_aocc() const
Definition madness/chem/molecularbasis.h:369
Tensor< double > avec
Definition madness/chem/molecularbasis.h:221
void set_guess_info(const Tensor< double > &dmat, const Tensor< double > &dmatpsp, const Tensor< double > &avec, const Tensor< double > &bvec, const Tensor< double > &aocc, const Tensor< double > &bocc, const Tensor< double > &aeps, const Tensor< double > &beps, const Tensor< double > &aoccpsp, const Tensor< double > &boccpsp)
Definition madness/chem/molecularbasis.h:236
void set_bocc(Tensor< double > &occ)
Definition madness/chem/molecularbasis.h:381
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 madness/chem/molecularbasis.h:273
Tensor< double > dmatpsp
Definition madness/chem/molecularbasis.h:221
void set_dmat(Tensor< double > &mat)
Definition madness/chem/molecularbasis.h:337
Tensor< double > aeps
Definition madness/chem/molecularbasis.h:221
AtomicBasis(const std::vector< ContractedGaussianShell > &g)
Definition madness/chem/molecularbasis.h:226
AtomicBasis()
Definition madness/chem/molecularbasis.h:224
const Tensor< double > & get_aeps() const
Definition madness/chem/molecularbasis.h:361
int nbf() const
Returns the number of basis functions on the center.
Definition madness/chem/molecularbasis.h:254
Tensor< double > boccpsp
Definition madness/chem/molecularbasis.h:221
const Tensor< double > & get_dmat() const
Definition madness/chem/molecularbasis.h:333
void serialize(Archive &ar)
Definition madness/chem/molecularbasis.h:402
bool has_guess_info() const
Definition madness/chem/molecularbasis.h:329
const Tensor< double > & get_avec() const
Definition madness/chem/molecularbasis.h:353
Tensor< double > bvec
Definition madness/chem/molecularbasis.h:221
void set_boccpsp(Tensor< double > &occ)
Definition madness/chem/molecularbasis.h:397
bool has_guesspsp_info() const
Definition madness/chem/molecularbasis.h:341
const Tensor< double > & get_beps() const
Definition madness/chem/molecularbasis.h:365
const Tensor< double > & get_bocc() const
Definition madness/chem/molecularbasis.h:373
Tensor< double > aoccpsp
Definition madness/chem/molecularbasis.h:221
double eval_guess_density(double x, double y, double z, bool pspat) const
Evaluates the guess atomic density at point x, y, z relative to atomic center.
Definition madness/chem/molecularbasis.h:290
Tensor< double > dmat
Definition madness/chem/molecularbasis.h:221
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 madness/chem/molecularbasis.h:314
const Tensor< double > & get_boccpsp() const
Definition madness/chem/molecularbasis.h:389
const Tensor< double > & get_bvec() const
Definition madness/chem/molecularbasis.h:357
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
Represents a single shell of contracted, Cartesian, Gaussian primitives.
Definition madness/chem/molecularbasis.h:51
void serialize(Archive &ar)
Definition madness/chem/molecularbasis.h:211
static const int maxbf
Maximum number of basis functions in a shell.
Definition madness/chem/molecularbasis.h:54
int nbf() const
Returns the number of basis functions in the shell.
Definition madness/chem/molecularbasis.h:179
ContractedGaussianShell()
Definition madness/chem/molecularbasis.h:90
int numbf
Number of basis functions in shell (type+1)*(type+2)/2.
Definition madness/chem/molecularbasis.h:60
ContractedGaussianShell(int type, const std::vector< double > &coeff, const std::vector< double > &expnt, bool donorm=true)
Definition madness/chem/molecularbasis.h:93
double eval_radial(double rsq) const
Evaluates the radial part of the contracted function.
Definition madness/chem/molecularbasis.h:113
std::vector< double > coeff
Definition madness/chem/molecularbasis.h:57
double * eval(double rsq, double x, double y, double z, double *bf) const
Evaluates the entire shell returning the incremented result pointer.
Definition madness/chem/molecularbasis.h:125
const std::vector< double > & get_expnt() const
Returns a const reference to the exponents.
Definition madness/chem/molecularbasis.h:194
std::vector< double > expnt
Definition madness/chem/molecularbasis.h:58
const std::vector< double > & get_coeff() const
Returns a const reference to the coefficients.
Definition madness/chem/molecularbasis.h:189
int type
Angular momentum = 0, 1, 2, ...
Definition madness/chem/molecularbasis.h:56
int nprim() const
Returns the number of primitives in the contraction.
Definition madness/chem/molecularbasis.h:184
double rangesq() const
Returns square of the distance beyond which function is less than 1e-8.
Definition madness/chem/molecularbasis.h:107
int angular_momentum() const
Returns the shell angular momentum.
Definition madness/chem/molecularbasis.h:174
const char * get_desc(int ibf) const
Returns a string description of the basis function type.
Definition madness/chem/molecularbasis.h:199
static const int maxtype
Maximum angular momentum supported.
Definition madness/chem/molecularbasis.h:53
double rsqmax
Definition madness/chem/molecularbasis.h:59
void normalize()
Definition madness/chem/molecularbasis.h:62
Definition molecule.h:129
A tensor is a multidimensional array.
Definition tensor.h:317
T * ptr()
Returns a pointer to the internal data.
Definition tensor.h:1825
A simple, fixed dimension vector.
Definition vector.h:64
Defines common mathematical and physical constants.
static const double R
Definition csqrt.cc:46
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
static const double v
Definition hatom_sf_dirac.cc:20
static double pow(const double *a, const double *b)
Definition lda.h:74
Macros and tools pertaining to the configuration of MADNESS.
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:182
#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
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
const AtomicData & get_atomic_data(unsigned int atomic_number)
Definition atomutil.cc:167
NDIM & f
Definition mra.h:2481
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition mra.h:2066
static long abs(long a)
Definition tensor.h:218
static std::vector< int > at_to_bf
Definition preal.cc:19
static std::vector< int > at_nbf
Definition preal.cc:19
static const double m
Definition relops.cc:9
static const double thresh
Definition rk.cc:45
const char *const symbol
Definition atomutil.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
static Molecule molecule
Definition testperiodicdft.cc:39