36 PROFILE(
"sirius::DFT_ground_state::initial_state");
40 if (!
ctx_.full_potential()) {
41 if (
ctx_.cfg().parameters().precision_wf() ==
"fp32") {
42#if defined(SIRIUS_USE_FP32)
46 RTE_THROW(
"not compiled with FP32 support");
57DFT_ground_state::create_H0()
59 PROFILE(
"sirius::DFT_ground_state::create_H0");
67 PROFILE(
"sirius::DFT_ground_state::update");
74 if (!
ctx_.full_potential()) {
80DFT_ground_state::energy_kin_sum_pw()
const
84 for (
auto it :
kset_.spl_num_kpoints()) {
85 auto kp =
kset_.get<
double>(it.i);
87 #pragma omp parallel for schedule(static) reduction(+:ekin)
88 for (
int igloc = 0; igloc < kp->num_gkvec_loc(); igloc++) {
93 for (
int i = 0; i < kp->num_occupied_bands(ispin); i++) {
94 double f = kp->band_occupancy(i, ispin);
95 auto z = kp->spinor_wave_functions().pw_coeffs(igloc, wf::spin_index(ispin), wf::band_index(i));
96 d += f * (std::pow(z.real(), 2) + std::pow(z.imag(), 2));
99 if (kp->gkvec().reduced()) {
102 ekin += 0.5 * d * kp->weight() * Gk.length2();
110DFT_ground_state::total_energy()
const
116DFT_ground_state::serialize()
125 if (
ctx_.full_potential()) {
134 bool transform_to_rg{
true};
137 bool precompute_lapw{
true};
142 ::sirius::diagonalize<double, double>(H0,
kset_,
ctx_.cfg().settings().itsol_tol_min());
157 for (
int ig = 0; ig <
ctx_.
gvec().count(); ig++) {
158 rms += std::pow(std::abs(a.component(j).rg().f_pw_local(ig) - b.component(j).rg().f_pw_local(ig)), 2);
162 return std::sqrt(rms /
ctx_.
gvec().num_gvec());
165 double rms = calc_rms(
density_, rho);
170 dict[
"detot"] = gs0[
"energy"][
"total"].get<
double>() - gs1[
"energy"][
"total"].get<double>();
172 rms_veff = std::sqrt(rms_veff /
ctx_.
gvec().num_gvec());
174 if (
ctx_.verbosity() >= 1) {
175 RTE_OUT(
ctx_.
out()) <<
"RMS_rho: " << dict[
"rms"].get<
double>() << std::endl
176 <<
"RMS_veff: " << rms_veff << std::endl
177 <<
"Eold: " << gs0[
"energy"][
"total"].get<double>()
178 <<
" Enew: " << gs1[
"energy"][
"total"].get<
double>() << std::endl;
180 std::vector<std::string> labels({
"total",
"vha",
"vxc",
"exc",
"bxc",
"veff",
"eval_sum",
"kin",
"ewald",
181 "vloc",
"scf_correction",
"entropy_sum"});
183 for (
auto e: labels) {
184 RTE_OUT(
ctx_.
out()) <<
"energy component: " << e <<
", diff: " <<
185 std::abs(gs0[
"energy"][e].get<double>() - gs1[
"energy"][e].get<double>()) << std::endl;
196 PROFILE(
"sirius::DFT_ground_state::scf_loop");
198 auto tstart = std::chrono::high_resolution_clock::now();
200 double eold{0}, rms{0};
205 std::vector<double> rms_hist;
206 std::vector<double> etot_hist;
211 s <<
"density_tol : " << density_tol__ << std::endl
212 <<
"energy_tol : " << energy_tol__ << std::endl
213 <<
"iter_solver_tol (initial) : " << iter_solver_tol__ << std::endl
214 <<
"iter_solver_tol (target) : " <<
ctx_.cfg().settings().itsol_tol_min() << std::endl
215 <<
"num_dft_iter : " << num_dft_iter__;
218 for (
int iter = 0; iter < num_dft_iter__; iter++) {
219 PROFILE(
"sirius::DFT_ground_state::scf_loop|iteration");
222 s <<
"+------------------------------+" << std::endl
223 <<
"| SCF iteration " << std::setw(3) << iter <<
" out of " << std::setw(3) << num_dft_iter__ <<
'|' << std::endl
224 <<
"+------------------------------+" << std::endl;
229 if (
ctx_.cfg().parameters().precision_wf() ==
"fp32") {
230#if defined(SIRIUS_USE_FP32)
233 if (
ctx_.cfg().parameters().precision_hs() ==
"fp32") {
234 result = sirius::diagonalize<float, float>(H0,
kset_, iter_solver_tol__);
236 result = sirius::diagonalize<float, double>(H0,
kset_, iter_solver_tol__);
243 RTE_THROW(
"not compiled with FP32 support");
248 result = sirius::diagonalize<double, double>(H0,
kset_, iter_solver_tol__);
261 double eha_res = density_residual_hartree_energy(
density_, rho1);
265 if (
ctx_.cfg().mixer().use_hartree()) {
269 tol = std::min(
ctx_.cfg().settings().itsol_tol_scale()[0] * tol,
270 ctx_.cfg().settings().itsol_tol_scale()[1] * iter_solver_tol__);
272 iter_solver_tol__ = std::max(
ctx_.cfg().settings().itsol_tol_min(), tol);
274 bool iter_solver_converged{
true};
275 if (
ctx_.cfg().iterative_solver().type() !=
"exact") {
276 iter_solver_converged = (tol <=
ctx_.cfg().settings().itsol_tol_min());
279#if defined(SIRIUS_USE_FP32)
280 if (
ctx_.cfg().parameters().precision_gs() !=
"auto") {
282 if (
ctx_.cfg().parameters().precision_gs() ==
"fp64" &&
ctx_.cfg().parameters().precision_wf() ==
"fp32") {
284 if ((
ctx_.cfg().settings().fp32_to_fp64_rms() == 0 && iter_solver_tol__ <=
ctx_.cfg().settings().itsol_tol_min()) ||
285 (rms <
ctx_.cfg().settings().fp32_to_fp64_rms())) {
286 std::cout <<
"switching to FP64" << std::endl;
288 ctx_.cfg().settings().itsol_tol_min(std::numeric_limits<double>::epsilon() * 10);
289 ctx_.cfg().parameters().precision_wf(
"fp64");
290 ctx_.cfg().parameters().precision_hs(
"fp64");
293 for (
auto it :
kset_.spl_num_kpoints()) {
295 wf::copy(sddk::memory_t::host,
kset_.get<
float>(it.i)->spinor_wave_functions(),
304 kset_.get<
double>(ik)->band_energy(j, ispn,
kset_.get<
float>(ik)->band_energy(j, ispn));
305 kset_.get<
double>(ik)->band_occupancy(j, ispn,
kset_.get<
float>(ik)->band_occupancy(j, ispn));
313 if (
ctx_.cfg().control().verification() >= 1) {
321 if (!
ctx_.full_potential() &&
ctx_.cfg().control().verification() >= 2) {
322 if (
ctx_.verbosity() >= 1) {
323 RTE_OUT(
ctx_.
out()) <<
"checking functional derivative of Exc\n";
325 sirius::check_xc_potential(
density_);
328 if (
ctx_.cfg().parameters().use_scf_correction()) {
329 double e2 = energy_potential(rho1,
potential_);
334 double etot = total_energy();
336 etot_hist.push_back(etot);
338 rms_hist.push_back(rms);
341 std::stringstream out;
345 out <<
"iteration : " << iter <<
", RMS : " << std::setprecision(12) << std::scientific << rms
346 <<
", energy difference : " << std::setprecision(12) << std::scientific << etot - eold;
347 if (!
ctx_.full_potential()) {
349 <<
"Hartree energy of density residual : " << eha_res << std::endl
350 <<
"bands are converged : " << boolstr(result.converged);
352 if (
ctx_.cfg().iterative_solver().type() !=
"exact") {
354 <<
"iterative solver converged : " << boolstr(iter_solver_converged);
359 bool converged{
true};
360 converged = (std::abs(eold - etot) < energy_tol__) && result.converged && iter_solver_converged;
361 if (
ctx_.cfg().mixer().use_hartree()) {
362 converged = converged && (eha_res < density_tol__);
364 converged = converged && (rms < density_tol__);
367 std::stringstream out;
369 out <<
"converged after " << iter + 1 <<
" SCF iterations!" << std::endl;
378 std::stringstream out;
384 ctx_.create_storage_file(storage_file_name);
385 if (
ctx_.full_potential()) {
388 density_.mag(j).rg().fft_transform(-1);
396 auto tstop = std::chrono::high_resolution_clock::now();
398 auto dict = serialize();
403 double rho_min{1e100};
404 for (
int ir = 0; ir <
density_.
rho().rg().spfft().local_slice_size(); ir++) {
405 rho_min = std::min(rho_min,
density_.
rho().rg().value(ir));
407 dict[
"rho_min"] = rho_min;
411 dict[
"scf_time"] = std::chrono::duration_cast<std::chrono::duration<double>>(tstop - tstart).count();
412 dict[
"etot_history"] = etot_hist;
414 dict[
"converged"] =
true;
415 dict[
"num_scf_iterations"] = num_iter;
416 dict[
"rms_history"] = rms_hist;
418 dict[
"converged"] =
false;
421 if (env::check_scf_density()) {
445 double etot = total_energy();
447 double ef =
kset_.energy_fermi();
450 double one_elec_en = evalsum1 - (evxc + evha + ebxc);
454 one_elec_en -= hub_one_elec;
460 out__ <<
"Energy" << std::endl
461 <<
hbar(80,
'-') << std::endl;
463 auto write_energy = [&](std::string label__,
double value__) {
464 out__ << std::left << std::setw(30) << label__ <<
" : " << std::right << std::setw(16) << std::setprecision(8)
465 << std::fixed << value__ << std::endl;
468 auto write_energy2 = [&](std::string label__,
double value__)
470 out__ << std::left << std::setw(30) << label__ <<
" : " << std::right
471 << std::setw(16) << std::setprecision(8) << std::fixed << value__ <<
" (Ha), "
472 << std::setw(16) << std::setprecision(8) << std::fixed << value__ * 2 <<
" (Ry)" << std::endl;
475 write_energy(
"valence_eval_sum", evalsum1);
476 if (
ctx_.full_potential()) {
477 write_energy(
"core_eval_sum", evalsum2);
478 write_energy(
"kinetic energy", ekin);
479 write_energy(
"enuc", enuc);
481 write_energy(
"<rho|V^{XC}>", evxc);
482 write_energy(
"<rho|E^{XC}>", eexc);
483 write_energy(
"<mag|B^{XC}>", ebxc);
484 write_energy(
"<rho|V^{H}>", evha);
485 if (!
ctx_.full_potential()) {
486 write_energy2(
"one-electron contribution", one_elec_en);
487 write_energy(
"hartree contribution", 0.5 * evha);
488 write_energy(
"xc contribution", eexc);
492 write_energy(
"smearing (-TS)", s_sum);
494 if (
ctx_.hubbard_correction()) {
495 auto e = ::sirius::energy(
density_.occupation_matrix());
496 write_energy2(
"Hubbard energy", e);
497 write_energy2(
"Hubbard one-el contribution", hub_one_elec);
499 write_energy2(
"Total energy", etot);
501 write_energy(
"band gap (eV)", gap);
502 write_energy(
"Efermi", ef);
namespace for Niels Lohmann
Potential potential_
Instance of the Potential class.
double scf_correction_energy_
Correction to total energy from the SCF density minimisation.
void initial_state()
Generate initial density, potential and a subspace of wave-functions.
void print_info(std::ostream &out__) const
Print the basic information (total energy, charges, moments, etc.).
Density density_
Instance of the Density class.
Simulation_context & ctx_
Context of simulation.
K_point_set & kset_
Set of k-points that are used to generate density.
double ewald_energy_
Store Ewald energy which is computed once and which doesn't change during the run.
json find(double density_tol, double energy_tol, double initial_tolerance, int num_dft_iter, bool write_state)
Run the SCF ground state calculation and find a total energy minimum.
std::shared_ptr< Hamiltonian0< double > > H0_
k-point independent part of the Hamiltonian.
void update()
Update the parameters after the change of lattice vectors or atomic positions.
Unit_cell & unit_cell_
Alias of the unit cell.
json check_scf_density()
A quick check of self-constent density in case of pseudopotential.
Generate charge density and magnetization from occupied spinor wave-functions.
void mixer_init(config_t::mixer_t const &mixer_cfg__)
Initialize density mixer.
void update()
Update internal parameters after a change of lattice vectors or atomic coordinates.
bool check_num_electrons() const
Check total density for the correct number of electrons.
double mix()
Mix new density.
auto const & rho() const
Return const reference to charge density (scalar functions).
void initial_density()
Generate initial charge density and magnetization.
void generate(K_point_set const &ks__, bool symmetrize__, bool add_core__, bool transform_to_rg__)
Generate full charge density (valence + core) and magnetization from the wave functions.
Four-component function consisting of scalar and vector parts.
Represent the k-point independent part of Hamiltonian.
int num_kpoints() const
Return total number of k-points.
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>.
void find_band_occupancies()
Find Fermi energy and band occupation numbers.
void update()
Update k-points after moving atoms or changing the lattice vectors.
Generate effective potential from charge density and magnetization.
void generate(Density const &density__, bool use_sym__, bool transform_to_rg__)
Generate effective potential and magnetic field from charge density and magnetization.
void update()
Recompute some variables that depend on atomic positions or the muffin-tin radius.
auto const & gvec() const
Return const reference to Gvec object.
void update()
Update context after setting new lattice vectors or atomic coordinates.
std::ostream & out() const
Return output stream.
mpi::Communicator const & comm() const
Total communicator of the simulation.
void message(int level__, char const *label__, std::stringstream const &s) const
Print message from the stringstream.
int num_spins() const
Number of spin components.
int num_spinors() const
Number of spinor wave-functions labeled by a sinlge band index.
int num_bands(int num_bands__)
Set the number of bands.
int num_mag_dims() const
Number of dimensions in the magnetization vector.
bool use_symmetry() const
Get a use_symmetry flag.
double num_electrons() const
Total number of electrons (core + valence).
void allreduce(T *buffer__, int count__) const
Perform the in-place (the output buffer is used as the input buffer) all-to-all reduction.
Describe a range of bands.
Contains definition and partial implementation of sirius::DFT_ground_state class.
Entry point for Hamiltonain diagonalization.
Create intial subspace from atomic-like wave-functions.
void copy(sddk::memory_t mem__, Wave_functions< T > const &in__, wf::spin_index s_in__, wf::band_range br_in__, Wave_functions< F > &out__, wf::spin_index s_out__, wf::band_range br_out__)
Copy wave-functions.
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.
@ pseudopotential
Pseudopotential (ultrasoft, norm-conserving, PAW).
double energy_vxc(Density const &density, Potential const &potential)
Returns exchange correlation potential.
double energy_vha(Potential const &potential)
Returns Hatree potential.
const double ha2ev
Hartree in electron-volt units.
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.
void initialize_subspace(Hamiltonian_k< T > const &Hk__, K_point< T > &kp__, int num_ao__)
Initialize the wave-functions subspace at a given k-point.
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.