comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Add INTENSITY_LIST, START_FRAME and END_FRAME keywords for periodic efield

Browse source 
Co-authored-by: Leo Decking <leo@decking.dev>
This commit is contained in:
Leo Decking 2025-02-26 09:44:50 +01:00 • committed by GitHub
parent 71e52877b8
commit 60cf13d4c5
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 240 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/cp_control_types.F
View file

@ -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
! **************************************************************************************************

64
src/cp_control_utils.F
View file

@ -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")

16
src/eeq_method.F
View file

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

36
src/efield_tb_methods.F
View file

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

35
src/input_cp2k_field.F
View file

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

41
src/qs_efield_berry.F
View file

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

70
tests/QS/regtest-p-efield/H2O-field-list.inp Normal file
View file

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

1
tests/QS/regtest-p-efield/TEST_FILES
View file

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