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 9 * Copyright (c) 2010, 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 18 #include "global.h" 19 20 #ifndef ERKALE_BASISLIB 21 #define ERKALE_BASISLIB 22 23 #include <armadillo> 24 #include <vector> 25 #include <string> 26 27 #include "basis.h" 28 29 /// Compute overlap of normalized Gaussian primitives 30 arma::mat overlap(const arma::vec & z, const arma::vec & zp, int am); 31 32 /// Find angular momentum 33 int find_am(char am); 34 35 /// Find basis set file 36 std::string find_basis(const std::string & filename, bool verbose=true); 37 38 /// System default location for basis sets 39 #ifndef ERKALE_SYSTEM_LIBRARY 40 #define ERKALE_SYSTEM_LIBRARY "/usr/share/erkale/basis" 41 #endif 42 43 /** 44 * \class FunctionShell 45 * 46 * \brief A shell of functions 47 * 48 * This class defines a shell of functions of the same angular 49 * momentum, used in the ElementBasisSet class. 50 * 51 * \author Susi Lehtola 52 * \date 2011/05/08 17:10 53 */ 54 55 class FunctionShell { 56 /// Angular momentum 57 int am; 58 /// Exponential contraction 59 std::vector<contr_t> C; 60 61 public: 62 /// Construct a shell with angular momentum am 63 FunctionShell(int am=-1); 64 /// Construct a shell with angular momentum am and given contraction 65 FunctionShell(int am, const std::vector<contr_t> & c); 66 /// Destructor 67 ~FunctionShell(); 68 69 /// Add exponent into contraction 70 void add_exponent(double C, double z); 71 /// Comparison operator for ordering in decreasing angular momentum and exponent 72 bool operator<(const FunctionShell &rhs) const; 73 /// Comparison operator for identity 74 bool operator==(const FunctionShell &rhs) const; 75 76 /// Get angular momentum 77 int get_am() const; 78 /// Get contraction length 79 size_t get_Ncontr() const; 80 81 /// Get contraction coefficients 82 std::vector<contr_t> get_contr() const; 83 84 /// Sort exponents in decreasing order 85 void sort(); 86 /// Normalize coefficients 87 void normalize(); 88 /// Print out info 89 void print() const; 90 91 friend class BasisSetLibrary; 92 }; 93 94 /** 95 * \class ElementBasisSet 96 * 97 * \brief Basis set for an element 98 * 99 * This class defines a basis set for an element, used by the 100 * BasisSetLibrary class. 101 * 102 * \author Susi Lehtola 103 * \date 2011/05/08 17:10 104 */ 105 106 class ElementBasisSet { 107 108 /// Symbol of element 109 std::string symbol; 110 /// Atom id for which the basis is for (0 for all atoms) 111 size_t number; 112 /// List of shells 113 std::vector<FunctionShell> bf; 114 115 public: 116 /// Dummy constructor 117 ElementBasisSet(); 118 /// Constructor 119 ElementBasisSet(std::string sym, size_t number=0); 120 /// Destructor 121 ~ElementBasisSet(); 122 123 /// Add a shell of functions to the basis 124 void add_function(FunctionShell f); 125 /// Sort the shells in decreasing angular momentum 126 void sort(); 127 /// Print out basis set 128 void print() const; 129 130 /// Get the symbol of the element 131 std::string get_symbol() const; 132 /// Get the number 133 size_t get_number() const; 134 /// Set the number 135 void set_number(size_t num); 136 137 /// Comparison operator for sorting 138 bool operator<(const ElementBasisSet &rhs) const; 139 140 /// Get the shells 141 std::vector<FunctionShell> get_shells() const; 142 /// Get the shells for am value 143 std::vector<FunctionShell> get_shells(int am) const; 144 145 /// Get exponents and contraction coefficients of angular momentum shell am 146 void get_primitives(arma::vec & exps, arma::mat & coeffs, int am) const; 147 /// Get free primitives and generally contracted part 148 void get_primitives(arma::vec & zfree, arma::vec & zgen, arma::mat & cgen, int am) const; 149 150 /// Orthonormalize generally contracted functions. NB! This can screw up your computational efficiency! 151 void orthonormalize(); 152 153 /** 154 * P-orthogonalization of basis set. 155 * A refinement on the Davidson refinement (a.k.a. "Davidson rotation"). 156 * 157 * The procedure is based on the article 158 * 159 * F. Jensen, "Unifying General and Segmented Contracted Basis 160 * Sets. Segmented Polarization Consistent Basis Sets", 161 * J. Chem. Theory Comput. 10, 1074 (2014). 162 */ 163 void P_orthogonalize(double cutoff=1e-8, double Cortho=1e-4); 164 165 /// Prune primitives with large overlap, replacing them with the geometric average 166 void prune(double cutoff=0.9, bool coulomb=false); 167 168 /// Merge primitives with large overlap (also decontracts basis). This routine is intended to merging core-correlation functions 169 void merge(double cutoff=0.9, bool verbose=true, bool coulomb=false); 170 171 /// Get maximum angular momentum used in the shells 172 int get_max_am() const; 173 /// Get angular momentum of i:th shell 174 int get_am(size_t ind) const; 175 /// Get amount of functions 176 int get_Nbf() const; 177 /// Get maximum contraction depth 178 size_t get_max_Ncontr() const; 179 180 /// Normalize coefficients 181 void normalize(); 182 183 /// Decontract set 184 void decontract(); 185 /** 186 * Generate density fitting basis 187 * 188 * The procedure has been documented in the article 189 * 190 * R. Yang, A. P. Rendell and M. J. Frisch, "Automatically generated 191 * Coulomb fitting basis sets: Design and accuracy for systems 192 * containing H to Kr", J. Chem. Phys. 127 (2007), 074102. 193 */ 194 ElementBasisSet density_fitting(int lmaxinc, double fsam) const; 195 /** 196 * Generate product basis set. 197 * 198 * Brute-force generation of orbital products, followed by merging 199 * product exponents that deviate by fsam at maximum. 200 */ 201 ElementBasisSet product_set(int lmaxinc, double fsam) const; 202 203 /// Form compact Cholesky set 204 ElementBasisSet cholesky_set(double thr, int maxam, double ovlthr) const; 205 206 /// Augment the basis with steep functions 207 void augment_steep(int naug); 208 /// Augment the basis with diffuse functions 209 void augment_diffuse(int naug); 210 211 friend class BasisSetLibrary; 212 }; 213 214 /// Helper for P-orthogonalization - compute inner product for inside out purification 215 double P_innerprod_inout(const arma::vec & ai, const arma::mat & S, const arma::vec & aj, size_t P); 216 /// Helper for P-orthogonalization - compute inner product for outside in purification 217 double P_innerprod_outin(const arma::vec & ai, const arma::mat & S, const arma::vec & aj, size_t P); 218 219 /// Helper for P-orthogonalization - count amount of shared exponents 220 size_t count_shared(const arma::vec & ai, const arma::vec & aj); 221 /// Helper for P-orthogonalization: check if functions have already been treated 222 bool treated_inout(const arma::mat & c, size_t i, size_t j); 223 /// Helper for P-orthogonalization: check if functions have already been treated 224 bool treated_outin(const arma::mat & c, size_t i, size_t j); 225 226 /** 227 * \class BasisSetLibrary 228 * 229 * \brief Basis set library class 230 * 231 * This class defines a basis set library class that can be used 232 * e.g. to read in basis sets from files, or to save basis sets to 233 * files. 234 * 235 * \author Susi Lehtola 236 * \date 2011/05/08 17:10 237 */ 238 239 class BasisSetLibrary { 240 241 /// Name of basis set 242 std::string name; 243 /// List of elements included in basis set 244 std::vector<ElementBasisSet> elements; 245 246 public: 247 /// Constructor 248 BasisSetLibrary(); 249 /// Destructor 250 ~BasisSetLibrary(); 251 252 /// Load basis set from file in Gaussian'94 format. Handles parsing for special Pople sets 253 void load_basis(const std::string & filename, bool verbose=true); 254 255 /// Load basis set from file in Gaussian'94 format 256 void load_gaussian94(const std::string & filename, bool verbose=true); 257 /// Save basis set to file in Gaussian'94 format 258 void save_gaussian94(const std::string & filename, bool append=false) const; 259 260 /// Save basis set to file in CFOUR format 261 void save_cfour(const std::string & filename, const std::string & basisname, bool newformat=true, bool append=false) const; 262 /// Save basis set to file in Dalton format 263 void save_dalton(const std::string & filename, bool append=false) const; 264 /// Save basis set to file in Molpro format 265 void save_molpro(const std::string & filename, bool append=false) const; 266 267 /// Add element to basis set 268 void add_element(const ElementBasisSet & el); 269 270 /// Sort library 271 void sort(); 272 273 /// Get number of elements 274 size_t get_Nel() const; 275 /// Get symbol of ind'th element 276 std::string get_symbol(size_t ind) const; 277 /// Get elements 278 std::vector<ElementBasisSet> get_elements() const; 279 280 /// Get maximum angular momentum used in basis set 281 int get_max_am() const; 282 /// Get maximum contraction depth 283 size_t get_max_Ncontr() const; 284 285 /// Get basis set for wanted element 286 ElementBasisSet get_element(std::string el, size_t number=0) const; 287 288 /// Normalize coefficients 289 void normalize(); 290 291 /// Orthonormalize generally contracted functions. NB! This can screw up your computational efficiency! 292 void orthonormalize(); 293 294 /// Decontract basis set 295 void decontract(); 296 297 /// Generate density fitting set 298 BasisSetLibrary density_fitting(int lvalinc, double fsam) const; 299 /// Generate product set 300 BasisSetLibrary product_set(int lvalinc, double fsam) const; 301 /// Generate compact Cholesky set 302 BasisSetLibrary cholesky_set(double thr, int maxam, double ovlthr) const; 303 304 /** 305 * P-orthogonalization [F. Jensen, JCTC 10, 1074 (2014)]. 306 * 307 * The cutoff drops out primitives with coefficients smaller than 308 * the threshold in the intermediate normalization (largest 309 * coefficient is set to unity). 310 */ 311 void P_orthogonalize(double cutoff=1e-8, double Cortho=1e-4); 312 313 /// Merge primitives with large overlap (also decontracts basis) 314 void merge(double cutoff=0.9, bool verbose=true); 315 316 /// Augment the basis with steep functions 317 void augment_steep(int naug); 318 /// Augment the basis with diffuse functions 319 void augment_diffuse(int naug); 320 321 /// Print out library 322 void print() const; 323 }; 324 325 326 #endif 327