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_MOLECULE_H
33 #define MADNESS_MOLECULE_H
34 
35 /// \file moldft/molecule.h
36 /// \brief Declaration of molecule related classes and functions
37 
38 #include <polar/corepotential.h>
39 #include <polar/atomutil.h>
40 #include <madness/world/vector.h>
41 #include <vector>
42 #include <string>
43 #include <iostream>
44 #include <fstream>
45 #include <sstream>
46 #include <algorithm>
47 #include <ctype.h>
48 #include <cmath>
49 #include <madness/tensor/tensor.h>
50 #include <madness/misc/misc.h>
51 
52 
53 class Atom {
54 public:
55     double x, y, z, q;          ///< Coordinates and charge in atomic units
56     unsigned int atomic_number; ///< Atomic number
57 
Atom(double x,double y,double z,double q,unsigned int atomic_number)58     explicit Atom(double x, double y, double z, double q, unsigned int atomic_number)
59             : x(x), y(y), z(z), q(q), atomic_number(atomic_number) {}
60 
Atom(const Atom & a)61     Atom(const Atom& a)
62             : x(a.x), y(a.y), z(a.z), q(a.q), atomic_number(a.atomic_number) {}
63 
64     /// Default construct makes a zero charge ghost atom at origin
Atom()65     Atom()
66             : x(0), y(0), z(0), q(0), atomic_number(0) {}
67 
68 
get_coords()69     madness::Vector<double,3> get_coords() const {
70         return madness::Vector<double,3>{x, y, z};
71     }
72 
73     template <typename Archive>
serialize(Archive & ar)74     void serialize(Archive& ar) {
75         ar & x & y & z & q & atomic_number;
76     }
77 };
78 
79 std::ostream& operator<<(std::ostream& s, const Atom& atom);
80 
81 class Molecule {
82 private:
83     // If you add more fields don't forget to serialize them
84     std::vector<Atom> atoms;
85     std::vector<double> rcut;  // Reciprocal of the smoothing radius
86     double eprec;              // Error in energy/atom due to smoothing
87     CorePotentialManager core_pot;
88     madness::Tensor<double> field;
89 
90     void swapaxes(int ix, int iy);
91 
92     template <typename opT>
93     bool test_for_op(double xaxis, double yaxis, double zaxis, opT op) const;
94 
95     bool test_for_c2(double xaxis, double yaxis, double zaxis) const;
96 
97     bool test_for_sigma(double xaxis, double yaxis, double zaxis) const;
98 
99     bool test_for_inverse() const;
100 
101 public:
102     /// Makes a molecule with zero atoms
Molecule()103     Molecule() : atoms(), rcut(), eprec(1e-4), core_pot(), field(3L) {};
104 
105     Molecule(const std::string& filename);
106 
107     void read_file(const std::string& filename);
108 
109     void read_core_file(const std::string& filename);
110 
guess_file()111     std::string guess_file() const { return core_pot.guess_file(); };
112 
113     unsigned int n_core_orb_all() const ;
114 
n_core_orb(unsigned int atn)115     unsigned int n_core_orb(unsigned int atn) const {
116         if (core_pot.is_defined(atn))
117             return core_pot.n_core_orb_base(atn);
118         else
119             return 0;
120     };
121 
get_core_l(unsigned int atn,unsigned int c)122     unsigned int get_core_l(unsigned int atn, unsigned int c) const {
123         return core_pot.get_core_l(atn, c);
124     }
125 
get_core_bc(unsigned int atn,unsigned int c)126     double get_core_bc(unsigned int atn, unsigned int c) const {
127         return core_pot.get_core_bc(atn, c);
128     }
129 
130     double core_eval(int atom, unsigned int core, int m, double x, double y, double z) const;
131 
132     double core_derivative(int atom, int axis, unsigned int core, int m, double x, double y, double z) const;
133 
is_potential_defined(unsigned int atn)134     bool is_potential_defined(unsigned int atn) const { return core_pot.is_defined(atn); };
135 
is_potential_defined_atom(int i)136     bool is_potential_defined_atom(int i) const { return core_pot.is_defined(atoms[i].atomic_number); };
137 
138     void add_atom(double x, double y, double z,  double q, int atn);
139 
natom()140     int natom() const {
141         return atoms.size();
142     };
143 
144     void set_atom_coords(unsigned int i, double x, double y, double z);
145 
146     madness::Tensor<double> get_all_coords() const;
147 
148     std::vector< madness::Vector<double,3> > get_all_coords_vec() const;
149 
150     std::vector<double> atomic_radii;
151 
152     void set_all_coords(const madness::Tensor<double>& newcoords);
153 
154     void set_eprec(double value);
155 
156     void set_rcut(double value);
157 
set_core_eprec(double value)158     void set_core_eprec(double value) {
159         core_pot.set_eprec(value);
160     }
161 
set_core_rcut(double value)162     void set_core_rcut(double value) {
163         core_pot.set_rcut(value);
164     }
165 
get_eprec()166     double get_eprec() const {
167         return eprec;
168     }
169 
170     double bounding_cube() const;
171 
172     const Atom& get_atom(unsigned int i) const;
173 
174     void print() const;
175 
176     double inter_atomic_distance(unsigned int i,unsigned int j) const;
177 
178     double nuclear_repulsion_energy() const;
179 
180     double nuclear_repulsion_derivative(int i, int j) const;
181 
182     double nuclear_dipole(int axis) const;
183 
184     double nuclear_charge_density(double x, double y, double z) const;
185 
186     double mol_nuclear_charge_density(double x, double y, double z) const;
187 
188     double smallest_length_scale() const;
189 
190     void identify_point_group();
191 
192     void center();
193 
194     void orient();
195 
196     double total_nuclear_charge() const;
197 
198     double nuclear_attraction_potential(double x, double y, double z) const;
199 
200     double molecular_core_potential(double x, double y, double z) const;
201 
202     double core_potential_derivative(int atom, int axis, double x, double y, double z) const;
203 
204     double nuclear_attraction_potential_derivative(int atom, int axis, double x, double y, double z) const;
205 
206     template <typename Archive>
serialize(Archive & ar)207     void serialize(Archive& ar) {
208         ar & atoms & rcut & eprec & core_pot;
209     }
210 };
211 
212 
213 #endif
214