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