From a46081d96a03ccbdfab4e16c06047c69628c5253 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Tue, 31 Dec 2024 15:57:06 -0800 Subject: [PATCH 01/12] Correct the force reference. --- 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 48ad2bd9..7d766a6f 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -117,7 +117,7 @@ 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 + // - Force: Page 41, Eq. (2.41) of https://digital.csic.es/bitstream/10261/44512/1/xandrade_phd.pdf // - Stress formulas: Eq. (33) of https://arxiv.org/pdf/1809.08157 From 2aae5551731626d377548c7cfebe64538ed887ac Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 6 Jan 2025 00:27:33 -0800 Subject: [PATCH 02/12] Return the xc potential in xc_term. --- src/hamiltonian/self_consistency.hpp | 8 +++++++- src/hamiltonian/xc_term.hpp | 29 +++++++++------------------- 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index c7dbc824..6e996ae8 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -133,7 +133,13 @@ class self_consistency { // XC double exc, nvxc; - xc_(spin_density, core_density_, vks, exc, nvxc); + auto vxc = xc_(spin_density, core_density_, exc, nvxc); + + gpu::run(vxc.local_set_size(), vxc.basis().local_size(), + [vx = begin(vxc.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip){ + vk[ip][is] += vx[ip][is]; + }); + energy.xc(exc); energy.nvxc(nvxc); diff --git a/src/hamiltonian/xc_term.hpp b/src/hamiltonian/xc_term.hpp index e3b65ba6..2f2f8802 100644 --- a/src/hamiltonian/xc_term.hpp +++ b/src/hamiltonian/xc_term.hpp @@ -87,17 +87,6 @@ class xc_term { return full_density; } - //////////////////////////////////////////////////////////////////////////////////////////// - - template - void process_potential(VXC const & vxc, VKS & vks) const { - - gpu::run(vxc.local_set_size(), vxc.basis().local_size(), - [vx = begin(vxc.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip){ - vk[ip][is] += vx[ip][is]; - }); - } - /////////////////////////////////////////////////////////////////////////////////////////// template @@ -122,20 +111,19 @@ class xc_term { //////////////////////////////////////////////////////////////////////////////////////////// - template - void operator()(SpinDensityType const & spin_density, CoreDensityType const & core_density, VKSType & vks, double & exc, double & nvxc) const { - - exc = 0.0; + template + auto operator()(SpinDensityType const & spin_density, CoreDensityType const & core_density, double & exc, double & nvxc) const { + + basis::field_set vxc(spin_density.skeleton()); + vxc.fill(0.0); + exc = 0.0; nvxc = 0.0; - if(not any_true_functional()) return; + if(not any_true_functional()) return vxc; auto full_density = process_density(spin_density, core_density); double efunc = 0.0; - basis::field_set vxc(spin_density.skeleton()); - vxc.fill(0.0); - basis::field_set vfunc(full_density.skeleton()); auto density_gradient = std::optional{}; if(any_requires_gradient()) density_gradient.emplace(operations::gradient(full_density)); @@ -149,8 +137,9 @@ class xc_term { exc += efunc; } - process_potential(vxc, vks); nvxc += compute_nvxc(spin_density, vxc); + + return vxc; } //////////////////////////////////////////////////////////////////////////////////////////// From 313e10dfc49752eb0fddec41e138c26a47a6e211 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 6 Jan 2025 00:36:49 -0800 Subject: [PATCH 03/12] Save the xc potential in hamiltonian. --- src/hamiltonian/ks_hamiltonian.hpp | 2 ++ src/hamiltonian/self_consistency.hpp | 6 +++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/hamiltonian/ks_hamiltonian.hpp b/src/hamiltonian/ks_hamiltonian.hpp index 887241b9..4ecc7f6e 100644 --- a/src/hamiltonian/ks_hamiltonian.hpp +++ b/src/hamiltonian/ks_hamiltonian.hpp @@ -42,6 +42,7 @@ class ks_hamiltonian { private: exchange_operator exchange_; + basis::field_set vxc_; basis::field_set scalar_potential_; vector3 uniform_vector_potential_; projector_all projectors_all_; @@ -91,6 +92,7 @@ class ks_hamiltonian { ks_hamiltonian(const basis::real_space & basis, ionic::brillouin const & bzone, states::ks_states const & states, atomic_potential const & pot, systems::ions const & ions, const double exchange_coefficient, bool use_ace = false): exchange_(basis.cell(), bzone, exchange_coefficient, use_ace), + vxc_(basis, states.num_density_components()), scalar_potential_(basis, states.num_density_components()), uniform_vector_potential_({0.0, 0.0, 0.0}), states_(states) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 6e996ae8..708a18f3 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -133,10 +133,10 @@ class self_consistency { // XC double exc, nvxc; - auto vxc = xc_(spin_density, core_density_, exc, nvxc); + hamiltonian.vxc_ = xc_(spin_density, core_density_, exc, nvxc); - gpu::run(vxc.local_set_size(), vxc.basis().local_size(), - [vx = begin(vxc.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip){ + gpu::run(hamiltonian.vxc_.local_set_size(), hamiltonian.vxc_.basis().local_size(), + [vx = begin(hamiltonian.vxc_.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip){ vk[ip][is] += vx[ip][is]; }); From c2c689081e72c88ca027f5d0d870ad8a681d7871 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 6 Jan 2025 01:06:23 -0800 Subject: [PATCH 04/12] Add a (dummy) function to calculate the nlcc forces. --- external_libs/pseudopod | 2 +- src/hamiltonian/atomic_potential.hpp | 13 +++++++++++++ src/hamiltonian/ks_hamiltonian.hpp | 6 ++++++ src/observables/forces_stress.hpp | 9 ++++----- 4 files changed, 24 insertions(+), 6 deletions(-) diff --git a/external_libs/pseudopod b/external_libs/pseudopod index 8aa6a94c..b21353a5 160000 --- a/external_libs/pseudopod +++ b/external_libs/pseudopod @@ -1 +1 @@ -Subproject commit 8aa6a94ca2c40ccb9bff82d626ffc3635d4af7b0 +Subproject commit b21353a5374989f49d24fe28c335c585a3420e76 diff --git a/src/hamiltonian/atomic_potential.hpp b/src/hamiltonian/atomic_potential.hpp index 9bf87e2b..724bdd4e 100644 --- a/src/hamiltonian/atomic_potential.hpp +++ b/src/hamiltonian/atomic_potential.hpp @@ -300,6 +300,19 @@ namespace hamiltonian { density.all_reduce(comm); return density; } + + //////////////////////////////////////////////////////////////////////////////////// + + template + auto nlcc_forces(Comm & comm, basis::field_set const & vxc, Ions const & ions) const { + CALI_CXX_MARK_FUNCTION; + + gpu::array, 1> forces(ions.size(), {0.0, 0.0, 0.0}); + + return forces; + } + + //////////////////////////////////////////////////////////////////////////////////// template void info(output_stream & out) const { diff --git a/src/hamiltonian/ks_hamiltonian.hpp b/src/hamiltonian/ks_hamiltonian.hpp index 4ecc7f6e..482c2d7e 100644 --- a/src/hamiltonian/ks_hamiltonian.hpp +++ b/src/hamiltonian/ks_hamiltonian.hpp @@ -230,6 +230,12 @@ class ks_hamiltonian { return uniform_vector_potential_; } + //////////////////////////////////////////////////////////////////////////////////////////// + + auto & vxc() const { + return vxc_; + } + //////////////////////////////////////////////////////////////////////////////////////////// template diff --git a/src/observables/forces_stress.hpp b/src/observables/forces_stress.hpp index 7d766a6f..79b2eec8 100644 --- a/src/observables/forces_stress.hpp +++ b/src/observables/forces_stress.hpp @@ -180,13 +180,12 @@ struct forces_stress { electrons.density_basis().comm().all_reduce_n(reinterpret_cast(raw_pointer_cast(forces_local.data_elements())), 3*forces_local.size()); } } - - + + auto forces_nlcc = electrons.atomic_pot().nlcc_forces(electrons.states_comm(), ham.vxc(), ions); + for(int iatom = 0; iatom < ions.size(); iatom++){ - forces[iatom] = ionic_forces[iatom] + forces_local[iatom] + forces_non_local[iatom]; + forces[iatom] = ionic_forces[iatom] + forces_local[iatom] + forces_non_local[iatom] + forces_nlcc[iatom]; } - - // MISSING: the non-linear core correction term to the force // THE XC CONTRIBUTION TO THE STRESS for(int alpha = 0; alpha < 3; alpha++) { From 72880bfc7feb4a694760f6e94b89241f574d4804 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 6 Jan 2025 02:30:58 -0800 Subject: [PATCH 05/12] Implementation of the NLCC forces --- external_libs/pseudopod | 2 +- src/hamiltonian/atomic_potential.hpp | 36 ++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/external_libs/pseudopod b/external_libs/pseudopod index b21353a5..8aa6a94c 160000 --- a/external_libs/pseudopod +++ b/external_libs/pseudopod @@ -1 +1 @@ -Subproject commit b21353a5374989f49d24fe28c335c585a3420e76 +Subproject commit 8aa6a94ca2c40ccb9bff82d626ffc3635d4af7b0 diff --git a/src/hamiltonian/atomic_potential.hpp b/src/hamiltonian/atomic_potential.hpp index 724bdd4e..60d575e4 100644 --- a/src/hamiltonian/atomic_potential.hpp +++ b/src/hamiltonian/atomic_potential.hpp @@ -307,7 +307,43 @@ namespace hamiltonian { auto nlcc_forces(Comm & comm, basis::field_set const & vxc, Ions const & ions) const { CALI_CXX_MARK_FUNCTION; + // NLCC Forces from Kronik et al., J. Chem. Phys. 115, 4322 (2001): Eq. 9 + gpu::array, 1> forces(ions.size(), {0.0, 0.0, 0.0}); + + auto && basis = vxc.basis(); + + parallel::partition part(ions.size(), comm); + + for(auto iatom = part.start(); iatom < part.end(); iatom++){ + + auto atom_position = ions.positions()[iatom]; + + auto & ps = pseudo_for_element(ions.species(iatom)); + + if(not ps.has_nlcc_density()) continue; + + assert(has_nlcc()); + + basis::spherical_grid sphere(basis, atom_position, ps.nlcc_density_radius()); + + auto ff = gpu::run(gpu::reduce(sphere.size()), zero>(), + [vx = begin(vxc.hypercubic()), nspin = std::min(2l, vxc.local_set_size()), sph = sphere.ref(), spline = ps.nlcc_density().function()] + GPU_LAMBDA (auto ipoint){ + auto vv = 0.0; + for(int ispin = 0; ispin < nspin; ispin++) vv += vx[sph.grid_point(ipoint)[0]][sph.grid_point(ipoint)[1]][sph.grid_point(ipoint)[2]][ispin]; + auto rr = sph.distance(ipoint); + auto grad = spline.derivative(rr)*sph.point_pos(ipoint)/rr; + return -vv*grad; + }); + + forces[iatom] = basis.volume_element()*basis.cell().metric().to_cartesian(ff); + + } + + if(comm.size() > 1) { + comm.all_reduce_n(reinterpret_cast(raw_pointer_cast(forces.data_elements())), 3*forces.size()); + } return forces; } From 5c3765719cb3bd227895c83f352c489effb707b5 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 6 Jan 2025 11:43:44 -0800 Subject: [PATCH 06/12] Fix cuda compilation. --- src/hamiltonian/atomic_potential.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hamiltonian/atomic_potential.hpp b/src/hamiltonian/atomic_potential.hpp index 60d575e4..7f2fd0bc 100644 --- a/src/hamiltonian/atomic_potential.hpp +++ b/src/hamiltonian/atomic_potential.hpp @@ -304,7 +304,7 @@ namespace hamiltonian { //////////////////////////////////////////////////////////////////////////////////// template - auto nlcc_forces(Comm & comm, basis::field_set const & vxc, Ions const & ions) const { + gpu::array, 1> nlcc_forces(Comm & comm, basis::field_set const & vxc, Ions const & ions) const { CALI_CXX_MARK_FUNCTION; // NLCC Forces from Kronik et al., J. Chem. Phys. 115, 4322 (2001): Eq. 9 From 7c46bd2394c3b0412ef1571e4b40f75de11e73d4 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 6 Jan 2025 14:08:07 -0800 Subject: [PATCH 07/12] Corect NLCC forces sign. --- src/hamiltonian/atomic_potential.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hamiltonian/atomic_potential.hpp b/src/hamiltonian/atomic_potential.hpp index 7f2fd0bc..389f874a 100644 --- a/src/hamiltonian/atomic_potential.hpp +++ b/src/hamiltonian/atomic_potential.hpp @@ -334,7 +334,7 @@ namespace hamiltonian { for(int ispin = 0; ispin < nspin; ispin++) vv += vx[sph.grid_point(ipoint)[0]][sph.grid_point(ipoint)[1]][sph.grid_point(ipoint)[2]][ispin]; auto rr = sph.distance(ipoint); auto grad = spline.derivative(rr)*sph.point_pos(ipoint)/rr; - return -vv*grad; + return vv*grad; }); forces[iatom] = basis.volume_element()*basis.cell().metric().to_cartesian(ff); From 4cbc01183c46a0b1a1cda62f34094a7c8102c957 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Mon, 6 Jan 2025 19:25:30 -0800 Subject: [PATCH 08/12] Increase the integration sphere a bit to account for the derivative. --- src/hamiltonian/atomic_potential.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hamiltonian/atomic_potential.hpp b/src/hamiltonian/atomic_potential.hpp index 389f874a..3dbbfeff 100644 --- a/src/hamiltonian/atomic_potential.hpp +++ b/src/hamiltonian/atomic_potential.hpp @@ -325,7 +325,7 @@ namespace hamiltonian { assert(has_nlcc()); - basis::spherical_grid sphere(basis, atom_position, ps.nlcc_density_radius()); + basis::spherical_grid sphere(basis, atom_position, 1.1*ps.nlcc_density_radius()); auto ff = gpu::run(gpu::reduce(sphere.size()), zero>(), [vx = begin(vxc.hypercubic()), nspin = std::min(2l, vxc.local_set_size()), sph = sphere.ref(), spline = ps.nlcc_density().function()] From 30e0d0127cb3c3aa2d16711f0d8182585a2a9746 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Tue, 7 Jan 2025 14:01:48 -0800 Subject: [PATCH 09/12] Update the tests of forces to account for the nlcc. --- tests/al4h1.sh | 11 ++++++----- tests/na2+.py | 4 ++-- tests/nitrogen.sh | 6 +++--- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/tests/al4h1.sh b/tests/al4h1.sh index 423aed92..b96d9f65 100755 --- a/tests/al4h1.sh +++ b/tests/al4h1.sh @@ -32,11 +32,11 @@ inq util match `inq results ground-state energy nvxc` -4.9007783475 inq util match `inq results ground-state energy exact-exchange` 0.000000000000 3e-5 inq util match `inq results ground-state energy ion` -10.318372113231 3e-5 -inq util match `inq results ground-state forces 0` -0.022483431037 -0.041215997171 -0.052723786483 3e-5 -inq util match `inq results ground-state forces 1` -0.022476660700 0.052697035680 0.041207478998 3e-5 -inq util match `inq results ground-state forces 2` 0.005730135670 -0.012778476335 0.012775275108 3e-5 -inq util match `inq results ground-state forces 3` 0.007076613283 0.012276399154 -0.012280307956 3e-5 -inq util match `inq results ground-state forces 4` 0.027652090218 -0.010193515961 0.010356483661 3e-5 +inq util match `inq results ground-state forces 0` -2.09878683861885427520e-02 -3.86439030449507309184e-02 -4.92550174029726950398e-02 3e-5 +inq util match `inq results ground-state forces 1` -2.09841584974472528113e-02 4.92322743523460754078e-02 3.86387752796131661581e-02 3e-5 +inq util match `inq results ground-state forces 2` 6.67547441607678308795e-03 -1.28302905479014666551e-02 1.28262908762605074703e-02 3e-5 +inq util match `inq results ground-state forces 3` 7.67798179337165348501e-03 1.25067841618366794310e-02 -1.25106042421354955035e-02 3e-5 +inq util match `inq results ground-state forces 4` 2.76562191520605543671e-02 -1.01874661145976881660e-02 1.03629838578964024454e-02 3e-5 inq electrons extra states 0 inq real-time num-steps 30 @@ -50,3 +50,4 @@ inq util match `inq results real-time total-energy 0` -9.8023385899 inq util match `inq results real-time total-energy 10` -9.80233858990918349718e+00 3e-5 inq util match `inq results real-time total-energy 20` -9.80233858990920126075e+00 3e-5 inq util match `inq results real-time total-energy 30` -9.80233858990928652588e+00 3e-5 + diff --git a/tests/na2+.py b/tests/na2+.py index d64062cd..297b475b 100755 --- a/tests/na2+.py +++ b/tests/na2+.py @@ -30,8 +30,8 @@ assert pinq.util.match(pinq.results.ground_state.energy.exact_exchange(), 0.000000000000, 3e-5) assert pinq.util.match(pinq.results.ground_state.energy.ion(), -0.365304170454, 3e-5) -assert pinq.util.match(pinq.results.ground_state.forces()[0], [7.26217835896163880644e-11, -8.66073022870166591849e-11, 1.46550554264721099619e-03], 3e-5) -assert pinq.util.match(pinq.results.ground_state.forces()[1], [8.64681672462082104618e-11, -5.46191861868080879087e-11, -1.46550559849585658283e-03], 3e-5) +assert pinq.util.match(pinq.results.ground_state.forces()[0], [3.74786824263793630366e-11, -1.05853471242772412887e-10, 1.31008348141167883100e-03], 3e-5) +assert pinq.util.match(pinq.results.ground_state.forces()[1], [6.54160820496192356550e-11, -6.38554194489870691095e-11, -1.31008351907754708524e-03], 3e-5) pinq.real_time.num_steps(3000) pinq.real_time.time_step(0.075, "atu") diff --git a/tests/nitrogen.sh b/tests/nitrogen.sh index 0f73cb0b..fff8142f 100755 --- a/tests/nitrogen.sh +++ b/tests/nitrogen.sh @@ -25,10 +25,10 @@ inq util match `inq results ground-state energy nvxc` -6.353978631002 inq util match `inq results ground-state energy exact-exchange` 0.000000000000 1e-5 inq util match `inq results ground-state energy ion` -1.517434464849 1e-5 -inq util match `inq results ground-state forces 0` -0.00000000250910 -0.00000000185763 0.14573092924241 3e-5 -inq util match `inq results ground-state forces 1` 0.00000000341961 0.00000000029088 -0.14573050416125 3e-5 +inq util match `inq results ground-state forces 0` -1.17619241701800515970e-08 3.74561683428810830648e-09 1.38894117877148809415e-01 3e-5 +inq util match `inq results ground-state forces 1` 1.47940666596688138219e-08 -2.21046899539461590873e-09 -1.38893781930214815779e-01 3e-5 #extra check to verify the indiviual component query -inq util match `inq results ground-state forces 1 z` -0.14573050416125 3e-5 +inq util match `inq results ground-state forces 1 z` -1.38893781930214815779e-01 3e-5 inq real-time time-step 0.025 atu inq real-time num-steps 10 From b223d3923c0df7bc33e377a66be0dd1b9130c516 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Tue, 7 Jan 2025 14:53:32 -0800 Subject: [PATCH 10/12] Avoid division by zero when an atom is on top of a grid point. --- src/hamiltonian/atomic_potential.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hamiltonian/atomic_potential.hpp b/src/hamiltonian/atomic_potential.hpp index 3dbbfeff..0cf8c370 100644 --- a/src/hamiltonian/atomic_potential.hpp +++ b/src/hamiltonian/atomic_potential.hpp @@ -333,12 +333,12 @@ namespace hamiltonian { auto vv = 0.0; for(int ispin = 0; ispin < nspin; ispin++) vv += vx[sph.grid_point(ipoint)[0]][sph.grid_point(ipoint)[1]][sph.grid_point(ipoint)[2]][ispin]; auto rr = sph.distance(ipoint); + if(rr < 1e-15) return zero>(); auto grad = spline.derivative(rr)*sph.point_pos(ipoint)/rr; return vv*grad; }); forces[iatom] = basis.volume_element()*basis.cell().metric().to_cartesian(ff); - } if(comm.size() > 1) { From 5d54a56c5361a43bfed0216026347f1810586371 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Tue, 7 Jan 2025 20:19:40 -0800 Subject: [PATCH 11/12] Add an assertion. --- src/hamiltonian/self_consistency.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 708a18f3..e31f7824 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -134,6 +134,8 @@ class self_consistency { // XC double exc, nvxc; hamiltonian.vxc_ = xc_(spin_density, core_density_, exc, nvxc); + + assert(hamiltonian.vxc_.set_size() == vks.set_size()); gpu::run(hamiltonian.vxc_.local_set_size(), hamiltonian.vxc_.basis().local_size(), [vx = begin(hamiltonian.vxc_.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip){ From d7424aa052ecbec6bebc34b7c325ac510b65b1f3 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Tue, 7 Jan 2025 21:58:29 -0800 Subject: [PATCH 12/12] Missing reduction. --- src/hamiltonian/atomic_potential.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/hamiltonian/atomic_potential.hpp b/src/hamiltonian/atomic_potential.hpp index 0cf8c370..f1214c2c 100644 --- a/src/hamiltonian/atomic_potential.hpp +++ b/src/hamiltonian/atomic_potential.hpp @@ -341,6 +341,11 @@ namespace hamiltonian { forces[iatom] = basis.volume_element()*basis.cell().metric().to_cartesian(ff); } + if(basis.comm().size() > 1) { + basis.comm().all_reduce_n(reinterpret_cast(raw_pointer_cast(forces.data_elements())), 3*forces.size()); + } + + //we should use allgather here if(comm.size() > 1) { comm.all_reduce_n(reinterpret_cast(raw_pointer_cast(forces.data_elements())), 3*forces.size()); }