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

EC| Add print section for EC properties

Browse source
This commit is contained in:
Fabian Belleflamme 2021-12-03 11:58:20 +01:00 • committed by Ole Schütt
parent 1e94a9dfe3
commit ffa6531213
19 changed files with 516 additions and 96 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/ec_env_types.F
View file

@ -83,6 +83,7 @@ MODULE ec_env_types
TYPE(qs_dispersion_type), POINTER :: dispersion_env
! matrices in complete basis
! KS: Kohn-Sham; H: Core; S: overlap; T: kinetic energy;
! P: Harris density, W: Harris energy weighted density
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, &
matrix_h, &
matrix_s, &
@ -95,6 +96,8 @@ MODULE ec_env_types
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: cpmos
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_z, matrix_wz, z_admm
! Harris (rhoout), and response density (rhoz) on grid
TYPE(pw_p_type), DIMENSION(:), POINTER :: rhoout_r, rhoz_r
! potentials from input density
TYPE(pw_p_type), POINTER :: vh_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace

9
src/ec_environment.F
View file

@ -107,7 +107,9 @@ CONTAINS
TYPE(section_vals_type), POINTER :: dft_section
TYPE(section_vals_type), OPTIONAL, POINTER :: ec_section
INTEGER :: ikind, maxlgto, nkind, unit_nr
CHARACTER(LEN=*), PARAMETER :: routineN = 'init_ec_env', routineP = moduleN//':'//routineN
INTEGER :: handle, ikind, maxlgto, nkind, unit_nr
LOGICAL :: explicit
REAL(KIND=dp) :: eps_pgf_orb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@ -121,6 +123,8 @@ CONTAINS
TYPE(section_vals_type), POINTER :: nl_section, pp_section, section1, &
section2, xc_section
CALL timeset(routineN, handle)
NULLIFY (atomic_kind_set, dispersion_env, ec_env%ls_env, para_env)
NULLIFY (ec_env%sab_orb, ec_env%sac_ppl, ec_env%sap_ppnl)
NULLIFY (ec_env%matrix_ks, ec_env%matrix_h, ec_env%matrix_s)
@ -140,6 +144,7 @@ CONTAINS
NULLIFY (ec_env%vxc_rspace)
NULLIFY (ec_env%vtau_rspace)
NULLIFY (ec_env%vadmm_rspace)
NULLIFY (ec_env%rhoout_r, ec_env%rhoz_r)
ec_env%should_update = .TRUE.
ec_env%mao = .FALSE.
@ -288,6 +293,8 @@ CONTAINS
END IF
CALL timestop(handle)
END SUBROUTINE init_ec_env
! **************************************************************************************************

162
src/energy_corrections.F
View file

@ -118,12 +118,14 @@ MODULE energy_corrections
USE pw_poisson_types, ONLY: pw_poisson_type
USE pw_pool_types, ONLY: pw_pool_create_pw,&
pw_pool_give_back_pw,&
pw_pool_p_type,&
pw_pool_type
USE pw_types, ONLY: COMPLEXDATA1D,&
REALDATA3D,&
REALSPACE,&
RECIPROCALSPACE,&
pw_p_type
pw_p_type,&
pw_type
USE qs_collocate_density, ONLY: calculate_rho_elec
USE qs_core_energies, ONLY: calculate_ecore_overlap,&
calculate_ptrace
@ -178,6 +180,7 @@ MODULE energy_corrections
virial_create,&
virial_release,&
virial_type
USE voronoi_interface, ONLY: entry_voronoi_or_bqb
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -351,7 +354,11 @@ CONTAINS
LOGICAL, INTENT(IN) :: calculate_forces
INTEGER, INTENT(IN) :: unit_nr
INTEGER :: ispin, nspins
REAL(KIND=dp) :: exc
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
IF (ec_env%should_update) THEN
CALL ec_build_neighborlist(qs_env, ec_env)
@ -385,6 +392,19 @@ CONTAINS
CALL response_calculation(qs_env, ec_env)
! Allocate response density on real space grid for use in properties
! Calculated in response_force
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
nspins = dft_control%nspins
CPASSERT(ASSOCIATED(pw_env))
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
ALLOCATE (ec_env%rhoz_r(nspins))
DO ispin = 1, nspins
NULLIFY (ec_env%rhoz_r(nspins)%pw)
CALL pw_pool_create_pw(auxbas_pw_pool, ec_env%rhoz_r(ispin)%pw, &
use_data=REALDATA3D, in_space=REALSPACE)
END DO
CALL response_force(qs_env, &
vh_rspace=ec_env%vh_rspace, &
vxc_rspace=ec_env%vxc_rspace, &
@ -394,12 +414,20 @@ CONTAINS
matrix_pz=ec_env%matrix_z, &
matrix_pz_admm=ec_env%z_admm, &
matrix_wz=ec_env%matrix_wz, &
rhopz_r=ec_env%rhoz_r, &
zehartree=ec_env%ehartree, &
p_env=ec_env%p_env, &
zexc=ec_env%exc)
CALL ec_properties(qs_env, ec_env)
! Deallocate Harris density and response density on grid
DO ispin = 1, nspins
CALL pw_pool_give_back_pw(auxbas_pw_pool, ec_env%rhoout_r(ispin)%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, ec_env%rhoz_r(ispin)%pw)
END DO
DEALLOCATE (ec_env%rhoout_r, ec_env%rhoz_r)
! Deallocate matrices
IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
@ -1016,6 +1044,15 @@ CONTAINS
task_list_external=ec_env%task_list)
END DO
! Save Harris on real space grid for use in properties
ALLOCATE (ec_env%rhoout_r(nspins))
DO ispin = 1, nspins
NULLIFY (ec_env%rhoout_r(ispin)%pw)
CALL pw_pool_create_pw(auxbas_pw_pool, ec_env%rhoout_r(ispin)%pw, &
use_data=REALDATA3D, in_space=REALSPACE)
CALL pw_copy(rhoout_r(ispin)%pw, ec_env%rhoout_r(ispin)%pw)
END DO
IF (use_virial) THEN
! Calculate the Hartree potential
@ -1985,12 +2022,12 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_properties'
CHARACTER(LEN=8), DIMENSION(3) :: rlab
CHARACTER(LEN=default_path_length) :: filename
CHARACTER(LEN=default_path_length) :: filename, my_pos_voro
CHARACTER(LEN=default_string_length) :: description
INTEGER :: akind, handle, i, ia, iatom, idir, &
ikind, iounit, ispin, maxmom, nspins, &
reference, unit_nr
LOGICAL :: magnetic, periodic
INTEGER :: akind, handle, i, ia, iatom, idir, ikind, iounit, ispin, maxmom, nspins, &
reference, should_print_bqb, should_print_voro, unit_nr, unit_nr_voro
LOGICAL :: append_voro, magnetic, periodic, &
voro_print_txt
REAL(KIND=dp) :: charge, dd, focc, tmp
REAL(KIND=dp), DIMENSION(3) :: cdip, pdip, rcc, rdip, ria, tdip
REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
@ -2000,10 +2037,17 @@ CONTAINS
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(cp_result_type), POINTER :: results
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
TYPE(dft_control_type), POINTER :: dft_control
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: rho_elec_rspace
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(pw_type), POINTER :: mb_rho
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(section_vals_type), POINTER :: print_key
TYPE(section_vals_type), POINTER :: ec_section, print_key, print_key_bqb, &
print_key_voro
CALL timeset(routineN, handle)
@ -2018,23 +2062,27 @@ CONTAINS
iounit = -1
END IF
print_key => section_vals_get_subs_vals(section_vals=qs_env%input, &
subsection_name="DFT%PRINT%MOMENTS")
NULLIFY (dft_control)
CALL get_qs_env(qs_env, dft_control=dft_control)
ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
print_key => section_vals_get_subs_vals(section_vals=ec_section, &
subsection_name="PRINT%MOMENTS")
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
maxmom = section_get_ival(section_vals=qs_env%input, &
keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT")
periodic = section_get_lval(section_vals=qs_env%input, &
keyword_name="DFT%PRINT%MOMENTS%PERIODIC")
reference = section_get_ival(section_vals=qs_env%input, &
keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
magnetic = section_get_lval(section_vals=qs_env%input, &
keyword_name="DFT%PRINT%MOMENTS%MAGNETIC")
maxmom = section_get_ival(section_vals=ec_section, &
keyword_name="PRINT%MOMENTS%MAX_MOMENT")
periodic = section_get_lval(section_vals=ec_section, &
keyword_name="PRINT%MOMENTS%PERIODIC")
reference = section_get_ival(section_vals=ec_section, &
keyword_name="PRINT%MOMENTS%REFERENCE")
magnetic = section_get_lval(section_vals=ec_section, &
keyword_name="PRINT%MOMENTS%MAGNETIC")
NULLIFY (ref_point)
CALL section_vals_val_get(qs_env%input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=qs_env%input, &
print_key_path="DFT%PRINT%MOMENTS", extension=".dat", &
CALL section_vals_val_get(ec_section, "PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=ec_section, &
print_key_path="PRINT%MOMENTS", extension=".dat", &
middle_name="moments", log_filename=.FALSE.)
IF (iounit > 0) THEN
@ -2122,13 +2170,85 @@ CONTAINS
END IF
CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
basis_section=qs_env%input, print_key_path="DFT%PRINT%MOMENTS")
basis_section=ec_section, print_key_path="PRINT%MOMENTS")
CALL get_qs_env(qs_env=qs_env, results=results)
description = "[DIPOLE]"
CALL cp_results_erase(results=results, description=description)
CALL put_results(results=results, description=description, values=tdip(1:3))
END IF
! Do a Voronoi Integration or write a compressed BQB File
print_key_voro => section_vals_get_subs_vals(ec_section, "PRINT%VORONOI")
print_key_bqb => section_vals_get_subs_vals(ec_section, "PRINT%E_DENSITY_BQB")
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key_voro), cp_p_file)) THEN
should_print_voro = 1
ELSE
should_print_voro = 0
END IF
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key_bqb), cp_p_file)) THEN
should_print_bqb = 1
ELSE
should_print_bqb = 0
END IF
IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0)) THEN
CALL get_qs_env(qs_env=qs_env, &
pw_env=pw_env)
CALL pw_env_get(pw_env=pw_env, &
auxbas_pw_pool=auxbas_pw_pool, &
pw_pools=pw_pools)
CALL pw_pool_create_pw(pool=auxbas_pw_pool, &
pw=rho_elec_rspace%pw, &
use_data=REALDATA3D, &
in_space=REALSPACE)
IF (dft_control%nspins > 1) THEN
! add Pout and Pz
CALL pw_copy(ec_env%rhoout_r(1)%pw, rho_elec_rspace%pw)
CALL pw_axpy(ec_env%rhoout_r(2)%pw, rho_elec_rspace%pw)
CALL pw_axpy(ec_env%rhoz_r(1)%pw, rho_elec_rspace%pw)
CALL pw_axpy(ec_env%rhoz_r(2)%pw, rho_elec_rspace%pw)
mb_rho => rho_elec_rspace%pw
ELSE
! add Pout and Pz
CALL pw_copy(ec_env%rhoout_r(1)%pw, rho_elec_rspace%pw)
CALL pw_axpy(ec_env%rhoz_r(1)%pw, rho_elec_rspace%pw)
mb_rho => rho_elec_rspace%pw
END IF ! nspins
IF (should_print_voro /= 0) THEN
CALL section_vals_val_get(print_key_voro, "OUTPUT_TEXT", l_val=voro_print_txt)
IF (voro_print_txt) THEN
append_voro = section_get_lval(ec_section, "PRINT%VORONOI%APPEND")
my_pos_voro = "REWIND"
IF (append_voro) THEN
my_pos_voro = "APPEND"
END IF
unit_nr_voro = cp_print_key_unit_nr(logger, ec_section, "PRINT%VORONOI", extension=".voronoi", &
file_position=my_pos_voro, log_filename=.FALSE.)
ELSE
unit_nr_voro = 0
END IF
ELSE
unit_nr_voro = 0
END IF
CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
unit_nr_voro, qs_env, mb_rho)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_elec_rspace%pw)
IF (unit_nr_voro > 0) THEN
CALL cp_print_key_finished_output(unit_nr_voro, logger, ec_section, "PRINT%VORONOI")
END IF
END IF
CALL timestop(handle)
END SUBROUTINE ec_properties

222
src/input_cp2k_ec.F
View file

@ -13,15 +13,26 @@
! **************************************************************************************************
MODULE input_cp2k_ec
USE bibliography, ONLY: Niklasson2003,&
VandeVondele2012
VandeVondele2012,&
brehm2018,&
brehm2020,&
brehm2021,&
rycroft2009,&
thomas2015
USE cp_output_handling, ONLY: cp_print_key_section_create,&
debug_print_level,&
high_print_level
USE input_constants, ONLY: &
bqb_opt_exhaustive, bqb_opt_normal, bqb_opt_off, bqb_opt_patient, bqb_opt_quick, &
ec_diagonalization, ec_functional_harris, ec_ls_solver, ec_matrix_sign, ec_matrix_tc2, &
ec_matrix_trs4, ec_mo_solver, ec_ot_atomic, ec_ot_diag, ec_ot_gs, kg_cholesky, &
ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, &
ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, &
ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_sign_ns, &
ls_scf_sign_proot, ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, precond_mlp
ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, precond_mlp, &
voro_radii_cov, voro_radii_unity, voro_radii_user, voro_radii_vdw
USE input_cp2k_mm, ONLY: create_dipoles_section
USE input_cp2k_xc, ONLY: create_xc_section
USE input_keyword_types, ONLY: keyword_create,&
keyword_release,&
@ -31,7 +42,9 @@ MODULE input_cp2k_ec
section_create,&
section_release,&
section_type
USE input_val_types, ONLY: char_t
USE input_val_types, ONLY: char_t,&
integer_t,&
real_t
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
#include "./base/base_uses.f90"
@ -83,6 +96,12 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
! add a print section for properties
NULLIFY (subsection)
CALL create_ec_print_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_FUNCTIONAL", &
description="Functional used in energy correction", &
usage="ENERGY_FUNCTIONAL HARRIS", &
@ -490,4 +509,201 @@ CONTAINS
END SUBROUTINE create_ec_solver_section
! **************************************************************************************************
!> \brief Create the print dft section
!> \param section the section to create
!> \author fbelle - from create_print_dft_section
! **************************************************************************************************
SUBROUTINE create_ec_print_section(section)
TYPE(section_type), POINTER :: section
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: print_key
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="PRINT", &
description="Section of possible print options in EC code.", &
n_keywords=0, n_subsections=1, repeats=.FALSE.)
NULLIFY (print_key, keyword)
! Output of BQB volumetric files
CALL cp_print_key_section_create(print_key, __LOCATION__, name="E_DENSITY_BQB", &
description="Controls the output of the electron density to the losslessly"// &
" compressed BQB file format, see [Brehm2018] (via LibBQB"// &
" see <a href=""https://brehm-research.de/bqb"" target=""_blank"">"// &
" https://brehm-research.de/bqb</a> ). Currently does not work"// &
" with changing cell vector (NpT ensemble).", &
print_level=debug_print_level + 1, filename="", &
citations=(/Brehm2018/))
CALL keyword_create(keyword, __LOCATION__, name="SKIP_FIRST", &
description="Skips the first step of a MD run (avoids duplicate step if restarted).", &
usage="SKIP_FIRST T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="STORE_STEP_NUMBER", &
description="Stores the step number and simulation time in the comment line of each BQB"// &
" frame. Switch it off for binary compatibility with original CP2k CUBE files.", &
usage="STORE_STEP_NUMBER F", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CHECK", &
description="Performs an on-the-fly decompression of each compressed BQB frame to check"// &
" whether the volumetric data exactly matches, and aborts the run if not so.", &
usage="CHECK T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OVERWRITE", &
description="Specify this keyword to overwrite the output BQB file if"// &
" it already exists. By default, the data is appended to an existing file.", &
usage="OVERWRITE T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="HISTORY", &
description="Controls how many previous steps are taken into account for extrapolation in"// &
" compression. Use a value of 1 to compress the frames independently.", &
usage="HISTORY 10", n_var=1, default_i_val=10, type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="PARAMETER_KEY", &
description="Allows to supply previously optimized compression parameters via a"// &
" parameter key (alphanumeric character sequence starting with 'at')."// &
" Just leave away the 'at' sign here, because CP2k will otherwise"// &
" assume it is a variable name in the input", &
usage="PARAMETER_KEY <KEY>", n_var=1, default_c_val="", type_of_var=char_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OPTIMIZE", &
description="Controls the time spent to optimize the parameters for compression efficiency.", &
usage="OPTIMIZE {OFF,QUICK,NORMAL,PATIENT,EXHAUSTIVE}", repeats=.FALSE., n_var=1, &
default_i_val=bqb_opt_quick, &
enum_c_vals=s2a("OFF", "QUICK", "NORMAL", "PATIENT", "EXHAUSTIVE"), &
enum_desc=s2a("No optimization (use defaults)", "Quick optimization", &
"Standard optimization", "Precise optimization", "Exhaustive optimization"), &
enum_i_vals=(/bqb_opt_off, bqb_opt_quick, bqb_opt_normal, bqb_opt_patient, bqb_opt_exhaustive/))
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Voronoi Integration via LibVori
CALL cp_print_key_section_create(print_key, __LOCATION__, name="VORONOI", &
description="Controls the Voronoi integration of the total electron density"// &
" for the computation of electromagnetic moments, see [Thomas2015],"// &
" [Brehm2020], and [Brehm2021] (via LibVori,"// &
" see <a href=""https://brehm-research.de/voronoi"" target=""_blank"">"// &
" https://brehm-research.de/voronoi</a> ).", &
print_level=debug_print_level + 1, filename="", &
citations=(/Rycroft2009, Thomas2015, Brehm2018, Brehm2020, Brehm2021/))
CALL keyword_create(keyword, __LOCATION__, name="APPEND", &
description="Appends frames to already existing .voronoi file.", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SANITY_CHECK", &
description="Performs a sanity check before each Voronoi integration, i.e.,"// &
" checks if every grid point is located in exactly one Voronoi cell.", &
usage="SANITY_CHECK T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OVERWRITE", &
description="Specify this keyword to overwrite any existing ""properties.emp"" file if"// &
" it already exists. By default, the data is appended to an existing .emp file.", &
usage="OVERWRITE T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SKIP_FIRST", &
description="Skips the first step of a MD run (avoids duplicate step if restarted).", &
usage="SKIP_FIRST T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VERBOSE", &
description="Switches on verbose screen output of the Voronoi integration.", &
usage="VERBOSE T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OUTPUT_EMP", &
description="Writes the resulting electromagnetic moments to a binary file ""properties.emp""."// &
" The file name cannot be changed.", &
usage="OUTPUT_EMP T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OUTPUT_TEXT", &
description="Writes the resulting electromagnetic moments to text files (*.voronoi)."// &
" The file name is specified via FILENAME.", &
usage="OUTPUT_TEXT T", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="REFINEMENT_FACTOR", &
description="Sets the refinement factor for the Voronoi integration.", &
usage="REFINEMENT 2", n_var=1, default_i_val=1, type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VORONOI_RADII", &
description="Which atomic radii to use for the radical Voronoi tessellation.", &
usage="VORONOI_RADII {UNITY,VDW,COVALENT,USER}", repeats=.FALSE., n_var=1, &
default_i_val=voro_radii_vdw, &
enum_c_vals=s2a("UNITY", "VDW", "COVALENT", "USER"), &
enum_desc=s2a("Use unity radii (non-radical Voronoi tessellation)", "Use VdW atom radii", &
"Use covalent atom radii", "Use user-specified atom radii"), &
enum_i_vals=(/voro_radii_unity, voro_radii_vdw, voro_radii_cov, voro_radii_user/))
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="USER_RADII", &
description="Defines user atom radii for the radical Voronoi tessellation (one per atom).", &
usage="USER_RADII {real} {real} {real}", repeats=.FALSE., &
unit_str="angstrom", &
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
!Printing of Moments
CALL create_dipoles_section(print_key, "MOMENTS", high_print_level)
CALL keyword_create( &
keyword, __LOCATION__, &
name="MAX_MOMENT", &
description="Maximum moment to be calculated. Values higher than 1 not implemented under periodic boundaries.", &
usage="MAX_MOMENT {integer}", &
repeats=.FALSE., &
n_var=1, &
type_of_var=integer_t, &
default_i_val=1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, &
name="MAGNETIC", &
description="Calculate also magnetic moments, only implemented without periodic boundaries", &
usage="MAGNETIC yes", &
repeats=.FALSE., &
n_var=1, &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
END SUBROUTINE create_ec_print_section
END MODULE input_cp2k_ec

12
src/response_solver.F
View file

@ -80,6 +80,7 @@ MODULE response_solver
pw_env_type
USE pw_grid_types, ONLY: pw_grid_type
USE pw_methods, ONLY: pw_axpy,&
pw_copy,&
pw_integral_ab,&
pw_scale,&
pw_transfer,&
@ -785,18 +786,20 @@ CONTAINS
!> \param matrix_wz Energy-weighted linear response density
!> \param zehartree Hartree volume response contribution to stress tensor
!> \param zexc XC volume response contribution to stress tensor
!> \param rhopz_r Response density on real space grid
!> \param p_env ...
!> \param ex_env ...
! **************************************************************************************************
SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, &
zehartree, zexc, p_env, ex_env)
zehartree, zexc, rhopz_r, p_env, ex_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(pw_p_type), POINTER :: vh_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_pz, matrix_pz_admm, &
matrix_wz
REAL(KIND=dp), OPTIONAL :: zehartree, zexc
TYPE(pw_p_type), DIMENSION(:), OPTIONAL, POINTER :: rhopz_r
TYPE(qs_p_env_type), OPTIONAL, POINTER :: p_env
TYPE(excited_energy_type), OPTIONAL, POINTER :: ex_env
@ -1199,6 +1202,13 @@ CONTAINS
CALL pw_axpy(rhoz_g(ispin)%pw, rhoz_tot_gspace%pw)
END DO
! Save response density on real space grid for properties
IF (PRESENT(rhopz_r)) THEN
DO ispin = 1, nspins
CALL pw_copy(rhoz_r(ispin)%pw, rhopz_r(ispin)%pw)
END DO
END IF
! Stress-tensor contribution second derivative
! Volume : int v_H[n^z]*n_in
! Volume : int epsilon_xc*n_z

12
tests/QS/regtest-debug-5/h2o_dip26.inp
View file

@ -16,6 +16,12 @@
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -34,12 +40,6 @@
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip27.inp
View file

@ -19,6 +19,12 @@
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -37,12 +43,6 @@
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip28.inp
View file

@ -16,6 +16,12 @@
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -36,12 +42,6 @@
&END
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip29.inp
View file

@ -22,6 +22,12 @@
EPS 1.E-5
EPS_FILTER 1.E-6
&END
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -40,12 +46,6 @@
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip30.inp
View file

@ -22,6 +22,12 @@
EPS 1.E-5
EPS_FILTER 1.E-6
&END
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -40,12 +46,6 @@
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip31.inp
View file

@ -22,6 +22,12 @@
EPS 1.E-5
EPS_FILTER 1.E-6
&END
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -45,12 +51,6 @@
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip32.inp
View file

@ -19,6 +19,12 @@
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -37,12 +43,6 @@
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip33.inp
View file

@ -22,6 +22,12 @@
EPS 1.E-5
EPS_FILTER 1.E-6
&END
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -40,12 +46,6 @@
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-debug-5/h2o_dip34.inp
View file

@ -30,6 +30,12 @@
EPS 1.E-5
EPS_FILTER 1.E-6
&END
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
SCF_GUESS ATOMIC
@ -48,12 +54,6 @@
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL

12
tests/QS/regtest-ec/HF_ec_dipole.inp
View file

@ -1,12 +1,6 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
@ -24,6 +18,12 @@
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-7

12
tests/QS/regtest-ec/HF_ec_field.inp
View file

@ -5,12 +5,6 @@
INTENSITY 0.001
POLARISATION 0.0 0.0 1.0
&END
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
@ -28,6 +22,12 @@
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-7

66
tests/QS/regtest-ec/HF_ec_voronoi.inp Normal file
View file

@ -0,0 +1,66 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 400
&END MGRID
&QS
EPS_DEFAULT 1.E-14
&END QS
&ENERGY_CORRECTION
ENERGY_FUNCTIONAL HARRIS
HARRIS_BASIS ORBITAL
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&VORONOI
SANITY_CHECK T
VERBOSE T
REFINEMENT_FACTOR 1
&END
&END
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-7
SCF_GUESS ATOMIC
&END
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
H 0.000000 0.000000 0.000000
F 0.000000 0.000000 1.050000
&END COORD
&KIND H
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q1
&END KIND
&KIND F
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q7
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT HF
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL

4
tests/QS/regtest-ec/N2_ec_tc2.inp
View file

@ -1,10 +1,6 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
&PRINT
&DERIVATIVES
&END
&END
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID

2
tests/QS/regtest-ec/TEST_FILES
View file

@ -13,6 +13,8 @@ H2_H2O_ecprim.inp 11 1e-10
2H2O_ecmao2.inp 11 1e-08 -34.4437114917
# ENERGY : LS + KG embed + Harris (Diag) dipole
HF_ec_dipole.inp 11 1e-08 -24.8917060368
# ENERGY : LS + KG embed + Harris (Diag) dipole
HF_ec_voronoi.inp 11 1e-08 -24.8917060368
# ENERGY : LS + KG embed + Harris (Diag) efield
HF_ec_field.inp 11 1e-08 -24.8908706031
# ENERGY: Test Harris functional solvers