From f1eb7d55b906f1efd150ecc101bf7623c6b37740 Mon Sep 17 00:00:00 2001 From: Samuel Andermatt Date: Thu, 21 Apr 2016 14:18:13 +0000 Subject: [PATCH] first implementation to add custom efields svn-origin-rev: 16834 --- src/cp_control_utils.F | 36 ++++++-- src/efield_utils.F | 14 ++- src/input_constants.F | 3 +- src/input_cp2k_dft.F | 52 +++++++++-- .../QS/regtest-rtp-2/H2-emd-efield-custom.inp | 86 +++++++++++++++++++ tests/QS/regtest-rtp-2/TEST_FILES | 1 + tests/QS/regtest-rtp-2/field.txt | 5 ++ 7 files changed, 176 insertions(+), 21 deletions(-) create mode 100644 tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp create mode 100644 tests/QS/regtest-rtp-2/field.txt diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 3de59d60c1..0bf683319a 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -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 diff --git a/src/efield_utils.F b/src/efield_utils.F index 21a2ddfc07..aa89e6f6f5 100644 --- a/src/efield_utils.F +++ b/src/efield_utils.F @@ -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 diff --git a/src/input_constants.F b/src/input_constants.F index fd8f39afe7..477856d393 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -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, & diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index e812abf23e..9bcf37480e 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -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 ... diff --git a/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp b/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp new file mode 100644 index 0000000000..7f0adf1c5c --- /dev/null +++ b/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp @@ -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 diff --git a/tests/QS/regtest-rtp-2/TEST_FILES b/tests/QS/regtest-rtp-2/TEST_FILES index 800afc8c2d..5920dfd4ab 100644 --- a/tests/QS/regtest-rtp-2/TEST_FILES +++ b/tests/QS/regtest-rtp-2/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-rtp-2/field.txt b/tests/QS/regtest-rtp-2/field.txt new file mode 100644 index 0000000000..24bade207f --- /dev/null +++ b/tests/QS/regtest-rtp-2/field.txt @@ -0,0 +1,5 @@ +0 +1000000000 +2000000000 +1000000000 +0