From 1e7ec3eaf9706db6cc0af7aebba764eaf718ec61 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 13:37:00 -0700 Subject: [PATCH 01/27] Implementation of the kinetic energy part of the stress sensor. --- src/observables/forces_stress.hpp | 36 ++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 366f7895..2c667055 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -36,6 +36,23 @@ struct forces_stress { #ifndef ENABLE_GPU private: #endif + + GPU_FUNCTION static void stress_component(int const index, int & alpha, int & beta) { + alpha = index; + beta = index; + if(index == 3) { + alpha = 0; + beta = 1; + } + if(index == 4) { + alpha = 1; + beta = 2; + } + if(index == 5) { + alpha = 0; + beta = 2; + } + } template void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){ @@ -46,7 +63,7 @@ struct forces_stress { gdensity.fill(vector3{0.0, 0.0, 0.0}); gpu::array, 1> forces_non_local(ions.size(), {0.0, 0.0, 0.0}); - + auto iphi = 0; for(auto & phi : electrons.kpin()){ @@ -54,6 +71,23 @@ struct forces_stress { observables::density::calculate_gradient_add(electrons.occupations()[iphi], phi, gphi, gdensity); ham.projectors_all().force(phi, gphi, ions.cell().metric(), electrons.occupations()[iphi], ham.uniform_vector_potential(), forces_non_local); + + auto stress_kinetic = gpu::run(6, gpu::reduce(gphi.local_set_size()), gpu::reduce(gphi.basis().local_size()), + [gph = begin(gphi.matrix()), occ = begin(electrons.occupations()[iphi])] GPU_LAMBDA (auto index, auto ist, auto ip) { + int alpha, beta; + stress_component(index, alpha, beta); + return occ[ist]*conj(gph[ip][ist][alpha])*gph[ip][ist][beta]; + }); + + if(gphi.full_comm().size() > 1){ + gphi.full_comm().all_reduce_n(raw_pointer_cast(stress_kinetic.data_elements()), 6);; + } + + for(auto index = 0; index < 6; index++) { + int alpha, beta; + stress_component(index, alpha, beta); + stress[alpha][beta] += -2.0/gphi.basis().cell().volume()*real(stress_kinetic[index]); + } iphi++; } From a45000759efb617c9b87d6a748c2ddcf0b9f5111 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 14:36:01 -0700 Subject: [PATCH 02/27] Pass the energy to calculate the stress. --- src/ground_state/calculator.hpp | 2 +- src/observables/forces_stress.hpp | 11 ++++++----- src/real_time/propagate.hpp | 6 +++--- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/ground_state/calculator.hpp b/src/ground_state/calculator.hpp index 6394246a..9a445eb8 100644 --- a/src/ground_state/calculator.hpp +++ b/src/ground_state/calculator.hpp @@ -238,7 +238,7 @@ class calculator { auto normres = res.energy.calculate(ham_, electrons); if(solver_.calc_forces() and electrons.states().spinor_dim() == 1) { - res.forces = observables::forces_stress{ions_, electrons, ham_}.forces; + res.forces = observables::forces_stress{ions_, electrons, ham_, res.energy}.forces; } if(solver_.calc_forces() and electrons.states().spinor_dim() == 2) { diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 2c667055..1c01848c 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -25,12 +25,12 @@ struct forces_stress { forces_stress() = default; - template - forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham): + template + forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham, Energy const & energy): forces(ions.size()), stress({3, 3}, 0.0) { - calculate(ions, electrons, ham); + calculate(ions, electrons, ham, energy); } #ifndef ENABLE_GPU @@ -54,8 +54,8 @@ struct forces_stress { } } - template - void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){ + template + void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham, Energy const & energy){ CALI_CXX_MARK_FUNCTION; @@ -72,6 +72,7 @@ struct forces_stress { ham.projectors_all().force(phi, gphi, ions.cell().metric(), electrons.occupations()[iphi], ham.uniform_vector_potential(), forces_non_local); + //STRESS KINETIC auto stress_kinetic = gpu::run(6, gpu::reduce(gphi.local_set_size()), gpu::reduce(gphi.basis().local_size()), [gph = begin(gphi.matrix()), occ = begin(electrons.occupations()[iphi])] GPU_LAMBDA (auto index, auto ist, auto ip) { int alpha, beta; diff --git a/src/real_time/propagate.hpp b/src/real_time/propagate.hpp index 638a801f..ae55ad34 100644 --- a/src/real_time/propagate.hpp +++ b/src/real_time/propagate.hpp @@ -72,8 +72,8 @@ void propagate(systems::ions & ions, systems::electrons & electrons, ProcessFunc energy.calculate(ham, electrons); energy.ion(ionic::interaction_energy(ions.cell(), ions, electrons.atomic_pot())); - auto forces = decltype(observables::forces_stress{ions, electrons, ham}.forces){}; - if(ion_propagator.needs_force()) forces = observables::forces_stress{ions, electrons, ham}.forces; + auto forces = decltype(observables::forces_stress{ions, electrons, ham, energy}.forces){}; + if(ion_propagator.needs_force()) forces = observables::forces_stress{ions, electrons, ham, energy}.forces; auto current = vector3{0.0, 0.0, 0.0}; if(sc.has_induced_vector_potential()) current = observables::current(ions, electrons, ham); @@ -98,7 +98,7 @@ void propagate(systems::ions & ions, systems::electrons & electrons, ProcessFunc energy.calculate(ham, electrons); - if(ion_propagator.needs_force()) forces = observables::forces_stress{ions, electrons, ham}.forces; + if(ion_propagator.needs_force()) forces = observables::forces_stress{ions, electrons, ham, energy}.forces; //propagate ionic velocities to t + dt ion_propagator.propagate_velocities(dt, ions, forces); From caaab6ecc3dad6bf2dbd042efeb16aad52939138 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 14:51:45 -0700 Subject: [PATCH 03/27] The diagonal XC term of the stress. --- src/observables/forces_stress.hpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 1c01848c..c22d4265 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -58,7 +58,7 @@ struct forces_stress { void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham, Energy const & energy){ CALI_CXX_MARK_FUNCTION; - + basis::field> gdensity(electrons.density_basis()); gdensity.fill(vector3{0.0, 0.0, 0.0}); @@ -130,7 +130,15 @@ struct forces_stress { forces[iatom] = ionic_forces[iatom] + forces_local[iatom] + forces_non_local[iatom]; } - // MISSING: the non-linear core correction term + // MISSING: the non-linear core correction term to the force + + // THE XC CONTRIBUTION TO THE STRESS + for(int alpha = 0; alpha < 3; alpha++) { + stress(alpha, alpha) += energy.xc() - energy.nvxc(); + } + + //missing: the XC gradient term + } }; From eab298407ae4778002530b7b266a597753568d1a Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 14:58:21 -0700 Subject: [PATCH 04/27] Add links with the formulas for forces and stress. --- src/observables/forces_stress.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index c22d4265..c6fdd1b1 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -56,7 +56,11 @@ struct forces_stress { template void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham, Energy const & energy){ - + // This function calculates the force and the stress. Sources: + // - Force: Eq. (2.40) of https://digital.csic.es/bitstream/10261/44512/1/xandrade_phd.pdf + // - Stress formulas: Eq. (33) of https://arxiv.org/pdf/1809.08157 + + CALI_CXX_MARK_FUNCTION; basis::field> gdensity(electrons.density_basis()); From af7a6dffabd7b10d31d38125adf3000700302c8f Mon Sep 17 00:00:00 2001 From: Alfredo Correa Date: Mon, 7 Oct 2024 15:19:58 -0700 Subject: [PATCH 05/27] reintroduce electron gas test to validate stress from finite differences --- tests/electron_gas_stress.cpp | 54 +++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 tests/electron_gas_stress.cpp diff --git a/tests/electron_gas_stress.cpp b/tests/electron_gas_stress.cpp new file mode 100644 index 00000000..796feb96 --- /dev/null +++ b/tests/electron_gas_stress.cpp @@ -0,0 +1,54 @@ +/* -*- indent-tabs-mode: t -*- */ + +/* + Copyright (C) 2019 Xavier Andrade + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA +*/ + +#include +#include +#include +#include +#include + +#include + +int main(int argc, char ** argv){ + + using namespace inq; + using namespace inq::magnitude; + + input::environment env; + + utils::match energy_match(1.0e-5); + + { + auto box = systems::cell::cubic(10.0_b); + + systems::ions ions(box); + + systems::electrons electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).temperature(300.0_K).extra_electrons(14.0).extra_states(2), input::kpoints::grid({1, 1, 1}, false)); + + ground_state::initial_guess(ions, electrons); + + auto result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.energy_tolerance(1e-9_Ha)); + + energy_match.check("total energy", result.energy.total() , -0.567967321401); + energy_match.check("kinetic energy", result.energy.kinetic() , 2.485678165423); +// energy_match.check("hartree", result.energy.hartree , 0.000000732036); +// energy_match.check("XC energy", result.energy.xc , -3.053646218860); + } +} From 7cb683c17759b858073d2a45a5875e6f81f04de0 Mon Sep 17 00:00:00 2001 From: Alfredo Correa Date: Mon, 7 Oct 2024 15:55:49 -0700 Subject: [PATCH 06/27] convert energy calculation into a function --- tests/electron_gas_stress.cpp | 42 +++++++++++++++++++++-------------- 1 file changed, 25 insertions(+), 17 deletions(-) diff --git a/tests/electron_gas_stress.cpp b/tests/electron_gas_stress.cpp index 796feb96..20ce6f7a 100644 --- a/tests/electron_gas_stress.cpp +++ b/tests/electron_gas_stress.cpp @@ -26,7 +26,23 @@ #include -int main(int argc, char ** argv){ +template +auto total_energy(inq::input::environment& env, LatticeParameters lat) { + using namespace inq; + using namespace inq::magnitude; + + auto const box = systems::cell::cubic(lat); + + systems::ions const ions(box); + + systems::electrons electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).temperature(300.0_K).extra_electrons(14.0).extra_states(2), input::kpoints::grid({1, 1, 1}, false)); + + ground_state::initial_guess(ions, electrons); + + return ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.energy_tolerance(1e-9_Ha)); +} + +int main(int argc, char** argv){ using namespace inq; using namespace inq::magnitude; @@ -35,20 +51,12 @@ int main(int argc, char ** argv){ utils::match energy_match(1.0e-5); - { - auto box = systems::cell::cubic(10.0_b); - - systems::ions ions(box); - - systems::electrons electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).temperature(300.0_K).extra_electrons(14.0).extra_states(2), input::kpoints::grid({1, 1, 1}, false)); - - ground_state::initial_guess(ions, electrons); - - auto result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.energy_tolerance(1e-9_Ha)); - - energy_match.check("total energy", result.energy.total() , -0.567967321401); - energy_match.check("kinetic energy", result.energy.kinetic() , 2.485678165423); -// energy_match.check("hartree", result.energy.hartree , 0.000000732036); -// energy_match.check("XC energy", result.energy.xc , -3.053646218860); - } + auto const result = total_energy(env, 10.0_b); + + energy_match.check("total energy", result.energy.total() , -0.567967321401); + energy_match.check("kinetic energy", result.energy.kinetic() , 2.485678165423); + + // energy_match.check("hartree", result.energy.hartree , 0.000000732036); + // energy_match.check("XC energy", result.energy.xc , -3.053646218860); + } From 8620eb6d7dc6389fab702d0ff6d678a86e8e8d40 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 15:16:49 -0700 Subject: [PATCH 07/27] Divide by the cell volume at the end, and symmetrize. --- src/observables/forces_stress.hpp | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index c6fdd1b1..ede00977 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -91,7 +91,7 @@ struct forces_stress { for(auto index = 0; index < 6; index++) { int alpha, beta; stress_component(index, alpha, beta); - stress[alpha][beta] += -2.0/gphi.basis().cell().volume()*real(stress_kinetic[index]); + stress[alpha][beta] += -2.0*real(stress_kinetic[index]); } iphi++; @@ -142,6 +142,22 @@ struct forces_stress { } //missing: the XC gradient term + + //MISSING: stress electrostatic term + + // Divide by the cell volume + for(auto alpha = 0; alpha < 3; alpha++){ + for(auto beta = 0; beta < 3; beta++){ + stress[alpha][beta] /= electrons.density_basis().cell().volume(); + } + } + + //symmetrize + for(auto alpha = 0; alpha < 3; alpha++){ + for(auto beta = alpha + 1; beta < 3; beta++){ + stress[alpha][beta] = stress[beta][alpha]; + } + } } From e04b6bc769592e93afd77b8e1af4b18d0b9faa91 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 15:21:10 -0700 Subject: [PATCH 08/27] Add the stress to ground_state::results. --- src/ground_state/calculator.hpp | 4 +++- src/ground_state/results.hpp | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/ground_state/calculator.hpp b/src/ground_state/calculator.hpp index 9a445eb8..2b8c5d29 100644 --- a/src/ground_state/calculator.hpp +++ b/src/ground_state/calculator.hpp @@ -238,7 +238,9 @@ class calculator { auto normres = res.energy.calculate(ham_, electrons); if(solver_.calc_forces() and electrons.states().spinor_dim() == 1) { - res.forces = observables::forces_stress{ions_, electrons, ham_, res.energy}.forces; + auto fas = observables::forces_stress{ions_, electrons, ham_, res.energy}; + res.forces = fas.forces; + res.stress = fas.stress; } if(solver_.calc_forces() and electrons.states().spinor_dim() == 2) { diff --git a/src/ground_state/results.hpp b/src/ground_state/results.hpp index b69c4fbc..49de47b5 100644 --- a/src/ground_state/results.hpp +++ b/src/ground_state/results.hpp @@ -25,6 +25,7 @@ struct results { vector3 magnetization; energy_type energy; forces_type forces; + gpu::array stress; void save(parallel::communicator & comm, std::string const & dirname) const { auto error_message = "INQ error: Cannot save the ground_state::results to directory '" + dirname + "'."; From fccbd35cd55e1a2a6d4292ecefce86bc8f9ef1eb Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 15:41:15 -0700 Subject: [PATCH 09/27] Use a vector3 of vector3s for the stress. --- src/ground_state/results.hpp | 2 +- src/observables/forces_stress.hpp | 14 ++++++++++---- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/src/ground_state/results.hpp b/src/ground_state/results.hpp index 49de47b5..2fe99903 100644 --- a/src/ground_state/results.hpp +++ b/src/ground_state/results.hpp @@ -25,7 +25,7 @@ struct results { vector3 magnetization; energy_type energy; forces_type forces; - gpu::array stress; + vector3> stress; void save(parallel::communicator & comm, std::string const & dirname) const { auto error_message = "INQ error: Cannot save the ground_state::results to directory '" + dirname + "'."; diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index ede00977..7b17b80b 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -21,14 +21,13 @@ namespace observables { struct forces_stress { gpu::array, 1> forces; - gpu::array stress; + vector3> stress; forces_stress() = default; template forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham, Energy const & energy): - forces(ions.size()), - stress({3, 3}, 0.0) + forces(ions.size()) { calculate(ions, electrons, ham, energy); } @@ -63,6 +62,13 @@ struct forces_stress { CALI_CXX_MARK_FUNCTION; + // SET THE STRESS TO ZERO + for(auto alpha = 0; alpha < 3; alpha++){ + for(auto beta = 0; beta < 3; beta++){ + stress[alpha][beta] = 0.0; + } + } + basis::field> gdensity(electrons.density_basis()); gdensity.fill(vector3{0.0, 0.0, 0.0}); @@ -138,7 +144,7 @@ struct forces_stress { // THE XC CONTRIBUTION TO THE STRESS for(int alpha = 0; alpha < 3; alpha++) { - stress(alpha, alpha) += energy.xc() - energy.nvxc(); + stress[alpha][alpha] += energy.xc() - energy.nvxc(); } //missing: the XC gradient term From 210a3b5a099e51623e4f3730fc71977798eb5eae Mon Sep 17 00:00:00 2001 From: Alfredo Correa Date: Mon, 7 Oct 2024 16:12:30 -0700 Subject: [PATCH 10/27] calculate two volumes --- tests/electron_gas_stress.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/tests/electron_gas_stress.cpp b/tests/electron_gas_stress.cpp index 20ce6f7a..87dc1649 100644 --- a/tests/electron_gas_stress.cpp +++ b/tests/electron_gas_stress.cpp @@ -51,12 +51,16 @@ int main(int argc, char** argv){ utils::match energy_match(1.0e-5); - auto const result = total_energy(env, 10.0_b); + auto const lat1_in_bohr = 10.0; + auto const vol1_in_bohr3 = std::pow(lat1_in_bohr, 3); - energy_match.check("total energy", result.energy.total() , -0.567967321401); - energy_match.check("kinetic energy", result.energy.kinetic() , 2.485678165423); + auto const result1 = total_energy(env, lat1_in_bohr*1.0_b); - // energy_match.check("hartree", result.energy.hartree , 0.000000732036); - // energy_match.check("XC energy", result.energy.xc , -3.053646218860); + energy_match.check("total energy", result1.energy.total() , -0.567967321401); + energy_match.check("kinetic energy", result1.energy.kinetic() , 2.485678165423); + auto const lat2_in_bohr = 10.01; + auto const vol2_in_bohr3 = std::pow(lat2_in_bohr, 3); + + auto const result2 = total_energy(env, lat2_in_bohr*1.0_b); } From 5c75a1d7c27abbcefbfd6551fbb382d9cb36436a Mon Sep 17 00:00:00 2001 From: Alfredo Correa Date: Mon, 7 Oct 2024 16:35:11 -0700 Subject: [PATCH 11/27] check result from finite difference, electrostatic term seems to be missing --- tests/electron_gas_stress.cpp | 87 +++++++++++++++++++++++------------ 1 file changed, 58 insertions(+), 29 deletions(-) diff --git a/tests/electron_gas_stress.cpp b/tests/electron_gas_stress.cpp index 87dc1649..06799e55 100644 --- a/tests/electron_gas_stress.cpp +++ b/tests/electron_gas_stress.cpp @@ -7,60 +7,89 @@ it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. - + This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. - + You should have received a copy of the GNU Lesser General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */ -#include +#include +#include #include +#include #include -#include -#include #include -template -auto total_energy(inq::input::environment& env, LatticeParameters lat) { - using namespace inq; - using namespace inq::magnitude; +template +auto total_energy(inq::input::environment &env, LatticeParameters lat) { + using namespace inq; + using namespace inq::magnitude; - auto const box = systems::cell::cubic(lat); + auto const box = systems::cell::cubic(lat); - systems::ions const ions(box); + systems::ions const ions(box); - systems::electrons electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).temperature(300.0_K).extra_electrons(14.0).extra_states(2), input::kpoints::grid({1, 1, 1}, false)); + systems::electrons electrons(env.par(), ions, + options::electrons{} + .cutoff(30.0_Ha) + .temperature(300.0_K) + .extra_electrons(14.0) + .extra_states(2), + input::kpoints::grid({1, 1, 1}, false)); - ground_state::initial_guess(ions, electrons); + ground_state::initial_guess(ions, electrons); - return ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.energy_tolerance(1e-9_Ha)); + return ground_state::calculate(ions, electrons, options::theory{}.lda(), + inq::options::ground_state{} + .energy_tolerance(1e-9_Ha) + .calculate_forces()); } -int main(int argc, char** argv){ +int main(int argc, char **argv) { + + using namespace inq; + using namespace inq::magnitude; + + input::environment env; + + utils::match energy_match(1.0e-5); + + auto const lat1_in_bohr = 10.0; + auto const vol1_in_bohr3 = std::pow(lat1_in_bohr, 3); + + auto const result1 = total_energy(env, lat1_in_bohr * 1.0_b); - using namespace inq; - using namespace inq::magnitude; - - input::environment env; - - utils::match energy_match(1.0e-5); + energy_match.check("total energy", result1.energy.total(), -0.567967321401); + energy_match.check("kinetic energy", result1.energy.kinetic(), + 2.485678165423); - auto const lat1_in_bohr = 10.0; - auto const vol1_in_bohr3 = std::pow(lat1_in_bohr, 3); + auto const lat2_in_bohr = 10.01; + auto const vol2_in_bohr3 = std::pow(lat2_in_bohr, 3); - auto const result1 = total_energy(env, lat1_in_bohr*1.0_b); + auto const result2 = total_energy(env, lat2_in_bohr * 1.0_b); - energy_match.check("total energy", result1.energy.total() , -0.567967321401); - energy_match.check("kinetic energy", result1.energy.kinetic() , 2.485678165423); + auto const stress = result1.stress; + std::cout << "stress" << std::endl; + std::cout << stress[0][0] << " " << stress[0][1] << " " << stress[0][2] + << '\n'; + std::cout << stress[1][0] << " " << stress[1][1] << " " << stress[1][2] + << '\n'; + std::cout << stress[2][0] << " " << stress[2][1] << " " << stress[2][2] + << '\n'; - auto const lat2_in_bohr = 10.01; - auto const vol2_in_bohr3 = std::pow(lat2_in_bohr, 3); + std::cout << "energies (1, 2) = " << result1.energy.total() << " " + << result2.energy.total() << std::endl; + std::cout << "volumes (1, 2) = " << vol1_in_bohr3 << " " << vol2_in_bohr3 + << std::endl; - auto const result2 = total_energy(env, lat2_in_bohr*1.0_b); + std::cout << "(energy(2) - energy(1) ) / (vo2 - vol1) = " + << (result2.energy.total() - result1.energy.total()) / + (vol2_in_bohr3 - vol1_in_bohr3) + << std::endl; } From 44049e3ebb5e26284e4ebe71f049294f3d94d4fc Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 16:22:34 -0700 Subject: [PATCH 12/27] Print the stress and pressure in results. --- src/ground_state/results.hpp | 7 ++++++- tests/diamond.cpp | 4 +++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/ground_state/results.hpp b/src/ground_state/results.hpp index 2fe99903..32ed8eb8 100644 --- a/src/ground_state/results.hpp +++ b/src/ground_state/results.hpp @@ -65,8 +65,13 @@ struct results { std::cout << "Ground-state results:\n"; std::cout << " iterations = " << self.total_iter << '\n'; - std::cout << " dipole = " << self.dipole << '\n'; + std::cout << " dipole [a.u.] = " << self.dipole << '\n'; std::cout << " magnetization = " << self.magnetization << '\n'; + std::cout << " stress [a.u.] = " << self.stress[0] << '\n'; + std::cout << " = " << self.stress[1] << '\n'; + std::cout << " = " << self.stress[2] << '\n'; + auto pressure = -(self.stress[0][0] + self.stress[1][1] + self.stress[2][2])/3.0; + std::cout << " pressure = " << pressure << " Ha/b^3 | " << pressure*29421.016 << " GPa \n"; std::cout << " total energy = " << utils::num_to_str("%.8f", self.energy.total()) << " Ha | " << utils::num_to_str("%.8f", self.energy.total()/in_atomic_units(1.0_eV)) << " eV \n"; diff --git a/tests/diamond.cpp b/tests/diamond.cpp index 7a7a8c61..15b1236a 100644 --- a/tests/diamond.cpp +++ b/tests/diamond.cpp @@ -95,8 +95,10 @@ int main(int argc, char ** argv){ if(energy_match.fail()) return energy_match.fail(); { - auto result = ground_state::calculate(ions, electrons, options::theory{}.hartree_fock(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1e-8_Ha)); + auto result = ground_state::calculate(ions, electrons, options::theory{}.hartree_fock(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1e-8_Ha).calculate_forces()); + std::cout << result << std::endl; + energy_match.check("total energy", result.energy.total(), -9.788709725748); energy_match.check("kinetic energy", result.energy.kinetic(), 8.151819376871); energy_match.check("eigenvalues", result.energy.eigenvalues(), -0.249826944617); From 56478529d874648dc2988d8424d808f6fdb6ff7e Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 16:41:55 -0700 Subject: [PATCH 13/27] Save and load the stress in results. --- src/ground_state/results.hpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/ground_state/results.hpp b/src/ground_state/results.hpp index 32ed8eb8..123ba73c 100644 --- a/src/ground_state/results.hpp +++ b/src/ground_state/results.hpp @@ -36,7 +36,8 @@ struct results { utils::save_value(comm, dirname + "/dipole", dipole, error_message); utils::save_value(comm, dirname + "/magnetization", magnetization, error_message); utils::save_value(comm, dirname + "/num_atoms", forces.size(), error_message); - utils::save_container(comm, dirname + "/forces", forces, error_message); + utils::save_value(comm, dirname + "/stress", stress, error_message); + utils::save_container(comm, dirname + "/forces", forces, error_message); } static auto load(std::string const & dirname) { @@ -48,7 +49,8 @@ struct results { utils::load_value(dirname + "/total_iter", res.total_iter, error_message); utils::load_value(dirname + "/dipole", res.dipole, error_message); utils::load_value(dirname + "/magnetization", res.magnetization, error_message); - + utils::load_value(dirname + "/stress", res.stress, error_message); + int num_atoms; utils::load_value(dirname + "/num_atoms", num_atoms, error_message); @@ -111,6 +113,9 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { res.energy.xc(7.55); res.energy.nvxc(8.55); res.energy.exact_exchange(10.55); + res.stress[0] = vector3{1.28, 7.62, 5.56}; + res.stress[1] = vector3{7.62, 3.03, 0.57}; + res.stress[2] = vector3{5.56, 0.57, 8.88};; res.forces = ground_state::results::forces_type{65, vector3{3.55, 4.55, 5.55}}; CHECK(res.total_iter == 333); @@ -125,6 +130,9 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { CHECK(res.energy.xc() == 7.55); CHECK(res.energy.nvxc() == 8.55); CHECK(res.energy.exact_exchange() == 10.55); + CHECK(res.stress[0] == vector3{1.28, 7.62, 5.56}); + CHECK(res.stress[1] == vector3{7.62, 3.03, 0.57}); + CHECK(res.stress[2] == vector3{5.56, 0.57, 8.88}); CHECK(res.forces == ground_state::results::forces_type{65, vector3{3.55, 4.55, 5.55}}); res.save(comm, "save_results"); @@ -142,6 +150,9 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { CHECK(read_res.energy.xc() == 7.55); CHECK(read_res.energy.nvxc() == 8.55); CHECK(read_res.energy.exact_exchange() == 10.55); + CHECK(read_res.stress[0] == vector3{1.28, 7.62, 5.56}); + CHECK(read_res.stress[1] == vector3{7.62, 3.03, 0.57}); + CHECK(read_res.stress[2] == vector3{5.56, 0.57, 8.88}); CHECK(read_res.forces == ground_state::results::forces_type{65, vector3{3.55, 4.55, 5.55}}); } From 90ff9a510578701e724106d2f91ed0b31633d94c Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 7 Oct 2024 23:16:20 -0700 Subject: [PATCH 14/27] Convert the gradient to cartesian to calculate the stress. --- src/observables/forces_stress.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 7b17b80b..18d83cf7 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -84,10 +84,11 @@ struct forces_stress { //STRESS KINETIC auto stress_kinetic = gpu::run(6, gpu::reduce(gphi.local_set_size()), gpu::reduce(gphi.basis().local_size()), - [gph = begin(gphi.matrix()), occ = begin(electrons.occupations()[iphi])] GPU_LAMBDA (auto index, auto ist, auto ip) { + [metric = ions.cell().metric(), gph = begin(gphi.matrix()), occ = begin(electrons.occupations()[iphi])] GPU_LAMBDA (auto index, auto ist, auto ip) { int alpha, beta; stress_component(index, alpha, beta); - return occ[ist]*conj(gph[ip][ist][alpha])*gph[ip][ist][beta]; + auto grad_cart = metric.to_cartesian(gph[ip][ist]); + return occ[ist]*conj(grad_cart[alpha])*grad_cart[beta]; }); if(gphi.full_comm().size() > 1){ From 14ea9daccbcee9fb3f54c98aa68d7c4cf579ae86 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 14 Oct 2024 15:08:59 -0700 Subject: [PATCH 15/27] Add a (non-compiling) implementation of the electrostatic stress. --- src/observables/forces_stress.hpp | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 18d83cf7..d2f69163 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -52,6 +52,32 @@ struct forces_stress { beta = 2; } } + + template + auto stress_electrostatic(Density const & density) { + + auto potential = solvers::poisson::solve(density); + auto efield = operations::gradient(potential); + + auto stress1d = gpu::run(6, gpu::reduce(efield.basis().local_size()), + [metric = efield.basis().cell().metric(), ef = begin(efield.linear())] GPU_LAMBDA (auto index, auto ip) { + int alpha, beta; + stress_component(index, alpha, beta); + auto ef_cart = metric.to_cartesian(ef[ip]); + return ef_cart[alpha]*ef_cart[beta]; + }); + + vector3> stress; + + for(auto index = 0; index < 6; index++) { + int alpha, beta; + stress_component(index, alpha, beta); + stress[alpha][beta] = stress1d[index]/(4.0*M_PI); + if(beta /= alpha) stress[beta][alpha] = stress1d[index]/(4.0*M_PI); + } + + return stress; + } template void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham, Energy const & energy){ @@ -150,7 +176,7 @@ struct forces_stress { //missing: the XC gradient term - //MISSING: stress electrostatic term + stress += stress_electrostatic(electrons.density()); // Divide by the cell volume for(auto alpha = 0; alpha < 3; alpha++){ From ccfa92477cdf19d7efef99a393c4caa32290e405 Mon Sep 17 00:00:00 2001 From: Alfredo Correa Date: Tue, 15 Oct 2024 14:50:13 -0700 Subject: [PATCH 16/27] hypercubic is adding too much constness to the subarray reference --- src/basis/field.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/basis/field.hpp b/src/basis/field.hpp index 5e860569..4645977f 100644 --- a/src/basis/field.hpp +++ b/src/basis/field.hpp @@ -77,7 +77,7 @@ namespace basis { linear_ = parallel::get_remote_points(old, rem_points); } - explicit field(const field & coeff) = default; //avoid unadverted copies + explicit field(const field & coeff) = default; //avoid unadverted copies field(field && coeff) = default; field & operator=(const field & coeff) = default; field & operator=(field && coeff) = default; @@ -134,7 +134,7 @@ namespace basis { // emulate a field_set auto hypercubic() const { - return cubic().template reinterpret_array_cast(1); + return cubic().template reinterpret_array_cast(1); } auto hypercubic() { @@ -193,7 +193,7 @@ namespace basis { field complex_field(field const & rfield) { - field cfield(rfield.skeleton()); + field cfield(rfield.skeleton()); gpu::run(rfield.basis().part().local_size(), [cp = begin(cfield.linear()), rp = begin(rfield.linear())] GPU_LAMBDA (auto ip){ @@ -216,14 +216,14 @@ field> complex_field(field } field real_field(field const & cfield) { - field rfield(cfield.skeleton()); + field rfield(cfield.skeleton()); rfield.linear() = boost::multi::blas::real(cfield.linear()); return rfield; } template field> real_field(field> const & cfield) { - field> rfield(cfield.skeleton()); + field> rfield(cfield.skeleton()); gpu::run(3, cfield.basis().part().local_size(), [rp = begin(rfield.linear()), cp = begin(cfield.linear())] GPU_LAMBDA (auto idir, auto ip){ @@ -248,7 +248,7 @@ field> real_field(field ff_copy(ff.skeleton()); @@ -299,7 +299,7 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){ CHECK(std::get<1>(sizes(ff.hypercubic())) == 11); CHECK(std::get<2>(sizes(ff.hypercubic())) == 20); - CHECK(std::get<3>(sizes(ff.hypercubic())) == 1); + CHECK(std::get<3>(sizes(ff.hypercubic())) == 1); //Make sure the hypercubic array is correctly ordered, so it can be flattened auto strd = strides(ff.hypercubic()); From 39827e49b7428e92ef22ea7f4b79ee6bb2eb01db Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Tue, 15 Oct 2024 19:09:56 -0700 Subject: [PATCH 17/27] Fix compilation of the gradient for field. --- src/operations/gradient.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operations/gradient.hpp b/src/operations/gradient.hpp index 08dd40b5..b778991b 100644 --- a/src/operations/gradient.hpp +++ b/src/operations/gradient.hpp @@ -40,7 +40,7 @@ ResultType gradient(FieldSetType const & ff, double factor = 1.0, vector3 Date: Wed, 16 Oct 2024 02:02:01 -0700 Subject: [PATCH 18/27] Fix cuda compilation. --- src/observables/forces_stress.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index d2f69163..2289ebc2 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -54,7 +54,7 @@ struct forces_stress { } template - auto stress_electrostatic(Density const & density) { + vector3> stress_electrostatic(Density const & density) { auto potential = solvers::poisson::solve(density); auto efield = operations::gradient(potential); From e21c63788ef861550115de7fe7004ba829ab7a87 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 16 Oct 2024 12:47:24 -0700 Subject: [PATCH 19/27] Corrected not equal operator. This is not Fortran! --- src/observables/forces_stress.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 2289ebc2..5b859b2d 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -73,7 +73,7 @@ struct forces_stress { int alpha, beta; stress_component(index, alpha, beta); stress[alpha][beta] = stress1d[index]/(4.0*M_PI); - if(beta /= alpha) stress[beta][alpha] = stress1d[index]/(4.0*M_PI); + if(beta != alpha) stress[beta][alpha] = stress1d[index]/(4.0*M_PI); } return stress; From 4913050432fa8c00d3cf198f515c3d168a27b147 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 16 Oct 2024 16:44:15 -0700 Subject: [PATCH 20/27] Move the kinetic stress to its own function. --- src/observables/forces_stress.hpp | 46 ++++++++++++++++++------------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 5b859b2d..7b1d6251 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -53,6 +53,31 @@ struct forces_stress { } } + template + vector3> stress_kinetic(GPhi const & gphi, Occupations const & occupations) { + + auto stress1d = gpu::run(6, gpu::reduce(gphi.local_set_size()), gpu::reduce(gphi.basis().local_size()), + [metric = gphi.basis().cell().metric(), gph = begin(gphi.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto index, auto ist, auto ip) { + int alpha, beta; + stress_component(index, alpha, beta); + auto grad_cart = metric.to_cartesian(gph[ip][ist]); + return occ[ist]*real(conj(grad_cart[alpha])*grad_cart[beta]); + }); + + if(gphi.full_comm().size() > 1) gphi.full_comm().all_reduce_n(raw_pointer_cast(stress1d.data_elements()), 6);; + + vector3> stress; + + for(auto index = 0; index < 6; index++) { + int alpha, beta; + stress_component(index, alpha, beta); + stress[alpha][beta] = stress1d[index]; + if(beta != alpha) stress[beta][alpha] = stress1d[index]; + } + + return -2.0*stress; + } + template vector3> stress_electrostatic(Density const & density) { @@ -108,25 +133,8 @@ struct forces_stress { ham.projectors_all().force(phi, gphi, ions.cell().metric(), electrons.occupations()[iphi], ham.uniform_vector_potential(), forces_non_local); - //STRESS KINETIC - auto stress_kinetic = gpu::run(6, gpu::reduce(gphi.local_set_size()), gpu::reduce(gphi.basis().local_size()), - [metric = ions.cell().metric(), gph = begin(gphi.matrix()), occ = begin(electrons.occupations()[iphi])] GPU_LAMBDA (auto index, auto ist, auto ip) { - int alpha, beta; - stress_component(index, alpha, beta); - auto grad_cart = metric.to_cartesian(gph[ip][ist]); - return occ[ist]*conj(grad_cart[alpha])*grad_cart[beta]; - }); - - if(gphi.full_comm().size() > 1){ - gphi.full_comm().all_reduce_n(raw_pointer_cast(stress_kinetic.data_elements()), 6);; - } - - for(auto index = 0; index < 6; index++) { - int alpha, beta; - stress_component(index, alpha, beta); - stress[alpha][beta] += -2.0*real(stress_kinetic[index]); - } - + stress += stress_kinetic(gphi, electrons.occupations()[iphi]); + iphi++; } From 87c89f64775d8cef590a31148129808dacd18ef5 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 16 Oct 2024 16:48:53 -0700 Subject: [PATCH 21/27] Define a function to build the tensor from a 1D representation. --- src/observables/forces_stress.hpp | 36 ++++++++++++++----------------- 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 7b1d6251..cbacd9b8 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -53,6 +53,20 @@ struct forces_stress { } } + template + auto tensor(Stress1D const & stress1d) { + vector3> stress; + + for(auto index = 0; index < 6; index++) { + int alpha, beta; + stress_component(index, alpha, beta); + stress[alpha][beta] = stress1d[index]; + if(beta != alpha) stress[beta][alpha] = stress1d[index]; + } + + return stress; + } + template vector3> stress_kinetic(GPhi const & gphi, Occupations const & occupations) { @@ -66,16 +80,7 @@ struct forces_stress { if(gphi.full_comm().size() > 1) gphi.full_comm().all_reduce_n(raw_pointer_cast(stress1d.data_elements()), 6);; - vector3> stress; - - for(auto index = 0; index < 6; index++) { - int alpha, beta; - stress_component(index, alpha, beta); - stress[alpha][beta] = stress1d[index]; - if(beta != alpha) stress[beta][alpha] = stress1d[index]; - } - - return -2.0*stress; + return -2.0*tensor(stress1d); } template @@ -92,16 +97,7 @@ struct forces_stress { return ef_cart[alpha]*ef_cart[beta]; }); - vector3> stress; - - for(auto index = 0; index < 6; index++) { - int alpha, beta; - stress_component(index, alpha, beta); - stress[alpha][beta] = stress1d[index]/(4.0*M_PI); - if(beta != alpha) stress[beta][alpha] = stress1d[index]/(4.0*M_PI); - } - - return stress; + return tensor(stress1d)/(4.0*M_PI); } template From 71c179989b6a514d0fbbe8f94deb83dd96037eb3 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 16 Oct 2024 17:00:53 -0700 Subject: [PATCH 22/27] Do not symmetrize the stress tensor. --- src/observables/forces_stress.hpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index cbacd9b8..2c7a3510 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -188,13 +188,6 @@ struct forces_stress { stress[alpha][beta] /= electrons.density_basis().cell().volume(); } } - - //symmetrize - for(auto alpha = 0; alpha < 3; alpha++){ - for(auto beta = alpha + 1; beta < 3; beta++){ - stress[alpha][beta] = stress[beta][alpha]; - } - } } From 1295b2d3e12beda286715894d53b434b7b1a3279 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 16 Oct 2024 17:08:32 -0700 Subject: [PATCH 23/27] The volume element was missing from the integrals. --- src/observables/forces_stress.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 2c7a3510..30343729 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -80,7 +80,7 @@ struct forces_stress { if(gphi.full_comm().size() > 1) gphi.full_comm().all_reduce_n(raw_pointer_cast(stress1d.data_elements()), 6);; - return -2.0*tensor(stress1d); + return -2.0*gphi.basis().volume_element()*tensor(stress1d); } template @@ -97,7 +97,7 @@ struct forces_stress { return ef_cart[alpha]*ef_cart[beta]; }); - return tensor(stress1d)/(4.0*M_PI); + return density.basis().volume_element()/(4.0*M_PI)*tensor(stress1d); } template From a06f0a143903aa0a0a3ae6c9328d3fcc913205e5 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 16 Oct 2024 17:34:31 -0700 Subject: [PATCH 24/27] Remove the factor 2 that comes from occupations. --- src/observables/forces_stress.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 30343729..3367a9cf 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -80,7 +80,7 @@ struct forces_stress { if(gphi.full_comm().size() > 1) gphi.full_comm().all_reduce_n(raw_pointer_cast(stress1d.data_elements()), 6);; - return -2.0*gphi.basis().volume_element()*tensor(stress1d); + return -gphi.basis().volume_element()*tensor(stress1d); } template From d0b3bf07e7a87e912d2ef4bb9a28f3346f5a7a02 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 21 Oct 2024 13:23:35 -0700 Subject: [PATCH 25/27] Check the calculation of the stress tensor for the electron gas. This is working now. --- tests/electron_gas_stress.cpp | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/tests/electron_gas_stress.cpp b/tests/electron_gas_stress.cpp index 06799e55..127036ba 100644 --- a/tests/electron_gas_stress.cpp +++ b/tests/electron_gas_stress.cpp @@ -58,22 +58,18 @@ int main(int argc, char **argv) { input::environment env; - utils::match energy_match(1.0e-5); + utils::match match(1.0e-5); auto const lat1_in_bohr = 10.0; auto const vol1_in_bohr3 = std::pow(lat1_in_bohr, 3); - auto const result1 = total_energy(env, lat1_in_bohr * 1.0_b); - - energy_match.check("total energy", result1.energy.total(), -0.567967321401); - energy_match.check("kinetic energy", result1.energy.kinetic(), - 2.485678165423); + auto const result1 = total_energy(env, lat1_in_bohr*1.0_b); auto const lat2_in_bohr = 10.01; auto const vol2_in_bohr3 = std::pow(lat2_in_bohr, 3); - auto const result2 = total_energy(env, lat2_in_bohr * 1.0_b); - + auto const result2 = total_energy(env, lat2_in_bohr*1.0_b); + auto const stress = result1.stress; std::cout << "stress" << std::endl; std::cout << stress[0][0] << " " << stress[0][1] << " " << stress[0][2] @@ -89,7 +85,25 @@ int main(int argc, char **argv) { << std::endl; std::cout << "(energy(2) - energy(1) ) / (vo2 - vol1) = " - << (result2.energy.total() - result1.energy.total()) / - (vol2_in_bohr3 - vol1_in_bohr3) + << (result2.energy.total() - result1.energy.total())/(vol2_in_bohr3 - vol1_in_bohr3) << std::endl; + + + match.check("total energy 1", result1.energy.total(), -0.684940553454); + match.check("total energy 2", result2.energy.total(), -0.686911961588); + match.check("finite diff pressure", (result2.energy.total() - result1.energy.total())/(vol2_in_bohr3 - vol1_in_bohr3), -0.000656479347, 1e-7); + + match.check("stress xx", stress[0][0], -0.000658598015, 1e-8); + match.check("stress yy", stress[1][1], -0.000658598015, 1e-8); + match.check("stress zz", stress[2][2], -0.000658598015, 1e-8); + + match.check("stress xy", stress[0][1], 0.0, 1e-8); + match.check("stress xz", stress[0][2], 0.0, 1e-8); + match.check("stress yz", stress[1][2], 0.0, 1e-8); + + match.check("stress yx", stress[1][0], 0.0, 1e-8); + match.check("stress zx", stress[2][0], 0.0, 1e-8); + match.check("stress zy", stress[2][1], 0.0, 1e-8); + + return match.fail(); } From 457de92721651d9c7ed98d066d88d555b330920a Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Tue, 22 Oct 2024 23:06:33 -0700 Subject: [PATCH 26/27] Add initial values for the stress tensor reductions. --- src/observables/forces_stress.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 3367a9cf..4ca2829b 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -70,7 +70,7 @@ struct forces_stress { template vector3> stress_kinetic(GPhi const & gphi, Occupations const & occupations) { - auto stress1d = gpu::run(6, gpu::reduce(gphi.local_set_size()), gpu::reduce(gphi.basis().local_size()), + auto stress1d = gpu::run(6, gpu::reduce(gphi.local_set_size()), gpu::reduce(gphi.basis().local_size()), 0.0, [metric = gphi.basis().cell().metric(), gph = begin(gphi.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto index, auto ist, auto ip) { int alpha, beta; stress_component(index, alpha, beta); @@ -89,7 +89,7 @@ struct forces_stress { auto potential = solvers::poisson::solve(density); auto efield = operations::gradient(potential); - auto stress1d = gpu::run(6, gpu::reduce(efield.basis().local_size()), + auto stress1d = gpu::run(6, gpu::reduce(efield.basis().local_size()), 0.0, [metric = efield.basis().cell().metric(), ef = begin(efield.linear())] GPU_LAMBDA (auto index, auto ip) { int alpha, beta; stress_component(index, alpha, beta); From c27c3b81067f20493ee2a16279f6a00bd28ba551 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 23 Oct 2024 00:39:49 -0700 Subject: [PATCH 27/27] Missing reduction for the electrostatic stress term. --- src/observables/forces_stress.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 4ca2829b..c4bf8eea 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -78,7 +78,7 @@ struct forces_stress { return occ[ist]*real(conj(grad_cart[alpha])*grad_cart[beta]); }); - if(gphi.full_comm().size() > 1) gphi.full_comm().all_reduce_n(raw_pointer_cast(stress1d.data_elements()), 6);; + if(gphi.full_comm().size() > 1) gphi.full_comm().all_reduce_n(raw_pointer_cast(stress1d.data_elements()), 6); return -gphi.basis().volume_element()*tensor(stress1d); } @@ -97,6 +97,8 @@ struct forces_stress { return ef_cart[alpha]*ef_cart[beta]; }); + if(efield.basis().comm().size() > 1) efield.basis().comm().all_reduce_n(raw_pointer_cast(stress1d.data_elements()), 6); + return density.basis().volume_element()/(4.0*M_PI)*tensor(stress1d); }