1 /****************************************************************************** 2 3 This source file is part of the Avogadro project. 4 5 Copyright 2008-2009 Marcus D. Hanwell 6 Copyright 2008 Albert De Fusco 7 Copyright 2010-2013 Kitware, Inc. 8 9 This source code is released under the New BSD License, (the "License"). 10 11 Unless required by applicable law or agreed to in writing, software 12 distributed under the License is distributed on an "AS IS" BASIS, 13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 See the License for the specific language governing permissions and 15 limitations under the License. 16 17 ******************************************************************************/ 18 19 #ifndef AVOGADRO_CORE_GAUSSIANSET_H 20 #define AVOGADRO_CORE_GAUSSIANSET_H 21 22 #include "basisset.h" 23 24 #include <avogadro/core/matrix.h> 25 #include <avogadro/core/vector.h> 26 27 #include <vector> 28 29 namespace Avogadro { 30 namespace Core { 31 32 /** 33 * Enumeration of the SCF type. 34 */ 35 enum ScfType 36 { 37 Rhf, 38 Uhf, 39 Rohf, 40 Unknown 41 }; 42 43 /** 44 * @class GaussianSet gaussianset.h <avogadro/core/gaussianset.h> 45 * @brief A container for Gaussian type outputs from QM codes. 46 * @author Marcus D. Hanwell 47 * 48 * The GaussianSet class has a transparent data structure for storing the basis 49 * sets output by many quantum mechanical codes. It has a certain hierarchy 50 * where shells are built up from n primitives, in this case Gaussian Type 51 * Orbitals (GTOs). Each shell has a type (S, P, D, F, etc) and is composed of 52 * one or more GTOs. Each GTO has a contraction coefficient, c, and an exponent, 53 * a. 54 * 55 * When calculating Molecular Orbitals (MOs) each orthogonal shell has an 56 * independent coefficient. That is the S type orbitals have one coefficient, 57 * the P type orbitals have three coefficients (Px, Py and Pz), the D type 58 * orbitals have five (or six if cartesian types) coefficients, and so on. 59 */ 60 61 class AVOGADROCORE_EXPORT GaussianSet : public BasisSet 62 { 63 public: 64 /** 65 * Constructor. 66 */ 67 GaussianSet(); 68 69 /** 70 * Destructor. 71 */ 72 ~GaussianSet() override; 73 74 /** 75 * Clone. 76 */ clone()77 GaussianSet* clone() const override { return new GaussianSet(*this); } 78 79 /** 80 * Enumeration of the Gaussian type orbitals. 81 */ 82 enum orbital 83 { 84 S, 85 SP, 86 P, 87 D, 88 D5, 89 F, 90 F7, 91 G, 92 G9, 93 H, 94 H11, 95 I, 96 I13, 97 UU 98 }; 99 100 /** 101 * Add a basis to the basis set. 102 * @param atom Index of the atom to add the Basis to. 103 * @param type The type of the Basis being added. 104 * @return The index of the added Basis. 105 */ 106 unsigned int addBasis(unsigned int atom, orbital type); 107 108 /** 109 * Add a GTO to the supplied basis. 110 * @param basis The index of the Basis to add the GTO to. 111 * @param c The contraction coefficient of the GTO. 112 * @param a The exponent of the GTO. 113 * @return The index of the added GTO. 114 */ 115 unsigned int addGto(unsigned int basis, double c, double a); 116 117 /** 118 * Set the molecular orbital (MO) coefficients to the GaussianSet. 119 * @param MOs Vector containing the MO coefficients for the GaussianSet. 120 * @param type The type of the MOs (Paired, Alpha, Beta). 121 */ 122 void setMolecularOrbitals(const std::vector<double>& MOs, 123 ElectronType type = Paired); 124 125 /** 126 * Set the molecular orbital (MO) coefficients for a given index. Note 127 * that this must be used with coordinate sets to work correctly. 128 * @param MOs Vector containing the MO coefficients for the GaussianSet. 129 * @param type The type of the MOs (Paired, Alpha, Beta). 130 * @param index The index of the MO in the sequence. 131 */ 132 void setMolecularOrbitals(const std::vector<double>& MOs, ElectronType type, 133 Index index); 134 135 /** 136 * Get the number of elements in the set. 137 */ setCount()138 int setCount() { return static_cast<int>(m_moMatrixSet[0].size()); } 139 140 /** 141 * Set the active element in the set, this expects a corresponding 142 * coordinate set element, and will change the active MO matrix. 143 */ 144 bool setActiveSetStep(int index); 145 146 /** 147 * @brief Set the molecular orbital energies, expected in Hartrees. 148 * @param energies The vector containing energies for the MOs of type 149 * @param type The type of the electrons being set. 150 */ 151 void setMolecularOrbitalEnergy(const std::vector<double>& energies, 152 ElectronType type = Paired); 153 154 /** 155 * @brief Set the molecular orbital occupancies. 156 * @param occ The occupancies for the MOs of type. 157 * @param type The type of the electrons being set. 158 */ 159 void setMolecularOrbitalOccupancy(const std::vector<unsigned char>& occ, 160 ElectronType type = Paired); 161 162 /** 163 * @brief This enables support of sparse orbital sets, and provides a mapping 164 * from the index in memory to the actual molecular orbital number. 165 * @param nums The MO numbers (starting with an index of 1 for the first one). 166 * @param type The MO type (Paired, Alpha, Beta). 167 */ 168 void setMolecularOrbitalNumber(const std::vector<unsigned int>& nums, 169 ElectronType type = Paired); 170 171 /** 172 * Set the SCF density matrix for the GaussianSet. 173 */ 174 bool setDensityMatrix(const MatrixX& m); 175 176 /** 177 * Set the spin density matrix for the GaussianSet. 178 */ 179 bool setSpinDensityMatrix(const MatrixX& m); 180 181 /** 182 * @brief Generate the density matrix if we have the required information. 183 * @return True on success, false on failure. 184 */ 185 bool generateDensityMatrix(); 186 187 /** 188 * @return The number of molecular orbitals in the GaussianSet. 189 */ 190 unsigned int molecularOrbitalCount(ElectronType type = Paired) override; 191 192 /** 193 * Debug routine, outputs all of the data in the GaussianSet. 194 * @param type The electrons to output the information for. 195 */ 196 void outputAll(ElectronType type = Paired); 197 198 /** 199 * @return True of the basis set is valid, false otherwise. 200 * Default is true, if false then the basis set is likely unusable. 201 */ 202 bool isValid() override; 203 204 /** 205 * Set the SCF type for the object. 206 */ setScfType(ScfType type)207 void setScfType(ScfType type) { m_scfType = type; } 208 209 /** 210 * Get the SCF type for the object. 211 */ scfType()212 ScfType scfType() const { return m_scfType; } 213 214 /** 215 * Set the functional name (if applicable). 216 */ setFunctionalName(const std::string & name)217 void setFunctionalName(const std::string& name) { m_functionalName = name; } 218 219 /** 220 * Get the functional name (empty if none used). 221 */ functionalName()222 std::string functionalName() const { return m_functionalName; } 223 224 /** 225 * Initialize the calculation, this must normally be done before anything. 226 */ 227 void initCalculation(); 228 229 /** 230 * Accessors for the various properties of the GaussianSet. 231 */ symmetry()232 std::vector<int>& symmetry() { return m_symmetry; } symmetry()233 std::vector<int> symmetry() const { return m_symmetry; } atomIndices()234 std::vector<unsigned int>& atomIndices() { return m_atomIndices; } atomIndices()235 std::vector<unsigned int> atomIndices() const { return m_atomIndices; } moIndices()236 std::vector<unsigned int>& moIndices() { return m_moIndices; } moIndices()237 std::vector<unsigned int> moIndices() const { return m_moIndices; } gtoIndices()238 std::vector<unsigned int>& gtoIndices() { return m_gtoIndices; } gtoIndices()239 std::vector<unsigned int> gtoIndices() const { return m_gtoIndices; } cIndices()240 std::vector<unsigned int>& cIndices() { return m_cIndices; } cIndices()241 std::vector<unsigned int> cIndices() const { return m_cIndices; } gtoA()242 std::vector<double>& gtoA() { return m_gtoA; } gtoA()243 std::vector<double> gtoA() const { return m_gtoA; } gtoC()244 std::vector<double>& gtoC() { return m_gtoC; } gtoC()245 std::vector<double> gtoC() const { return m_gtoC; } gtoCN()246 std::vector<double>& gtoCN() 247 { 248 initCalculation(); 249 return m_gtoCN; 250 } 251 252 MatrixX& moMatrix(ElectronType type = Paired) 253 { 254 if (type == Paired || type == Alpha) 255 return m_moMatrix[0]; 256 else 257 return m_moMatrix[1]; 258 } 259 260 MatrixX moMatrix(ElectronType type = Paired) const 261 { 262 if (type == Paired || type == Alpha) 263 return m_moMatrix[0]; 264 else 265 return m_moMatrix[1]; 266 } 267 268 std::vector<double>& moEnergy(ElectronType type = Paired) 269 { 270 if (type == Paired || type == Alpha) 271 return m_moEnergy[0]; 272 else 273 return m_moEnergy[1]; 274 } 275 276 std::vector<double> moEnergy(ElectronType type = Paired) const 277 { 278 if (type == Paired || type == Alpha) 279 return m_moEnergy[0]; 280 else 281 return m_moEnergy[1]; 282 } 283 284 std::vector<unsigned char>& moOccupancy(ElectronType type = Paired) 285 { 286 if (type == Paired || type == Alpha) 287 return m_moOccupancy[0]; 288 else 289 return m_moOccupancy[1]; 290 } 291 292 std::vector<unsigned char> moOccupancy(ElectronType type = Paired) const 293 { 294 if (type == Paired || type == Alpha) 295 return m_moOccupancy[0]; 296 else 297 return m_moOccupancy[1]; 298 } 299 300 std::vector<unsigned int>& moNumber(ElectronType type = Paired) 301 { 302 if (type == Paired || type == Alpha) 303 return m_moNumber[0]; 304 else 305 return m_moNumber[1]; 306 } 307 308 std::vector<unsigned int> moNumber(ElectronType type = Paired) const 309 { 310 if (type == Paired || type == Alpha) 311 return m_moNumber[0]; 312 else 313 return m_moNumber[1]; 314 } 315 densityMatrix()316 MatrixX& densityMatrix() { return m_density; } spinDensityMatrix()317 MatrixX& spinDensityMatrix() { return m_spinDensity; } 318 319 private: 320 /** 321 * @brief This group is used once, and refers to the entire molecule. 322 */ 323 std::vector<int> m_symmetry; //! Symmetry of the basis, S, P... 324 std::vector<unsigned int> m_atomIndices; //! Indices into the atomPos vector 325 std::vector<unsigned int> m_moIndices; //! Indices into the MO/density matrix 326 std::vector<unsigned int> m_gtoIndices; //! Indices into the GTO vector 327 std::vector<unsigned int> m_cIndices; //! Indices into m_gtoCN 328 std::vector<double> m_gtoA; //! The GTO exponent 329 std::vector<double> m_gtoC; //! The GTO contraction coefficient 330 std::vector<double> m_gtoCN; //! The GTO contraction coefficient (normalized) 331 332 /** 333 * @brief This block can be once (doubly) or in two parts (alpha and beta) for 334 * open shell calculations. 335 */ 336 MatrixX m_moMatrix[2]; //! MO coefficient matrix 337 338 /** 339 * @brief If there are a sequence of related MOs, they are stored here, and 340 * set as the active MOs upon demand. Alpha will store Paired or the Alpha, 341 * Beta will store Beta coefficients for the appropriate calculation types. 342 */ 343 std::vector<MatrixX> m_moMatrixSet[2]; 344 345 /** 346 * @brief This block stores energies for the molecular orbitals (same 347 * convention as the molecular orbital coefficients). 348 */ 349 std::vector<double> m_moEnergy[2]; 350 351 /** 352 * @brief The occupancy of the molecular orbitals. 353 */ 354 std::vector<unsigned char> m_moOccupancy[2]; 355 356 /** 357 * @brief This stores the molecular orbital number (when they are sparse). It 358 * is used to lookup the actual index of the molecular orbital data. 359 */ 360 std::vector<unsigned int> m_moNumber[2]; 361 362 MatrixX m_density; //! Density matrix 363 MatrixX m_spinDensity; //! Spin Density matrix 364 365 unsigned int m_numMOs; //! The number of GTOs (not always!) 366 bool m_init; //! Has the calculation been initialised? 367 368 ScfType m_scfType; 369 370 std::string m_functionalName; 371 372 /** 373 * @brief Generate the density matrix if we have the required information. 374 * @return True on success, false on failure. 375 */ 376 bool generateDensity(); 377 378 /** 379 * @brief Generate the spin density matrix if we have the required 380 * information. 381 * @return True on success, false on failure. 382 */ 383 bool generateSpinDensity(); 384 }; 385 386 } // End Core namespace 387 } // End Avogadro namespace 388 389 #endif 390