diff --git a/src/ec_env_types.F b/src/ec_env_types.F
index 11f0a5042e..02ad3a4507 100644
--- a/src/ec_env_types.F
+++ b/src/ec_env_types.F
@@ -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
diff --git a/src/ec_environment.F b/src/ec_environment.F
index ee575e6a03..04233c3d9f 100644
--- a/src/ec_environment.F
+++ b/src/ec_environment.F
@@ -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
! **************************************************************************************************
diff --git a/src/energy_corrections.F b/src/energy_corrections.F
index 4b17a6a101..e76a6b5f85 100644
--- a/src/energy_corrections.F
+++ b/src/energy_corrections.F
@@ -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
diff --git a/src/input_cp2k_ec.F b/src/input_cp2k_ec.F
index 457d8f069d..52a013a23c 100644
--- a/src/input_cp2k_ec.F
+++ b/src/input_cp2k_ec.F
@@ -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 "// &
+ " https://brehm-research.de/bqb ). 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 ", 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 "// &
+ " https://brehm-research.de/voronoi ).", &
+ 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
diff --git a/src/response_solver.F b/src/response_solver.F
index 17e6a01831..d5a30817e9 100644
--- a/src/response_solver.F
+++ b/src/response_solver.F
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip26.inp b/tests/QS/regtest-debug-5/h2o_dip26.inp
index 3167542d4f..4ea46531f5 100644
--- a/tests/QS/regtest-debug-5/h2o_dip26.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip26.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip27.inp b/tests/QS/regtest-debug-5/h2o_dip27.inp
index f86d9080c4..606fb7fb14 100644
--- a/tests/QS/regtest-debug-5/h2o_dip27.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip27.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip28.inp b/tests/QS/regtest-debug-5/h2o_dip28.inp
index 6eb2e1501b..864941c176 100644
--- a/tests/QS/regtest-debug-5/h2o_dip28.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip28.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip29.inp b/tests/QS/regtest-debug-5/h2o_dip29.inp
index ad05f2f1b2..6f42d56698 100644
--- a/tests/QS/regtest-debug-5/h2o_dip29.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip29.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip30.inp b/tests/QS/regtest-debug-5/h2o_dip30.inp
index 2cebe7e102..8499144e83 100644
--- a/tests/QS/regtest-debug-5/h2o_dip30.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip30.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip31.inp b/tests/QS/regtest-debug-5/h2o_dip31.inp
index 68ed125d49..589fafe124 100644
--- a/tests/QS/regtest-debug-5/h2o_dip31.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip31.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip32.inp b/tests/QS/regtest-debug-5/h2o_dip32.inp
index 13ac4d159e..506487d615 100644
--- a/tests/QS/regtest-debug-5/h2o_dip32.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip32.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip33.inp b/tests/QS/regtest-debug-5/h2o_dip33.inp
index 74f7863aea..e5d4874665 100644
--- a/tests/QS/regtest-debug-5/h2o_dip33.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip33.inp
@@ -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
diff --git a/tests/QS/regtest-debug-5/h2o_dip34.inp b/tests/QS/regtest-debug-5/h2o_dip34.inp
index 7e4ec08700..4d47b13971 100644
--- a/tests/QS/regtest-debug-5/h2o_dip34.inp
+++ b/tests/QS/regtest-debug-5/h2o_dip34.inp
@@ -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
diff --git a/tests/QS/regtest-ec/HF_ec_dipole.inp b/tests/QS/regtest-ec/HF_ec_dipole.inp
index f30048edea..e3af5627be 100644
--- a/tests/QS/regtest-ec/HF_ec_dipole.inp
+++ b/tests/QS/regtest-ec/HF_ec_dipole.inp
@@ -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
diff --git a/tests/QS/regtest-ec/HF_ec_field.inp b/tests/QS/regtest-ec/HF_ec_field.inp
index fd85fde093..48f09a711d 100644
--- a/tests/QS/regtest-ec/HF_ec_field.inp
+++ b/tests/QS/regtest-ec/HF_ec_field.inp
@@ -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
diff --git a/tests/QS/regtest-ec/HF_ec_voronoi.inp b/tests/QS/regtest-ec/HF_ec_voronoi.inp
new file mode 100644
index 0000000000..a602e83192
--- /dev/null
+++ b/tests/QS/regtest-ec/HF_ec_voronoi.inp
@@ -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
diff --git a/tests/QS/regtest-ec/N2_ec_tc2.inp b/tests/QS/regtest-ec/N2_ec_tc2.inp
index e924760feb..5d4629fede 100644
--- a/tests/QS/regtest-ec/N2_ec_tc2.inp
+++ b/tests/QS/regtest-ec/N2_ec_tc2.inp
@@ -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
diff --git a/tests/QS/regtest-ec/TEST_FILES b/tests/QS/regtest-ec/TEST_FILES
index 078808a9ef..cc70815a42 100644
--- a/tests/QS/regtest-ec/TEST_FILES
+++ b/tests/QS/regtest-ec/TEST_FILES
@@ -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