fixing the core H for EMD

cp2k · Oct 10, 2023 · 6a02f61 · 6a02f61
1 parent 188e445
commit 6a02f61
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 5 deletions.
diff --git a/src/qs_force.F b/src/qs_force.F
@@ -28,7 +28,8 @@ MODULE qs_force
                                               dbcsr_set
    USE dft_plus_u,                      ONLY: plus_u
    USE ec_env_types,                    ONLY: energy_correction_type
-   USE efield_utils,                    ONLY: calculate_ecore_efield
+   USE efield_utils,                    ONLY: calculate_ecore_efield,&
+                                              efield_potential_lengh_gauge
    USE energy_corrections,              ONLY: energy_correction
    USE excited_states,                  ONLY: excited_state_energy
    USE hfx_exx,                         ONLY: calculate_exx
@@ -290,6 +291,11 @@ SUBROUTINE qs_forces(qs_env)
       ELSE
          ! Dispersion energy and forces are calculated in qs_energy?
          CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.)
+         ! The above line reset the core H, which should be re-updated in case a TD field is applied:
+         IF (qs_env%run_rtp) THEN
+            IF (dft_control%apply_efield_field) &
+               CALL efield_potential_lengh_gauge(qs_env)
+         END IF
          CALL calculate_ecore_self(qs_env)
          CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.TRUE.)
          CALL calculate_ecore_efield(qs_env, calculate_forces=.TRUE.)
@@ -587,6 +593,7 @@ SUBROUTINE write_forces(qs_force, atomic_kind_set, ftype, output_unit, &
                   iatom, ikind, "    dispersion", qs_force(ikind)%dispersion(1:3, i), &
                   iatom, ikind, "           gCP", qs_force(ikind)%gcp(1:3, i), &
                   iatom, ikind, "       fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
+                  iatom, ikind, "     ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
                   iatom, ikind, "        efield", qs_force(ikind)%efield(1:3, i), &
                   iatom, ikind, "           eev", qs_force(ikind)%eev(1:3, i), &
                   iatom, ikind, "   mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &

diff --git a/tests/QS/regtest-rtp-2/TEST_FILES b/tests/QS/regtest-rtp-2/TEST_FILES
@@ -10,9 +10,9 @@ H2-rtp_restart-1.inp                                   1      3e-13
 H2-rtp-efield.inp                                      1      2e-10              -0.82613324527108
 H2-rtp-efield-vg.inp                                   1      2e-10              -0.83550227784861
 H2-rtp-efield-vg-restart.inp                           1      2e-10              -0.79589429986106
-H2-emd-efield.inp                                      2      4e-11                -0.902242711172
-H2-emd-efield-ramp.inp                                 2      5e-12                -0.902242710949
-H2-emd-efield-custom.inp                               2      1e-14                -0.902242709393
+H2-emd-efield.inp                                      2      4e-11                -0.894818188949
+H2-emd-efield-ramp.inp                                 2      5e-12                -0.885822019185
+H2-emd-efield-custom.inp                               2      1e-14                -0.902238043295
 H2-rtp_ETRS_ARNOLDI.inp                                1      3e-13              -0.90223968349550
-H2-emd_ETRS_ARNOLDI.inp                                1      1e-10             -17.08557183219799
+H2-emd_ETRS_ARNOLDI.inp                                1      1e-10             -17.08556436304139
 #EOF