1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2011
9 * Copyright (c) 2010-2011, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16
17 #ifndef ERKALE_BASIS
18 #define ERKALE_BASIS
19
20 #include "global.h"
21
22 #include <armadillo>
23 #include <vector>
24 #include <string>
25 #include <cfloat>
26
27 /// Forward declaration
28 class Settings;
29
30 #include "xyzutils.h"
31
32 /// Angular momentum notation for shells
33 const char shell_types[]={'S','P','D','F','G','H','I','J','K','L','M','N','O','Q','R'};
34 /// Maximum angular momentum supported in current version of ERKALE
35 const int max_am=sizeof(shell_types)/sizeof(shell_types[0])-1;
36
37 /// Structure for defining shells of functions
38 typedef struct {
39 /// Angular momentum in x direction
40 int l;
41 /// Angular momentum in y direction
42 int m;
43 /// Angular momentum in z direction
44 int n;
45 /// Relative normalization coefficient
46 double relnorm;
47 } shellf_t;
48
49 /// Coordinates structure
50 typedef struct {
51 /// x coordinate
52 double x;
53 /// y coordinate
54 double y;
55 /// z coordinate
56 double z;
57 } coords_t;
58
59 arma::vec coords_to_vec(const coords_t & c);
60 coords_t vec_to_coords(const arma::vec & v);
61
62 // Forward declaration
63 class GaussianShell;
64
65 /// Nucleus structure
66 typedef struct {
67 /// Index of nucleus
68 size_t ind;
69 /// Location of nucleus
70 coords_t r;
71 /// Counterpoise nucleus..?
72 bool bsse;
73
74 /// Type of nucleus
75 std::string symbol;
76 /// Nuclear charge
77 int Z;
78 /// Net charge in system (used for atomic guess)
79 int Q;
80
81 /// List of shells located on nucleus
82 std::vector<const GaussianShell *> shells;
83 } nucleus_t;
84
85 /// Comparison operator
86 bool operator==(const nucleus_t & lhs, const nucleus_t & rhs);
87
88 /// Structure for unique shell pairs
89 typedef struct {
90 /// Index of first shell
91 size_t is;
92 /// Index of first functions on first shell
93 size_t i0;
94 /// Angular momentum of first shell
95 int li;
96
97 /// Index of second shell
98 size_t js;
99 /// Index of first function on second shell
100 size_t j0;
101 /// Angular momentum of second shell
102 int lj;
103 } shellpair_t;
104
105 /// Comparison operator for shellpairs for ordering into libint order
106 bool operator<(const shellpair_t & lhs, const shellpair_t & rhs);
107
108 /// Helper for integral sorts
109 typedef struct {
110 /// First shell
111 size_t is;
112 /// First function on shell
113 size_t i0;
114 /// Amount of functions on shell
115 size_t Ni;
116
117 /// Second shell
118 size_t js;
119 /// First function on shell
120 size_t j0;
121 /// Amount of functions on shell
122 size_t Nj;
123
124 /// Maximum (uv|uv)^1/2 on shell
125 double eri;
126 } eripair_t;
127
128 /// Comparison operator
129 bool operator<(const eripair_t & lhs, const eripair_t & rhs);
130
131 // Forward declaration
132 class BasisSetLibrary;
133 class ElementBasisSet;
134
135 /// Order shells solely on merit of exponents (for forming density fitting basis)
136 bool exponent_compare(const GaussianShell & lhs, const GaussianShell & rhs);
137
138 /// Comparison operator
139 bool operator==(const coords_t & lhs, const coords_t & rhs);
140 /// Compute displacement
141 coords_t operator-(const coords_t & lhs, const coords_t & rhs);
142 /// Compute sum
143 coords_t operator+(const coords_t & lhs, const coords_t & rhs);
144 /// Compute scaling by division
145 coords_t operator/(const coords_t & lhs, double fac);
146 /// Compute scaling by multiplication
147 coords_t operator*(const coords_t & lhs, double fac);
148
149 /// Compute squared norm
normsq(const coords_t & r)150 inline double normsq(const coords_t & r) {
151 return r.x*r.x + r.y*r.y + r.z*r.z;
152 }
153
154 /// Compute norm
norm(const coords_t & r)155 inline double norm(const coords_t & r) {
156 return sqrt(normsq(r));
157 }
158
159 /// Structure for contractions
160 typedef struct {
161 /// Coefficient
162 double c;
163 /// Exponent
164 double z;
165 } contr_t;
166 /// Comparison for contractions
167 bool operator<(const contr_t & lhs, const contr_t & rhs);
168 /// Identity for contractions
169 bool operator==(const contr_t & lhs, const contr_t & rhs);
170
171 #include "basislibrary.h"
172
173
174 /**
175 * \class BasisSet
176 *
177 * \brief Class for a Gaussian basis set
178 *
179 * This class contains the data structures necessary for Gaussian
180 * basis sets, and functions for evaluating integrals over them.
181 *
182 * \author Susi Lehtola
183 * \date 2011/05/05 20:17
184 */
185
186 /// Basis set
187 class BasisSet {
188 /// Nuclei
189 std::vector<nucleus_t> nuclei;
190 /// Basis functions
191 std::vector<GaussianShell> shells;
192
193 /// Use spherical harmonics by default as basis?
194 bool uselm;
195 /// Use cartesian s and p functions if spherical harmonics are used?
196 bool optlm;
197
198 /// Internuclear distances
199 arma::mat nucleardist;
200 /// List of unique shell pairs
201 std::vector<shellpair_t> shellpairs;
202
203 /// Ranges of shells
204 std::vector<double> shell_ranges;
205
206 /// Check for same geometry
207 bool same_geometry(const BasisSet & rhs) const;
208 /// Check for same shells
209 bool same_shells(const BasisSet & rhs) const;
210
211 public:
212 /// Dummy constructor
213 BasisSet();
214 /// Construct basis set with Nat atoms, using given settings
215 BasisSet(size_t Nat);
216 /// Destructor
217 ~BasisSet();
218
219 /**
220 * Generate density fitting basis
221 *
222 * The procedure has been documented in the article
223 *
224 * R. Yang, A. P. Rendell and M. J. Frisch, "Automatically generated
225 * Coulomb fitting basis sets: Design and accuracy for systems
226 * containing H to Kr", J. Chem. Phys. 127 (2007), 074102.
227 */
228 BasisSet density_fitting(double fsam=1.5, int lmaxinc=1) const;
229
230 /**
231 * Generate Coulomb and exchange fitting basis
232 *
233 * The procedure has been documented in the article
234 *
235 * O. Vahtras, J. Almlöf and M. W. Feyereisen, "Integral approximations
236 * for LCAO-SCF calculations", Chem. Phys. Lett. 213, p. 514 - 518 (1993).
237 */
238 BasisSet exchange_fitting() const;
239
240 /// Decontract basis set, m gives mapping from old functions to new ones
241 BasisSet decontract(arma::mat & m) const;
242
243 /// Add nucleus
244 void add_nucleus(const nucleus_t & nuc);
245 /// Add a shell to a nucleus and sort functions if wanted
246 void add_shell(size_t nucind, const GaussianShell & sh, bool sort=true);
247 /// Add a shell to a nucleus and sort functions if wanted
248 void add_shell(size_t nucind, int am, bool uselm, const std::vector<contr_t> & C, bool sort=true);
249 /// Add all shells to a nucleus and sort functions if wanted
250 void add_shells(size_t nucind, ElementBasisSet el, bool sort=true);
251
252 /// Sort shells in nuclear order, then by angular momentum, then by exponents
253 void sort();
254 /// Check numbering of basis functions
255 void check_numbering();
256 /// Update the nuclear list of shells
257 void update_nuclear_shell_list();
258
259 /* Finalization routines */
260
261 /// Compute nuclear distance table
262 void compute_nuclear_distances();
263
264 /// Form list of unique shell pairs
265 void form_unique_shellpairs();
266 /// Get list of unique shell pairs
267 std::vector<shellpair_t> get_unique_shellpairs() const;
268
269 /// Get list of ERI pairs
270 std::vector<eripair_t> get_eripairs(arma::mat & Q, arma::mat & M, double thr, double omega=0.0, double alpha=1.0, double beta=0.0, bool verbose=false) const;
271
272 /// Convert contractions from normalized primitives to unnormalized primitives
273 void convert_contractions();
274 /// Convert contraction on given shell
275 void convert_contraction(size_t ish);
276
277 /// Normalize contractions. If !coeffs, only cartesian factors are calculated.
278 void normalize(bool coeffs=true);
279 /// Normalize contractions in Coulomb norm (for density fitting)
280 void coulomb_normalize();
281
282 /// Do all of the above
283 void finalize(bool convert=false, bool normalize=true);
284
285 /// Get distance of nuclei
286 double nuclear_distance(size_t i, size_t j) const;
287 /// Get nuclear distances
288 arma::mat nuclear_distances() const;
289
290 /// Get angular momentum of shell
291 int get_am(size_t shind) const;
292 /// Get maximum angular momentum in basis set
293 int get_max_am() const;
294 /// Get maximum number of contractions
295 size_t get_max_Ncontr() const;
296
297 /// Get index of last function, throws an exception if no functions exist
298 size_t get_last_ind() const;
299 /// Get index of first function on shell
300 size_t get_first_ind(size_t shind) const;
301 /// Get index of last function on shell
302 size_t get_last_ind(size_t shind) const;
303
304 /// Get R^2 expectation value (measure of basis function extent)
305 arma::vec get_bf_Rsquared() const;
306 /// Get shell indices of basis functions
307 arma::uvec shell_indices() const;
308 /// Find shell index of basis function
309 size_t find_shell_ind(size_t find) const;
310
311 /// Get shells in basis set
312 std::vector<GaussianShell> get_shells() const;
313 /// Get ind:th shell
314 GaussianShell get_shell(size_t shind) const;
315 /// Get index of the center of the ind'th shell
316 size_t get_shell_center_ind(size_t shind) const;
317 /// Get coordinates of center of ind'th shell
318 coords_t get_shell_center(size_t shind) const;
319
320 /// Get exponential contraction of the ind:th shell
321 std::vector<contr_t> get_contr(size_t ind) const;
322 /// Get normalized exponential contraction of the ind:th shell
323 std::vector<contr_t> get_contr_normalized(size_t ind) const;
324 /// Get the cartesian functions on the ind:th shell
325 std::vector<shellf_t> get_cart(size_t ind) const;
326
327 /// Are spherical harmonics the default for new shells?
328 bool is_lm_default() const;
329 /// Is shell ind using spherical harmonics?
330 bool lm_in_use(size_t ind) const;
331 /// Toggle the use of spherical harmonics on shell ind
332 void set_lm(size_t ind, bool lm);
333
334 /// Get m values of basis functions
335 arma::ivec get_m_values() const;
336 /// Unique m values in basis set
337 arma::ivec unique_m_values() const;
338 /// Mapping from m value to unique m value
339 std::map<int, arma::uword> unique_m_map() const;
340 /// Count occupied orbitals
341 arma::imat count_m_occupied(const arma::mat & C) const;
342 /// Count occupied orbitals
343 arma::imat count_m_occupied(const arma::mat & Ca, const arma::mat & Cb) const;
344
345 /// Get indices of basis functions with wanted m value
346 arma::uvec m_indices(int m) const;
347
348 /// Get transformation matrix
349 arma::mat get_trans(size_t ind) const;
350
351 /// Fill spherical harmonics transformation table
352 void fill_sph_transmat();
353
354 /// Get size of basis set
355 size_t get_Nbf() const;
356 /// Get amount of cartesian functions on shell
357 size_t get_Ncart() const;
358 /// Get amount of spherical harmonics on shell
359 size_t get_Nlm() const;
360
361 /// Get number of shells
362 size_t get_Nshells() const;
363 /// Get number of basis functions on shell
364 size_t get_Nbf(size_t ind) const;
365 /// Get number of cartesians on shell
366 size_t get_Ncart(size_t ind) const;
367
368 /**
369 * Get range of shell (distance at which functions have dropped below epsilon)
370 *
371 * The default parameter is from the reference R. E. Stratmann,
372 * G. E. Scuseria, and M. J. Frisch, "Achieving linear scaling in
373 * exchage-correlation density functional quadratures",
374 * Chem. Phys. Lett. 257 (1993), pp. 213-223.
375 */
376 void compute_shell_ranges(double eps=1e-10);
377 /// Get precomputed ranges of shells
378 std::vector<double> get_shell_ranges() const;
379 /// Get range of shells with given value of epsilon
380 std::vector<double> get_shell_ranges(double eps) const;
381
382 /// Get distances to other nuclei
383 std::vector<double> get_nuclear_distances(size_t inuc) const;
384
385 /// Get number of nuclei
386 size_t get_Nnuc() const;
387 /// Get nucleus
388 nucleus_t get_nucleus(size_t inuc) const;
389 /// Get nuclei
390 std::vector<nucleus_t> get_nuclei() const;
391
392 /// Get coordinates of all nuclei
393 arma::mat get_nuclear_coords() const;
394 /// Set coordinates of all nuclei
395 void set_nuclear_coords(const arma::mat & coords);
396
397 /// Get coordinates of nucleus
398 coords_t get_nuclear_coords(size_t inuc) const;
399 /// Get charge of nucleus
400 int get_Z(size_t inuc) const;
401 /// Get symbol of nucleus
402 std::string get_symbol(size_t inuc) const;
403 /// Get human readable symbol of nucleus (-Bq)
404 std::string get_symbol_hr(size_t inuc) const;
405
406 /// Get basis functions centered on a given atom
407 std::vector<GaussianShell> get_funcs(size_t inuc) const;
408 /// Get indices of shells centered on a given atom
409 std::vector<size_t> get_shell_inds(size_t inuc) const;
410
411 /// Evaluate functions at (x,y,z)
412 arma::vec eval_func(double x, double y, double z) const;
413 /// Evaluate gradient at (x,y,z)
414 arma::mat eval_grad(double x, double y, double z) const;
415 /// Evaluate Hessian at (x,y,z)
416 arma::mat eval_hess(double x, double y, double z) const;
417
418 /// Evaluate functions of shell ish at (x,y,z)
419 arma::vec eval_func(size_t ish, double x, double y, double z) const;
420 /// Evaluate gradients of shell ish at (x,y,z)
421 arma::mat eval_grad(size_t ish, double x, double y, double z) const;
422 /// Evaluate laplacian of shell ish at (x,y,z)
423 arma::vec eval_lapl(size_t ish, double x, double y, double z) const;
424 /// Evaluate Hessian of shell ish at (x,y,z)
425 arma::mat eval_hess(size_t ish, double x, double y, double z) const;
426 /// Evaluate gradient of laplacian of shell ish at (x,y,z)
427 arma::mat eval_laplgrad(size_t ish, double x, double y, double z) const;
428
429 /// Print out basis set
430 void print(bool verbose=false) const;
431
432 /// Calculate transformation matrix from cartesians to spherical harmonics
433 arma::mat cart_to_sph_trans() const;
434 /// Calculate transfomration matrix from spherical harmonics to cartesians
435 arma::mat sph_to_cart_trans() const;
436
437 /// Calculate overlap matrix
438 arma::mat overlap() const;
439 /// Calculate overlap matrix in Coulomb metric
440 arma::mat coulomb_overlap() const;
441 /// Calculate overlap with another basis set
442 arma::mat overlap(const BasisSet & rhs) const;
443 /// Calculate overlap in Coulomb metric with another basis set
444 arma::mat coulomb_overlap(const BasisSet & rhs) const;
445 /// Calculate kinetic energy matrix
446 arma::mat kinetic() const;
447 /// Calculate nuclear repulsion matrix
448 arma::mat nuclear() const;
449 /// Calculate electric potential matrix
450 arma::mat potential(coords_t r) const;
451
452 /**
453 Calculates the ERI screening matrices
454 \f$ Q_{\mu \nu} = (\mu \nu | \mu \nu)^{1/2} \f$
455 and
456 \f$ M_{\mu \nu} = (\mu \mu | \nu \nu)^{1/2} \f$
457 as described in J. Chem. Phys. 147, 144101 (2017).
458 */
459 void eri_screening(arma::mat & Q, arma::mat & M, double omega=0.0, double alpha=1.0, double beta=0.0) const;
460
461 /// Calculate nuclear Pulay forces
462 arma::vec nuclear_pulay(const arma::mat & P) const;
463 /// Calculate nuclear Hellman-Feynman force
464 arma::vec nuclear_der(const arma::mat & P) const;
465 /// Calculate kinetic Pulay forces
466 arma::vec kinetic_pulay(const arma::mat & P) const;
467 /// Calculate overlap derivative force
468 arma::vec overlap_der(const arma::mat & W) const;
469 /// Calculate nuclear repulsion force
470 arma::vec nuclear_force() const;
471
472 /// Compute moment integral around (x,y,z)
473 std::vector<arma::mat> moment(int mom, double x=0.0, double y=0.0, double z=0.0) const;
474
475 /// Compute integrals of basis functions (used in xc-fitting)
476 arma::vec integral() const;
477
478 /// Calculate nuclear charge
479 int Ztot() const;
480
481 /// Calculate nuclear repulsion energy
482 double Enuc() const;
483
484 /// Project MOs to new basis set
485 void projectMOs(const BasisSet & old, const arma::colvec & oldE, const arma::mat & oldMOs, arma::colvec & E, arma::mat & MOs, size_t nocc) const;
486 /// Translate OMOs to new geometry, assuming same basis set. Virtuals are discarded and regenerated.
487 void projectOMOs(const BasisSet & old, const arma::cx_mat & oldOMOs, arma::cx_mat & OMOs, size_t nocc) const;
488
489 /// Are the basis sets the same?
490 bool operator==(const BasisSet & rhs) const;
491
492 /// Find "identical" shells in basis set.
493 std::vector< std::vector<size_t> > find_identical_shells() const;
494 };
495
496 /// Compute three-center overlap integral
497 arma::cube three_overlap(const GaussianShell *is, const GaussianShell *js, const GaussianShell *ks);
498
499 /**
500 * \class GaussianShell
501 *
502 * \brief A shell of Gaussian basis functions of a given angular
503 * momentum, sharing the same exponential contraction.
504 *
505 * \author Susi Lehtola
506 * \date 2011/05/05 20:17
507 */
508 class GaussianShell {
509 /// Number of first function on shell
510 size_t indstart;
511
512 /// Coordinates of center
513 coords_t cen;
514 /// Index of center
515 size_t cenind;
516
517 /// Use spherical harmonics?
518 bool uselm;
519 /// Transformation matrix to spherical basis
520 arma::mat transmat;
521
522 /**
523 * Contraction of unnormalized primitives.
524 * N.B. Normalization is wrt first function of shell.
525 */
526 std::vector<contr_t> c;
527
528 /// Angular momentum of shell
529 int am;
530
531 /**
532 * Table of cartesians, containing am indices
533 * and relative normalization factors.
534 */
535 std::vector<shellf_t> cart;
536
537 public:
538 /// Dummy constructor
539 GaussianShell();
540 /// Constructor, need also to set index of first function and nucleus (see below)
541 GaussianShell(int am, bool lm, const std::vector<contr_t> & C);
542 /// Destructor
543 ~GaussianShell();
544
545 /// Set index of first basis function
546 void set_first_ind(size_t ind);
547 /// Set center
548 void set_center(const coords_t & cenv, size_t cenindv);
549
550 /// Sort exponents in decreasing order
551 void sort();
552
553 /**
554 * Convert contraction from coefficients of normalized primitives to
555 * those of unnormalized primitives.
556 */
557 void convert_contraction();
558
559 /**
560 * Convert contraction from coefficients of normalized density
561 * primitives to those of unnormalized primitives.
562 */
563 void convert_sap_contraction();
564
565 /// Normalize contractions
566 void normalize(bool coeffs=true);
567 /// Normalize contractions in Coulomb norm (for density fitting)
568 void coulomb_normalize();
569
570 /// Get the exponential contraction
571 std::vector<contr_t> get_contr() const;
572 /// Get cartesians
573 std::vector<shellf_t> get_cart() const;
574
575 /**
576 * Get contraction coefficients of normalized primitives. For some
577 * reason these are at variance with what is put in, e.g. the
578 * cc-pVXZ basis set data from the ESML basis set exchange. Maybe
579 * the input isn't really normalized..?
580 */
581 std::vector<contr_t> get_contr_normalized() const;
582
583 /// Number of functions on shell
584 size_t get_Nbf() const;
585 /// Number of cartesians on shell
586 size_t get_Ncart() const;
587 /// Number of spherical harmonics on shell
588 size_t get_Nlm() const;
589
590 /**
591 * Compute range of shell - how far must one go for absolute value
592 * of functions to drop below eps?
593 */
594 double range(double eps) const;
595
596 /// Are spherical harmonics in use?
597 bool lm_in_use() const;
598 /// Toggle use of spherical harmonics
599 void set_lm(bool lm);
600 /// Get transformation matrix to spherical harmonics basis
601 arma::mat get_trans() const;
602
603 /// Get number of contractions
604 size_t get_Ncontr() const;
605 /// Get angular momentum
606 int get_am() const;
607 /// Get nucleus index
608 size_t get_center_ind() const;
609 /// Get coordinates
610 coords_t get_center() const;
611
612 /// Comparison operator for angular momentum ordering
613 bool operator<(const GaussianShell & rhs) const;
614 /// Are the two the same?
615 bool operator==(const GaussianShell & rhs) const;
616
617 /// Get index of first function on shell
618 size_t get_first_ind() const;
619 /// Get index of last function on shell
620 size_t get_last_ind() const;
621
622 /// Print out information about shell
623 void print() const;
624
625 /// Evaluate functions at (x,y,z)
626 arma::vec eval_func(double x, double y, double z) const;
627 /// Evaluate gradients at (x,y,z)
628 arma::mat eval_grad(double x, double y, double z) const;
629 /// Evaluate laplacian at (x,y,z)
630 arma::vec eval_lapl(double x, double y, double z) const;
631 /// Evaluate Hessian at (x,y,z)
632 arma::mat eval_hess(double x, double y, double z) const;
633 /// Evaluate gradient of laplacian at (x,y,z)
634 arma::mat eval_laplgrad(double x, double y, double z) const;
635
636 /// Calculate block overlap matrix between shells
637 arma::mat overlap(const GaussianShell & rhs) const;
638 /// Calculate block Coulomb overlap matrix between shells
639 arma::mat coulomb_overlap(const GaussianShell & rhs) const;
640 /// Calculate kinetic energy matrix between shells
641 arma::mat kinetic(const GaussianShell & rhs) const;
642 /// Calculate nuclear repulsion matrix between shells
643 arma::mat nuclear(double cx, double cy, double cz, const GaussianShell & rhs) const;
644
645 /// Calculate nuclear Pulay forces
646 arma::vec nuclear_pulay(double cx, double cy, double cz, const arma::mat & P, const GaussianShell & rhs) const;
647 /// Calculate nuclear Hellman-Feynman force
648 arma::vec nuclear_der(double cx, double cy, double cz, const arma::mat & P, const GaussianShell & rhs) const;
649 /// Calculate kinetic Pulay forces
650 arma::vec kinetic_pulay(const arma::mat & P, const GaussianShell & rhs) const;
651 /// Calculate overlap derivative
652 arma::vec overlap_der(const arma::mat & W, const GaussianShell & rhs) const;
653
654 /// Calculate integral over function (used in xc-fitting)
655 arma::vec integral() const;
656
657 /// Calculate moment integrals around (x,y,z) between shells
658 std::vector<arma::mat> moment(int mom, double x, double y, double z, const GaussianShell & rhs) const;
659 };
660
661 /// Get dummy shell
662 GaussianShell dummyshell();
663
664 /// Form index helper table: i*(i+1)/2
665 std::vector<size_t> i_idx(size_t N);
666
667 /// Construct basis set from input
668 void construct_basis(BasisSet & basis, const std::vector<atom_t> & atoms, const BasisSetLibrary & baslib);
669 /// Constuct basis set from input
670 void construct_basis(BasisSet & basis, const std::vector<nucleus_t> & nuclei, const BasisSetLibrary & baslib);
671
672 /// Compute values of orbitals at given point
673 arma::vec compute_orbitals(const arma::mat & C, const BasisSet & bas, const coords_t & r);
674 /// Compute Fermi-Löwdin orbitals
675 arma::mat fermi_lowdin_orbitals(const arma::mat & C, const BasisSet & bas, const arma::mat & r);
676
677 /// Compute density at given point
678 double compute_density(const arma::mat & P, const BasisSet & bas, const coords_t & r);
679 /// Compute density and gradient at a given point
680 void compute_density_gradient(const arma::mat & P, const BasisSet & bas, const coords_t & r, double & d, arma::vec & g);
681 /// Compute density, gradient and hessian at a given point
682 void compute_density_gradient_hessian(const arma::mat & P, const BasisSet & bas, const coords_t & r, double & d, arma::vec & g, arma::mat & h);
683
684 /// Compute electrostatic potential at given point
685 double compute_potential(const arma::mat & P, const BasisSet & bas, const coords_t & r);
686
687 /**
688 * Compute electron localization function
689 *
690 * A. D. Becke and K. E. Edgecombe, "A simple measure of electron
691 * localization in atomic and molecular systems", J. Chem. Phys. 92,
692 * 5397 (1990).
693 */
694 double compute_elf(const arma::mat & P, const BasisSet & bas, const coords_t & r);
695
696 /// Calculate difference from orthogonality
697 double orth_diff(const arma::mat & C, const arma::mat & S);
698 /// Check orthonormality of complex molecular orbitals
699 double orth_diff(const arma::cx_mat & C, const arma::mat & S);
700
701 /// Check orthonormality of real molecular orbitals
702 void check_orth(const arma::mat & C, const arma::mat & S, bool verbose, double thr=sqrt(DBL_EPSILON));
703 /// Check orthonormality of complex molecular orbitals
704 void check_orth(const arma::cx_mat & C, const arma::mat & S, bool verbose, double thr=sqrt(DBL_EPSILON));
705
706 /**
707 * Construct intrinsic atomic orbitals.
708 *
709 * G. Knizia, "Intrinsic Atomic Orbitals: An Unbiased Bridge between
710 * Quantum Theory and Chemical Concepts", J. Chem. Theory
711 * Comput. 9, 4834 (2013). doi: 10.1021/ct400687b.
712 *
713 * The algorithm returns the IAO matrix, and stores the atomic indices
714 * in idx.
715 */
716 arma::mat construct_IAO(const BasisSet & basis, const arma::mat & C, std::vector< std::vector<size_t> > & idx, bool verbose=true, std::string minbas="MINAO.gbs");
717 /// Same, but using complex coefficients
718 arma::cx_mat construct_IAO(const BasisSet & basis, const arma::cx_mat & C, std::vector< std::vector<size_t> > & idx, bool verbose=true, std::string minbas="MINAO.gbs");
719
720 /// Block matrix by m value
721 arma::mat block_m(const arma::mat & F, const arma::ivec & mv);
722 /// Get norm by m value block
723 arma::mat m_norm(const arma::mat & C, const arma::ivec & mv);
724 /// Classify by m value
725 arma::ivec m_classify(const arma::mat & C, const arma::ivec & mv);
726
727 #endif
728