comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 13:45:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

E-field restart and output (#3354)

Browse source
This commit is contained in:
marcella 2024-04-16 12:53:47 +02:00 • committed by GitHub
parent e5fdd8129e
commit 1a29073fbc
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
6 changed files with 76 additions and 33 deletions
Download patch file Download diff file Expand all files Collapse all files

15
src/colvar_methods.F
View file

@ -3899,7 +3899,8 @@ CONTAINS
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 @@ CONTAINS
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

1
src/cp_control_types.F
View file

@ -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.

2
src/efield_utils.F
View file

@ -234,10 +234,12 @@ CONTAINS
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

57
src/emd/rt_propagation_output.F
View file

@ -826,49 +826,56 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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

33
src/input_restart_force_eval.F
View file

@ -89,19 +89,22 @@ CONTAINS
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 @@ CONTAINS
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

1
src/rt_propagation_velocity_gauge.F
View file

@ -267,6 +267,7 @@ CONTAINS
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