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