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

first implementation to add custom efields

Browse source 
svn-origin-rev: 16834
This commit is contained in:
Samuel Andermatt 2016-04-21 14:18:13 +00:00
parent 428da2a35b
commit f1eb7d55b9
7 changed files with 176 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

36
src/cp_control_utils.F
View file

@ -16,13 +16,16 @@ MODULE cp_control_utils
admm_control_create, admm_control_type, ddapc_control_create, ddapc_restraint_type, &
dft_control_create, dft_control_type, efield_type, qs_control_type, sccs_control_create, &
tddfpt_control_create, tddfpt_control_type
USE cp_files, ONLY: close_file,&
open_file
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE cp_units, ONLY: cp_unit_from_cp2k
USE cp_units, ONLY: cp_unit_from_cp2k,&
cp_unit_to_cp2k
USE input_constants, ONLY: &
constant_env, do_admm_basis_projection, do_admm_blocked_projection, &
constant_env, custom_env, do_admm_basis_projection, do_admm_blocked_projection, &
do_admm_blocking_purify_full, do_admm_charge_constrained_projection, &
do_admm_exch_scaling_merlot, do_admm_purify_mcweeny, do_admm_purify_mo_diag, &
do_admm_purify_mo_no_diag, do_admm_purify_none, do_admm_purify_none_dm, &
@ -1633,7 +1636,8 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'read_efield_sections', &
routineP = moduleN//':'//routineN
INTEGER :: i
CHARACTER(len=default_path_length) :: file_name
INTEGER :: i, io, j, n, unit_nr
REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_vals
TYPE(efield_type), POINTER :: efield
TYPE(section_vals_type), POINTER :: tmp_section
@ -1656,6 +1660,7 @@ CONTAINS
i_val=efield%envelop_id)
CALL section_vals_val_get(efield_section, "WAVELENGTH", i_rep_section=i, &
r_val=efield%wavelength)
IF (efield%envelop_id == constant_env) THEN
ALLOCATE (efield%envelop_i_vars(2))
tmp_section => section_vals_get_subs_vals(efield_section, "CONSTANT_ENV", i_rep_section=i)
@ -1663,16 +1668,14 @@ CONTAINS
i_val=efield%envelop_i_vars(1))
CALL section_vals_val_get(tmp_section, "END_STEP", &
i_val=efield%envelop_i_vars(2))
END IF
IF (efield%envelop_id == gaussian_env) THEN
ELSE IF (efield%envelop_id == gaussian_env) THEN
ALLOCATE (efield%envelop_r_vars(2))
tmp_section => section_vals_get_subs_vals(efield_section, "GAUSSIAN_ENV", i_rep_section=i)
CALL section_vals_val_get(tmp_section, "T0", &
r_val=efield%envelop_r_vars(1))
CALL section_vals_val_get(tmp_section, "SIGMA", &
r_val=efield%envelop_r_vars(2))
END IF
IF (efield%envelop_id == ramp_env) THEN
ELSE IF (efield%envelop_id == ramp_env) THEN
ALLOCATE (efield%envelop_i_vars(4))
tmp_section => section_vals_get_subs_vals(efield_section, "RAMP_ENV", i_rep_section=i)
CALL section_vals_val_get(tmp_section, "START_STEP_IN", &
@ -1683,6 +1686,25 @@ CONTAINS
i_val=efield%envelop_i_vars(3))
CALL section_vals_val_get(tmp_section, "END_STEP_OUT", &
i_val=efield%envelop_i_vars(4))
ELSE IF (efield%envelop_id == custom_env) THEN
tmp_section => section_vals_get_subs_vals(efield_section, "CUSTOM_ENV", i_rep_section=i)
CALL section_vals_val_get(tmp_section, "EFIELD_FILE_NAME", c_val=file_name)
CALL open_file(file_name=TRIM(file_name), file_action="READ", file_status="OLD", unit_number=unit_nr)
!Determine the number of lines in file
n = -1
io = 0
DO WHILE (io == 0)
READ (unit_nr, *, iostat=io)
n = n+1
END DO
!Read the file
REWIND (unit_nr)
ALLOCATE (efield%envelop_r_vars(n))
DO j = 1, n
READ (unit_nr, *) efield%envelop_r_vars(j)
efield%envelop_r_vars(j) = cp_unit_to_cp2k(cp_unit_to_cp2k(efield%envelop_r_vars(j), "hartree"), "m", power=-1)
END DO
CALL close_file(unit_nr)
END IF
END DO
END SUBROUTINE read_efield_sections

14
src/efield_utils.F
View file

@ -16,6 +16,7 @@ MODULE efield_utils
efield_type
USE cp_para_types, ONLY: cp_para_env_type
USE input_constants, ONLY: constant_env,&
custom_env,&
gaussian_env,&
ramp_env
USE kahan_sum, ONLY: accurate_sum
@ -150,8 +151,7 @@ CONTAINS
field = field+strength*COS(sim_time*nu*2.0_dp*pi+ &
efield%phase_offset*pi)*pol(:)
END IF
END IF
IF (efield%envelop_id == ramp_env) THEN
ELSE IF (efield%envelop_id == ramp_env) THEN
IF (sim_step .GE. efield%envelop_i_vars(1) .AND. sim_step .LE. efield%envelop_i_vars(2)) &
strength = strength*(sim_step-efield%envelop_i_vars(1))/(efield%envelop_i_vars(2)-efield%envelop_i_vars(1))
IF (sim_step .GE. efield%envelop_i_vars(3) .AND. sim_step .LE. efield%envelop_i_vars(4)) &
@ -160,11 +160,17 @@ CONTAINS
IF (sim_step .LE. efield%envelop_i_vars(1)) strength = 0.0_dp
field = field+strength*COS(sim_time*nu*2.0_dp*pi+ &
efield%phase_offset*pi)*pol(:)
END IF
IF (efield%envelop_id == gaussian_env) THEN
ELSE IF (efield%envelop_id == gaussian_env) THEN
env = EXP(-0.5_dp*((sim_time-efield%envelop_r_vars(1))/efield%envelop_r_vars(2))**2.0_dp)
field = field+strength*env*COS(sim_time*nu*2.0_dp*pi+ &
efield%phase_offset*pi)*pol(:)
ELSE IF (efield%envelop_id == custom_env) THEN
IF (sim_step .LE. SIZE(efield%envelop_r_vars)) THEN
strength = efield%envelop_r_vars(sim_step)
ELSE
strength = 0.0_dp
ENDIF
field = field+strength*pol(:)
END IF
END DO

3
src/input_constants.F
View file

@ -790,7 +790,8 @@ MODULE input_constants
INTEGER, PARAMETER, PUBLIC :: constant_env = 1, &
gaussian_env = 2, &
ramp_env = 3
ramp_env = 3, &
custom_env = 4
! how to solve polarizable force fields
INTEGER, PARAMETER, PUBLIC :: do_fist_pol_none = 1, &

52
src/input_cp2k_dft.F
View file

@ -29,11 +29,11 @@ MODULE input_cp2k_dft
USE cp_units, ONLY: cp_unit_to_cp2k
USE input_constants, ONLY: &
atomic_guess, cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_reduce, &
cholesky_restore, constant_env, core_guess, diag_block_davidson, diag_block_krylov, &
diag_filter_matrix, diag_ot, diag_standard, dispersion_d3, dispersion_uff, &
do_admm_aux_exch_func_bee, do_admm_aux_exch_func_default, do_admm_aux_exch_func_none, &
do_admm_aux_exch_func_opt, do_admm_aux_exch_func_pbex, do_admm_basis_projection, &
do_admm_blocked_projection, do_admm_blocking_purify_full, &
cholesky_restore, constant_env, core_guess, custom_env, diag_block_davidson, &
diag_block_krylov, diag_filter_matrix, diag_ot, diag_standard, dispersion_d3, &
dispersion_uff, do_admm_aux_exch_func_bee, do_admm_aux_exch_func_default, &
do_admm_aux_exch_func_none, do_admm_aux_exch_func_opt, do_admm_aux_exch_func_pbex, &
do_admm_basis_projection, do_admm_blocked_projection, do_admm_blocking_purify_full, &
do_admm_charge_constrained_projection, do_admm_exch_scaling_merlot, &
do_admm_exch_scaling_none, do_admm_purify_cauchy, do_admm_purify_cauchy_subspace, &
do_admm_purify_mcweeny, do_admm_purify_mo_diag, do_admm_purify_mo_no_diag, &
@ -5569,11 +5569,12 @@ CONTAINS
description="Shape of the efield pulse used in real-time propagation (RTP) calculations.", &
usage="ENVELOP CONSTANT", &
default_i_val=constant_env, &
enum_c_vals=s2a("CONSTANT", "GAUSSIAN", "RAMP"), &
enum_c_vals=s2a("CONSTANT", "GAUSSIAN", "RAMP", "CUSTOM"), &
enum_desc=s2a("No envelop function is applied to the strength", &
"A Gaussian function is used as envelop ", &
"Linear tune in/out of the field"), &
enum_i_vals=(/constant_env, gaussian_env, ramp_env/))
"Linear tune in/out of the field", &
"A custom field read from a file"), &
enum_i_vals=(/constant_env, gaussian_env, ramp_env, custom_env/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
@ -5581,11 +5582,15 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_gaussian_env_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_ramp_env_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_gaussian_env_section(subsection)
CALL create_custom_env_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
@ -7014,6 +7019,7 @@ CONTAINS
CALL keyword_release(keyword)
END SUBROUTINE create_gaussian_env_section
! **************************************************************************************************
!> \brief ...
!> \param section ...
@ -7063,6 +7069,34 @@ CONTAINS
END SUBROUTINE create_ramp_env_section
! **************************************************************************************************
!> \brief ...
!> \param section ...
! **************************************************************************************************
SUBROUTINE create_custom_env_section(section)
TYPE(section_type), POINTER :: section
CHARACTER(len=*), PARAMETER :: routineN = 'create_custom_env_section', &
routineP = moduleN//':'//routineN
TYPE(keyword_type), POINTER :: keyword
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, "CUSTOM_ENV", &
description="Parameters for a custom efield", &
n_keywords=6, n_subsections=1, repeats=.TRUE.)
NULLIFY (keyword)
CALL keyword_create(keyword, name="EFIELD_FILE_NAME", &
description="Specify file that contains the electric field [V/m]", &
usage="EFIELD_FILE_NAME filename", &
n_var=1, type_of_var=char_t, default_c_val="")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_custom_env_section
! **************************************************************************************************
!> \brief Create CP2K input section for the SCCS model
!> \param section ...

86
tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp Normal file
View file

@ -0,0 +1,86 @@
&FORCE_EVAL
METHOD QUICKSTEP
&DFT
&REAL_TIME_PROPAGATION
MAX_ITER 5
MAT_EXP TAYLOR
EPS_ITER 1.0E-9
INITIAL_WFN SCF_WFN
&END
&EFIELD
POLARISATION 0 1 0
ENVELOP CUSTOM
&CUSTOM_ENV
EFIELD_FILE_NAME field.txt
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 100
&END MGRID
&QS
EPS_DEFAULT 1.0E-10
MAP_CONSISTENT
&END QS
&SCF
MAX_SCF 60
EPS_SCF 1.0E-3
SCF_GUESS ATOMIC
&OUTER_SCF
EPS_SCF 1.0E-3
&END
&OT ON
ROTATION
&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
H 6.000000 4.5000000 6.000000
H 6.000000 1.5000000 6.000000
&END COORD
&KIND H
BASIS_SET DZVP-GTH-PBE
POTENTIAL GTH-PBE-q1
&END KIND
&KIND HX
ELEMENT H
GHOST T
BASIS_SET DZVP-GTH-PBE
POTENTIAL GTH-PBE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PBE
POTENTIAL GTH-PBE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2-emd-efield-custom
RUN_TYPE EHRENFEST_DYN
! RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
STEPS 6
TIMESTEP [au_t] 0.25
TEMPERATURE 300.0
&END MD
&PRINT
&TRAJECTORY
&EACH
MD 1
&END
&END
&END
&END MOTION

1
tests/QS/regtest-rtp-2/TEST_FILES
View file

@ -10,6 +10,7 @@ H2-rtp_restart-1.inp 1 3e-13 -0.90223968361437
H2-rtp-efield.inp 1 2e-10 -0.8296879654
H2-emd-efield.inp 2 4e-11 -0.894866854822E+00
H2-emd-efield-ramp.inp 2 1e-14 -0.886674608725
H2-emd-efield-custom.inp 2 1e-14 -0.896362421571
H2-rtp_ETRS_ARNOLDI.inp 1 3e-13 -0.90223968349550
H2-emd_ETRS_ARNOLDI.inp 1 1e-10 -17.085427015
#EOF

5
tests/QS/regtest-rtp-2/field.txt Normal file
View file

@ -0,0 +1,5 @@
0
1000000000
2000000000
1000000000
0