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