From 1a29073fbcb476dc3ab79655744d12eef97a2b59 Mon Sep 17 00:00:00 2001 From: marcella <40494362+marci73@users.noreply.github.com> Date: Tue, 16 Apr 2024 12:53:47 +0200 Subject: [PATCH] E-field restart and output (#3354) --- src/colvar_methods.F | 15 +++++--- src/cp_control_types.F | 1 + src/efield_utils.F | 2 + src/emd/rt_propagation_output.F | 57 ++++++++++++++++++----------- src/input_restart_force_eval.F | 33 ++++++++++++++--- src/rt_propagation_velocity_gauge.F | 1 + 6 files changed, 76 insertions(+), 33 deletions(-) diff --git a/src/colvar_methods.F b/src/colvar_methods.F index cf196f1f0d..585d642953 100644 --- a/src/colvar_methods.F +++ b/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 diff --git a/src/cp_control_types.F b/src/cp_control_types.F index 6a23fa3fc3..7953a6ca28 100644 --- a/src/cp_control_types.F +++ b/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. diff --git a/src/efield_utils.F b/src/efield_utils.F index 5959812ebf..0ea2e8d2f4 100644 --- a/src/efield_utils.F +++ b/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 diff --git a/src/emd/rt_propagation_output.F b/src/emd/rt_propagation_output.F index 845e824b6f..09dfc03ee3 100644 --- a/src/emd/rt_propagation_output.F +++ b/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 diff --git a/src/input_restart_force_eval.F b/src/input_restart_force_eval.F index 33380944e9..a9c7bafe6a 100644 --- a/src/input_restart_force_eval.F +++ b/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 diff --git a/src/rt_propagation_velocity_gauge.F b/src/rt_propagation_velocity_gauge.F index a0910d8e22..db1e58a93e 100644 --- a/src/rt_propagation_velocity_gauge.F +++ b/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