From 60cf13d4c5115fdbc5fba5faab15483a475585c2 Mon Sep 17 00:00:00 2001 From: Leo Decking Date: Wed, 26 Feb 2025 09:44:50 +0100 Subject: [PATCH] Add INTENSITY_LIST, START_FRAME and END_FRAME keywords for periodic efield Co-authored-by: Leo Decking --- src/cp_control_types.F | 3 + src/cp_control_utils.F | 64 ++++++++++++++++-- src/eeq_method.F | 16 +++++ src/efield_tb_methods.F | 36 +++++++--- src/input_cp2k_field.F | 35 +++++++++- src/qs_efield_berry.F | 41 ++++++++---- tests/QS/regtest-p-efield/H2O-field-list.inp | 70 ++++++++++++++++++++ tests/QS/regtest-p-efield/TEST_FILES | 1 + 8 files changed, 240 insertions(+), 26 deletions(-) create mode 100644 tests/QS/regtest-p-efield/H2O-field-list.inp diff --git a/src/cp_control_types.F b/src/cp_control_types.F index b12bbaf160..837de0879f 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -270,6 +270,9 @@ MODULE cp_control_types REAL(KIND=dp), DIMENSION(3) :: polarisation = 0.0_dp REAL(KIND=dp), DIMENSION(3) :: d_filter = 0.0_dp REAL(KIND=dp) :: strength = 0.0_dp + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: strength_list + INTEGER :: start_frame = 0 + INTEGER :: end_frame = -1 END TYPE period_efield_type ! ************************************************************************************************** diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 70aa70b360..fb77ed3bdb 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -513,6 +513,39 @@ CONTAINS dft_control%period_efield%displacement_field = .FALSE. CALL section_vals_val_get(tmp_section, "DISPLACEMENT_FIELD", & l_val=dft_control%period_efield%displacement_field) + + CALL section_vals_val_get(tmp_section, "INTENSITY_LIST", r_vals=pol) + + IF (SIZE(pol) > 1 .OR. pol(1) /= 0.0_dp) THEN + ! if INTENSITY_LIST is present, INTENSITY must not be present + IF (dft_control%period_efield%strength /= 0.0_dp) THEN + CPABORT("[PERIODIC FIELD] Only one of INTENSITY or INTENSITY_LIST can be specified.") + END IF + + ALLOCATE (dft_control%period_efield%strength_list(SIZE(pol))) + dft_control%period_efield%strength_list(1:SIZE(pol)) = pol(1:SIZE(pol)) + END IF + + CALL section_vals_val_get(tmp_section, "START_FRAME", & + i_val=dft_control%period_efield%start_frame) + CALL section_vals_val_get(tmp_section, "END_FRAME", & + i_val=dft_control%period_efield%end_frame) + + IF (dft_control%period_efield%end_frame /= -1) THEN + ! check if valid bounds are given + ! if an end frame is given, the number of active frames must be a + ! multiple of the number of intensities + IF (dft_control%period_efield%start_frame > dft_control%period_efield%end_frame) THEN + CPABORT("[PERIODIC FIELD] START_FRAME > END_FRAME") + ELSE IF (dft_control%period_efield%start_frame < 1) THEN + CPABORT("[PERIODIC FIELD] START_FRAME < 1") + ELSE IF (MOD(dft_control%period_efield%end_frame - & + dft_control%period_efield%start_frame + 1, SIZE(pol)) /= 0) THEN + CALL cp_abort(__LOCATION__, & + "[PERIODIC FIELD] Number of active frames must be a multiple of the number of intensities") + END IF + END IF + ! periodic fields don't work with RTP CPASSERT(.NOT. do_rtp) IF (dft_control%period_efield%displacement_field) THEN @@ -1804,7 +1837,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'write_dft_control' CHARACTER(LEN=20) :: tmpStr - INTEGER :: handle, output_unit + INTEGER :: handle, i, output_unit REAL(kind=dp) :: density_cut, density_smooth_cut_range, & gradient_cut, tau_cut TYPE(cp_logger_type), POINTER :: logger @@ -1939,9 +1972,32 @@ CONTAINS "PERIODIC_EFIELD| y", & dft_control%period_efield%polarisation(2), & "PERIODIC_EFIELD| z", & - dft_control%period_efield%polarisation(3), & - "PERIODIC_EFIELD| Intensity [a.u.]:", & - dft_control%period_efield%strength + dft_control%period_efield%polarisation(3) + + WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,I14)") & + "PERIODIC_EFIELD| Start Frame:", & + dft_control%period_efield%start_frame, & + "PERIODIC_EFIELD| End Frame:", & + dft_control%period_efield%end_frame + + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,I14)") & + "PERIODIC_EFIELD| Number of Intensities:", & + SIZE(dft_control%period_efield%strength_list) + WRITE (UNIT=output_unit, FMT="(T2,A,I10,T66,1X,ES14.6)") & + "PERIODIC_EFIELD| Intensity List [a.u.] ", & + 1, dft_control%period_efield%strength_list(1) + DO i = 2, SIZE(dft_control%period_efield%strength_list) + WRITE (UNIT=output_unit, FMT="(T2,A,I10,T66,1X,ES14.6)") & + "PERIODIC_EFIELD| ", & + i, dft_control%period_efield%strength_list(i) + END DO + ELSE + WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,ES14.6)") & + "PERIODIC_EFIELD| Intensity [a.u.]:", & + dft_control%period_efield%strength + END IF + IF (SQRT(DOT_PRODUCT(dft_control%period_efield%polarisation, & dft_control%period_efield%polarisation)) < EPSILON(0.0_dp)) THEN CPABORT("Invalid (too small) polarisation vector specified for PERIODIC_EFIELD") diff --git a/src/eeq_method.F b/src/eeq_method.F index b03cf7401d..6bb811a1c9 100644 --- a/src/eeq_method.F +++ b/src/eeq_method.F @@ -1529,6 +1529,10 @@ CONTAINS IF (dft_control%apply_period_efield) THEN dfield = dft_control%period_efield%displacement_field + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + CPABORT("use of strength_list not implemented for eeq_efield_energy") + END IF + fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) fieldpol = -fieldpol*dft_control%period_efield%strength @@ -1607,6 +1611,10 @@ CONTAINS IF (dft_control%apply_period_efield) THEN dfield = dft_control%period_efield%displacement_field + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + CPABORT("use of strength_list not implemented for eeq_efield_pot") + END IF + fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) fieldpol = -fieldpol*dft_control%period_efield%strength @@ -1673,6 +1681,10 @@ CONTAINS natom = SIZE(particle_set) + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + CPABORT("use of strength_list not implemented for eeq_dfield_pot") + END IF + dfilter(1:3) = dft_control%period_efield%d_filter(1:3) fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) @@ -1788,6 +1800,10 @@ CONTAINS dfield = dft_control%period_efield%displacement_field CPASSERT(.NOT. dfield) + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + CPABORT("use of strength_list not implemented for eeq_efield_force_periodic") + END IF + fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) fieldpol = -fieldpol*dft_control%period_efield%strength diff --git a/src/efield_tb_methods.F b/src/efield_tb_methods.F index da2a9b28f1..15ce8fddf5 100644 --- a/src/efield_tb_methods.F +++ b/src/efield_tb_methods.F @@ -92,10 +92,15 @@ CONTAINS CALL get_qs_env(qs_env, dft_control=dft_control) IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN IF (dft_control%apply_period_efield) THEN - IF (dft_control%period_efield%displacement_field) THEN - CALL dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) - ELSE - CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) + ! check if the periodic efield should be applied in the current step + IF (dft_control%period_efield%start_frame <= qs_env%sim_step .AND. & + (dft_control%period_efield%end_frame == -1 .OR. dft_control%period_efield%end_frame >= qs_env%sim_step)) THEN + + IF (dft_control%period_efield%displacement_field) THEN + CALL dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) + ELSE + CALL efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) + END IF END IF ELSE IF (dft_control%apply_efield) THEN CALL efield_tb_local(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy) @@ -285,7 +290,7 @@ CONTAINS INTEGER, DIMENSION(3) :: cellind INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: found, use_virial - REAL(KIND=dp) :: charge, dd, dr, fdir, fi + REAL(KIND=dp) :: charge, dd, dr, fdir, fi, strength REAL(KIND=dp), DIMENSION(3) :: fieldpol, fij, forcea, fpolvec, kvec, & qi, rab, ria, rib, rij REAL(KIND=dp), DIMENSION(3, 3) :: hmat @@ -321,9 +326,16 @@ CONTAINS use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) use_virial = use_virial .AND. calculate_forces + ! if an intensities list is given, select the value for the current step + strength = dft_control%period_efield%strength + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + strength = dft_control%period_efield%strength_list(MOD(qs_env%sim_step & + - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1) + END IF + fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) - fieldpol = -fieldpol*dft_control%period_efield%strength + fieldpol = -fieldpol*strength hmat = cell%hmat(:, :)/twopi DO idir = 1, 3 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir) @@ -608,7 +620,8 @@ CONTAINS INTEGER, DIMENSION(3) :: cellind INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: found, use_virial - REAL(KIND=dp) :: charge, dd, ener_field, fdir, omega + REAL(KIND=dp) :: charge, dd, ener_field, fdir, omega, & + strength REAL(KIND=dp), DIMENSION(3) :: ci, cqi, dfilter, di, fieldpol, fij, & hdi, kvec, qi, rab, ria, rib REAL(KIND=dp), DIMENSION(3, 3) :: hmat @@ -653,9 +666,16 @@ CONTAINS dfilter(1:3) = dft_control%period_efield%d_filter(1:3) + ! if an intensities list is given, select the value for the current step + strength = dft_control%period_efield%strength + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + strength = dft_control%period_efield%strength_list(MOD(qs_env%sim_step & + - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1) + END IF + fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) - fieldpol = fieldpol*dft_control%period_efield%strength + fieldpol = fieldpol*strength omega = cell%deth hmat = cell%hmat(:, :)/twopi diff --git a/src/input_cp2k_field.F b/src/input_cp2k_field.F index 77fe3fc88d..f3c4b4c362 100644 --- a/src/input_cp2k_field.F +++ b/src/input_cp2k_field.F @@ -65,7 +65,8 @@ CONTAINS NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="INTENSITY", & - description="Intensity of the electric field in a.u", & + description="Intensity of the electric field in a.u, "// & + "not allowed together with INTENSITY_LIST", & usage="INTENSITY 0.001", & default_r_val=0._dp) CALL section_add_keyword(section, keyword) @@ -89,13 +90,43 @@ CONTAINS CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="D_FILTER", & - description="Filter for displacement field (x,y,z-dirction)", & + description="Filter for displacement field (x,y,z-direction)", & usage="D_FILTER 1.0 0.0 0.0", & repeats=.FALSE., n_var=3, & type_of_var=real_t, default_r_vals=(/1.0_dp, 1.0_dp, 1.0_dp/)) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="INTENSITY_LIST", & + description="Intensities of the electric field in a.u. "// & + "They are applied sequentially, one per frame. "// & + "If the number of frames exceeds the number of values, "// & + "the list is cyclically repeated. Attention: not implemented for eeq.", & + usage="INTENSITY_LIST {real} {real} .. {real}", & + n_var=-1, type_of_var=real_t, default_r_vals=(/0.0_dp/)) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="START_FRAME", & + description="First frame the field is applied. "// & + "(0: first frame) "// & + "Attention: ignored for eeq", & + usage="START_FRAME 0", & + default_i_val=0) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="END_FRAME", & + description="Last frame the field is applied. "// & + "If an end frame is specified, the number of active frames "// & + "must be a multiple of the number of "// & + "the given intensity values. (-1: no end) "// & + "Attention: ignored for eeq", & + usage="END_FRAME -1", & + default_i_val=-1) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_per_efield_section ! ************************************************************************************************** !> \brief creates the section for time dependent nonperiodic fields diff --git a/src/qs_efield_berry.F b/src/qs_efield_berry.F index 23e5cae0d2..644141b609 100644 --- a/src/qs_efield_berry.F +++ b/src/qs_efield_berry.F @@ -121,11 +121,16 @@ CONTAINS dft_control=dft_control) IF (dft_control%apply_period_efield) THEN - IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env) - IF (dft_control%period_efield%displacement_field) THEN - CALL qs_dispfield_derivatives(qs_env, just_energy, calculate_forces) - ELSE - CALL qs_efield_derivatives(qs_env, just_energy, calculate_forces) + ! check if the periodic efield should be applied in the current step + IF (dft_control%period_efield%start_frame <= qs_env%sim_step .AND. & + (dft_control%period_efield%end_frame == -1 .OR. dft_control%period_efield%end_frame >= qs_env%sim_step)) THEN + + IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env) + IF (dft_control%period_efield%displacement_field) THEN + CALL qs_dispfield_derivatives(qs_env, just_energy, calculate_forces) + ELSE + CALL qs_efield_derivatives(qs_env, just_energy, calculate_forces) + END IF END IF END IF @@ -199,7 +204,7 @@ CONTAINS REAL(dp) :: charge, ci(3), cqi(3), dab, dd, & ener_field, f0, fab, fieldpol(3), & focc, fpolvec(3), hmat(3, 3), occ, & - qi(3), ti(3) + qi(3), strength, ti(3) REAL(dp), DIMENSION(3) :: forcea, forceb, kvec, ra, rab, rb, ria REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab @@ -250,9 +255,16 @@ CONTAINS CPABORT("Stress tensor for periodic E-field not implemented") END IF + ! if an intensities list is given, select the value for the current step + strength = dft_control%period_efield%strength + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + strength = dft_control%period_efield%strength_list(MOD(qs_env%sim_step & + - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1) + END IF + fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) - fieldpol = -fieldpol*dft_control%period_efield%strength + fieldpol = -fieldpol*strength hmat = cell%hmat(:, :)/twopi DO idir = 1, 3 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir) @@ -677,10 +689,8 @@ CONTAINS npgfb, nsgfa, nsgfb INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb LOGICAL :: found, uniform, use_virial - REAL(dp) :: charge, ci(3), cqi(3), dab, dd, di(3), & - ener_field, fab, fieldpol(3), focc, & - hmat(3, 3), occ, omega, qi(3), & - rlog(3), zlog(3) + REAL(dp) :: charge, ci(3), cqi(3), dab, dd, di(3), ener_field, fab, fieldpol(3), focc, & + hmat(3, 3), occ, omega, qi(3), rlog(3), strength, zlog(3) REAL(dp), DIMENSION(3) :: dfilter, forcea, forceb, kvec, ra, rab, & rb, ria REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work @@ -734,9 +744,16 @@ CONTAINS dfilter(1:3) = dft_control%period_efield%d_filter(1:3) + ! if an intensities list is given, select the value for the current step + strength = dft_control%period_efield%strength + IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN + strength = dft_control%period_efield%strength_list(MOD(qs_env%sim_step & + - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1) + END IF + fieldpol = dft_control%period_efield%polarisation fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) - fieldpol = fieldpol*dft_control%period_efield%strength + fieldpol = fieldpol*strength omega = cell%deth hmat = cell%hmat(:, :)/(twopi*omega) diff --git a/tests/QS/regtest-p-efield/H2O-field-list.inp b/tests/QS/regtest-p-efield/H2O-field-list.inp new file mode 100644 index 0000000000..1b567095fb --- /dev/null +++ b/tests/QS/regtest-p-efield/H2O-field-list.inp @@ -0,0 +1,70 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT H2O-field-list + RUN_TYPE MD +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 32 + TEMPERATURE 298 + TIMESTEP 0.4 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 200 + &END MGRID + &PERIODIC_EFIELD + END_FRAME 27 + INTENSITY_LIST 0.0002 0.0005 0.001 0.0005 0.0002 0.0 -0.0002 -0.0005 -0.001 -0.0005 -0.0002 0.0 + POLARISATION 1 0 0 + START_FRAME 4 + &END PERIODIC_EFIELD + &PRINT + &MOMENTS + &END MOMENTS + &END PRINT + &QS + EPS_DEFAULT 1.0E-12 + EXTRAPOLATION PS + EXTRAPOLATION_ORDER 2 + &END QS + &SCF + EPS_SCF 1.0E-6 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 6 + SCF_GUESS ATOMIC + &OT + &END OT + &END SCF + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 H2O + H 0.000000 -0.757136 0.520545 H2O + H 0.000000 0.757136 0.520545 H2O + &END COORD + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-p-efield/TEST_FILES b/tests/QS/regtest-p-efield/TEST_FILES index 4dbbd5bc37..05b0005922 100644 --- a/tests/QS/regtest-p-efield/TEST_FILES +++ b/tests/QS/regtest-p-efield/TEST_FILES @@ -8,6 +8,7 @@ H2O-field-gopt.inp 1 3e-11 H2O-field-gopt-lsd.inp 1 4e-12 -16.84348812829990 H2O-field.inp 1 4e-13 -15.90320150373305 H2O-field-lsd.inp 1 4e-14 -15.92914916395207 +H2O-field-list.inp 2 1e-14 -0.171095788220e+02 HF-field.inp 1 1e-12 -24.70767796237102 HF-field-gopt.inp 1 1e-12 -24.70292863028889 HF-field-debug.inp 0 1e-12