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