E-field restart and output (cp2k#3354)

src/colvar_methods.F

@@ -3899,7 +3899,8 @@ SUBROUTINE acid_hyd_dist_colvar(colvar, cell, subsys, particles)
       fcut = num_cut*invden_cut
 
       !Final distance acid - hydronium
-      fbrace = rion_num/rion_den/2.0_dp
+!      fbrace = rion_num/rion_den/2.0_dp
+      fbrace = rion_num/rion_den/n_oxygens_acid
       rion = fcut*fbrace
       colvar%ss = rion
 
@@ -3922,14 +3923,18 @@ SUBROUTINE acid_hyd_dist_colvar(colvar, cell, subsys, particles)
       DO kk = 1, n_oxygens_acid
          DO ii = 1, n_oxygens_water
             colvar%dsdr(1:3, offsetO + kk) = colvar%dsdr(1:3, offsetO + kk) &
-                                             + fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/2.0_dp
+                                             + fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
+!                                             + fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/2.0_dp
             colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
-                                   - fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/2.0_dp
+                                   - fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
+!                                   - fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/2.0_dp
             DO jj = 1, n_hydrogens
                colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
-                                      + fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/2.0_dp
+                                      + fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
+!                                      + fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/2.0_dp
                colvar%dsdr(1:3, offsetH + jj) = colvar%dsdr(1:3, offsetH + jj) &
-                                                - fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/2.0_dp
+                                                - fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
+!                                                - fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/2.0_dp
             END DO
          END DO
       END DO

src/cp_control_types.F

@@ -87,6 +87,7 @@ MODULE cp_control_types
       INTEGER, DIMENSION(3)                :: delta_pulse_direction = 0
       REAL(KIND=dp)                        :: delta_pulse_scale = 0.0_dp
       LOGICAL                              :: velocity_gauge = .FALSE.
+      REAL(KIND=dp), DIMENSION(3)          :: field = 0.0_dp
       REAL(KIND=dp), DIMENSION(3)          :: vec_pot = 0.0_dp
       LOGICAL                              :: nl_gauge_transform = .FALSE.
       LOGICAL                              :: is_proj_mo = .FALSE.

src/efield_utils.F

@@ -234,10 +234,12 @@ SUBROUTINE calculate_ecore_efield(qs_env, calculate_forces)
                   CALL get_qs_env(qs_env=qs_env, force=force)
                   force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) - field*zeff
                END IF
+!               END IF
             END DO
 
          END DO
          IF (dft_control%apply_efield_field) energy%efield_core = efield_ener
+!          energy%efield_core = efield_ener
       END IF
       CALL timestop(handle)
    END SUBROUTINE calculate_ecore_efield

src/emd/rt_propagation_output.F

@@ -826,49 +826,56 @@ SUBROUTINE print_field_applied(qs_env, dft_section)
 
       CHARACTER(LEN=3), DIMENSION(3)                     :: rlab
       CHARACTER(LEN=default_path_length)                 :: filename
-      INTEGER                                            :: i, output_unit, unit_nr
-      REAL(kind=dp)                                      :: field(3), to_write(3)
+      INTEGER                                            :: i, i_step, output_unit, unit_nr
+      LOGICAL                                            :: new_file
+      REAL(kind=dp)                                      :: field(3)
       TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(dft_control_type), POINTER                    :: dft_control
+      TYPE(rt_prop_type), POINTER                        :: rtp
 
       NULLIFY (dft_control)
 
       logger => cp_get_default_logger()
       output_unit = cp_logger_get_default_io_unit(logger)
 
-      CALL get_qs_env(qs_env, dft_control=dft_control)
+      CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp)
+
+      i_step = rtp%istep
 
       unit_nr = cp_print_key_unit_nr(logger, dft_section, &
-                                     "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat")
+                                     "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat", is_new_file=new_file)
 
       IF (output_unit > 0) THEN
+         rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
          IF (unit_nr /= output_unit) THEN
             INQUIRE (UNIT=unit_nr, NAME=filename)
             WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
                "FIELD", "The field applied is written to the file:", &
                TRIM(filename)
-            WRITE (UNIT=unit_nr, FMT="(/,(T2,A,T40,I6))") &
-               "Real time propagation step:", qs_env%sim_step
          ELSE
-            WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED"
+            WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED [a.u.]"
+            WRITE (UNIT=output_unit, FMT="(T5,3(A,A,E16.8,1X))") &
+               (TRIM(rlab(i)), "=", dft_control%rtp_control%field(i), i=1, 3)
          END IF
 
-         rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
+         IF (new_file) THEN
+            IF (dft_control%apply_efield_field) THEN
+               WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,3(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z"
+            ELSE IF (dft_control%apply_vector_potential) THEN
+               WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,6(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z", &
+                  "Vec. Pot. X", "Vec. Pot. Y", "Vec. Pot. Z"
+            END IF
+         END IF
 
          IF (dft_control%apply_efield_field) THEN
-            WRITE (unit_nr, "(T3,A)") "Electric Field (LG) in atomic units:"
             CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
-            to_write = field
+            WRITE (UNIT=unit_nr, FMT="(I10,F20.6,3(E16.8,1X))") i_step, i_step*rtp%dt*femtoseconds, field(1), field(2), field(3)
          ELSE IF (dft_control%apply_vector_potential) THEN
-            WRITE (unit_nr, "(T3,A)") "Vector potential (VG) in atomic units:"
-            to_write = dft_control%rtp_control%vec_pot
-         ELSE
-            WRITE (unit_nr, "(T3,A)") "No electric field applied"
-            to_write = 0._dp
+            WRITE (UNIT=unit_nr, FMT="(I10,F20.6,6(E16.8,1X))") i_step, i_step*rtp%dt*femtoseconds, &
+               dft_control%rtp_control%field(1), dft_control%rtp_control%field(2), dft_control%rtp_control%field(3), &
+               dft_control%rtp_control%vec_pot(1), dft_control%rtp_control%vec_pot(2), dft_control%rtp_control%vec_pot(3)
          END IF
 
-         WRITE (unit_nr, "(T5,3(A,A,E16.8,1X))") &
-            (TRIM(rlab(i)), "=", to_write(i), i=1, 3)
       END IF
 
       CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
@@ -891,13 +898,16 @@ SUBROUTINE print_rtp_energy_components(qs_env, dft_section)
       INTEGER                                            :: i_step, output_unit, unit_nr
       LOGICAL                                            :: new_file
       TYPE(cp_logger_type), POINTER                      :: logger
+      TYPE(dft_control_type), POINTER                    :: dft_control
       TYPE(qs_energy_type), POINTER                      :: energy
       TYPE(rt_prop_type), POINTER                        :: rtp
 
+      NULLIFY (dft_control, energy, rtp)
+
       logger => cp_get_default_logger()
       output_unit = cp_logger_get_default_io_unit(logger)
 
-      CALL get_qs_env(qs_env, rtp=rtp, energy=energy)
+      CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
       i_step = rtp%istep
 
       unit_nr = cp_print_key_unit_nr(logger, dft_section, &
@@ -917,11 +927,14 @@ SUBROUTINE print_rtp_energy_components(qs_env, dft_section)
          IF (new_file) THEN
             ! NOTE that these are not all terms contributing to the total energy for RTP, only a selection of those
             ! most significant / impactful. Therefore the printed components likely will not add up to the total energy.
-            WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,6X,A,6X,A,6X,A,6X,A,6X,A)') "Step Nr.", "Time[fs]", &
-               "Total ener.[a.u.]", "core.[a.u.]", "core_self.[a.u.]", "hartree.[a.u.]", "exc.[a.u.]"
+            WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,10(6X,A))') "Step Nr.", "Time[fs]", &
+               "Total ener.[a.u.]", "core[a.u.]   ", " overlap [a.u.]", "hartree[a.u.]", " exc. [a.u.] ", &
+               " hartree 1c[a.u.]", "exc 1c[a.u.] ", "exc admm[a.u.]", "exc 1c admm[a.u.]", "efield LG"
+
          END IF
-         WRITE (UNIT=unit_nr, FMT="(I10,F20.6,F20.9,F20.9,F20.9,F20.9,F20.9)") &
-            i_step, i_step*rtp%dt*femtoseconds, energy%total, energy%core, energy%core_self, energy%hartree, energy%exc
+         WRITE (UNIT=unit_nr, FMT="(I10,F20.6,10(F20.9))") &
+            i_step, i_step*rtp%dt*femtoseconds, energy%total, energy%core, energy%core_overlap, energy%hartree, energy%exc, &
+            energy%hartree_1c, energy%exc1, energy%exc_aux_fit, energy%exc1_aux_fit, energy%efield_core
 
       END IF
 

src/input_restart_force_eval.F

@@ -89,19 +89,22 @@ SUBROUTINE update_force_eval(force_env, root_section, &
       LOGICAL, INTENT(IN)                                :: write_binary_restart_file
       LOGICAL, OPTIONAL                                  :: respa
 
-      INTEGER                                            :: i_subsys, iforce_eval, myid, nforce_eval
+      INTEGER                                            :: i, i_subsys, iforce_eval, myid, &
+                                                            nforce_eval
       INTEGER, DIMENSION(:), POINTER                     :: i_force_eval
-      LOGICAL                                            :: multiple_subsys, skip_vel_section
+      LOGICAL                                            :: is_present, multiple_subsys, &
+                                                            skip_vel_section
       REAL(KIND=dp), DIMENSION(:), POINTER               :: work
       TYPE(cell_type), POINTER                           :: cell
       TYPE(cp_subsys_type), POINTER                      :: subsys
       TYPE(dft_control_type), POINTER                    :: dft_control
       TYPE(section_vals_type), POINTER                   :: cell_section, dft_section, &
                                                             force_env_sections, qmmm_section, &
-                                                            rng_section, subsys_section
+                                                            rng_section, subsys_section, &
+                                                            tmp_section
       TYPE(virial_type), POINTER                         :: virial
 
-      NULLIFY (rng_section, subsys_section, cell_section, virial, subsys, cell, dft_control, work)
+      NULLIFY (rng_section, subsys_section, cell_section, virial, subsys, cell, dft_control, work, tmp_section)
       ! If it's not a dynamical run don't update the velocity section
       CALL section_vals_val_get(root_section, "GLOBAL%RUN_TYPE", i_val=myid)
       skip_vel_section = ((myid /= mol_dyn_run) .AND. &
@@ -187,8 +190,26 @@ SUBROUTINE update_force_eval(force_env, root_section, &
          END IF
       END IF
 
-      IF (myid == ehrenfest) CALL section_vals_val_set(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION%INITIAL_WFN", &
-                                                       i_val=use_rt_restart)
+      IF (ASSOCIATED(force_env%qs_env)) THEN
+         CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
+         dft_section => section_vals_get_subs_vals(force_env_sections, "DFT", &
+                                                   i_rep_section=i_force_eval(iforce_eval))
+         tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION")
+         CALL section_vals_get(tmp_section, explicit=is_present)
+         IF (is_present) THEN
+            CALL section_vals_val_set(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION%INITIAL_WFN", &
+                                      i_val=use_rt_restart)
+            IF (ASSOCIATED(dft_control%efield_fields)) THEN
+               tmp_section => section_vals_get_subs_vals(dft_section, "EFIELD")
+               ALLOCATE (work(3))
+               DO i = 1, SIZE(dft_control%efield_fields)
+                  work = dft_control%efield_fields(1)%efield%vec_pot_initial
+                  CALL section_vals_val_set(tmp_section, "VEC_POT_INITIAL", i_rep_section=i, &
+                                            r_vals_ptr=work)
+               END DO
+            END IF
+         END IF
+      END IF
       DEALLOCATE (i_force_eval)
 
    END SUBROUTINE update_force_eval

src/rt_propagation_velocity_gauge.F

@@ -267,6 +267,7 @@ SUBROUTINE update_vector_potential(qs_env, dft_control)
       REAL(kind=dp)                                      :: field(3)
 
       CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
+      dft_control%rtp_control%field = field
       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