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 9247 : 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 9247 : CPASSERT(.NOT. ASSOCIATED(ec_env))
113 120211 : ALLOCATE (ec_env)
114 9247 : CALL init_ec_env(qs_env, ec_env, dft_section, ec_section)
115 :
116 9247 : 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 9247 : 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 9247 : 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 9247 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_set
151 : TYPE(mp_para_env_type), POINTER :: para_env
152 9247 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
153 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
154 9247 : 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 9247 : CALL timeset(routineN, handle)
162 :
163 9247 : NULLIFY (atomic_kind_set, dispersion_env, ec_env%ls_env, para_env)
164 9247 : NULLIFY (ec_env%sab_orb, ec_env%sac_ae, ec_env%sac_ppl, ec_env%sap_ppnl)
165 9247 : NULLIFY (ec_env%matrix_ks, ec_env%matrix_h, ec_env%matrix_s)
166 9247 : NULLIFY (ec_env%matrix_t, ec_env%matrix_p, ec_env%matrix_w)
167 9247 : NULLIFY (ec_env%task_list)
168 9247 : NULLIFY (ec_env%mao_coef)
169 9247 : NULLIFY (ec_env%force)
170 9247 : NULLIFY (ec_env%dispersion_env)
171 9247 : NULLIFY (ec_env%xc_section)
172 9247 : NULLIFY (ec_env%matrix_z)
173 9247 : NULLIFY (ec_env%matrix_hz)
174 9247 : NULLIFY (ec_env%matrix_wz)
175 9247 : NULLIFY (ec_env%z_admm)
176 9247 : NULLIFY (ec_env%p_env)
177 9247 : NULLIFY (ec_env%vxc_rspace)
178 9247 : NULLIFY (ec_env%vtau_rspace)
179 9247 : NULLIFY (ec_env%vadmm_rspace)
180 9247 : NULLIFY (ec_env%vadmm_tau_rspace)
181 9247 : NULLIFY (ec_env%rhoout_r, ec_env%rhoz_r)
182 9247 : NULLIFY (ec_env%x_data)
183 9247 : ec_env%should_update = .TRUE.
184 9247 : ec_env%mao = .FALSE.
185 9247 : ec_env%do_ec_admm = .FALSE.
186 9247 : ec_env%do_kpoints = .FALSE.
187 9247 : ec_env%do_ec_hfx = .FALSE.
188 9247 : ec_env%reuse_hfx = .FALSE.
189 :
190 9247 : IF (qs_env%energy_correction) THEN
191 :
192 306 : CPASSERT(PRESENT(ec_section))
193 :
194 : ! get a useful output_unit
195 306 : logger => cp_get_default_logger()
196 306 : IF (logger%para_env%is_source()) THEN
197 153 : 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 306 : i_val=ec_env%ks_solver)
204 : CALL section_vals_val_get(ec_section, "ENERGY_FUNCTIONAL", &
205 306 : i_val=ec_env%energy_functional)
206 : CALL section_vals_val_get(ec_section, "FACTORIZATION", &
207 306 : i_val=ec_env%factorization)
208 : CALL section_vals_val_get(ec_section, "OT_INITIAL_GUESS", &
209 306 : i_val=ec_env%ec_initial_guess)
210 : CALL section_vals_val_get(ec_section, "EPS_DEFAULT", &
211 306 : r_val=ec_env%eps_default)
212 : CALL section_vals_val_get(ec_section, "HARRIS_BASIS", &
213 306 : c_val=ec_env%basis)
214 : CALL section_vals_val_get(ec_section, "ELECTRONIC_TEMPERATURE", &
215 306 : r_val=etemp)
216 : CALL section_vals_val_get(ec_section, "MAO", &
217 306 : l_val=ec_env%mao)
218 : CALL section_vals_val_get(ec_section, "MAO_MAX_ITER", &
219 306 : i_val=ec_env%mao_max_iter)
220 : CALL section_vals_val_get(ec_section, "MAO_EPS_GRAD", &
221 306 : r_val=ec_env%mao_eps_grad)
222 : CALL section_vals_val_get(ec_section, "MAO_EPS1", &
223 306 : r_val=ec_env%mao_eps1)
224 : CALL section_vals_val_get(ec_section, "MAO_IOLEVEL", &
225 306 : 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 306 : l_val=ec_env%skip_ec)
229 : ! Debug output
230 : CALL section_vals_val_get(ec_section, "DEBUG_FORCES", &
231 306 : l_val=ec_env%debug_forces)
232 : CALL section_vals_val_get(ec_section, "DEBUG_STRESS", &
233 306 : l_val=ec_env%debug_stress)
234 : CALL section_vals_val_get(ec_section, "DEBUG_EXTERNAL_METHOD", &
235 306 : l_val=ec_env%debug_external)
236 : ! WFN output
237 306 : section1 => section_vals_get_subs_vals(ec_section, "PRINT%HARRIS_OUTPUT_WFN")
238 306 : CALL section_vals_get(section1, explicit=ec_env%write_harris_wfn)
239 : ! ADMM
240 306 : 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 306 : c_val=ec_env%exresp_fn)
244 : CALL section_vals_val_get(ec_section, "EXTERNAL_RESPONSE_ERROR_FILENAME", &
245 306 : c_val=ec_env%exresperr_fn)
246 : CALL section_vals_val_get(ec_section, "EXTERNAL_RESULT_FILENAME", &
247 306 : c_val=ec_env%exresult_fn)
248 : CALL section_vals_val_get(ec_section, "ERROR_ESTIMATION", &
249 306 : l_val=ec_env%do_error)
250 : CALL section_vals_val_get(ec_section, "ERROR_ESTIMATION_METHOD", &
251 306 : c_val=ec_env%error_method)
252 306 : CALL uppercase(ec_env%error_method)
253 : CALL section_vals_val_get(ec_section, "ERROR_CUTOFF", &
254 306 : r_val=ec_env%error_cutoff)
255 : CALL section_vals_val_get(ec_section, "ERROR_SUBSPACE_SIZE", &
256 306 : i_val=ec_env%error_subspace)
257 :
258 306 : ec_env%do_skip = .FALSE.
259 :
260 : ! smearing
261 306 : 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 306 : ec_env%smear%do_smear = .FALSE.
269 306 : 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 306 : CALL get_qs_env(qs_env, do_kpoints=gs_kpoints)
279 306 : kp_section => section_vals_get_subs_vals(ec_section, "KPOINTS")
280 306 : CALL section_vals_get(kp_section, explicit=explicit)
281 306 : ec_env%do_kpoints = gs_kpoints .OR. explicit
282 306 : IF (ec_env%do_kpoints) THEN
283 14 : IF (.NOT. explicit) THEN
284 0 : kp_section => section_vals_get_subs_vals(dft_section, "KPOINTS")
285 : END IF
286 14 : CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
287 14 : CALL kpoint_create(ec_env%kpoints)
288 14 : CALL read_kpoint_section(ec_env%kpoints, kp_section, cell%hmat, cell)
289 14 : CALL kpoint_initialize(ec_env%kpoints, particle_set, cell)
290 14 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
291 14 : CALL kpoint_env_initialize(ec_env%kpoints, para_env, blacs_env)
292 : ELSE
293 292 : NULLIFY (ec_env%kpoints)
294 : END IF
295 :
296 : ! set basis
297 306 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
298 306 : CALL uppercase(ec_env%basis)
299 518 : SELECT CASE (ec_env%basis)
300 : CASE ("ORBITAL")
301 470 : DO ikind = 1, nkind
302 258 : qs_kind => qs_kind_set(ikind)
303 258 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
304 470 : IF (ASSOCIATED(basis_set)) THEN
305 258 : NULLIFY (harris_basis)
306 258 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
307 258 : IF (ASSOCIATED(harris_basis)) THEN
308 6 : CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
309 : END IF
310 258 : NULLIFY (harris_basis)
311 258 : CALL copy_gto_basis_set(basis_set, harris_basis)
312 258 : 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 224 : DO ikind = 1, nkind
336 132 : qs_kind => qs_kind_set(ikind)
337 132 : NULLIFY (harris_basis)
338 132 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
339 224 : 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 306 : CPABORT("Unknown basis set for energy correction (Harris functional)")
345 : END SELECT
346 : !
347 306 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="HARRIS")
348 306 : CALL init_orbital_pointers(maxlgto + 1)
349 : ! GAPW: Generate soft version of Harris basis
350 306 : CALL get_qs_env(qs_env, dft_control=dft_control)
351 306 : 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 510 : CALL add_basis_set_to_container(qs_kind%basis_sets, harris_soft_basis, "HARRIS_SOFT")
366 : END DO
367 : END IF
368 : !
369 306 : CALL uppercase(ec_env%basis)
370 :
371 : ! Basis may only differ from ground-state if explicitly added
372 306 : ec_env%basis_inconsistent = .FALSE.
373 306 : IF (ec_env%basis == "HARRIS") THEN
374 224 : DO ikind = 1, nkind
375 132 : qs_kind => qs_kind_set(ikind)
376 : ! Basis sets of ground-state
377 132 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
378 : ! Basis sets of energy correction
379 132 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
380 :
381 224 : IF (basis_set%name /= harris_basis%name) THEN
382 70 : 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 306 : 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 306 : 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 468 : SELECT CASE (ec_env%energy_functional)
401 : CASE (ec_functional_harris)
402 162 : ec_env%ec_name = "Harris"
403 : CASE (ec_functional_dc)
404 128 : ec_env%ec_name = "DC-DFT"
405 : CASE (ec_functional_ext)
406 16 : ec_env%ec_name = "External Energy"
407 : CASE DEFAULT
408 306 : CPABORT("unknown energy correction")
409 : END SELECT
410 : ! select the XC section
411 306 : NULLIFY (xc_section)
412 306 : xc_section => section_vals_get_subs_vals(dft_section, "XC")
413 306 : section1 => section_vals_get_subs_vals(ec_section, "XC")
414 306 : section2 => section_vals_get_subs_vals(ec_section, "XC%XC_FUNCTIONAL")
415 306 : CALL section_vals_get(section2, explicit=explicit)
416 306 : IF (explicit) THEN
417 290 : CALL xc_functionals_expand(section2, section1)
418 290 : 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 306 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
424 306 : 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 306 : 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 306 : (section_get_ival(xc_section, "XC_GRID%XC_DERIV") == xc_deriv_collocate))
431 : ! dispersion
432 1530 : ALLOCATE (dispersion_env)
433 306 : NULLIFY (xc_section)
434 306 : xc_section => ec_env%xc_section
435 306 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, para_env=para_env)
436 306 : CALL qs_dispersion_env_set(dispersion_env, xc_section)
437 306 : 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 306 : ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
442 4 : NULLIFY (nl_section)
443 4 : nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
444 4 : CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
445 : END IF
446 306 : ec_env%dispersion_env => dispersion_env
447 :
448 : ! Check if hybrid functional are used
449 306 : ec_hfx_section => section_vals_get_subs_vals(ec_section, "XC%HF")
450 306 : CALL section_vals_get(ec_hfx_section, explicit=ec_env%do_ec_hfx)
451 :
452 : ! Initialize Harris LS solver environment
453 306 : ec_env%use_ls_solver = .FALSE.
454 : ec_env%use_ls_solver = (ec_env%ks_solver == ec_matrix_sign) &
455 : .OR. (ec_env%ks_solver == ec_matrix_trs4) &
456 306 : .OR. (ec_env%ks_solver == ec_matrix_tc2)
457 :
458 306 : IF (ec_env%use_ls_solver) THEN
459 22 : CALL ec_ls_create(qs_env, ec_env)
460 : END IF
461 :
462 : ! check that Harris functional with electronic temperature uses diagonalization
463 306 : IF (ec_env%energy_functional == ec_functional_harris) THEN
464 162 : IF (ec_env%smear%do_smear .AND. ec_env%ks_solver /= ec_diagonalization) THEN
465 0 : CPABORT("Harris functional with Fermi-Dirac smearing needs diagonalization solver.")
466 : END IF
467 162 : IF (ec_env%do_kpoints .AND. ec_env%ks_solver /= ec_diagonalization) THEN
468 0 : CPABORT("Harris functional with K-points needs diagonalization solver.")
469 : END IF
470 : END IF
471 :
472 : ! initialize Kpoint MOs
473 306 : IF (ec_env%do_kpoints) THEN
474 14 : CALL get_qs_env(qs_env, dft_control=dft_control)
475 14 : nspins = dft_control%nspins
476 14 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf, basis_type="HARRIS")
477 14 : focc = 2.0_dp + REAL(1 - nspins, dp)
478 14 : flexible_electron_count = dft_control%relax_multiplicity
479 14 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
480 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
481 14 : nrow_global=nsgf, ncol_global=nsgf)
482 14 : CALL get_qs_env(qs_env, nelectron_spin=nel)
483 56 : ALLOCATE (mo_set(nspins))
484 28 : DO ispin = 1, nspins
485 14 : n_el_f = nel(ispin)
486 : CALL allocate_mo_set(mo_set(ispin), nsgf, nsgf, nel(ispin), n_el_f, &
487 14 : focc, flexible_electron_count)
488 28 : CALL init_mo_set(mo_set(ispin), fm_struct=fm_struct, name="MO")
489 : END DO
490 14 : CALL kpoint_initialize_mos(ec_env%kpoints, mo_set)
491 14 : CALL kpoint_initialize_mo_set(ec_env%kpoints)
492 14 : CALL qs_basis_rotation(qs_env, ec_env%kpoints, basis_type="HARRIS")
493 28 : DO ispin = 1, nspins
494 28 : CALL deallocate_mo_set(mo_set(ispin))
495 : END DO
496 14 : DEALLOCATE (mo_set)
497 28 : CALL cp_fm_struct_release(fm_struct)
498 : END IF
499 :
500 : END IF
501 :
502 9247 : CALL timestop(handle)
503 :
504 9247 : END SUBROUTINE init_ec_env
505 :
506 : ! **************************************************************************************************
507 : !> \brief Initializes linear scaling environment for LS based solver of
508 : !> Harris energy functional and parses input section
509 : !> \param qs_env ...
510 : !> \param ec_env ...
511 : !> \par History
512 : !> 2020.10 created [Fabian Belleflamme]
513 : !> \author Fabian Belleflamme
514 : ! **************************************************************************************************
515 22 : SUBROUTINE ec_ls_create(qs_env, ec_env)
516 : TYPE(qs_environment_type), POINTER :: qs_env
517 : TYPE(energy_correction_type), POINTER :: ec_env
518 :
519 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ls_create'
520 :
521 : INTEGER :: handle
522 : REAL(KIND=dp) :: mu
523 : TYPE(dft_control_type), POINTER :: dft_control
524 : TYPE(ls_scf_env_type), POINTER :: ls_env
525 22 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
526 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
527 : TYPE(section_vals_type), POINTER :: ec_section, input
528 :
529 22 : CALL timeset(routineN, handle)
530 :
531 858 : ALLOCATE (ec_env%ls_env)
532 22 : ls_env => ec_env%ls_env
533 :
534 22 : NULLIFY (dft_control, input, ls_env%para_env)
535 :
536 : CALL get_qs_env(qs_env, &
537 : dft_control=dft_control, &
538 : input=input, &
539 : molecule_set=molecule_set, &
540 : particle_set=particle_set, &
541 : para_env=ls_env%para_env, &
542 22 : nelectron_spin=ls_env%nelectron_spin)
543 :
544 : ! copy some basic stuff
545 22 : ls_env%nspins = dft_control%nspins
546 22 : ls_env%natoms = SIZE(particle_set, 1)
547 22 : CALL ls_env%para_env%retain()
548 :
549 : ! initialize block to group to defined molecules
550 66 : ALLOCATE (ls_env%ls_mstruct%atom_to_molecule(ls_env%natoms))
551 22 : CALL molecule_of_atom(molecule_set, atom_to_mol=ls_env%ls_mstruct%atom_to_molecule)
552 :
553 22 : ls_env%do_transport = .FALSE.
554 22 : ls_env%do_pao = .FALSE.
555 22 : ls_env%ls_mstruct%do_pao = ls_env%do_pao
556 22 : ls_env%do_pexsi = .FALSE.
557 22 : ls_env%has_unit_metric = .FALSE.
558 :
559 22 : ec_section => section_vals_get_subs_vals(input, "DFT%ENERGY_CORRECTION")
560 22 : CALL section_vals_val_get(ec_section, "EPS_FILTER", r_val=ls_env%eps_filter)
561 22 : CALL section_vals_val_get(ec_section, "MU", r_val=mu)
562 22 : CALL section_vals_val_get(ec_section, "FIXED_MU", l_val=ls_env%fixed_mu)
563 66 : ls_env%mu_spin = mu
564 22 : CALL section_vals_val_get(ec_section, "S_PRECONDITIONER", i_val=ls_env%s_preconditioner_type)
565 22 : CALL section_vals_val_get(ec_section, "MATRIX_CLUSTER_TYPE", i_val=ls_env%ls_mstruct%cluster_type)
566 22 : CALL section_vals_val_get(ec_section, "S_INVERSION", i_val=ls_env%s_inversion_type)
567 22 : CALL section_vals_val_get(ec_section, "CHECK_S_INV", l_val=ls_env%check_s_inv)
568 22 : CALL section_vals_val_get(ec_section, "REPORT_ALL_SPARSITIES", l_val=ls_env%report_all_sparsities)
569 22 : CALL section_vals_val_get(ec_section, "SIGN_METHOD", i_val=ls_env%sign_method)
570 22 : CALL section_vals_val_get(ec_section, "SIGN_ORDER", i_val=ls_env%sign_order)
571 22 : CALL section_vals_val_get(ec_section, "DYNAMIC_THRESHOLD", l_val=ls_env%dynamic_threshold)
572 22 : CALL section_vals_val_get(ec_section, "NON_MONOTONIC", l_val=ls_env%non_monotonic)
573 22 : CALL section_vals_val_get(ec_section, "S_SQRT_METHOD", i_val=ls_env%s_sqrt_method)
574 22 : CALL section_vals_val_get(ec_section, "S_SQRT_ORDER", i_val=ls_env%s_sqrt_order)
575 22 : CALL section_vals_val_get(ec_section, "EPS_LANCZOS", r_val=ls_env%eps_lanczos)
576 22 : CALL section_vals_val_get(ec_section, "MAX_ITER_LANCZOS", i_val=ls_env%max_iter_lanczos)
577 :
578 24 : SELECT CASE (ec_env%ks_solver)
579 : CASE (ec_matrix_sign)
580 : ! S inverse required for Sign matrix algorithm,
581 : ! calculated either by Hotelling or multiplying S matrix sqrt inv
582 24 : SELECT CASE (ls_env%s_inversion_type)
583 : CASE (ls_s_inversion_sign_sqrt)
584 2 : ls_env%needs_s_inv = .TRUE.
585 2 : ls_env%use_s_sqrt = .TRUE.
586 : CASE (ls_s_inversion_hotelling)
587 0 : ls_env%needs_s_inv = .TRUE.
588 0 : ls_env%use_s_sqrt = .FALSE.
589 : CASE (ls_s_inversion_none)
590 0 : ls_env%needs_s_inv = .FALSE.
591 0 : ls_env%use_s_sqrt = .FALSE.
592 : CASE DEFAULT
593 2 : CPABORT("Unknown ls_env%s_inversion_type")
594 : END SELECT
595 : CASE (ec_matrix_trs4, ec_matrix_tc2)
596 20 : ls_env%needs_s_inv = .FALSE.
597 20 : ls_env%use_s_sqrt = .TRUE.
598 : CASE DEFAULT
599 22 : CPABORT("Unknown ec_env%ks_solver")
600 : END SELECT
601 :
602 22 : SELECT CASE (ls_env%s_preconditioner_type)
603 : CASE (ls_s_preconditioner_none)
604 0 : ls_env%has_s_preconditioner = .FALSE.
605 : CASE DEFAULT
606 22 : ls_env%has_s_preconditioner = .TRUE.
607 : END SELECT
608 :
609 : ! buffer for the history of matrices, not needed here
610 22 : ls_env%extrapolation_order = 0
611 22 : ls_env%scf_history%nstore = 0
612 22 : ls_env%scf_history%istore = 0
613 44 : ALLOCATE (ls_env%scf_history%matrix(ls_env%nspins, ls_env%scf_history%nstore))
614 :
615 22 : NULLIFY (ls_env%mixing_store)
616 :
617 22 : CALL timestop(handle)
618 :
619 44 : END SUBROUTINE ec_ls_create
620 :
621 : ! **************************************************************************************************
622 : !> \brief Print out the energy correction input section
623 : !>
624 : !> \param ec_env ...
625 : !> \par History
626 : !> 2020.10 created [Fabian Belleflamme]
627 : !> \author Fabian Belleflamme
628 : ! **************************************************************************************************
629 306 : SUBROUTINE ec_write_input(ec_env)
630 : TYPE(energy_correction_type), POINTER :: ec_env
631 :
632 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_write_input'
633 :
634 : INTEGER :: handle, unit_nr
635 : TYPE(cp_logger_type), POINTER :: logger
636 : TYPE(ls_scf_env_type), POINTER :: ls_env
637 :
638 306 : CALL timeset(routineN, handle)
639 :
640 306 : logger => cp_get_default_logger()
641 306 : IF (logger%para_env%is_source()) THEN
642 153 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
643 : ELSE
644 153 : unit_nr = -1
645 : END IF
646 :
647 306 : IF (unit_nr > 0) THEN
648 :
649 : WRITE (unit_nr, '(T2,A)') &
650 153 : "!"//REPEAT("-", 29)//" Energy Correction "//REPEAT("-", 29)//"!"
651 :
652 : ! Type of energy correction
653 234 : SELECT CASE (ec_env%energy_functional)
654 : CASE (ec_functional_harris)
655 81 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "HARRIS FUNCTIONAL"
656 : CASE (ec_functional_dc)
657 64 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "DC-DFT"
658 : CASE (ec_functional_ext)
659 153 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "External"
660 : END SELECT
661 153 : WRITE (unit_nr, '()')
662 :
663 : ! Energy correction parameters
664 153 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_default:", ec_env%eps_default
665 :
666 153 : CALL uppercase(ec_env%basis)
667 259 : SELECT CASE (ec_env%basis)
668 : CASE ("ORBITAL")
669 106 : WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "ORBITAL"
670 : CASE ("PRIMITIVE")
671 1 : WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "PRIMITIVE"
672 : CASE ("HARRIS")
673 153 : WRITE (unit_nr, '(T2,A,T61,A20)') "EC Basis: ", "HARRIS"
674 : END SELECT
675 :
676 : ! Info how HFX in energy correction is treated
677 153 : IF (ec_env%do_ec_hfx) THEN
678 :
679 14 : WRITE (unit_nr, '(T2,A,T61,L20)') "DC-DFT with HFX", ec_env%do_ec_hfx
680 14 : WRITE (unit_nr, '(T2,A,T61,L20)') "Reuse HFX integrals", ec_env%reuse_hfx
681 14 : WRITE (unit_nr, '(T2,A,T61,L20)') "DC-DFT HFX with ADMM", ec_env%do_ec_admm
682 :
683 : END IF ! ec_env%do_ec_hfx
684 :
685 153 : IF (ec_env%do_kpoints) THEN
686 7 : CALL write_kpoint_info(ec_env%kpoints, iounit=unit_nr)
687 : END IF
688 :
689 : ! Parameters for Harris functional solver
690 153 : IF (ec_env%energy_functional == ec_functional_harris) THEN
691 :
692 : ! Algorithm
693 149 : SELECT CASE (ec_env%ks_solver)
694 : CASE (ec_diagonalization)
695 68 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "DIAGONALIZATION"
696 : CASE (ec_ot_diag)
697 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "OT DIAGONALIZATION"
698 : CASE (ec_matrix_sign)
699 1 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "MATRIX_SIGN"
700 : CASE (ec_matrix_trs4)
701 9 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TRS4"
702 9 : CALL cite_reference(Niklasson2003)
703 : CASE (ec_matrix_tc2)
704 1 : WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TC2"
705 82 : CALL cite_reference(Niklasson2014)
706 : END SELECT
707 81 : WRITE (unit_nr, '()')
708 :
709 : ! MAO
710 81 : IF (ec_env%mao) THEN
711 2 : WRITE (unit_nr, '(T2,A,T61,L20)') "MAO:", ec_env%mao
712 2 : WRITE (unit_nr, '(T2,A,T61,L20)') "MAO_IOLEVEL:", ec_env%mao_iolevel
713 2 : WRITE (unit_nr, '(T2,A,T61,I20)') "MAO_MAX_ITER:", ec_env%mao_max_iter
714 2 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS_GRAD:", ec_env%mao_eps_grad
715 2 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS1:", ec_env%mao_eps1
716 2 : WRITE (unit_nr, '()')
717 : END IF
718 :
719 : ! Parameters for linear response solver
720 81 : IF (.NOT. ec_env%use_ls_solver) THEN
721 :
722 70 : WRITE (unit_nr, '(T2,A)') "MO Solver"
723 70 : WRITE (unit_nr, '()')
724 :
725 138 : SELECT CASE (ec_env%ks_solver)
726 : CASE (ec_diagonalization)
727 :
728 138 : SELECT CASE (ec_env%factorization)
729 : CASE (kg_cholesky)
730 68 : WRITE (unit_nr, '(T2,A,T61,A20)') "Factorization: ", "CHOLESKY"
731 : END SELECT
732 :
733 : CASE (ec_ot_diag)
734 :
735 : ! OT Diagonalization
736 : ! Initial guess : 1) block diagonal initial guess
737 : ! 2) GS-density matrix (might require trafo if basis diff)
738 :
739 3 : SELECT CASE (ec_env%ec_initial_guess)
740 : CASE (ec_ot_atomic)
741 1 : WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "ATOMIC"
742 : CASE (ec_ot_gs)
743 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "GROUND STATE DM"
744 : END SELECT
745 :
746 : CASE DEFAULT
747 70 : CPABORT("Unknown Diagonalization algorithm for Harris functional")
748 : END SELECT
749 :
750 : ELSE
751 :
752 11 : WRITE (unit_nr, '(T2,A)') "AO Solver"
753 11 : WRITE (unit_nr, '()')
754 :
755 11 : ls_env => ec_env%ls_env
756 11 : WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", ls_env%eps_filter
757 11 : WRITE (unit_nr, '(T2,A,T61,L20)') "fixed chemical potential (mu)", ls_env%fixed_mu
758 11 : WRITE (unit_nr, '(T2,A,T61,L20)') "Computing inv(S):", ls_env%needs_s_inv
759 11 : WRITE (unit_nr, '(T2,A,T61,L20)') "Computing sqrt(S):", ls_env%use_s_sqrt
760 11 : WRITE (unit_nr, '(T2,A,T61,L20)') "Computing S preconditioner ", ls_env%has_s_preconditioner
761 :
762 11 : IF (ls_env%use_s_sqrt) THEN
763 21 : SELECT CASE (ls_env%s_sqrt_method)
764 : CASE (ls_s_sqrt_ns)
765 10 : WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ"
766 : CASE (ls_s_sqrt_proot)
767 1 : WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT"
768 : CASE DEFAULT
769 11 : CPABORT("Unknown sqrt method.")
770 : END SELECT
771 11 : WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", ls_env%s_sqrt_order
772 : END IF
773 :
774 11 : SELECT CASE (ls_env%s_preconditioner_type)
775 : CASE (ls_s_preconditioner_none)
776 0 : WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "NONE"
777 : CASE (ls_s_preconditioner_atomic)
778 11 : WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "ATOMIC"
779 : CASE (ls_s_preconditioner_molecular)
780 11 : WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "MOLECULAR"
781 : END SELECT
782 :
783 22 : SELECT CASE (ls_env%ls_mstruct%cluster_type)
784 : CASE (ls_cluster_atomic)
785 11 : WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("ATOMIC")
786 : CASE (ls_cluster_molecular)
787 0 : WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("MOLECULAR")
788 : CASE DEFAULT
789 11 : CPABORT("Unknown cluster type")
790 : END SELECT
791 :
792 : END IF
793 :
794 : END IF ! if ec_functional_harris
795 :
796 153 : WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
797 153 : WRITE (unit_nr, '()')
798 :
799 : END IF ! unit_nr
800 :
801 306 : CALL timestop(handle)
802 :
803 306 : END SUBROUTINE ec_write_input
804 :
805 : END MODULE ec_environment
|
