35 #pragma omp parallel for reduction(+ : ewald_g)
36 for (
int igloc = gvec.
skip_g0(); igloc < gvec.
count(); igloc++) {
39 std::complex<double> rho(0, 0);
41 for (
int ia = 0; ia < unit_cell.
num_atoms(); ia++) {
43 static_cast<double>(unit_cell.
atom(ia).zn());
46 ewald_g += std::pow(std::abs(rho), 2) * std::exp(-g2 / 4 / alpha) / g2;
54 ewald_g -= std::pow(unit_cell.
num_electrons(), 2) / alpha / 4;
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);
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;
72 return (ewald_g + ewald_r);
90 return potential.energy_vha();
98 ebxc += sirius::inner(density.mag(j), potential.effective_magnetic_field(j));
106 auto& unit_cell = ctx.unit_cell();
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);
121 return sirius::inner(potential.local_potential(), density.
rho().rg());
143 return sirius::inner(density.
rho(), potential.effective_potential());
153ks_energy(Simulation_context
const& ctx,
const std::map<std::string, double>& energies)
157 switch (ctx.electronic_structure_method()) {
159 tot_en = energies.at(
"ekin") + energies.at(
"exc") + 0.5 * energies.at(
"vha") + energies.at(
"enuc");
165 energies.at(
"valence_eval_sum") - energies.at(
"vxc") - energies.at(
"bxc") - energies.at(
"PAW_one_elec");
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");
179ks_energy(Simulation_context
const& ctx, K_point_set
const& kset, Density
const& density, Potential
const& potential,
182 return ks_energy(ctx, total_energy_components(ctx, kset, density, potential,
ewald_energy));
190 double eks = ks_energy(ctx, kset, density, potential,
ewald_energy);
193 switch (ctx.electronic_structure_method()) {
204 RTE_THROW(
"invalid electronic_structure_method");
211std::map<std::string, double>
212total_energy_components(Simulation_context
const& ctx,
const K_point_set& kset, Density
const& density,
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);
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);
231 table[
"exc"] =
energy_exc(density, potential);
233 table[
"PAW_total_energy"] = potential.PAW_total_energy(density);
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);
243 table[
"entropy"] = kset.entropy_sum();
249hubbard_energy(Density
const& density)
251 if (density.ctx().hubbard_correction()) {
252 return energy(density.occupation_matrix());
259one_electron_energy(Density
const& density, Potential
const& potential)
262 potential.PAW_one_elec_energy(density) + one_electron_energy_hubbard(density, potential);
266one_electron_energy_hubbard(Density
const& density, Potential
const& potential)
268 auto& ctx = density.ctx();
269 if (ctx.hubbard_correction()) {
270 return one_electron_energy_hubbard(density.occupation_matrix(), potential.hubbard_potential());
276energy_potential(Density
const& density, Potential
const& potential)
279 potential.PAW_one_elec_energy(density) + ::sirius::hubbard_energy(density);
int num_atoms() const
Return number of atoms belonging to the current symmetry class.
Generate charge density and magnetization from occupied spinor wave-functions.
auto const & rho() const
Return const reference to charge density (scalar functions).
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.
auto energy_vxc(Density const &density__) const
Integral of .
auto energy_exc(Density const &density__) const
Integral of .
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.
double omega() const
Unit cell volume.
Atom const & atom(int id__) const
Return const atom instance by id.
int num_atom_symmetry_classes() const
Number of atom symmetry classes.
int num_atoms() const
Number of atoms in the unit cell.
double num_electrons() const
Total number of electrons (core + valence).
Atom_symmetry_class const & atom_symmetry_class(int id__) const
Return const symmetry class instance by class id.
A set of G-vectors for FFTs and G+k basis functions.
double gvec_len(int ig__) const
Return length of the G-vector.
int skip_g0() const
Local starting index of G-vectors if G=0 is not counted.
int count() const
Number of G-vectors for a fine-grained distribution for the current MPI rank.
r3::vector< int > gvec(int ig__) const
Return G vector in fractional coordinates.
void allreduce(T *buffer__, int count__) const
Perform the in-place (the output buffer is used as the input buffer) all-to-all reduction.
Namespace of the SIRIUS library.
double core_eval_sum(Unit_cell const &unit_cell)
Return eigen-value sum of core states.
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.
@ 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.
double eval_sum(Unit_cell const &unit_cell, K_point_set const &kset)
TODO doc.
double energy_vxc(Density const &density, Potential const &potential)
Returns exchange correlation potential.
double energy_vloc(Density const &density, Potential const &potential)
TODO doc.
double energy_vha(Potential const &potential)
Returns Hatree potential.
double energy_kin(Simulation_context const &ctx, K_point_set const &kset, Density const &density, Potential const &potential)
Return kinetic energy.
double energy_enuc(Simulation_context const &ctx, Potential const &potential)
Return nucleus energy in the electrostatic field.
double energy_bxc(const Density &density, const Potential &potential)
TODO doc.
double energy_exc(Density const &density, Potential const &potential)
Returns exchange correlation energy.
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.