Merge branch 'stress_tensor' into 'master'

Implementation of  the stress sensor for the electron gas (no ionic terms yet)

See merge request npneq/inq!1142

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<Type const>(1);
+			return cubic().template reinterpret_array_cast<Type>(1);
 		}
 
 		auto hypercubic() {
@@ -193,7 +193,7 @@ namespace basis {
 
 
 field<basis::real_space, complex> complex_field(field<basis::real_space, double> const & rfield) {
-	field<basis::real_space, complex> cfield(rfield.skeleton());		
+	field<basis::real_space, complex> 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<basis::real_space, vector3<inq::complex, VectorSpace>> complex_field(field
 }
 
 field<basis::real_space, double> real_field(field<basis::real_space, complex> const & cfield) {
-	field<basis::real_space, double> rfield(cfield.skeleton());		
+	field<basis::real_space, double> rfield(cfield.skeleton());     
 	rfield.linear() = boost::multi::blas::real(cfield.linear());
 	return rfield;
 }
 
 template <class VectorSpace>
 field<basis::real_space, vector3<double, VectorSpace>> real_field(field<basis::real_space, vector3<complex, VectorSpace>> const & cfield) {
-	field<basis::real_space, vector3<double, VectorSpace>> rfield(cfield.skeleton());		
+	field<basis::real_space, vector3<double, VectorSpace>> 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<basis::real_space, vector3<double, VectorSpace>> real_field(field<basis::r
 TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){
 
 	using namespace inq;
-	using namespace inq::magnitude;	
+	using namespace inq::magnitude; 
 	using namespace Catch::literals;
 
 	parallel::communicator comm{boost::mpi3::environment::get_world_instance()};
@@ -276,7 +276,7 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){
 
 	ff.fill(12.2244);
 
-	for(int ii = 0; ii < rs.part().local_size(); ii++) CHECK(ff.linear()[ii] == 12.2244_a);	
+	for(int ii = 0; ii < rs.part().local_size(); ii++) CHECK(ff.linear()[ii] == 12.2244_a); 
 
 	basis::field<basis::real_space, double> 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());

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_}.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) {

src/ground_state/results.hpp

@@ -25,6 +25,7 @@ struct results {
   vector3<double> magnetization;
   energy_type energy;
   forces_type forces;
+	vector3<vector3<double>> 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 + "'.";
@@ -35,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) {
@@ -47,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);
 
@@ -64,8 +67,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";
@@ -105,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<double>{1.28, 7.62, 5.56};
+	res.stress[1] = vector3<double>{7.62, 3.03, 0.57};
+	res.stress[2] = vector3<double>{5.56, 0.57, 8.88};;
   res.forces = ground_state::results::forces_type{65, vector3<double>{3.55, 4.55, 5.55}};
 
   CHECK(res.total_iter == 333);
@@ -119,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<double>{1.28, 7.62, 5.56});
+	CHECK(res.stress[1] == vector3<double>{7.62, 3.03, 0.57});
+	CHECK(res.stress[2] == vector3<double>{5.56, 0.57, 8.88});
   CHECK(res.forces == ground_state::results::forces_type{65, vector3<double>{3.55, 4.55, 5.55}});
 
   res.save(comm, "save_results");
@@ -136,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<double>{1.28, 7.62, 5.56});
+	CHECK(read_res.stress[1] == vector3<double>{7.62, 3.03, 0.57});
+	CHECK(read_res.stress[2] == vector3<double>{5.56, 0.57, 8.88});
   CHECK(read_res.forces == ground_state::results::forces_type{65, vector3<double>{3.55, 4.55, 5.55}});
 
 }

src/observables/forces_stress.hpp

@@ -21,40 +21,118 @@ namespace observables {
 
 struct forces_stress {
 	gpu::array<vector3<double>, 1> forces;
-	gpu::array<double, 2>          stress;
+	vector3<vector3<double>> stress;
 
 	forces_stress() = default;
 
-	template <typename HamiltonianType>
-	forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham):
-		forces(ions.size()),
-		stress({3, 3}, 0.0)
+	template <typename HamiltonianType, typename Energy>
+	forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham, Energy const & energy):
+		forces(ions.size())
 	{
-		calculate(ions, electrons, ham);
+		calculate(ions, electrons, ham, energy);
 	}
 
 #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 <typename Stress1D>
+	auto tensor(Stress1D const & stress1d) {
+		vector3<vector3<double>> 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 <typename HamiltonianType>
-	void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){
+	template <typename GPhi, typename Occupations>
+	vector3<vector3<double>> 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()), 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);
+															 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);
+
+		return -gphi.basis().volume_element()*tensor(stress1d);
+	}
+
+	template <typename Density>
+	vector3<vector3<double>> 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()), 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);
+															 auto ef_cart = metric.to_cartesian(ef[ip]);
+															 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);
+	}
 
+	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
+		//   - Stress formulas: Eq. (33) of https://arxiv.org/pdf/1809.08157
+
+
 		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<basis::real_space, vector3<double, covariant>> gdensity(electrons.density_basis());
 		gdensity.fill(vector3<double, covariant>{0.0, 0.0, 0.0});
 
 		gpu::array<vector3<double>, 1> forces_non_local(ions.size(), {0.0, 0.0, 0.0});
-		
+
 		auto iphi = 0;
 		for(auto & phi : electrons.kpin()){
 
 			auto gphi = operations::gradient(phi, /* factor = */ 1.0, /*shift = */ phi.kpoint());
 			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);
-
+
+			stress += stress_kinetic(gphi, electrons.occupations()[iphi]);
+
 			iphi++;
 		}
 
@@ -95,7 +173,24 @@ 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
+
+		stress += stress_electrostatic(electrons.density());
+
+		// 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();
+			}
+		}
+
 	}
 
 };

src/operations/gradient.hpp

@@ -40,7 +40,7 @@ ResultType gradient(FieldSetType const & ff, double factor = 1.0, vector3<double
 
 		ResultType gradff(ff.skeleton());
 
-		gpu::run(gradff.set_part().local_size(), gradff.basis().local_sizes()[2], gradff.basis().local_sizes()[1], gradff.basis().local_sizes()[0],
+		gpu::run(gradff.local_set_size(), gradff.basis().local_sizes()[2], gradff.basis().local_sizes()[1], gradff.basis().local_sizes()[0],
 						 [point_op = ff.basis().point_op(), gradffcub = begin(gradff.hypercubic()), ffcub = begin(ff.hypercubic()), factor, shift]
 						 GPU_LAMBDA (auto ist, auto iz, auto iy, auto ix){
 							 auto grad = factor*complex(0.0, 1.0)*(point_op.gvector(ix, iy, iz) + shift);

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<double, covariant>{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);

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