Fixed a test and cleaned up the LDOS

mala-project · Apr 11, 2024 · d671a29 · d671a29
1 parent bb2acd3
commit d671a29
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 6 deletions.
diff --git a/examples/ex19_atomic_forces.py b/examples/ex19_atomic_forces.py
@@ -136,22 +136,42 @@ def hartree_contribution():
     density = ldos_calculator.density
     # mala_forces = density_calculator.force_contributions
     ldos_calculator.debug_forces_flag = "hartree"
+
     mala_forces = ldos_calculator.atomic_forces.copy()
+    fermi_energy = ldos_calculator.fermi_energy.copy()
 
     delta_numerical = 1e-6
     points = [0, 2000, 4000]
 
     for point in points:
         numerical_forces = []
         for i in range(0, parameters.targets.ldos_gridsize):
+            # Chemical potential constant
+            # ldos_plus = ldos.copy()
+            # ldos_plus[point, i] += delta_numerical * 0.5
+            # ldos_calculator.read_from_array(ldos_plus)
+            # _, derivative_plus = ldos_calculator.get_total_energy(ldos_plus,
+            #                                                       fermi_energy=fermi_energy,
+            #                                                       return_energy_contributions=True)
+            # derivative_plus = derivative_plus["e_hartree"]
+            #
+            # ldos_minus = ldos.copy()
+            # ldos_minus[point, i] -= delta_numerical * 0.5
+            # ldos_calculator.read_from_array(ldos_minus)
+            # _, derivative_minus = ldos_calculator.get_total_energy(ldos_minus,
+            #                                                       fermi_energy=fermi_energy,
+            #                                                       return_energy_contributions=True)
+            # derivative_minus = derivative_minus["e_hartree"]
+
+            # Chemical potential NOT constant
             ldos_plus = ldos.copy()
-            ldos_plus[0, i] += delta_numerical * 0.5
+            ldos_plus[point, i] += delta_numerical * 0.5
             ldos_calculator.read_from_array(ldos_plus)
             _, derivative_plus = ldos_calculator.get_total_energy(return_energy_contributions=True)
             derivative_plus = derivative_plus["e_hartree"]
 
             ldos_minus = ldos.copy()
-            ldos_minus[0, i] -= delta_numerical * 0.5
+            ldos_minus[point, i] -= delta_numerical * 0.5
             ldos_calculator.read_from_array(ldos_minus)
             _, derivative_minus = ldos_calculator.get_total_energy(return_energy_contributions=True)
             derivative_minus = derivative_minus["e_hartree"]

diff --git a/mala/targets/ldos.py b/mala/targets/ldos.py
@@ -1350,6 +1350,7 @@ def get_atomic_forces(self, ldos_data=None, dos_data=None,
             if self.debug_forces_flag is None or self.debug_forces_flag == \
                 "hartree":
                 print(self._density_calculator.force_contributions.shape)
+                d_E_h_d_n = self._density_calculator.force_contributions
                 d_n_d_d = analytical_integration_weights("F0", "F1", fermi_energy,
                                                          energy_grid, temperature)
 
@@ -1359,14 +1360,14 @@ def get_atomic_forces(self, ldos_data=None, dos_data=None,
                                                                  energy_grid, temperature)) / \
                            (2.0*delta)
                 d_n_d_mu = np.dot(ldos_data, d_f01_d_mu)
-                d_E_h_d_mu = np.dot(self._density_calculator.force_contributions[:,0],
-                                    d_n_d_mu)
+                d_E_h_d_mu = np.dot(d_E_h_d_n[:,0], d_n_d_mu)
                 d_mu_d_d = np.zeros_like(ldos_data)
                 d_mu_d_d[:] = -(-1.0) * (d_n_d_d*self.voxel.volume)
                 d_mu_d_d /= (self._density_of_states_calculator.
                              d_number_of_electrons_d_mu)
-                d_E_d_n = self._density_calculator.force_contributions * \
-                          d_n_d_d + d_E_h_d_mu * d_mu_d_d
+                d_E_d_n = d_E_h_d_n * d_n_d_d + d_E_h_d_mu * d_mu_d_d
+
+
                 d_E_d_d += d_E_d_n
 
             if self.input_data_derivative is not None: