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