Loading [MathJax]/extensions/TeX/AMSsymbols.js
SIRIUS 7.5.0
Electronic structure library and applications
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
energy.cpp
Go to the documentation of this file.
1// Copyright (c) 2013-2020 Anton Kozhevnikov, Thomas Schulthess
2// All rights reserved.
3//
4// Redistribution and use in source and binary forms, with or without modification, are permitted provided that
5// the following conditions are met:
6//
7// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the
8// following disclaimer.
9// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
10// and the following disclaimer in the documentation and/or other materials provided with the distribution.
11//
12// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
13// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
14// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
15// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
16// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
17// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
18// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
19
20/** \file energy.cpp
21 *
22 * \brief Total energy terms.
23 */
24
25#include "energy.hpp"
26
27namespace sirius {
28
29double
30ewald_energy(const Simulation_context& ctx, const fft::Gvec& gvec, const Unit_cell& unit_cell)
31{
32 double alpha{ctx.ewald_lambda()};
33 double ewald_g{0};
34
35 #pragma omp parallel for reduction(+ : ewald_g)
36 for (int igloc = gvec.skip_g0(); igloc < gvec.count(); igloc++) {
37 double g2 = std::pow(gvec.gvec_len<index_domain_t::local>(igloc), 2);
38
39 std::complex<double> rho(0, 0);
40
41 for (int ia = 0; ia < unit_cell.num_atoms(); ia++) {
42 rho += ctx.gvec_phase_factor(gvec.gvec<index_domain_t::local>(igloc), ia) *
43 static_cast<double>(unit_cell.atom(ia).zn());
44 }
45
46 ewald_g += std::pow(std::abs(rho), 2) * std::exp(-g2 / 4 / alpha) / g2;
47 }
48
49 ctx.comm().allreduce(&ewald_g, 1);
50 if (gvec.reduced()) {
51 ewald_g *= 2;
52 }
53 /* remaining G=0 contribution */
54 ewald_g -= std::pow(unit_cell.num_electrons(), 2) / alpha / 4;
55 ewald_g *= (twopi / unit_cell.omega());
56
57 /* remove self-interaction */
58 for (int ia = 0; ia < unit_cell.num_atoms(); ia++) {
59 ewald_g -= std::sqrt(alpha / pi) * std::pow(unit_cell.atom(ia).zn(), 2);
60 }
61
62 double ewald_r{0};
63 #pragma omp parallel for reduction(+ : ewald_r)
64 for (int ia = 0; ia < unit_cell.num_atoms(); ia++) {
65 for (int i = 1; i < unit_cell.num_nearest_neighbours(ia); i++) {
66 int ja = unit_cell.nearest_neighbour(i, ia).atom_id;
67 double d = unit_cell.nearest_neighbour(i, ia).distance;
68 ewald_r += 0.5 * unit_cell.atom(ia).zn() * unit_cell.atom(ja).zn() * std::erfc(std::sqrt(alpha) * d) / d;
69 }
70 }
71
72 return (ewald_g + ewald_r);
73}
74
75double
76energy_vxc(Density const& density, Potential const& potential)
77{
78 return potential.energy_vxc(density);
79}
80
81double
82energy_exc(Density const& density, Potential const& potential)
83{
84 return potential.energy_exc(density);
85}
86
87double
88energy_vha(Potential const& potential)
89{
90 return potential.energy_vha();
91}
92
93double
94energy_bxc(const Density& density, const Potential& potential)
95{
96 double ebxc{0};
97 for (int j = 0; j < density.ctx().num_mag_dims(); j++) {
98 ebxc += sirius::inner(density.mag(j), potential.effective_magnetic_field(j));
99 }
100 return ebxc;
101}
102
103double
104energy_enuc(Simulation_context const& ctx, Potential const& potential)
105{
106 auto& unit_cell = ctx.unit_cell();
107 double enuc{0};
108 if (ctx.full_potential()) {
109 for (auto it : unit_cell.spl_num_atoms()) {
110 int zn = unit_cell.atom(it.i).zn();
111 enuc -= 0.5 * zn * potential.vh_el(it.i);
112 }
113 ctx.comm().allreduce(&enuc, 1);
114 }
115 return enuc;
116}
117
118double
119energy_vloc(Density const& density, Potential const& potential)
120{
121 return sirius::inner(potential.local_potential(), density.rho().rg());
122}
123
124double
125core_eval_sum(Unit_cell const& unit_cell)
126{
127 double sum{0};
128 for (int ic = 0; ic < unit_cell.num_atom_symmetry_classes(); ic++) {
129 sum += unit_cell.atom_symmetry_class(ic).core_eval_sum() * unit_cell.atom_symmetry_class(ic).num_atoms();
130 }
131 return sum;
132}
133
134double
135eval_sum(Unit_cell const& unit_cell, K_point_set const& kset)
136{
137 return core_eval_sum(unit_cell) + kset.valence_eval_sum();
138}
139
140double
141energy_veff(Density const& density, Potential const& potential)
142{
143 return sirius::inner(density.rho(), potential.effective_potential());
144}
145
146double
147energy_kin(Simulation_context const& ctx, K_point_set const& kset, Density const& density, Potential const& potential)
148{
149 return eval_sum(ctx.unit_cell(), kset) - energy_veff(density, potential) - energy_bxc(density, potential);
150}
151
152double
153ks_energy(Simulation_context const& ctx, const std::map<std::string, double>& energies)
154{
155 double tot_en{0};
156
157 switch (ctx.electronic_structure_method()) {
159 tot_en = energies.at("ekin") + energies.at("exc") + 0.5 * energies.at("vha") + energies.at("enuc");
160 break;
161 }
162
164 tot_en =
165 energies.at("valence_eval_sum") - energies.at("vxc") - energies.at("bxc") - energies.at("PAW_one_elec");
166 tot_en +=
167 -0.5 * energies.at("vha") + energies.at("exc") + energies.at("PAW_total_energy") + energies.at("ewald");
168 if (ctx.hubbard_correction()) {
169 tot_en += energies.at("hubbard_energy") - energies.at("hubbard_one_el_contribution");
170 }
171 break;
172 }
173 }
174
175 return tot_en;
176}
177
178double
179ks_energy(Simulation_context const& ctx, K_point_set const& kset, Density const& density, Potential const& potential,
180 double ewald_energy)
181{
182 return ks_energy(ctx, total_energy_components(ctx, kset, density, potential, ewald_energy));
183}
184
185double
186total_energy(Simulation_context const& ctx, K_point_set const& kset, Density const& density, Potential const& potential,
187 double ewald_energy)
188{
189
190 double eks = ks_energy(ctx, kset, density, potential, ewald_energy);
191 double tot_en{0};
192
193 switch (ctx.electronic_structure_method()) {
195 tot_en = eks;
196 break;
197 }
198
200 tot_en = eks + kset.entropy_sum();
201 break;
202 }
203 default: {
204 RTE_THROW("invalid electronic_structure_method");
205 }
206 }
207
208 return tot_en;
209}
210
211std::map<std::string, double>
212total_energy_components(Simulation_context const& ctx, const K_point_set& kset, Density const& density,
213 Potential const& potential, double ewald_energy)
214{
215 std::map<std::string, double> table;
216 switch (ctx.electronic_structure_method()) {
218 table["ekin"] = energy_kin(ctx, kset, density, potential);
219 table["exc"] = energy_exc(density, potential);
220 table["vha"] = energy_vha(potential);
221 table["enuc"] = energy_enuc(ctx, potential);
222 break;
223 }
224
226 table["valence_eval_sum"] = kset.valence_eval_sum();
227 table["vxc"] = energy_vxc(density, potential);
228 table["bxc"] = energy_bxc(density, potential);
229 table["PAW_one_elec"] = potential.PAW_one_elec_energy(density);
230 table["vha"] = energy_vha(potential);
231 table["exc"] = energy_exc(density, potential);
232 table["ewald"] = ewald_energy;
233 table["PAW_total_energy"] = potential.PAW_total_energy(density);
234 break;
235 }
236 }
237
238 if (ctx.hubbard_correction()) {
239 table["hubbard_one_el_contribution"] = one_electron_energy_hubbard(density, potential);
240 table["hubbard_energy"] = ::sirius::hubbard_energy(density);
241 }
242
243 table["entropy"] = kset.entropy_sum();
244
245 return table;
246}
247
248double
249hubbard_energy(Density const& density)
250{
251 if (density.ctx().hubbard_correction()) {
252 return energy(density.occupation_matrix());
253 } else {
254 return 0.0;
255 }
256}
257
258double
259one_electron_energy(Density const& density, Potential const& potential)
260{
261 return energy_vha(potential) + energy_vxc(density, potential) + energy_bxc(density, potential) +
262 potential.PAW_one_elec_energy(density) + one_electron_energy_hubbard(density, potential);
263}
264
265double
266one_electron_energy_hubbard(Density const& density, Potential const& potential)
267{
268 auto& ctx = density.ctx();
269 if (ctx.hubbard_correction()) {
270 return one_electron_energy_hubbard(density.occupation_matrix(), potential.hubbard_potential());
271 }
272 return 0.0;
273}
274
275double
276energy_potential(Density const& density, Potential const& potential)
277{
278 const double e = energy_veff(density, potential) + energy_bxc(density, potential) +
279 potential.PAW_one_elec_energy(density) + ::sirius::hubbard_energy(density);
280 return e;
281}
282
283} // namespace sirius
int num_atoms() const
Return number of atoms belonging to the current symmetry class.
Generate charge density and magnetization from occupied spinor wave-functions.
Definition: density.hpp:214
auto const & rho() const
Return const reference to charge density (scalar functions).
Definition: density.hpp:416
Set of k-points.
Definition: k_point_set.hpp:41
double entropy_sum() const
Return entropy contribution from smearing store in Kpoint<T>.
double valence_eval_sum() const
Return sum of valence eigen-values store in Kpoint<T>.
Generate effective potential from charge density and magnetization.
Definition: potential.hpp:46
auto energy_vxc(Density const &density__) const
Integral of .
Definition: potential.hpp:819
auto energy_exc(Density const &density__) const
Integral of .
Definition: potential.hpp:831
Simulation context is a set of parameters and objects describing a single simulation.
auto gvec_phase_factor(r3::vector< int > G__, int ia__) const
Phase factors .
double ewald_lambda() const
Find the lambda parameter used in the Ewald summation.
mpi::Communicator const & comm() const
Total communicator of the simulation.
Representation of a unit cell.
Definition: unit_cell.hpp:43
double omega() const
Unit cell volume.
Definition: unit_cell.hpp:275
Atom const & atom(int id__) const
Return const atom instance by id.
Definition: unit_cell.hpp:344
int num_atom_symmetry_classes() const
Number of atom symmetry classes.
Definition: unit_cell.hpp:320
int num_atoms() const
Number of atoms in the unit cell.
Definition: unit_cell.hpp:338
double num_electrons() const
Total number of electrons (core + valence).
Definition: unit_cell.hpp:369
Atom_symmetry_class const & atom_symmetry_class(int id__) const
Return const symmetry class instance by class id.
Definition: unit_cell.hpp:326
A set of G-vectors for FFTs and G+k basis functions.
Definition: gvec.hpp:130
double gvec_len(int ig__) const
Return length of the G-vector.
Definition: gvec.hpp:637
int skip_g0() const
Local starting index of G-vectors if G=0 is not counted.
Definition: gvec.hpp:526
int count() const
Number of G-vectors for a fine-grained distribution for the current MPI rank.
Definition: gvec.hpp:506
r3::vector< int > gvec(int ig__) const
Return G vector in fractional coordinates.
Definition: gvec.hpp:540
void allreduce(T *buffer__, int count__) const
Perform the in-place (the output buffer is used as the input buffer) all-to-all reduction.
Total energy terms.
Namespace of the SIRIUS library.
Definition: sirius.f90:5
double core_eval_sum(Unit_cell const &unit_cell)
Return eigen-value sum of core states.
Definition: energy.cpp:125
double total_energy(Simulation_context const &ctx, K_point_set const &kset, Density const &density, Potential const &potential, double ewald_energy)
Total energy of the electronic subsystem.
Definition: energy.cpp:186
@ full_potential_lapwlo
Full potential linearized augmented plane waves with local orbitals.
@ pseudopotential
Pseudopotential (ultrasoft, norm-conserving, PAW).
double energy_veff(Density const &density, Potential const &potential)
TODO doc.
Definition: energy.cpp:141
double eval_sum(Unit_cell const &unit_cell, K_point_set const &kset)
TODO doc.
Definition: energy.cpp:135
double energy_vxc(Density const &density, Potential const &potential)
Returns exchange correlation potential.
Definition: energy.cpp:76
const double twopi
Definition: constants.hpp:45
double energy_vloc(Density const &density, Potential const &potential)
TODO doc.
Definition: energy.cpp:119
const double pi
Definition: constants.hpp:42
double energy_vha(Potential const &potential)
Returns Hatree potential.
Definition: energy.cpp:88
double energy_kin(Simulation_context const &ctx, K_point_set const &kset, Density const &density, Potential const &potential)
Return kinetic energy.
Definition: energy.cpp:147
double energy_enuc(Simulation_context const &ctx, Potential const &potential)
Return nucleus energy in the electrostatic field.
Definition: energy.cpp:104
double energy_bxc(const Density &density, const Potential &potential)
TODO doc.
Definition: energy.cpp:94
double energy_exc(Density const &density, Potential const &potential)
Returns exchange correlation energy.
Definition: energy.cpp:82
double ewald_energy(const Simulation_context &ctx, const fft::Gvec &gvec, const Unit_cell &unit_cell)
Compute the ion-ion electrostatic energy using Ewald method.
Definition: energy.cpp:30