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, 2013
9  * Copyright (c) 2013, 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_BADERGRID
18 #define ERKALE_BADERGRID
19 
20 #include "global.h"
21 #include "basis.h"
22 #include "dftgrid.h"
23 
24 /**
25  * The algorithms here are original work, but actually are very much
26  * similar to the work in
27  *
28  * J. I. Rodríguez, A. M. Köster, P. W. Ayers, A. Santos-Valle,
29  * A. Vela, and G. Merino, "An efficient grid-based scheme to compute
30  * QTAIM atomic properties without explicit calculation of zero-flux
31  * surfaces", J. Comp. Chem. 30, 1082 (2009).
32  */
33 
34 /**
35  * Integration over Bader regions
36  */
37 class BaderGrid {
38   /// Basis set
39   const BasisSet *basp;
40   /// Grid worker
41   AngularGrid wrk;
42 
43   /// Locations of maxima
44   std::vector<coords_t> maxima;
45   /// Grid points corresponding to the regions
46   std::vector< std::vector<gridpoint_t> > reggrid;
47   /// Amount of nuclei
48   size_t Nnuc;
49 
50   /// Verbose operation?
51   bool verbose;
52   /// Print maxima
53   void print_maxima() const;
54 
55  public:
56   /// Constructor
57   BaderGrid();
58   /// Destructor
59   ~BaderGrid();
60 
61   /// Set parameters
62   void set(const BasisSet & basis, bool verbose=true, bool lobatto=false);
63 
64   /// Construct grid with AO overlap matrix threshold thr and classify points into regions with P
65   void construct_bader(const arma::mat & P, double thr);
66   /// Construct grid with AO overlap matrix threshold thr and classify points into Voronoi regions
67   void construct_voronoi(double tol);
68   /// Get amount of regions
69   size_t get_Nmax() const;
70 
71   /// Compute regional charges
72   arma::vec regional_charges(const arma::mat & P);
73   /// Compute nuclear charges
74   arma::vec nuclear_charges(const arma::mat & P);
75 
76   /// Compute regional overlap matrices
77   std::vector<arma::mat> regional_overlap();
78   /// Compute regional overlap matrix
79   arma::mat regional_overlap(size_t ireg);
80 };
81 
82 /// Track point to maximum
83 coords_t track_to_maximum(const BasisSet & basis, const arma::mat & P, const coords_t r0, size_t & nd, size_t & ng);
84 
85 
86 #endif
87