1 /*------------------------------------------------------------------- 2 Copyright 2011 Ravishankar Sundararaman 3 Copyright 1996-2003 Sohrab Ismail-Beigi 4 5 This file is part of JDFTx. 6 7 JDFTx is free software: you can redistribute it and/or modify 8 it under the terms of the GNU General Public License as published by 9 the Free Software Foundation, either version 3 of the License, or 10 (at your option) any later version. 11 12 JDFTx is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 GNU General Public License for more details. 16 17 You should have received a copy of the GNU General Public License 18 along with JDFTx. If not, see <http://www.gnu.org/licenses/>. 19 -------------------------------------------------------------------*/ 20 21 #ifndef JDFTX_ELECTRONIC_IONINFO_H 22 #define JDFTX_ELECTRONIC_IONINFO_H 23 24 #include <electronic/SpeciesInfo.h> 25 #include <electronic/IonicMinimizer.h> 26 #include <core/matrix.h> 27 #include <core/ScalarField.h> 28 #include <core/Thread.h> 29 30 //! @addtogroup IonicSystem 31 //! @{ 32 //! @file IonInfo.h Class IonInfo and related definitions 33 34 //! Coordinate system for ion positions 35 enum CoordsType {CoordsLattice, CoordsCartesian}; 36 37 //! Coordinate system for force output: 38 enum ForcesOutputCoords { ForcesCoordsPositions, ForcesCoordsLattice, ForcesCoordsCartesian, ForcesCoordsContravariant }; 39 static EnumStringMap<ForcesOutputCoords> forcesOutputCoordsMap( 40 ForcesCoordsPositions, "Positions", 41 ForcesCoordsLattice, "Lattice", 42 ForcesCoordsCartesian, "Cartesian", 43 ForcesCoordsContravariant, "Contravariant"); 44 45 //! Check method used for determining whether pseudopotential cores overlap 46 enum coreOverlapCheck { additive, vector, none }; 47 48 //! Container class for ionic system: collection of species, each with several atoms 49 class IonInfo 50 { 51 public: 52 std::vector< std::shared_ptr<SpeciesInfo> > species; //!< list of ionic species 53 std::vector<string> pspFilenamePatterns; //!< list of wildcards for pseudopotential sets 54 55 CoordsType coordsType; //!< coordinate system for ionic positions etc. 56 ForcesOutputCoords forcesOutputCoords; //!< coordinate system to print forces in 57 coreOverlapCheck coreOverlapCondition; //! Check method used for determining whether pseudopotential cores overlap 58 bool vdWenable; //!< whether vdW pair-potential corrections are enabled 59 double vdWscale; //!< If non-zero, override the default scale parameter 60 double ljOverride; //!< If non-zero, replace electronic DFT with LJ pair potential with rCut=ljOverride (for testing geometry optimization and dynamics algorithms only) 61 62 IonicGradient forces; //!< forces at current atomic positions in latice coordinates 63 matrix3<> stress; //!< stresses at current lattice geometry in Eh/a0^3 (only calculated if optimizing lattice or dumping stress) 64 bool computeStress; //!< flag to control whether ionicEnergyAndGrad() computes the stress tensor in addition to forces 65 diagMatrix thermostat; //!< optional thermostat internal degrees of freedom used for IonicDynamics 66 diagMatrix barostat; //!< optional barostat internal degrees of freedom used for IonicDynamics 67 68 ScalarFieldTilde Vlocps; //!< Net local pseudopotential 69 ScalarFieldTilde rhoIon; //!< Total ionic charge density (with width ionWidth, used for interactions with fluid) 70 ScalarFieldTilde nChargeball; //!< Extra electron density around ionic cores to keep fluid out (DEPRECATED) 71 ScalarField nCore; //!< Core electron density for partial (nonlinear) core correction 72 ScalarField tauCore; //!< Model for the KE density of the core (TF+vW applied to nCore) (used by meta-GGAs) 73 74 IonInfo(); 75 76 void setup(const Everything&); 77 void printPositions(FILE*) const; 78 bool checkPositions() const; //!< check for overlapping atoms, return true if okay 79 double getZtot() const; //!< get total Z of all species and atoms 80 81 //! Update Vlocps, rhoIon, nChargeball, nCore and the energies dependent only on ionic positions 82 void update(class Energies&); 83 84 //! Return the total (free) energy and update the ionic gradient in IonInfo::forces 85 //! If IonInfo::computeStress = true, also update the lattice gradient in IonInfo::stress 86 double ionicEnergyAndGrad(); 87 88 //! Return the non-local pseudopotential energy due to a single state. 89 //! Optionally accumulate the corresponding electronic gradient in HCq and ionic gradient in forces 90 double EnlAndGrad(const QuantumNumber& qnum, const diagMatrix& Fq, const std::vector<matrix>& VdagCq, std::vector<matrix>& HVdagCq) const; 91 92 //! Accumulate pseudopotential dependent contribution to the overlap in OCq 93 void augmentOverlap(const ColumnBundle& Cq, ColumnBundle& OCq, std::vector<matrix>* VdagCq=0) const; 94 95 //Multi-stage density augmentation and gradient propagation (see corresponding functions in SpeciesInfo) 96 void augmentDensityInit() const; //!< initialize density augmentation 97 void augmentDensitySpherical(const QuantumNumber& qnum, const diagMatrix& Fq, const std::vector<matrix>& VdagCq) const; //!< calculate density augmentation in spherical functions 98 void augmentDensityGrid(ScalarFieldArray& n) const; //!< propagate from spherical functions to grid 99 void augmentDensityGridGrad(const ScalarFieldArray& E_n, IonicGradient* forces=0, matrix3<>* Eaug_RRT=0) const; //!< propagate grid gradients to spherical functions 100 void augmentDensitySphericalGrad(const QuantumNumber& qnum, const std::vector<matrix>& VdagCq, std::vector<matrix>& HVdagCq) const; //!< propagate spherical function gradients to wavefunctions 101 102 void project(const ColumnBundle& Cq, std::vector<matrix>& VdagCq, matrix* rotExisting=0) const; //Update pseudopotential projections (optionally retain non-zero ones with specified rotation) 103 void projectGrad(const std::vector<matrix>& HVdagCq, const ColumnBundle& Cq, ColumnBundle& HCq) const; //Propagate projected gradient (HVdagCq) to full gradient (HCq) 104 105 //! Compute U corrections (DFT+U in the simplified rotationally-invariant scheme [Dudarev et al, Phys. Rev. B 57, 1505]) 106 //rhoAtom is a flat array of atomic density matrices per U type, with index order (outer to inner): species, Uparam(n,l), spin, atom 107 size_t rhoAtom_nMatrices() const; //!< number of matrices in rhoAtom array 108 void rhoAtom_initZero(std::vector<matrix>& rhoAtom) const; //!< initialize matrices of appropriate size to zero 109 void rhoAtom_calc(const std::vector<diagMatrix>& F, const std::vector<ColumnBundle>& C, std::vector<matrix>& rhoAtom) const; //!< compute atomic density matrices 110 double rhoAtom_computeU(const std::vector<matrix>& rhoAtom, std::vector<matrix>& U_rhoAtom) const; //!< compute U energy and gradient w.r.t atomic density matrices 111 void rhoAtom_grad(const ColumnBundle& Cq, const std::vector<matrix>& U_rhoAtom, ColumnBundle& HCq) const; //!< propagate U_rhoAtom to wavefunction gradient (per k-point to enable band structure) 112 void rhoAtom_forces(const std::vector<diagMatrix>& F, const std::vector<ColumnBundle>& C, const std::vector<matrix>& U_rhoAtom, 113 IonicGradient& forces, matrix3<>* EU_RRT) const; //!< propagate U_rhoAtom to forces (and optionally stress) 114 115 matrix rHcommutator(const ColumnBundle &Y, int iDir, const matrix& YdagHY) const; //!< Expectation value of [r_iDir,H] = D_iDir + nonlocal corrections 116 117 int nAtomicOrbitals() const; //!< Get total number of atomic orbitals 118 ColumnBundle getAtomicOrbitals(int q, bool applyO, int extraCols=0) const; //!< Get all atomic orbitals of a given state number q, optionally with operator O pre-applied (with room for extra columns if specified) 119 120 //! Method for determining ion charge width 121 enum IonWidthMethod 122 { IonWidthEcut, //!< determine ion width from Ecut 123 IonWidthFFTbox, //!< determine ion width from grid spacing 124 IonWidthManual //!< manually specify the ion width 125 } 126 ionWidthMethod; //!< method for determining ion charge width 127 double ionWidth; //!< width for gaussian representation of nuclei 128 bool shouldPrintForceComponents; 129 130 private: 131 const Everything* e; 132 ScalarFieldTilde rhoIonBare; //rhoIon without ionWidth required for stress calculation 133 134 //! Compute all pair-potential terms in the energy, forces or lattice derivative (E_RRT) (electrostatic, and optionally vdW) 135 void pairPotentialsAndGrad(class Energies* ener=0, IonicGradient* forces=0, matrix3<>* E_RRT=0) const; 136 137 //! Compute pulay contributions to energy and optionally stress 138 double calcEpulay(matrix3<>* E_RRT=0) const; 139 }; 140 141 //! @} 142 #endif // JDFTX_ELECTRONIC_IONINFO_H 143