MADNESS  0.10.1
apps/periodic_old/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 MOLECULAR_BASIS_H
33 #define MOLECULAR_BASIS_H
34 
35 #include <madness/madness_config.h>
36 #include <madness/constants.h>
37 #include "mentity.h"
38 #include <madness/external/tinyxml/tinyxml.h>
39 #include <madness/tensor/tensor.h>
40 
41 #include <vector>
42 #include <algorithm>
43 #include <iostream>
44 #include <sstream>
45 #include <iomanip>
46 #include <cstdio>
47 
48 /// Represents a single shell of contracted, Cartesian, Gaussian primitives
50  int type; ///< Angular momentum = 0, 1, 2, ...
51  std::vector<double> coeff;
52  std::vector<double> expnt;
53  double rsqmax;
54  int numbf; ///< Number of basis functions in shell (type+1)*(type+2)/2
55 
56  void normalize() {
57  // nwcchem cartesian normalization conventions
58  // translation of nmcoeff.F into python and thence to c++
59  int np = coeff.size();
60  if (np == 1) coeff[0] = 1.0e0;
61 
62  double pi32=pow(madness::constants::pi,1.5);
63  int l_lim = 2*type - 1;
64  double f = 1.0e00;
65  for (int n=l_lim; n>1; n-=2) f *= n;
66 
67  for (int n=0; n<np; n++)
68  coeff[n] *= pow(2.e0*expnt[n]/madness::constants::pi,0.75e0)*pow(4.e0*expnt[n],0.5E0*type)/sqrt(f);
69 
70  double sum = 0.0;
71  for (int n1=0; n1<np; n1++) {
72  for (int n2=0; n2<np; n2++) {
73  double S =pi32/pow(expnt[n1]+expnt[n2],1.5e0+type)/pow(2e0,type);
74  sum = sum + coeff[n1]*coeff[n2]*S;
75  }
76  }
77  sum *= f;
78 
79  f = 1e0/sqrt(sum);
80  for (int n=0; n<np; n++) coeff[n] *= f;
81  }
82 
83 public:
85  : type(-1), coeff(), expnt(), rsqmax(0.0), numbf(0) {};
86 
88  const std::vector<double>& coeff,
89  const std::vector<double>& expnt,
90  bool donorm=true)
91  : type(type), coeff(coeff), expnt(expnt), numbf((type+1)*(type+2)/2) {
92  if (donorm) normalize();
93  double minexpnt = expnt[0];
94  for (unsigned int i=1; i<expnt.size(); i++)
95  minexpnt = std::min(minexpnt,expnt[i]);
96  rsqmax = 18.4/minexpnt; // 18.4 = 8*ln(10)
97  }
98 
99 
100  /// Returns square of the distance beyond which function is less than 1e-8.
101  double rangesq() const {
102  return rsqmax;
103  }
104 
105 
106  /// Evaluates the radial part of the contracted function
107  double eval_radial(double rsq) const {
108  if (rsq > rsqmax) return 0.0;
109  double sum = 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);
113  }
114  return sum;
115  }
116 
117 
118  /// Evaluates the entire shell returning the incremented result pointer
119  double* eval(double rsq, double x, double y, double z, double* bf) const {
120  double R = eval_radial(rsq);
121  if (fabs(R) < 1e-8) {
122  for (int i=0; i<numbf; i++) bf[i] = 0.0;
123 
124  }
125  else {
126  switch (type) {
127  case 0:
128  bf[0] = R;
129  break;
130  case 1:
131  bf[0] = R*x;
132  bf[1] = R*y;
133  bf[2] = R*z;
134  break;
135  case 2:
136  bf[0] = R*x*x;
137  bf[1] = R*x*y;
138  bf[2] = R*x*z ;
139  bf[3] = R*y*y ;
140  bf[4] = R*y*z ;
141  bf[5] = R*z*z;
142  break;
143  case 3:
144  bf[0] = R*x*x*x;
145  bf[1] = R*x*x*y;
146  bf[2] = R*x*x*z;
147  bf[3] = R*x*y*y;
148  bf[4] = R*x*y*z;
149  bf[5] = R*x*z*z;
150  bf[6] = R*y*y*y;
151  bf[7] = R*y*y*z;
152  bf[8] = R*y*z*z;
153  bf[9] = R*z*z*z;
154  break;
155 
156  default:
157  throw "UNKNOWN ANGULAR MOMENTUM";
158  }
159  }
160  return bf+numbf;
161  }
162 
163 
164  /// Returns the shell angular momentum
165  int angular_momentum() const {
166  return type;
167  }
168 
169  /// Returns the number of basis functions in the shell
170  int nbf() const {
171  return numbf;
172  }
173 
174  /// Returns the number of primitives in the contraction
175  int nprim() const {
176  return coeff.size();
177  }
178 
179  /// Returns a const reference to the coefficients
180  const std::vector<double>& get_coeff() const {
181  return coeff;
182  }
183 
184  /// Returns a const reference to the exponents
185  const std::vector<double>& get_expnt() const {
186  return expnt;
187  }
188 
189  /// Returns a string description of the basis function type
190  const char* get_desc(int ibf) const {
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"}
196  };
197  MADNESS_ASSERT(ibf<numbf && ibf >= 0);
198  return tags[type][ibf];
199  }
200 
201  template <typename Archive>
202  void serialize(Archive& ar) {
203  ar & type & coeff & expnt & rsqmax & numbf;
204  }
205 };
206 
207 
208 std::ostream& operator<<(std::ostream& s, const ContractedGaussianShell& c) {
209  static const char* tag[] = {"s","p","d","f","g"};
210  char buf[32768];
211  char* p = buf;
212  const std::vector<double>& coeff = c.get_coeff();
213  const std::vector<double>& expnt = c.get_expnt();
214 
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, ", ");
219  }
220  p += sprintf(p, "]");
221  s << buf;
222  return s;
223 }
224 
225 
226 /// Represents multiple shells of contracted gaussians on a single center
227 class AtomicBasis {
228  std::vector<ContractedGaussianShell> g;
229  double rmaxsq;
230  int numbf;
231  Tensor<double> dmat, avec, bvec;
232 
233 public:
234  AtomicBasis() : g(), rmaxsq(0.0), numbf(0) {};
235 
236  AtomicBasis(const std::vector<ContractedGaussianShell>& g)
237  : g(g) {
238  rmaxsq = 0.0;
239  numbf = 0;
240  for (unsigned int i=0; i<g.size(); i++) {
241  rmaxsq = std::max(rmaxsq, g[i].rangesq());
242  numbf += g[i].nbf();
243  }
244  }
245 
246  void set_guess_info(const Tensor<double>& dmat,
247  const Tensor<double>& avec, const Tensor<double>& bvec) {
248  this->dmat = copy(dmat);
249  this->avec = copy(avec);
250  this->bvec = copy(bvec);
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_ASSERT(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) const {
292  double rsq = x*x + y*y + z*z;
293  if (rsq > rmaxsq) return 0.0;
294 
295  double bf[numbf];
296  eval(x, y, z, bf);
297  const double* p = dmat.ptr();
298  double sum = 0.0;
299  for (int i=0; i<numbf; i++, p+=numbf) {
300  double sumj = 0.0;
301  for (int j=0; j<numbf; ++j)
302  sumj += p[j]*bf[j];
303  sum += bf[i]*sumj;
304  }
305  return sum;
306  }
307 
308  /// Return shell that contains basis function ibf and also return index of function in the shell
309  const ContractedGaussianShell& get_shell_from_basis_function(int ibf, int& ibf_in_shell) const {
310  int n=0;
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;
315  return g[i];
316  }
317  else {
318  n += g[i].nbf();
319  }
320  }
321  MADNESS_EXCEPTION("AtomicBasis: get_shell_from_basis_function", ibf*100000 + nbf());
322  }
323 
324  bool has_guess_info() const {
325  return dmat.size()>0;
326  }
327 
328  const Tensor<double>& get_dmat() const {
329  return dmat;
330  };
331 
332  const Tensor<double>& get_avec() const {
333  return avec;
334  };
335 
336  const Tensor<double>& get_bvec() const {
337  return bvec;
338  };
339 
340  template <typename Archive>
341  void serialize(Archive& ar) {
342  ar & g & rmaxsq & numbf & dmat & avec & bvec;
343  }
344 
345 };
346 
347 std::ostream& operator<<(std::ostream& s, const AtomicBasis& c) {
348  const std::vector<ContractedGaussianShell>& shells = c.get_shells();
349  for (int i=0; i<c.nshell(); i++) {
350  s << " " << shells[i] << std::endl;
351  }
352  if (c.has_guess_info()) {
353  s << " " << "Guess density matrix" << std::endl;
354  s << c.get_dmat();
355  }
356 
357  return s;
358 }
359 
360 /// Used to represent one basis function from a shell on a specific center
362 private:
363  const double xx, yy, zz; // Coordinates of the center
364  const ContractedGaussianShell& shell; // Reference to the underlying atomic shell
365  const int ibf; // Index of basis function in the shell (0, 1, ...)
366  const int nbf; // Number of functions in the shell
367 
368 public:
369  AtomicBasisFunction(double x, double y, double z,
370  const ContractedGaussianShell& shell, int ibf)
371  : xx(x), yy(y), zz(z), shell(shell), ibf(ibf), nbf(shell.nbf()) {}
372 
373 
375  : xx(aofunc.xx)
376  , yy(aofunc.yy)
377  , zz(aofunc.zz)
378  , shell(aofunc.shell)
379  , ibf(aofunc.ibf)
380  , nbf(aofunc.nbf) {}
381 
382  double operator()(double x, double y, double z) const {
383  double bf[nbf];
384  x-=xx;
385  y-=yy;
386  z-=zz;
387  double rsq = x*x + y*y + z*z;
388  shell.eval(rsq, x, y, z, bf);
389  return bf[ibf];
390  }
391 
392  void print_me(std::ostream& s) const {
393  s << "atomic basis function: center " << xx << " " << yy << " " << zz << " : ibf " << ibf << " nbf " << nbf << " : shell " << shell << std::endl;
394  }
395 
397  return shell;
398  }
399 
400  int get_index() const {
401  return ibf;
402  }
403 
404  const char* get_desc() const {
405  return shell.get_desc(ibf);
406  }
407 
408  void get_coords(double& x, double& y, double& z) const {
409  x=xx; y=yy; z=zz;
410  return;
411  }
412 
413  double rangesq() const {
414  return shell.rangesq();
415  }
416 };
417 
418 std::ostream& operator<<(std::ostream& s, const AtomicBasisFunction& a) {
419  a.print_me(s);
420  return s;
421 };
422 
423 /// Contracted Gaussian basis
425  std::string name;
426  std::vector<AtomicBasis> ag; //< Basis associated by atomic number = 1, 2, ...; 0=Bq.
427 
428  template <typename T>
429  std::vector<T> load_tixml_vector(TiXmlElement* node, int n, const char* name) {
430  TiXmlElement* child = node->FirstChildElement(name);
431  MADNESS_ASSERT(child);
432  std::istringstream s(child->GetText());
433  std::vector<T> r(n);
434  for (int i=0; i<n; i++) {
435  MADNESS_ASSERT(s >> r[i]);
436  }
437  return r;
438  }
439 
440  template <typename T>
441  Tensor<T> load_tixml_matrix(TiXmlElement* node, int n, int m, const char* name) {
442  TiXmlElement* child = node->FirstChildElement(name);
443  MADNESS_ASSERT(child);
444  std::istringstream s(child->GetText());
445  Tensor<T> r(n,m);
446  for (int i=0; i<n; i++) {
447  for (int j=0; j<m; j++) {
448  MADNESS_ASSERT(s >> r(i,j));
449  }
450  }
451  return r;
452  }
453 
454 public:
455  AtomicBasisSet() : name("unknown"), ag(110) {}
456 
457 
458  AtomicBasisSet(std::string filename) : name(""), ag(110) {
460  }
461 
462  void read_file(std::string filename) {
463  static const bool debug = false;
464  TiXmlDocument doc(filename);
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;
470  MADNESS_EXCEPTION("AtomicBasisSet: Failed loading basis set",0);
471  }
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;
476  }
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;
480  int atn = symbol_to_atomic_number(symbol);
481  std::vector<ContractedGaussianShell> g;
482  for (TiXmlElement* shell=node->FirstChildElement(); shell; shell=shell->NextSiblingElement()) {
483  const char* type = shell->Attribute("type");
484  int nprim=-1;
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");
491  g.push_back(ContractedGaussianShell(0,scoeff,expnt));
492  g.push_back(ContractedGaussianShell(1,pcoeff,expnt));
493  }
494  else {
495  static const char* tag[] = {"S","P","D","F","G"};
496  int i;
497  for (i=0; i<5; i++) {
498  if (strcmp(type,tag[i]) == 0) goto foundit;
499  }
500  MADNESS_EXCEPTION("Loading atomic basis set: bad shell type?",0);
501 foundit:
502  std::vector<double> coeff = load_tixml_vector<double>(shell, nprim, "coefficients");
503  g.push_back(ContractedGaussianShell(i, coeff, expnt));
504  }
505  }
506  ag[atn] = AtomicBasis(g);
507  }
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;
511  int atn = symbol_to_atomic_number(symbol);
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);
518  }
519  else {
520  MADNESS_EXCEPTION("Loading atomic basis set: unexpected XML element", 0);
521  }
522  }
523 
524  }
525 
526 
527  /// Makes map from atoms to first basis function on atom and number of basis functions on atom
528  void atoms_to_bfn(const MolecularEntity& mentity, std::vector<int>& at_to_bf, std::vector<int>& at_nbf) {
529  at_to_bf = std::vector<int>(mentity.natom());
530  at_nbf = std::vector<int>(mentity.natom());
531 
532  int n = 0;
533  for (int i=0; i<mentity.natom(); i++) {
534  const Atom& atom = mentity.get_atom(i);
535  const int atn = atom.atomic_number;
537  at_to_bf[i] = n;
538  at_nbf[i] = ag[atn].nbf();
539  n += at_nbf[i];
540  }
541  }
542 
543 
544  /// Returns the number of the atom the ibf'th basis function is on
545  int basisfn_to_atom(const MolecularEntity& mentity, int ibf) const {
546  MADNESS_ASSERT(ibf >= 0);
547  int n = 0;
548  for (int i=0; i<mentity.natom(); i++) {
549  // Is the desired function on this atom?
550  const Atom& atom = mentity.get_atom(i);
551  const int atn = atom.atomic_number;
553  const int nbf_on_atom = ag[atn].nbf();
554  if (ibf >= n && (n+nbf_on_atom) > ibf) {
555  return i;
556  }
557  else {
558  n += nbf_on_atom;
559  }
560  }
561  MADNESS_EXCEPTION("AtomicBasisSet: get_atomic_basis_function: confused?", ibf);
562  }
563 
564  /// Returns the ibf'th atomic basis function
566  MADNESS_ASSERT(ibf >= 0);
567  int n = 0;
568  for (int i=0; i<mentity.natom(); i++) {
569  // Is the desired function on this atom?
570  const Atom& atom = mentity.get_atom(i);
571  const int atn = atom.atomic_number;
573  const int nbf_on_atom = ag[atn].nbf();
574  if (ibf >= n && (n+nbf_on_atom) > ibf) {
575  int index;
576  const ContractedGaussianShell& shell =
577  ag[atn].get_shell_from_basis_function(ibf-n, index);
578  return AtomicBasisFunction(atom.x, atom.y, atom.z, shell, index);
579  }
580  else {
581  n += nbf_on_atom;
582  }
583  }
584  MADNESS_EXCEPTION("AtomicBasisSet: get_atomic_basis_function: confused?", ibf);
585  }
586 
587 
588  /// Given a molecule count the number of basis functions
589  int nbf(const MolecularEntity& mentity) const {
590  int n = 0;
591  for (int i=0; i<mentity.natom(); i++) {
592  const Atom& atom = mentity.get_atom(i);
593  const int atn = atom.atomic_number;
595  n += ag[atn].nbf();
596  }
597  return n;
598  }
599 
600  /// Evaluates the basis functions
601  void eval(const MolecularEntity& mentity, double x, double y, double z, double *bf) const {
602  for (int i=0; i<mentity.natom(); i++) {
603  const Atom& atom = mentity.get_atom(i);
604  const int atn = atom.atomic_number;
605  bf = ag[atn].eval(x-atom.x, y-atom.y, z-atom.z, bf);
606  }
607  }
608 
609 
610  /// Evaluates the guess density
611  double eval_guess_density(const MolecularEntity& mentity, double x, double y, double z) const {
612  double sum = 0.0;
613  for (int i=0; i<mentity.natom(); i++) {
614  const Atom& atom = mentity.get_atom(i);
615  const int atn = atom.atomic_number;
616  sum += ag[atn].eval_guess_density(x-atom.x, y-atom.y, z-atom.z);
617  }
618  return sum;
619  }
620 
621  bool is_supported(int atomic_number) const {
622  return ag[atomic_number].nbf() > 0;
623  }
624 
625  /// Print basis info for atoms in the molecule (once for each unique atom type)
626  void print(const MolecularEntity& mentity) const {
627  std::cout << "\n " << name << " atomic basis set" << std::endl;
628  for (int i=0; i<mentity.natom(); i++) {
629  const Atom& atom = mentity.get_atom(i);
630  const unsigned int atn = atom.atomic_number;
631  for (int j=0; j<i; j++) {
632  if (mentity.get_atom(j).atomic_number == atn)
633  goto doneitalready;
634  }
635  std::cout << std::endl;
636  std::cout << " " << get_atomic_data(i).symbol << std::endl;
637  std::cout << ag[atn];
638 doneitalready:
639  ;
640  }
641  }
642 
643  template <typename T>
645  const Tensor<T> v;
646  public:
647  AnalysisSorter(const Tensor<T>& v) : v(v) {}
648  bool operator()(long i, long j) const {
649  return std::abs(v[i]) > std::abs(v[j]);
650  }
651  };
652 
653  /// Given a vector of AO coefficients prints an analysis
654 
655  /// For each significant coeff it prints
656  /// - atomic symbol
657  /// - atom number
658  /// - basis function type (e.g., dxy)
659  /// - basis function number
660  /// - MO coeff
661  template <typename T>
662  void print_anal(const MolecularEntity& mentity, const Tensor<T>& v) {
663  const double thresh = 0.2*v.normf();
664  if (thresh == 0.0) {
665  printf(" zero vector\n");
666  return;
667  }
668  long nbf = int(v.dim[0]);
669  long list[nbf];
670  long ngot=0;
671  for (long i=0; i<nbf; i++) {
672  if (std::abs(v(i)) > thresh) {
673  list[ngot++] = i;
674  }
675  }
676  std::sort(list,list+ngot,AnalysisSorter<T>(v));
677 
678  const char* format;
679  if (mentity.natom() < 10) {
680  format = " %2s(%1d)%4s(%2ld)%6.3f ";
681  }
682  else if (mentity.natom() < 100) {
683  format = " %2s(%2d)%4s(%3ld)%6.3f ";
684  }
685  else if (mentity.natom() < 1000) {
686  format = " %2s(%3d)%4s(%4ld)%6.3f ";
687  }
688  else {
689  format = " %2s(%4d)%4s(%5ld)%6.3f ";
690  }
691  for (long ii=0; ii<ngot; ii++) {
692  long ibf = list[ii];
693 
694  const int iat = basisfn_to_atom(mentity, ibf);
695  const Atom& atom = mentity.get_atom(iat);
696  const AtomicBasisFunction ao = get_atomic_basis_function(mentity, ibf);
697  const char* desc = ao.get_desc();
698  const char* element = get_atomic_data(atom.atomic_number).symbol;
699 
700  // This will need wrapping in a template for a complex MO vector
701  printf(format, element, iat, desc, ibf, v[ibf]);
702  }
703  printf("\n");
704  }
705 
706  /// Print basis info for all supported atoms
707  void print_all() const {
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) {
711  std::cout << " " << get_atomic_data(i).symbol << std::endl;
712  std::cout << ag[i];
713  }
714  }
715  }
716 
717  template <typename Archive>
718  void serialize(Archive& ar) {
719  ar & name & ag;
720  }
721 };
722 
723 
724 
725 #endif
std::ostream & operator<<(std::ostream &s, const ContractedGaussianShell &c)
Definition: apps/periodic_old/molecularbasis.h:208
Definition: mentity.h:71
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
Definition: mentity.h:95
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