RTP: Fix restart for velocity gauge

cp2k · Apr 5, 2023 · 6a8f4e9 · 6a8f4e9
1 parent ef71e33
commit 6a8f4e9
Show file tree

Hide file tree

Showing 5 changed files with 43 additions and 8 deletions.
diff --git a/src/cp_control_types.F b/src/cp_control_types.F
@@ -194,13 +194,14 @@ MODULE cp_control_types
 ! **************************************************************************************************
    TYPE efield_type
       REAL(KIND=dp)                        :: actual_time
-      REAL(KIND=dp), DIMENSION(:), POINTER :: polarisation
+      REAL(KIND=dp), DIMENSION(:), POINTER :: polarisation => NULL()
       INTEGER                              :: envelop_id
-      REAL(KIND=dp), DIMENSION(:), POINTER :: envelop_r_vars
+      REAL(KIND=dp), DIMENSION(:), POINTER :: envelop_r_vars => NULL()
       INTEGER, DIMENSION(:), POINTER       :: envelop_i_vars
       REAL(KIND=dp)                        :: strength
       REAL(KIND=dp)                        :: phase_offset
       REAL(KIND=dp)                        :: wavelength
+      REAL(KIND=dp), DIMENSION(3)          :: vec_pot_initial = 0.0_dp
    END TYPE efield_type
 
    TYPE efield_p_type
@@ -948,9 +949,8 @@ SUBROUTINE efield_fields_release(efield_fields)
                IF (ASSOCIATED(efield_fields(i)%efield%envelop_i_vars)) THEN
                   DEALLOCATE (efield_fields(i)%efield%envelop_i_vars)
                END IF
-               IF (ASSOCIATED(efield_fields(i)%efield%polarisation)) THEN
+               IF (ASSOCIATED(efield_fields(i)%efield%polarisation)) &
                   DEALLOCATE (efield_fields(i)%efield%polarisation)
-               END IF
                DEALLOCATE (efield_fields(i)%efield)
             END IF
          END DO

diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F
@@ -410,7 +410,8 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
                dft_control%apply_efield_field = .TRUE.
             ELSE
                dft_control%apply_vector_potential = .TRUE.
-               dft_control%rtp_control%vec_pot = 0.0_dp
+               ! Use this input value of vector potential to (re)start RTP
+               dft_control%rtp_control%vec_pot = dft_control%efield_fields(1)%efield%vec_pot_initial
             END IF
          ELSE
             dft_control%apply_efield = .TRUE.
@@ -2167,6 +2168,9 @@ SUBROUTINE read_efield_sections(dft_control, efield_section)
                                    i_val=efield%envelop_id)
          CALL section_vals_val_get(efield_section, "WAVELENGTH", i_rep_section=i, &
                                    r_val=efield%wavelength)
+         CALL section_vals_val_get(efield_section, "VEC_POT_INITIAL", i_rep_section=i, &
+                                   r_vals=tmp_vals)
+         efield%vec_pot_initial = tmp_vals
 
          IF (efield%envelop_id == constant_env) THEN
             ALLOCATE (efield%envelop_i_vars(2))

diff --git a/src/input_cp2k_field.F b/src/input_cp2k_field.F
@@ -162,6 +162,19 @@ SUBROUTINE create_efield_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="VEC_POT_INITIAL", &
+                          description="Initial value of the vector "// &
+                          "potential (for velocity gauge). This input is "// &
+                          "made especially for restarting RTP calculation. "// &
+                          "Unit is atomic unit. "// &
+                          "Note that if several field sections are defined, only the first one will be used.", &
+                          usage="vec_pot_initial  1.0E-2 0.0 0.0", &
+                          repeats=.FALSE., &
+                          n_var=3, type_of_var=real_t, &
+                          default_r_vals=(/0.0_dp, 0.0_dp, 0.0_dp/))
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL create_constant_env_section(subsection)
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)

diff --git a/src/motion/rt_propagation.F b/src/motion/rt_propagation.F
@@ -375,7 +375,7 @@ SUBROUTINE run_propagation(qs_env, force_env, globenv)
             used_time = time_iter_stop - time_iter_start
             time_iter_start = time_iter_stop
             CALL rt_prop_output(qs_env, real_time_propagation, delta_iter=rtp%delta_iter, used_time=used_time)
-            CALL rt_write_input_restart(force_env=force_env)
+            CALL rt_write_input_restart(force_env=force_env, qs_env=qs_env)
             IF (should_stop) EXIT
          ELSE
             EXIT
@@ -395,16 +395,22 @@ END SUBROUTINE run_propagation
 !> \brief overwrites some values in the input file such that the .restart
 !>        file will contain the appropriate information
 !> \param md_env ...
+!> \param qs_env ...
 !> \param force_env ...
 !> \author Florian Schiffmann (02.09)
 ! **************************************************************************************************
 
-   SUBROUTINE rt_write_input_restart(md_env, force_env)
+   SUBROUTINE rt_write_input_restart(md_env, qs_env, force_env)
       TYPE(md_environment_type), OPTIONAL, POINTER       :: md_env
+      TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
       TYPE(force_env_type), POINTER                      :: force_env
 
-      TYPE(section_vals_type), POINTER                   :: motion_section, root_section, rt_section
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_vals
+      TYPE(dft_control_type), POINTER                    :: dft_control
+      TYPE(section_vals_type), POINTER                   :: efield_section, motion_section, &
+                                                            root_section, rt_section
 
+      CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
       root_section => force_env%root_section
       motion_section => section_vals_get_subs_vals(root_section, "MOTION")
       rt_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION")
@@ -414,6 +420,16 @@ SUBROUTINE rt_write_input_restart(md_env, force_env)
          CALL section_vals_val_set(motion_section, "MD%STEP_START_VAL", i_val=force_env%qs_env%sim_step)
       END IF
 
+      IF (dft_control%apply_vector_potential) THEN
+         efield_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%DFT%EFIELD")
+         NULLIFY (tmp_vals)
+         ALLOCATE (tmp_vals(3))
+         tmp_vals = dft_control%efield_fields(1)%efield%vec_pot_initial
+         CALL section_vals_val_set(efield_section, "VEC_POT_INITIAL", &
+                                   r_vals_ptr=tmp_vals, &
+                                   i_rep_section=1)
+      END IF
+
       CALL write_restart(md_env=md_env, root_section=root_section)
 
    END SUBROUTINE rt_write_input_restart

diff --git a/src/rt_propagation_velocity_gauge.F b/src/rt_propagation_velocity_gauge.F
@@ -267,6 +267,8 @@ SUBROUTINE update_vector_potential(qs_env, dft_control)
 
       CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
       dft_control%rtp_control%vec_pot = dft_control%rtp_control%vec_pot - field*qs_env%rtp%dt
+      ! Update the vec_pot_initial value for RTP restart:
+      dft_control%efield_fields(1)%efield%vec_pot_initial = dft_control%rtp_control%vec_pot
 
    END SUBROUTINE update_vector_potential