Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Energy correction environment setup and handling
10 : !> \par History
11 : !> 2019.09 created
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE ec_environment
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE basis_set_container_types, ONLY: add_basis_set_to_container,&
17 : remove_basis_from_container
18 : USE basis_set_types, ONLY: allocate_gto_basis_set,&
19 : copy_gto_basis_set,&
20 : create_primitive_basis_set,&
21 : gto_basis_set_type
22 : USE bibliography, ONLY: Niklasson2003,&
23 : Niklasson2014,&
24 : cite_reference
25 : USE cell_types, ONLY: cell_type
26 : USE cp_blacs_env, ONLY: cp_blacs_env_type
27 : USE cp_control_types, ONLY: dft_control_type
28 : USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr
29 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
30 : USE ec_env_types, ONLY: energy_correction_type
31 : USE input_constants, ONLY: &
32 : ec_diagonalization, ec_functional_dc, ec_functional_ext, ec_functional_harris, &
33 : ec_matrix_sign, ec_matrix_tc2, ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, &
34 : kg_cholesky, ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, &
35 : ls_s_inversion_none, ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, &
36 : ls_s_preconditioner_molecular, ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, &
37 : smear_fermi_dirac, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
38 : USE input_cp2k_check, ONLY: xc_functionals_expand
39 : USE input_section_types, ONLY: section_get_ival,&
40 : section_vals_get,&
41 : section_vals_get_subs_vals,&
42 : section_vals_type,&
43 : section_vals_val_get
44 : USE kinds, ONLY: dp
45 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
46 : kpoint_initialize
47 : USE kpoint_types, ONLY: kpoint_create,&
48 : read_kpoint_section,&
49 : write_kpoint_info
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE molecule_types, ONLY: molecule_of_atom,&
52 : molecule_type
53 : USE orbital_pointers, ONLY: init_orbital_pointers
54 : USE particle_types, ONLY: particle_type
55 : USE qs_dispersion_nonloc, ONLY: qs_dispersion_nonloc_init
56 : USE qs_dispersion_pairpot, ONLY: qs_dispersion_pairpot_init
57 : USE qs_dispersion_types, ONLY: qs_dispersion_type
58 : USE qs_dispersion_utils, ONLY: qs_dispersion_env_set
59 : USE qs_environment_types, ONLY: get_qs_env,&
60 : qs_environment_type
61 : USE qs_interactions, ONLY: init_interaction_radii_orb_basis
62 : USE qs_kind_types, ONLY: get_qs_kind,&
63 : get_qs_kind_set,&
64 : qs_kind_type
65 : USE qs_rho_types, ONLY: qs_rho_type
66 : USE soft_basis_set, ONLY: create_soft_basis
67 : USE string_utilities, ONLY: uppercase
68 : USE xc, ONLY: xc_uses_kinetic_energy_density,&
69 : xc_uses_norm_drho
70 : USE xc_input_constants, ONLY: xc_deriv_collocate
71 : #include "./base/base_uses.f90"
72 :
73 : IMPLICIT NONE
74 :
75 : PRIVATE
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_environment'
78 :
79 : PUBLIC :: ec_env_create
80 : PUBLIC :: ec_write_input
81 :
82 : CONTAINS
83 :
84 : ! **************************************************************************************************
85 : !> \brief Allocates and intitializes ec_env
86 : !> \param qs_env The QS environment
87 : !> \param ec_env The energy correction environment (the object to create)
88 : !> \param dft_section The DFT section
89 : !> \param ec_section The energy correction input section
90 : !> \par History
91 : !> 2019.09 created
92 : !> \author JGH
93 : ! **************************************************************************************************
94 7704 : SUBROUTINE ec_env_create(qs_env, ec_env, dft_section, ec_section)
95 : TYPE(qs_environment_type), POINTER :: qs_env
96 : TYPE(energy_correction_type), POINTER :: ec_env
97 : TYPE(section_vals_type), POINTER :: dft_section
98 : TYPE(section_vals_type), OPTIONAL, POINTER :: ec_section
99 :
100 7704 : CPASSERT(.NOT. ASSOCIATED(ec_env))
101 100152 : ALLOCATE (ec_env)
102 7704 : CALL init_ec_env(qs_env, ec_env, dft_section, ec_section)
103 :
104 7704 : END SUBROUTINE ec_env_create
105 :
106 : ! **************************************************************************************************
107 : !> \brief Initializes energy correction environment
108 : !> \param qs_env The QS environment
109 : !> \param ec_env The energy correction environment
110 : !> \param dft_section The DFT section
111 : !> \param ec_section The energy correction input section
112 : !> \par History
113 : !> 2019.09 created
114 : !> \author JGH
115 : ! **************************************************************************************************
116 7704 : SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section)
117 : TYPE(qs_environment_type), POINTER :: qs_env
118 : TYPE(energy_correction_type), POINTER :: ec_env
119 : TYPE(section_vals_type), POINTER :: dft_section
120 : TYPE(section_vals_type), OPTIONAL, POINTER :: ec_section
121 :
122 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_ec_env'
123 :
124 : INTEGER :: handle, ikind, maxlgto, nkind, unit_nr
125 : LOGICAL :: explicit, gpw, gs_kpoints, paw_atom
126 : REAL(KIND=dp) :: eps_pgf_orb, etemp, rc
127 7704 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
128 : TYPE(cell_type), POINTER :: cell
129 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
130 : TYPE(dft_control_type), POINTER :: dft_control
131 : TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis, &
132 : harris_soft_basis
133 : TYPE(mp_para_env_type), POINTER :: para_env
134 7704 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
135 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
136 7704 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
137 : TYPE(qs_kind_type), POINTER :: qs_kind
138 : TYPE(qs_rho_type), POINTER :: rho
139 : TYPE(section_vals_type), POINTER :: ec_hfx_section, kp_section, nl_section, &
140 : pp_section, section1, section2, &
141 : xc_fun_section, xc_section
142 :
143 7704 : CALL timeset(routineN, handle)
144 :
145 7704 : NULLIFY (atomic_kind_set, dispersion_env, ec_env%ls_env, para_env)
146 7704 : NULLIFY (ec_env%sab_orb, ec_env%sac_ae, ec_env%sac_ppl, ec_env%sap_ppnl)
147 7704 : NULLIFY (ec_env%matrix_ks, ec_env%matrix_h, ec_env%matrix_s)
148 7704 : NULLIFY (ec_env%matrix_t, ec_env%matrix_p, ec_env%matrix_w)
149 7704 : NULLIFY (ec_env%task_list)
150 7704 : NULLIFY (ec_env%mao_coef)
151 7704 : NULLIFY (ec_env%force)
152 7704 : NULLIFY (ec_env%dispersion_env)
153 7704 : NULLIFY (ec_env%xc_section)
154 7704 : NULLIFY (ec_env%matrix_z)
155 7704 : NULLIFY (ec_env%matrix_hz)
156 7704 : NULLIFY (ec_env%matrix_wz)
157 7704 : NULLIFY (ec_env%z_admm)
158 7704 : NULLIFY (ec_env%p_env)
159 7704 : NULLIFY (ec_env%vxc_rspace)
160 7704 : NULLIFY (ec_env%vtau_rspace)
161 7704 : NULLIFY (ec_env%vadmm_rspace)
162 7704 : NULLIFY (ec_env%rhoout_r, ec_env%rhoz_r)
163 7704 : NULLIFY (ec_env%x_data)
164 7704 : ec_env%should_update = .TRUE.
165 7704 : ec_env%mao = .FALSE.
166 7704 : ec_env%do_ec_admm = .FALSE.
167 7704 : ec_env%do_kpoints = .FALSE.
168 7704 : ec_env%do_ec_hfx = .FALSE.
169 7704 : ec_env%reuse_hfx = .FALSE.
170 :
171 7704 : IF (qs_env%energy_correction) THEN
172 :
173 286 : CPASSERT(PRESENT(ec_section))
174 :
175 : ! get a useful output_unit
176 286 : unit_nr = cp_logger_get_default_unit_nr(local=.FALSE.)
177 :
178 : CALL section_vals_val_get(ec_section, "ALGORITHM", &
179 286 : i_val=ec_env%ks_solver)
180 : CALL section_vals_val_get(ec_section, "ENERGY_FUNCTIONAL", &
181 286 : i_val=ec_env%energy_functional)
182 : CALL section_vals_val_get(ec_section, "FACTORIZATION", &
183 286 : i_val=ec_env%factorization)
184 : CALL section_vals_val_get(ec_section, "OT_INITIAL_GUESS", &
185 286 : i_val=ec_env%ec_initial_guess)
186 : CALL section_vals_val_get(ec_section, "EPS_DEFAULT", &
187 286 : r_val=ec_env%eps_default)
188 : CALL section_vals_val_get(ec_section, "HARRIS_BASIS", &
189 286 : c_val=ec_env%basis)
190 : CALL section_vals_val_get(ec_section, "ELECTRONIC_TEMPERATURE", &
191 286 : r_val=etemp)
192 : CALL section_vals_val_get(ec_section, "MAO", &
193 286 : l_val=ec_env%mao)
194 : CALL section_vals_val_get(ec_section, "MAO_MAX_ITER", &
195 286 : i_val=ec_env%mao_max_iter)
196 : CALL section_vals_val_get(ec_section, "MAO_EPS_GRAD", &
197 286 : r_val=ec_env%mao_eps_grad)
198 : CALL section_vals_val_get(ec_section, "MAO_EPS1", &
199 286 : r_val=ec_env%mao_eps1)
200 : CALL section_vals_val_get(ec_section, "MAO_IOLEVEL", &
201 286 : i_val=ec_env%mao_iolevel)
202 : ! Skip EC calculation if ground-state calculation did not converge
203 : CALL section_vals_val_get(ec_section, "SKIP_EC", &
204 286 : l_val=ec_env%skip_ec)
205 : ! Debug output
206 : CALL section_vals_val_get(ec_section, "DEBUG_FORCES", &
207 286 : l_val=ec_env%debug_forces)
208 : CALL section_vals_val_get(ec_section, "DEBUG_STRESS", &
209 286 : l_val=ec_env%debug_stress)
210 : CALL section_vals_val_get(ec_section, "DEBUG_EXTERNAL_METHOD", &
211 286 : l_val=ec_env%debug_external)
212 : ! ADMM
213 286 : CALL section_vals_val_get(ec_section, "ADMM", l_val=ec_env%do_ec_admm)
214 : ! EXTERNAL
215 : CALL section_vals_val_get(ec_section, "EXTERNAL_RESPONSE_FILENAME", &
216 286 : c_val=ec_env%exresp_fn)
217 : CALL section_vals_val_get(ec_section, "EXTERNAL_RESPONSE_ERROR_FILENAME", &
218 286 : c_val=ec_env%exresperr_fn)
219 : CALL section_vals_val_get(ec_section, "EXTERNAL_RESULT_FILENAME", &
220 286 : c_val=ec_env%exresult_fn)
221 : CALL section_vals_val_get(ec_section, "ERROR_ESTIMATION", &
222 286 : l_val=ec_env%do_error)
223 : CALL section_vals_val_get(ec_section, "ERROR_ESTIMATION_METHOD", &
224 286 : c_val=ec_env%error_method)
225 286 : CALL uppercase(ec_env%error_method)
226 : CALL section_vals_val_get(ec_section, "ERROR_CUTOFF", &
227 286 : r_val=ec_env%error_cutoff)
228 : CALL section_vals_val_get(ec_section, "ERROR_SUBSPACE_SIZE", &
229 286 : i_val=ec_env%error_subspace)
230 :
231 286 : ec_env%do_skip = .FALSE.
232 :
233 : ! smearing
234 286 : IF (etemp > 0.0_dp) THEN
235 0 : ec_env%smear%do_smear = .TRUE.
236 0 : ec_env%smear%method = smear_fermi_dirac
237 0 : ec_env%smear%electronic_temperature = etemp
238 0 : ec_env%smear%eps_fermi_dirac = 1.0E-5_dp
239 0 : ec_env%smear%fixed_mag_mom = -100.0_dp
240 : ELSE
241 286 : ec_env%smear%do_smear = .FALSE.
242 286 : ec_env%smear%electronic_temperature = 0.0_dp
243 : END IF
244 :
245 : ! kpoints
246 : ! Options: gs ec
247 : ! gamma gamma
248 : ! gamma KPec
249 : ! KPgs KPgs
250 : ! KPgs KPec
251 286 : CALL get_qs_env(qs_env, do_kpoints=gs_kpoints)
252 286 : kp_section => section_vals_get_subs_vals(ec_section, "KPOINTS")
253 286 : CALL section_vals_get(kp_section, explicit=explicit)
254 286 : ec_env%do_kpoints = gs_kpoints .OR. explicit
255 286 : IF (ec_env%do_kpoints) THEN
256 0 : IF (.NOT. explicit) THEN
257 0 : kp_section => section_vals_get_subs_vals(dft_section, "KPOINTS")
258 : END IF
259 0 : CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
260 0 : CALL kpoint_create(ec_env%kpoints)
261 0 : CALL read_kpoint_section(ec_env%kpoints, kp_section, cell%hmat)
262 0 : CALL kpoint_initialize(ec_env%kpoints, particle_set, cell)
263 0 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
264 0 : CALL kpoint_env_initialize(ec_env%kpoints, para_env, blacs_env)
265 : ELSE
266 286 : NULLIFY (ec_env%kpoints)
267 : END IF
268 :
269 : ! set basis
270 286 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
271 286 : CALL uppercase(ec_env%basis)
272 484 : SELECT CASE (ec_env%basis)
273 : CASE ("ORBITAL")
274 434 : DO ikind = 1, nkind
275 236 : qs_kind => qs_kind_set(ikind)
276 236 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
277 434 : IF (ASSOCIATED(basis_set)) THEN
278 236 : NULLIFY (harris_basis)
279 236 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
280 236 : IF (ASSOCIATED(harris_basis)) THEN
281 6 : CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
282 : END IF
283 236 : NULLIFY (harris_basis)
284 236 : CALL copy_gto_basis_set(basis_set, harris_basis)
285 236 : CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
286 : END IF
287 : END DO
288 : CASE ("PRIMITIVE")
289 6 : DO ikind = 1, nkind
290 4 : qs_kind => qs_kind_set(ikind)
291 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
292 6 : IF (ASSOCIATED(basis_set)) THEN
293 4 : NULLIFY (harris_basis)
294 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
295 4 : IF (ASSOCIATED(harris_basis)) THEN
296 0 : CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
297 : END IF
298 4 : NULLIFY (harris_basis)
299 4 : CALL create_primitive_basis_set(basis_set, harris_basis)
300 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
301 4 : eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
302 4 : CALL init_interaction_radii_orb_basis(harris_basis, eps_pgf_orb)
303 4 : harris_basis%kind_radius = basis_set%kind_radius
304 4 : CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
305 : END IF
306 : END DO
307 : CASE ("HARRIS")
308 212 : DO ikind = 1, nkind
309 126 : qs_kind => qs_kind_set(ikind)
310 126 : NULLIFY (harris_basis)
311 126 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
312 212 : IF (.NOT. ASSOCIATED(harris_basis)) THEN
313 0 : CPWARN("Harris Basis not defined for all types of atoms.")
314 : END IF
315 : END DO
316 : CASE DEFAULT
317 286 : CPABORT("Unknown basis set for energy correction (Harris functional)")
318 : END SELECT
319 : !
320 286 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="HARRIS")
321 286 : CALL init_orbital_pointers(maxlgto + 1)
322 : ! GAPW: Generate soft version of Harris basis
323 286 : CALL get_qs_env(qs_env, dft_control=dft_control)
324 286 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
325 50 : eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
326 118 : DO ikind = 1, nkind
327 68 : qs_kind => qs_kind_set(ikind)
328 68 : NULLIFY (harris_basis)
329 68 : CALL get_qs_kind(qs_kind, basis_set=harris_basis, basis_type="HARRIS")
330 68 : CALL get_qs_kind(qs_kind, hard_radius=rc, gpw_type_forced=gpw)
331 68 : NULLIFY (harris_soft_basis)
332 68 : CALL allocate_gto_basis_set(harris_soft_basis)
333 : CALL create_soft_basis(harris_basis, harris_soft_basis, &
334 : dft_control%qs_control%gapw_control%eps_fit, &
335 : rc, paw_atom, &
336 68 : dft_control%qs_control%gapw_control%force_paw, gpw)
337 68 : CALL init_interaction_radii_orb_basis(harris_soft_basis, eps_pgf_orb)
338 490 : CALL add_basis_set_to_container(qs_kind%basis_sets, harris_soft_basis, "HARRIS_SOFT")
339 : END DO
340 : END IF
341 : !
342 286 : CALL uppercase(ec_env%basis)
343 :
344 : ! Basis may only differ from ground-state if explicitly added
345 286 : ec_env%basis_inconsistent = .FALSE.
346 286 : IF (ec_env%basis == "HARRIS") THEN
347 212 : DO ikind = 1, nkind
348 126 : qs_kind => qs_kind_set(ikind)
349 : ! Basis sets of ground-state
350 126 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
351 : ! Basis sets of energy correction
352 126 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
353 :
354 212 : IF (basis_set%name /= harris_basis%name) THEN
355 64 : ec_env%basis_inconsistent = .TRUE.
356 : END IF
357 : END DO
358 : END IF
359 :
360 : !Density-corrected DFT must be performed with the same basis as ground-state
361 286 : IF (ec_env%energy_functional == ec_functional_dc .AND. ec_env%basis_inconsistent) THEN
362 : CALL cp_abort(__LOCATION__, &
363 : "DC-DFT: Correction and ground state need to use the same basis. "// &
364 0 : "Checked by comparing basis set names only.")
365 : END IF
366 286 : IF (ec_env%energy_functional == ec_functional_ext .AND. ec_env%basis_inconsistent) THEN
367 : CALL cp_abort(__LOCATION__, &
368 : "Exteranl Energy: Correction and ground state need to use the same basis. "// &
369 0 : "Checked by comparing basis set names only.")
370 : END IF
371 : !
372 : ! set functional
373 432 : SELECT CASE (ec_env%energy_functional)
374 : CASE (ec_functional_harris)
375 146 : ec_env%ec_name = "Harris"
376 : CASE (ec_functional_dc)
377 124 : ec_env%ec_name = "DC-DFT"
378 : CASE (ec_functional_ext)
379 16 : ec_env%ec_name = "External Energy"
380 : CASE DEFAULT
381 286 : CPABORT("unknown energy correction")
382 : END SELECT
383 : ! select the XC section
384 286 : NULLIFY (xc_section)
385 286 : xc_section => section_vals_get_subs_vals(dft_section, "XC")
386 286 : section1 => section_vals_get_subs_vals(ec_section, "XC")
387 286 : section2 => section_vals_get_subs_vals(ec_section, "XC%XC_FUNCTIONAL")
388 286 : CALL section_vals_get(section2, explicit=explicit)
389 286 : IF (explicit) THEN
390 270 : CALL xc_functionals_expand(section2, section1)
391 270 : ec_env%xc_section => section1
392 : ELSE
393 16 : ec_env%xc_section => xc_section
394 : END IF
395 : ! Check whether energy correction requires the kinetic energy density and rebuild rho if necessary
396 286 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
397 286 : xc_fun_section => section_vals_get_subs_vals(ec_env%xc_section, "XC_FUNCTIONAL")
398 : dft_control%use_kinetic_energy_density = dft_control%use_kinetic_energy_density .OR. &
399 286 : xc_uses_kinetic_energy_density(xc_fun_section, dft_control%lsd)
400 : ! Same for density gradient
401 : dft_control%drho_by_collocation = dft_control%drho_by_collocation .OR. &
402 : (xc_uses_norm_drho(xc_fun_section, dft_control%lsd) .AND. &
403 286 : (section_get_ival(xc_section, "XC_GRID%XC_DERIV") == xc_deriv_collocate))
404 : ! dispersion
405 1430 : ALLOCATE (dispersion_env)
406 : NULLIFY (xc_section)
407 286 : xc_section => ec_env%xc_section
408 286 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, para_env=para_env)
409 286 : CALL qs_dispersion_env_set(dispersion_env, xc_section)
410 286 : IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
411 0 : NULLIFY (pp_section)
412 0 : pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL")
413 0 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
414 286 : ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
415 0 : CPABORT("nl-vdW functionals not available for EC calculations")
416 0 : NULLIFY (nl_section)
417 0 : nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
418 0 : CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
419 : END IF
420 286 : ec_env%dispersion_env => dispersion_env
421 :
422 : ! Check if hybrid functional are used
423 286 : ec_hfx_section => section_vals_get_subs_vals(ec_section, "XC%HF")
424 286 : CALL section_vals_get(ec_hfx_section, explicit=ec_env%do_ec_hfx)
425 :
426 : ! Initialize Harris LS solver environment
427 286 : ec_env%use_ls_solver = .FALSE.
428 : ec_env%use_ls_solver = (ec_env%ks_solver == ec_matrix_sign) &
429 : .OR. (ec_env%ks_solver == ec_matrix_trs4) &
430 286 : .OR. (ec_env%ks_solver == ec_matrix_tc2)
431 :
432 286 : IF (ec_env%use_ls_solver) THEN
433 22 : CALL ec_ls_create(qs_env, ec_env)
434 : END IF
435 :
436 : ! check that Harris functional with electronic temperature uses diagonalization
437 286 : IF (ec_env%energy_functional == ec_functional_harris) THEN
438 146 : IF (ec_env%smear%do_smear .AND. ec_env%ks_solver /= ec_diagonalization) THEN
439 0 : CPABORT("Harris functional with Fermi-Dirac smearing needs diagonalization solver.")
440 : END IF
441 : END IF
442 :
443 : END IF
444 :
445 7704 : CALL timestop(handle)
446 :
447 7704 : END SUBROUTINE init_ec_env
448 :
449 : ! **************************************************************************************************
450 : !> \brief Initializes linear scaling environment for LS based solver of
451 : !> Harris energy functional and parses input section
452 : !> \param qs_env ...
453 : !> \param ec_env ...
454 : !> \par History
455 : !> 2020.10 created [Fabian Belleflamme]
456 : !> \author Fabian Belleflamme
457 : ! **************************************************************************************************
458 22 : SUBROUTINE ec_ls_create(qs_env, ec_env)
459 : TYPE(qs_environment_type), POINTER :: qs_env
460 : TYPE(energy_correction_type), POINTER :: ec_env
461 :
462 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ls_create'
463 :
464 : INTEGER :: handle
465 : REAL(KIND=dp) :: mu
466 : TYPE(dft_control_type), POINTER :: dft_control
467 : TYPE(ls_scf_env_type), POINTER :: ls_env
468 22 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
469 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
470 : TYPE(section_vals_type), POINTER :: ec_section, input
471 :
472 22 : CALL timeset(routineN, handle)
473 :
474 858 : ALLOCATE (ec_env%ls_env)
475 22 : ls_env => ec_env%ls_env
476 :
477 22 : NULLIFY (dft_control, input, ls_env%para_env)
478 :
479 : CALL get_qs_env(qs_env, &
480 : dft_control=dft_control, &
481 : input=input, &
482 : molecule_set=molecule_set, &
483 : particle_set=particle_set, &
484 : para_env=ls_env%para_env, &
485 22 : nelectron_spin=ls_env%nelectron_spin)
486 :
487 : ! copy some basic stuff
488 22 : ls_env%nspins = dft_control%nspins
489 22 : ls_env%natoms = SIZE(particle_set, 1)
490 22 : CALL ls_env%para_env%retain()
491 :
492 : ! initialize block to group to defined molecules
493 66 : ALLOCATE (ls_env%ls_mstruct%atom_to_molecule(ls_env%natoms))
494 22 : CALL molecule_of_atom(molecule_set, atom_to_mol=ls_env%ls_mstruct%atom_to_molecule)
495 :
496 22 : ls_env%do_transport = .FALSE.
497 22 : ls_env%do_pao = .FALSE.
498 22 : ls_env%ls_mstruct%do_pao = ls_env%do_pao
499 22 : ls_env%do_pexsi = .FALSE.
500 22 : ls_env%has_unit_metric = .FALSE.
501 :
502 22 : ec_section => section_vals_get_subs_vals(input, "DFT%ENERGY_CORRECTION")
503 22 : CALL section_vals_val_get(ec_section, "EPS_FILTER", r_val=ls_env%eps_filter)
504 22 : CALL section_vals_val_get(ec_section, "MU", r_val=mu)
505 22 : CALL section_vals_val_get(ec_section, "FIXED_MU", l_val=ls_env%fixed_mu)
506 66 : ls_env%mu_spin = mu
507 22 : CALL section_vals_val_get(ec_section, "S_PRECONDITIONER", i_val=ls_env%s_preconditioner_type)
508 22 : CALL section_vals_val_get(ec_section, "MATRIX_CLUSTER_TYPE", i_val=ls_env%ls_mstruct%cluster_type)
509 22 : CALL section_vals_val_get(ec_section, "S_INVERSION", i_val=ls_env%s_inversion_type)
510 22 : CALL section_vals_val_get(ec_section, "CHECK_S_INV", l_val=ls_env%check_s_inv)
511 22 : CALL section_vals_val_get(ec_section, "REPORT_ALL_SPARSITIES", l_val=ls_env%report_all_sparsities)
512 22 : CALL section_vals_val_get(ec_section, "SIGN_METHOD", i_val=ls_env%sign_method)
513 22 : CALL section_vals_val_get(ec_section, "SIGN_ORDER", i_val=ls_env%sign_order)
514 22 : CALL section_vals_val_get(ec_section, "DYNAMIC_THRESHOLD", l_val=ls_env%dynamic_threshold)
515 22 : CALL section_vals_val_get(ec_section, "NON_MONOTONIC", l_val=ls_env%non_monotonic)
516 22 : CALL section_vals_val_get(ec_section, "S_SQRT_METHOD", i_val=ls_env%s_sqrt_method)
517 22 : CALL section_vals_val_get(ec_section, "S_SQRT_ORDER", i_val=ls_env%s_sqrt_order)
518 22 : CALL section_vals_val_get(ec_section, "EPS_LANCZOS", r_val=ls_env%eps_lanczos)
519 22 : CALL section_vals_val_get(ec_section, "MAX_ITER_LANCZOS", i_val=ls_env%max_iter_lanczos)
520 :
521 24 : SELECT CASE (ec_env%ks_solver)
522 : CASE (ec_matrix_sign)
523 : ! S inverse required for Sign matrix algorithm,
524 : ! calculated either by Hotelling or multiplying S matrix sqrt inv
525 24 : SELECT CASE (ls_env%s_inversion_type)
526 : CASE (ls_s_inversion_sign_sqrt)
527 2 : ls_env%needs_s_inv = .TRUE.
528 2 : ls_env%use_s_sqrt = .TRUE.
529 : CASE (ls_s_inversion_hotelling)
530 0 : ls_env%needs_s_inv = .TRUE.
531 0 : ls_env%use_s_sqrt = .FALSE.
532 : CASE (ls_s_inversion_none)
533 0 : ls_env%needs_s_inv = .FALSE.
534 0 : ls_env%use_s_sqrt = .FALSE.
535 : CASE DEFAULT
536 2 : CPABORT("")
537 : END SELECT
538 : CASE (ec_matrix_trs4, ec_matrix_tc2)
539 20 : ls_env%needs_s_inv = .FALSE.
540 20 : ls_env%use_s_sqrt = .TRUE.
541 : CASE DEFAULT
542 22 : CPABORT("")
543 : END SELECT
544 :
545 22 : SELECT CASE (ls_env%s_preconditioner_type)
546 : CASE (ls_s_preconditioner_none)
547 0 : ls_env%has_s_preconditioner = .FALSE.
548 : CASE DEFAULT
549 22 : ls_env%has_s_preconditioner = .TRUE.
550 : END SELECT
551 :
552 : ! buffer for the history of matrices, not needed here
553 22 : ls_env%extrapolation_order = 0
554 22 : ls_env%scf_history%nstore = 0
555 22 : ls_env%scf_history%istore = 0
556 44 : ALLOCATE (ls_env%scf_history%matrix(ls_env%nspins, ls_env%scf_history%nstore))
557 :
558 22 : NULLIFY (ls_env%mixing_store)
559 :
560 22 : CALL timestop(handle)
561 :
562 44 : END SUBROUTINE ec_ls_create
563 :
564 : ! **************************************************************************************************
565 : !> \brief Print out the energy correction input section
566 : !>
567 : !> \param ec_env ...
568 : !> \par History
569 : !> 2020.10 created [Fabian Belleflamme]
570 : !> \author Fabian Belleflamme
571 : ! **************************************************************************************************
572 286 : SUBROUTINE ec_write_input(ec_env)
573 : TYPE(energy_correction_type), POINTER :: ec_env
574 :
575 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_write_input'
576 :
577 : INTEGER :: handle, unit_nr
578 : TYPE(ls_scf_env_type), POINTER :: ls_env
579 :
580 286 : CALL timeset(routineN, handle)
581 :
582 286 : unit_nr = cp_logger_get_default_unit_nr(local=.FALSE.)
583 :
584 286 : IF (unit_nr > 0) THEN
585 :
586 : WRITE (unit_nr, '(T2,A)') &
587 286 : "!"//REPEAT("-", 29)//" Energy Correction "//REPEAT("-", 29)//"!"
588 :
589 : ! Type of energy correction
590 432 : SELECT CASE (ec_env%energy_functional)
591 : CASE (ec_functional_harris)
592 146 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "HARRIS FUNCTIONAL"
593 : CASE (ec_functional_dc)
594 124 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "DC-DFT"
595 : CASE (ec_functional_ext)
596 286 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "External"
597 : END SELECT
598 286 : WRITE (unit_nr, '()')
599 :
600 : ! Energy correction parameters
601 286 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_default:", ec_env%eps_default
602 :
603 286 : CALL uppercase(ec_env%basis)
604 484 : SELECT CASE (ec_env%basis)
605 : CASE ("ORBITAL")
606 198 : WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "ORBITAL"
607 : CASE ("PRIMITIVE")
608 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "PRIMITIVE"
609 : CASE ("HARRIS")
610 286 : WRITE (unit_nr, '(T2,A,T61,A20)') "EC Basis: ", "HARRIS"
611 : END SELECT
612 :
613 : ! Info how HFX in energy correction is treated
614 286 : IF (ec_env%do_ec_hfx) THEN
615 :
616 28 : WRITE (unit_nr, '(T2,A,T61,L20)') "DC-DFT with HFX", ec_env%do_ec_hfx
617 28 : WRITE (unit_nr, '(T2,A,T61,L20)') "Reuse HFX integrals", ec_env%reuse_hfx
618 28 : WRITE (unit_nr, '(T2,A,T61,L20)') "DC-DFT HFX with ADMM", ec_env%do_ec_admm
619 :
620 : END IF ! ec_env%do_ec_hfx
621 :
622 286 : IF (ec_env%do_kpoints) THEN
623 0 : CALL write_kpoint_info(ec_env%kpoints, iounit=unit_nr)
624 : END IF
625 :
626 : ! Parameters for Harris functional solver
627 286 : IF (ec_env%energy_functional == ec_functional_harris) THEN
628 :
629 : ! Algorithm
630 266 : SELECT CASE (ec_env%ks_solver)
631 : CASE (ec_diagonalization)
632 120 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "DIAGONALIZATION"
633 : CASE (ec_ot_diag)
634 4 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "OT DIAGONALIZATION"
635 : CASE (ec_matrix_sign)
636 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "MATRIX_SIGN"
637 : CASE (ec_matrix_trs4)
638 18 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TRS4"
639 18 : CALL cite_reference(Niklasson2003)
640 : CASE (ec_matrix_tc2)
641 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TC2"
642 148 : CALL cite_reference(Niklasson2014)
643 : END SELECT
644 146 : WRITE (unit_nr, '()')
645 :
646 : ! MAO
647 146 : IF (ec_env%mao) THEN
648 4 : WRITE (unit_nr, '(T2,A,T61,L20)') "MAO:", ec_env%mao
649 4 : WRITE (unit_nr, '(T2,A,T61,L20)') "MAO_IOLEVEL:", ec_env%mao_iolevel
650 4 : WRITE (unit_nr, '(T2,A,T61,I20)') "MAO_MAX_ITER:", ec_env%mao_max_iter
651 4 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS_GRAD:", ec_env%mao_eps_grad
652 4 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS1:", ec_env%mao_eps1
653 4 : WRITE (unit_nr, '()')
654 : END IF
655 :
656 : ! Parameters for linear response solver
657 146 : IF (.NOT. ec_env%use_ls_solver) THEN
658 :
659 124 : WRITE (unit_nr, '(T2,A)') "MO Solver"
660 124 : WRITE (unit_nr, '()')
661 :
662 244 : SELECT CASE (ec_env%ks_solver)
663 : CASE (ec_diagonalization)
664 :
665 244 : SELECT CASE (ec_env%factorization)
666 : CASE (kg_cholesky)
667 120 : WRITE (unit_nr, '(T2,A,T61,A20)') "Factorization: ", "CHOLESKY"
668 : END SELECT
669 :
670 : CASE (ec_ot_diag)
671 :
672 : ! OT Diagonalization
673 : ! Initial guess : 1) block diagonal initial guess
674 : ! 2) GS-density matrix (might require trafo if basis diff)
675 :
676 6 : SELECT CASE (ec_env%ec_initial_guess)
677 : CASE (ec_ot_atomic)
678 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "ATOMIC"
679 : CASE (ec_ot_gs)
680 4 : WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "GROUND STATE DM"
681 : END SELECT
682 :
683 : CASE DEFAULT
684 124 : CPABORT("Unknown Diagonalization algorithm for Harris functional")
685 : END SELECT
686 :
687 : ELSE
688 :
689 22 : WRITE (unit_nr, '(T2,A)') "AO Solver"
690 22 : WRITE (unit_nr, '()')
691 :
692 22 : ls_env => ec_env%ls_env
693 22 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", ls_env%eps_filter
694 22 : WRITE (unit_nr, '(T2,A,T61,L20)') "fixed chemical potential (mu)", ls_env%fixed_mu
695 22 : WRITE (unit_nr, '(T2,A,T61,L20)') "Computing inv(S):", ls_env%needs_s_inv
696 22 : WRITE (unit_nr, '(T2,A,T61,L20)') "Computing sqrt(S):", ls_env%use_s_sqrt
697 22 : WRITE (unit_nr, '(T2,A,T61,L20)') "Computing S preconditioner ", ls_env%has_s_preconditioner
698 :
699 22 : IF (ls_env%use_s_sqrt) THEN
700 42 : SELECT CASE (ls_env%s_sqrt_method)
701 : CASE (ls_s_sqrt_ns)
702 20 : WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ"
703 : CASE (ls_s_sqrt_proot)
704 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT"
705 : CASE DEFAULT
706 22 : CPABORT("Unknown sqrt method.")
707 : END SELECT
708 22 : WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", ls_env%s_sqrt_order
709 : END IF
710 :
711 22 : SELECT CASE (ls_env%s_preconditioner_type)
712 : CASE (ls_s_preconditioner_none)
713 0 : WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "NONE"
714 : CASE (ls_s_preconditioner_atomic)
715 22 : WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "ATOMIC"
716 : CASE (ls_s_preconditioner_molecular)
717 22 : WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "MOLECULAR"
718 : END SELECT
719 :
720 44 : SELECT CASE (ls_env%ls_mstruct%cluster_type)
721 : CASE (ls_cluster_atomic)
722 22 : WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("ATOMIC")
723 : CASE (ls_cluster_molecular)
724 0 : WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("MOLECULAR")
725 : CASE DEFAULT
726 22 : CPABORT("Unknown cluster type")
727 : END SELECT
728 :
729 : END IF
730 :
731 : END IF ! if ec_functional_harris
732 :
733 286 : WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
734 286 : WRITE (unit_nr, '()')
735 :
736 : END IF ! unit_nr
737 :
738 286 : CALL timestop(handle)
739 :
740 286 : END SUBROUTINE ec_write_input
741 :
742 : END MODULE ec_environment
|
