Merge branch 'nlcc_forces' into 'master'

Add the non-linear core correction force term

See merge request npneq/inq!1186

src/hamiltonian/atomic_potential.hpp

@@ -300,6 +300,60 @@ namespace hamiltonian {
 			density.all_reduce(comm);
 			return density;			
 		}
+
+		////////////////////////////////////////////////////////////////////////////////////
+
+		template <class Comm, class Ions>
+		gpu::array<vector3<double>, 1> nlcc_forces(Comm & comm, basis::field_set<basis::real_space, double> 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<vector3<double>, 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, 1.1*ps.nlcc_density_radius());
+
+				auto ff = gpu::run(gpu::reduce(sphere.size()), zero<vector3<double, contravariant>>(),
+													 [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);
+														 if(rr < 1e-15) return zero<vector3<double, contravariant>>();
+														 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(basis.comm().size() > 1) {
+				basis.comm().all_reduce_n(reinterpret_cast<double *>(raw_pointer_cast(forces.data_elements())), 3*forces.size());
+			}
+
+			//we should use allgather here
+			if(comm.size() > 1) {
+				comm.all_reduce_n(reinterpret_cast<double *>(raw_pointer_cast(forces.data_elements())), 3*forces.size());
+			}
+
+			return forces;
+		}
+
+		////////////////////////////////////////////////////////////////////////////////////
 
 		template <class output_stream>
 		void info(output_stream & out) const {

src/hamiltonian/ks_hamiltonian.hpp

@@ -42,6 +42,7 @@ class ks_hamiltonian {
 private:
 
 	exchange_operator exchange_;
+	basis::field_set<basis::real_space, double> vxc_;
 	basis::field_set<basis::real_space, PotentialType> scalar_potential_;
 	vector3<double, covariant> 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)
@@ -228,6 +230,12 @@ class ks_hamiltonian {
 		return uniform_vector_potential_;
 	}
 
+	////////////////////////////////////////////////////////////////////////////////////////////
+
+	auto & vxc() const {
+		return vxc_;
+	}
+
 	////////////////////////////////////////////////////////////////////////////////////////////
 
 	template <typename Perturbation>

src/hamiltonian/self_consistency.hpp

@@ -133,7 +133,15 @@ class self_consistency {
 
 		// XC
 		double exc, nvxc;
-		xc_(spin_density, core_density_, vks, 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){
+							 vk[ip][is] += vx[ip][is];
+						 });
+
 		energy.xc(exc);
 		energy.nvxc(nvxc);
 

src/hamiltonian/xc_term.hpp

@@ -87,17 +87,6 @@ class xc_term {
 		return full_density;
 	}
 
-	////////////////////////////////////////////////////////////////////////////////////////////
-
-	template <typename VXC, typename VKS>
-	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 <typename SpinDensityType, typename VXC>
@@ -122,20 +111,19 @@ class xc_term {
 
   ////////////////////////////////////////////////////////////////////////////////////////////
 
-  template <typename SpinDensityType, typename CoreDensityType, typename VKSType>
-  void operator()(SpinDensityType const & spin_density, CoreDensityType const & core_density, VKSType & vks, double & exc, double & nvxc) const {
-
-      exc = 0.0;
+  template <typename SpinDensityType, typename CoreDensityType>
+  auto operator()(SpinDensityType const & spin_density, CoreDensityType const & core_density, double & exc, double & nvxc) const {
+
+		basis::field_set<basis::real_space, double> 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<basis::real_space, double> vxc(spin_density.skeleton());
-		vxc.fill(0.0);
-
 		basis::field_set<basis::real_space, double> vfunc(full_density.skeleton());
 		auto density_gradient = std::optional<decltype(operations::gradient(full_density))>{};
 		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;
 	}
 
 	////////////////////////////////////////////////////////////////////////////////////////////

src/observables/forces_stress.hpp

@@ -117,7 +117,7 @@ struct forces_stress {
 	template <typename HamiltonianType, typename Energy>
 	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
 
 
@@ -180,13 +180,12 @@ struct forces_stress {
 				electrons.density_basis().comm().all_reduce_n(reinterpret_cast<double *>(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++) {

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
+

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")

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