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 : !> \par History
10 : !> - Merged with the Quickstep MODULE method_specification (17.01.2002,MK)
11 : !> - USE statements cleaned, added
12 : !> (25.09.2002,MK)
13 : !> - Added more LSD structure (01.2003,Joost VandeVondele)
14 : !> - New molecule data types introduced (Sep. 2003,MK)
15 : !> - Cleaning; getting rid of pnode (02.10.2003,MK)
16 : !> - Sub-system setup added (08.10.2003,MK)
17 : !> \author MK (18.05.2000)
18 : ! **************************************************************************************************
19 : MODULE qs_environment
20 : USE almo_scf_env_methods, ONLY: almo_scf_env_create
21 : USE atom_kind_orbitals, ONLY: calculate_atomic_relkin
22 : USE atomic_kind_types, ONLY: atomic_kind_type
23 : USE auto_basis, ONLY: create_lri_aux_basis_set,&
24 : create_ri_aux_basis_set
25 : USE basis_set_container_types, ONLY: add_basis_set_to_container
26 : USE basis_set_types, ONLY: basis_sort_zet,&
27 : create_primitive_basis_set,&
28 : deallocate_gto_basis_set,&
29 : gto_basis_set_type
30 : USE bibliography, ONLY: Iannuzzi2006,&
31 : Iannuzzi2007,&
32 : cite_reference,&
33 : cp2kqs2020
34 : USE cell_types, ONLY: cell_type
35 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
36 : cp_blacs_env_release,&
37 : cp_blacs_env_type
38 : USE cp_control_types, ONLY: dft_control_type,&
39 : dftb_control_type,&
40 : gapw_control_type,&
41 : qs_control_type,&
42 : semi_empirical_control_type,&
43 : xtb_control_type
44 : USE cp_control_utils, ONLY: &
45 : read_ddapc_section, read_dft_control, read_mgrid_section, read_qs_section, &
46 : read_rixs_control, read_tddfpt2_control, write_admm_control, write_dft_control, &
47 : write_qs_control
48 : USE cp_ddapc_types, ONLY: cp_ddapc_ewald_create
49 : USE cp_log_handling, ONLY: cp_get_default_logger,&
50 : cp_logger_get_default_io_unit,&
51 : cp_logger_type,&
52 : cp_to_string
53 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
54 : cp_print_key_unit_nr
55 : USE cp_subsys_types, ONLY: cp_subsys_type
56 : USE cp_symmetry, ONLY: write_symmetry
57 : USE distribution_1d_types, ONLY: distribution_1d_release,&
58 : distribution_1d_type
59 : USE distribution_methods, ONLY: distribute_molecules_1d
60 : USE ec_env_types, ONLY: energy_correction_type
61 : USE ec_environment, ONLY: ec_env_create,&
62 : ec_write_input
63 : USE et_coupling_types, ONLY: et_coupling_create
64 : USE ewald_environment_types, ONLY: ewald_env_create,&
65 : ewald_env_get,&
66 : ewald_env_set,&
67 : ewald_environment_type,&
68 : read_ewald_section,&
69 : read_ewald_section_tb
70 : USE ewald_pw_methods, ONLY: ewald_pw_grid_update
71 : USE ewald_pw_types, ONLY: ewald_pw_create,&
72 : ewald_pw_type
73 : USE exstates_types, ONLY: excited_energy_type,&
74 : exstate_create
75 : USE external_potential_types, ONLY: get_potential,&
76 : init_potential,&
77 : set_potential
78 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_create,&
79 : fist_nonbond_env_type
80 : USE gamma, ONLY: init_md_ftable
81 : USE global_types, ONLY: global_environment_type
82 : USE hartree_local_methods, ONLY: init_coulomb_local
83 : USE header, ONLY: dftb_header,&
84 : qs_header,&
85 : se_header,&
86 : tblite_header,&
87 : xtb_header
88 : USE hfx_types, ONLY: compare_hfx_sections,&
89 : hfx_create
90 : USE input_constants, ONLY: &
91 : atomic_guess, cdft_charge_constraint, debug_run, diag_ot, dispersion_d2, dispersion_d3, &
92 : dispersion_d3bj, do_et_ddapc, do_method_am1, do_method_dftb, do_method_gapw, &
93 : do_method_gapw_xc, do_method_gpw, do_method_lrigpw, do_method_mndo, do_method_mndod, &
94 : do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, do_method_pm6fm, &
95 : do_method_pnnl, do_method_rigpw, do_method_rm1, do_method_xtb, do_qmmm_gauss, &
96 : do_qmmm_swave, do_spin_density, general_roks, gfn1xtb, hden_atomic, kg_tnadd_embed_ri, &
97 : linear_response_run, mtlr_atomic_perturbations, mtlr_reference_from_atomic, &
98 : mtlr_reference_from_restart, rel_none, rel_trans_atom, restart_guess, sic_none, &
99 : smear_energy_window, smear_fermi_dirac, smear_list, tblite_scc_mixer_tblite, &
100 : tddfpt_kernel_none, vdw_pairpot_dftd2, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, &
101 : vdw_pairpot_dftd4, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, &
102 : wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, xc_vdw_fun_none, xc_vdw_fun_nonloc, &
103 : xc_vdw_fun_pairpot, xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none
104 : USE input_section_types, ONLY: section_get_ival,&
105 : section_get_ivals,&
106 : section_vals_get,&
107 : section_vals_get_subs_vals,&
108 : section_vals_type,&
109 : section_vals_val_get
110 : USE kg_environment, ONLY: kg_env_create
111 : USE kinds, ONLY: default_string_length,&
112 : dp
113 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
114 : kpoint_initialize,&
115 : kpoint_initialize_mos
116 : USE kpoint_types, ONLY: get_kpoint_info,&
117 : kpoint_create,&
118 : kpoint_reset_initialization,&
119 : kpoint_type,&
120 : read_kpoint_section,&
121 : set_kpoint_info,&
122 : write_kpoint_info
123 : USE lri_environment_init, ONLY: lri_env_basis,&
124 : lri_env_init
125 : USE lri_environment_types, ONLY: lri_environment_type
126 : USE machine, ONLY: m_flush
127 : USE mathconstants, ONLY: pi
128 : USE message_passing, ONLY: mp_para_env_type
129 : USE molecule_kind_types, ONLY: molecule_kind_type,&
130 : write_molecule_kind_set
131 : USE molecule_types, ONLY: molecule_type
132 : USE mp2_setup, ONLY: read_mp2_section
133 : USE mp2_types, ONLY: mp2_env_create,&
134 : mp2_type
135 : USE multipole_types, ONLY: do_multipole_none
136 : USE orbital_pointers, ONLY: init_orbital_pointers
137 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
138 : USE particle_methods, ONLY: write_particle_distances,&
139 : write_qs_particle_coordinates,&
140 : write_structure_data
141 : USE particle_types, ONLY: particle_type
142 : USE physcon, ONLY: kelvin
143 : USE pw_env_types, ONLY: pw_env_type
144 : USE qmmm_types_low, ONLY: qmmm_env_qm_type
145 : USE qs_basis_rotation_methods, ONLY: qs_basis_rotation
146 : USE qs_dftb_parameters, ONLY: qs_dftb_param_init
147 : USE qs_dftb_types, ONLY: qs_dftb_atom_type,&
148 : qs_dftb_pairpot_type
149 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
150 : USE qs_dispersion_nonloc, ONLY: qs_dispersion_nonloc_init
151 : USE qs_dispersion_pairpot, ONLY: qs_dispersion_pairpot_init
152 : USE qs_dispersion_types, ONLY: qs_dispersion_type
153 : USE qs_dispersion_utils, ONLY: qs_dispersion_env_set,&
154 : qs_write_dispersion
155 : USE qs_energy_types, ONLY: allocate_qs_energy,&
156 : qs_energy_type
157 : USE qs_environment_methods, ONLY: qs_env_setup
158 : USE qs_environment_types, ONLY: get_qs_env,&
159 : qs_environment_type,&
160 : set_qs_env
161 : USE qs_force_types, ONLY: qs_force_type
162 : USE qs_gcp_types, ONLY: qs_gcp_type
163 : USE qs_gcp_utils, ONLY: qs_gcp_env_set,&
164 : qs_gcp_init
165 : USE qs_harris_types, ONLY: harris_rhoin_init,&
166 : harris_type
167 : USE qs_harris_utils, ONLY: harris_env_create,&
168 : harris_write_input
169 : USE qs_interactions, ONLY: init_interaction_radii,&
170 : init_se_nlradius,&
171 : write_core_charge_radii,&
172 : write_paw_radii,&
173 : write_pgf_orb_radii,&
174 : write_ppl_radii,&
175 : write_ppnl_radii
176 : USE qs_kind_types, ONLY: &
177 : check_qs_kind_set, get_qs_kind, get_qs_kind_set, init_cneo_basis_set, init_gapw_basis_set, &
178 : init_gapw_nlcc, init_qs_kind_set, qs_kind_type, set_qs_kind, write_gto_basis_sets, &
179 : write_qs_kind_set
180 : USE qs_ks_types, ONLY: qs_ks_env_create,&
181 : qs_ks_env_type,&
182 : set_ks_env
183 : USE qs_local_rho_types, ONLY: local_rho_type
184 : USE qs_mo_types, ONLY: allocate_mo_set,&
185 : mo_set_type
186 : USE qs_rho0_ggrid, ONLY: rho0_s_grid_create
187 : USE qs_rho0_methods, ONLY: init_rho0
188 : USE qs_rho0_types, ONLY: rho0_mpole_type
189 : USE qs_rho_atom_methods, ONLY: init_rho_atom
190 : USE qs_rho_atom_types, ONLY: rho_atom_type
191 : USE qs_subsys_methods, ONLY: qs_subsys_create
192 : USE qs_subsys_types, ONLY: qs_subsys_get,&
193 : qs_subsys_set,&
194 : qs_subsys_type
195 : USE qs_wf_history_methods, ONLY: wfi_create,&
196 : wfi_create_for_kp
197 : USE qs_wf_history_types, ONLY: qs_wf_history_type,&
198 : wfi_release
199 : USE rel_control_types, ONLY: rel_c_create,&
200 : rel_c_read_parameters,&
201 : rel_control_type
202 : USE scf_control_types, ONLY: scf_c_create,&
203 : scf_c_read_parameters,&
204 : scf_c_write_parameters,&
205 : scf_control_type
206 : USE semi_empirical_expns3_methods, ONLY: semi_empirical_expns3_setup
207 : USE semi_empirical_int_arrays, ONLY: init_se_intd_array
208 : USE semi_empirical_mpole_methods, ONLY: nddo_mpole_setup
209 : USE semi_empirical_mpole_types, ONLY: nddo_mpole_type
210 : USE semi_empirical_store_int_types, ONLY: semi_empirical_si_create,&
211 : semi_empirical_si_type
212 : USE semi_empirical_types, ONLY: se_taper_create,&
213 : se_taper_type
214 : USE semi_empirical_utils, ONLY: se_cutoff_compatible
215 : USE skala_gpw_functional, ONLY: xc_section_uses_native_skala_evaluator
216 : USE tblite_interface, ONLY: tb_get_basis,&
217 : tb_init_geometry,&
218 : tb_init_wf,&
219 : tb_set_calculator
220 : USE transport, ONLY: transport_env_create
221 : USE xc_derivatives, ONLY: xc_functionals_get_needs
222 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
223 : USE xtb_parameters, ONLY: init_xtb_basis,&
224 : xtb_parameters_init,&
225 : xtb_parameters_set,&
226 : xtb_spinpol_ext,&
227 : xtb_spinpol_init
228 : USE xtb_potentials, ONLY: xtb_pp_radius
229 : USE xtb_types, ONLY: allocate_xtb_atom_param,&
230 : get_xtb_atom_param,&
231 : set_xtb_atom_param,&
232 : write_xtb_atom_param
233 : #include "./base/base_uses.f90"
234 :
235 : IMPLICIT NONE
236 :
237 : PRIVATE
238 :
239 : ! *** Global parameters ***
240 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_environment'
241 :
242 : ! *** Public subroutines ***
243 : PUBLIC :: qs_init
244 :
245 : CONTAINS
246 :
247 : ! **************************************************************************************************
248 : !> \brief Reject unsupported GAPW low-spin ROKS inputs before basis, grid and HFX setup.
249 : !> \param dft_control DFT and Quickstep controls read from input
250 : !> \param dft_section DFT input section
251 : !> \param subsys_section Subsystem input section, including all atomic kinds
252 : ! **************************************************************************************************
253 9407 : SUBROUTINE check_gapw_low_spin_roks(dft_control, dft_section, subsys_section)
254 : TYPE(dft_control_type), POINTER :: dft_control
255 : TYPE(section_vals_type), POINTER :: dft_section, subsys_section
256 :
257 : INTEGER :: irep, nrep, option
258 : LOGICAL :: active, explicit
259 : TYPE(section_vals_type), POINTER :: section, subsection, xc_section
260 : TYPE(xc_rho_cflags_type) :: needs
261 :
262 9215 : IF (.NOT. dft_control%low_spin_roks) RETURN
263 40 : IF (dft_control%qs_control%gapw_xc) THEN
264 0 : CALL cp_abort(__LOCATION__, "GAPW_XC is not compatible with low spin ROKS method.")
265 : END IF
266 40 : IF (.NOT. dft_control%qs_control%gapw) RETURN
267 :
268 32 : IF (dft_control%qs_control%gapw_control%accurate_xcint) THEN
269 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support GAPW_ACCURATE_XCINT.")
270 : END IF
271 32 : IF (dft_control%do_admm) THEN
272 0 : CALL cp_abort(__LOCATION__, "ADMM not compatible with low spin ROKS method.")
273 : END IF
274 :
275 : ! Only XC and full-AO exchange are reevaluated for the spin configurations.
276 : ! Other spin-dependent corrections would retain the high-spin contribution.
277 32 : IF (dft_control%sic_method_id /= sic_none) THEN
278 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support SIC.")
279 : END IF
280 32 : IF (dft_control%qs_control%s2_restraint) THEN
281 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support S2 restraints.")
282 : END IF
283 32 : section => section_vals_get_subs_vals(subsys_section, "KIND")
284 32 : CALL section_vals_get(section, n_repetition=nrep)
285 64 : DO irep = 1, nrep
286 32 : subsection => section_vals_get_subs_vals(section, "DFT_PLUS_U", i_rep_section=irep)
287 32 : CALL section_vals_val_get(subsection, "_SECTION_PARAMETERS_", l_val=active)
288 64 : IF (active) CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support DFT+U.")
289 : END DO
290 32 : section => section_vals_get_subs_vals(dft_section, "QS%DDAPC_RESTRAINT")
291 32 : CALL section_vals_get(section, explicit=explicit, n_repetition=nrep)
292 32 : IF (explicit) THEN
293 0 : DO irep = 1, nrep
294 0 : CALL section_vals_val_get(section, "TYPE_OF_DENSITY", i_rep_section=irep, i_val=option)
295 0 : IF (option == do_spin_density) THEN
296 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support spin DDAPC restraints.")
297 : END IF
298 : END DO
299 : END IF
300 32 : section => section_vals_get_subs_vals(dft_section, "QS%CDFT%ATOM_GROUP")
301 32 : CALL section_vals_get(section, explicit=explicit, n_repetition=nrep)
302 32 : IF (explicit) THEN
303 0 : DO irep = 1, nrep
304 0 : CALL section_vals_val_get(section, "CONSTRAINT_TYPE", i_rep_section=irep, i_val=option)
305 0 : IF (option /= cdft_charge_constraint) THEN
306 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support spin CDFT constraints.")
307 : END IF
308 : END DO
309 : END IF
310 :
311 32 : xc_section => section_vals_get_subs_vals(dft_section, "XC")
312 32 : section => section_vals_get_subs_vals(xc_section, "HF")
313 32 : CALL section_vals_get(section, explicit=explicit, n_repetition=nrep)
314 32 : IF (explicit) THEN
315 6 : IF (nrep /= 1) CALL cp_abort(__LOCATION__, "GAPW low spin ROKS requires a single HF section.")
316 6 : CALL section_vals_val_get(section, "RI%_SECTION_PARAMETERS_", l_val=active)
317 6 : IF (active) CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support RI-HFX.")
318 6 : CALL section_vals_val_get(section, "ACE%ACTIVE", l_val=active)
319 6 : IF (active) CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support ACE-HFX.")
320 6 : CALL section_vals_val_get(section, "SCREENING%SCREEN_ON_INITIAL_P", l_val=active)
321 6 : IF (active) CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support SCREEN_ON_INITIAL_P.")
322 : END IF
323 32 : CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=option)
324 32 : IF (option == xc_vdw_fun_nonloc) THEN
325 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support nonlocal vdW functionals.")
326 : END IF
327 32 : IF (dft_control%use_gauxc) THEN
328 8 : IF (xc_section_uses_native_skala_evaluator(xc_section)) THEN
329 : CALL cp_abort(__LOCATION__, &
330 0 : "GAPW low spin ROKS requires the molecular GauXC grid, not native Skala grids.")
331 : END IF
332 : ELSE
333 24 : section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
334 24 : needs = xc_functionals_get_needs(section, lsd=.TRUE., calc_potential=.TRUE.)
335 24 : IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
336 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS does not support density-Laplacian functionals.")
337 : END IF
338 : END IF
339 9247 : END SUBROUTINE check_gapw_low_spin_roks
340 :
341 : ! **************************************************************************************************
342 : !> \brief Read the input and the database files for the setup of the
343 : !> QUICKSTEP environment.
344 : !> \param qs_env ...
345 : !> \param para_env ...
346 : !> \param root_section ...
347 : !> \param globenv ...
348 : !> \param cp_subsys ...
349 : !> \param kpoint_env ...
350 : !> \param qmmm ...
351 : !> \param qmmm_env_qm ...
352 : !> \param force_env_section ...
353 : !> \param subsys_section ...
354 : !> \param use_motion_section ...
355 : !> \param silent ...
356 : !> \param multip ...
357 : !> \param charge ...
358 : !> \author Creation (22.05.2000,MK)
359 : ! **************************************************************************************************
360 64729 : SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, &
361 : qmmm, qmmm_env_qm, force_env_section, subsys_section, &
362 : use_motion_section, silent, multip, charge)
363 :
364 : TYPE(qs_environment_type), POINTER :: qs_env
365 : TYPE(mp_para_env_type), POINTER :: para_env
366 : TYPE(section_vals_type), OPTIONAL, POINTER :: root_section
367 : TYPE(global_environment_type), OPTIONAL, POINTER :: globenv
368 : TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys
369 : TYPE(kpoint_type), OPTIONAL, POINTER :: kpoint_env
370 : LOGICAL, INTENT(IN), OPTIONAL :: qmmm
371 : TYPE(qmmm_env_qm_type), OPTIONAL, POINTER :: qmmm_env_qm
372 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
373 : LOGICAL, INTENT(IN) :: use_motion_section
374 : LOGICAL, INTENT(IN), OPTIONAL :: silent
375 : INTEGER, INTENT(IN), OPTIONAL :: multip, charge
376 :
377 : CHARACTER(LEN=default_string_length) :: basis_type
378 : INTEGER :: ikind, method_id, nelectron_total, &
379 : nkind, nkp_grid(3), tddfpt_kernel
380 : LOGICAL :: dftb_kpoint_sym_restricted, do_active_space, do_admm, do_admm_rpa, do_bse, &
381 : do_debug_fdiff, do_debug_forces, do_debug_stress_tensor, do_dftb_scc, do_dftb_scc_high_l, &
382 : do_ec_hfx, do_et, do_exx, do_gw, do_hfx, do_kpoints, do_linear_response, do_mp2, &
383 : do_ri_mp2, do_ri_rpa, do_ri_sos_mp2, do_tddfpt, do_tddfpt_unsupported_kpoints, &
384 : do_wfc_low_scaling, do_wfc_low_scaling_kpoints, do_xtb_tblite, final_kpoint_reinit, &
385 : is_identical, is_semi, kpoint_explicit, kpoint_verbose, mp2_present, my_qmmm, &
386 : ot_energies, owned_kpoints, qmmm_decoupl, same_except_frac, use_real_wfn, use_ref_cell
387 9247 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rtmat
388 9247 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
389 : TYPE(cell_type), POINTER :: my_cell, my_cell_ref
390 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
391 : TYPE(dft_control_type), POINTER :: dft_control
392 : TYPE(distribution_1d_type), POINTER :: local_particles
393 : TYPE(energy_correction_type), POINTER :: ec_env
394 : TYPE(excited_energy_type), POINTER :: exstate_env
395 : TYPE(harris_type), POINTER :: harris_env
396 : TYPE(kpoint_type), POINTER :: kpoints
397 : TYPE(lri_environment_type), POINTER :: lri_env
398 9247 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
399 9247 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
400 : TYPE(qs_ks_env_type), POINTER :: ks_env
401 : TYPE(qs_subsys_type), POINTER :: subsys
402 : TYPE(qs_wf_history_type), POINTER :: wf_history
403 : TYPE(rel_control_type), POINTER :: rel_control
404 : TYPE(scf_control_type), POINTER :: scf_control
405 : TYPE(section_vals_type), POINTER :: active_space_section, admm_section, dft_section, &
406 : ec_hfx_section, ec_section, et_coupling_section, gw_section, hfx_section, kpoint_section, &
407 : mp2_section, rpa_hfx_section, tddfpt_section, transport_section
408 :
409 9247 : NULLIFY (my_cell, my_cell_ref, atomic_kind_set, particle_set, &
410 9247 : qs_kind_set, kpoint_section, dft_section, ec_section, &
411 9247 : subsys, ks_env, dft_control, blacs_env)
412 :
413 9247 : CALL set_qs_env(qs_env, input=force_env_section)
414 9247 : IF (.NOT. ASSOCIATED(subsys_section)) THEN
415 118 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
416 : END IF
417 9247 : CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=method_id)
418 :
419 : ! QMMM
420 9247 : my_qmmm = .FALSE.
421 9247 : IF (PRESENT(qmmm)) my_qmmm = qmmm
422 9247 : qmmm_decoupl = .FALSE.
423 9247 : IF (PRESENT(qmmm_env_qm)) THEN
424 396 : IF (qmmm_env_qm%qmmm_coupl_type == do_qmmm_gauss .OR. &
425 : qmmm_env_qm%qmmm_coupl_type == do_qmmm_swave) THEN
426 : ! For GAUSS/SWAVE methods there could be a DDAPC decoupling requested
427 : qmmm_decoupl = my_qmmm .AND. qmmm_env_qm%periodic .AND. qmmm_env_qm%multipole .AND. &
428 466 : method_id /= do_method_dftb .AND. method_id /= do_method_xtb
429 : END IF
430 396 : qs_env%qmmm_env_qm => qmmm_env_qm
431 : END IF
432 9247 : CALL set_qs_env(qs_env=qs_env, qmmm=my_qmmm)
433 :
434 : ! Possibly initialize arrays for SE
435 998 : SELECT CASE (method_id)
436 : CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, &
437 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
438 998 : CALL init_se_intd_array()
439 998 : is_semi = .TRUE.
440 : CASE (do_method_xtb, do_method_dftb)
441 1530 : is_semi = .TRUE.
442 : CASE DEFAULT
443 9247 : is_semi = .FALSE.
444 : END SELECT
445 :
446 36988 : ALLOCATE (subsys)
447 : CALL qs_subsys_create(subsys, para_env, &
448 : force_env_section=force_env_section, &
449 : subsys_section=subsys_section, &
450 : use_motion_section=use_motion_section, &
451 : root_section=root_section, &
452 : cp_subsys=cp_subsys, &
453 9247 : elkind=is_semi, silent=silent)
454 :
455 9247 : ALLOCATE (ks_env)
456 9247 : CALL qs_ks_env_create(ks_env)
457 9247 : CALL set_ks_env(ks_env, subsys=subsys)
458 9247 : CALL set_qs_env(qs_env, ks_env=ks_env)
459 :
460 : CALL qs_subsys_get(subsys, &
461 : cell=my_cell, &
462 : cell_ref=my_cell_ref, &
463 : use_ref_cell=use_ref_cell, &
464 : atomic_kind_set=atomic_kind_set, &
465 : qs_kind_set=qs_kind_set, &
466 9247 : particle_set=particle_set)
467 :
468 9247 : CALL set_ks_env(ks_env, para_env=para_env)
469 9247 : IF (PRESENT(globenv)) THEN
470 : CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
471 9231 : globenv%blacs_repeatable)
472 : ELSE
473 16 : CALL cp_blacs_env_create(blacs_env, para_env)
474 : END IF
475 9247 : CALL set_ks_env(ks_env, blacs_env=blacs_env)
476 9247 : CALL cp_blacs_env_release(blacs_env)
477 :
478 : ! *** Setup the grids for the G-space Interpolation if any
479 : CALL cp_ddapc_ewald_create(qs_env%cp_ddapc_ewald, qmmm_decoupl, my_cell, &
480 9247 : force_env_section, subsys_section, para_env)
481 :
482 : ! kpoints
483 9247 : IF (PRESENT(kpoint_env)) THEN
484 12 : owned_kpoints = .FALSE.
485 12 : kpoints => kpoint_env
486 12 : CALL set_qs_env(qs_env=qs_env, kpoints=kpoints)
487 12 : CALL kpoint_initialize(kpoints, particle_set, my_cell)
488 : ELSE
489 9235 : owned_kpoints = .TRUE.
490 9235 : NULLIFY (kpoints)
491 9235 : CALL kpoint_create(kpoints)
492 9235 : CALL set_qs_env(qs_env=qs_env, kpoints=kpoints)
493 9235 : kpoint_section => section_vals_get_subs_vals(qs_env%input, "DFT%KPOINTS")
494 9235 : CALL read_kpoint_section(kpoints, kpoint_section, my_cell%hmat, my_cell)
495 9235 : CALL section_vals_get(kpoint_section, explicit=kpoint_explicit)
496 : ot_energies = .FALSE.
497 9235 : CALL section_vals_val_get(qs_env%input, "DFT%SCF%OT%ENERGIES", l_val=ot_energies)
498 9235 : IF (.NOT. kpoint_explicit .AND. ot_energies) THEN
499 : ! Use the production Mermin driver instead of maintaining a second sparse Gamma path.
500 2 : CALL set_kpoint_info(kpoints, kp_scheme="GAMMA")
501 : END IF
502 9235 : CALL get_kpoint_info(kpoints, verbose=kpoint_verbose)
503 9235 : IF (kpoint_verbose) CALL set_kpoint_info(kpoints, verbose=.FALSE.)
504 : do_hfx = .FALSE.
505 9235 : hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
506 9235 : CALL section_vals_get(hfx_section, explicit=do_hfx)
507 : do_exx = .FALSE.
508 9235 : rpa_hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
509 9235 : CALL section_vals_get(rpa_hfx_section, explicit=do_exx)
510 : do_admm = .FALSE.
511 9235 : admm_section => section_vals_get_subs_vals(qs_env%input, "DFT%AUXILIARY_DENSITY_MATRIX_METHOD")
512 9235 : CALL section_vals_get(admm_section, explicit=do_admm)
513 : do_gw = .FALSE.
514 9235 : gw_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%GW")
515 9235 : CALL section_vals_get(gw_section, explicit=do_gw)
516 9235 : IF (.NOT. do_gw) THEN
517 9101 : gw_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%BANDSTRUCTURE%GW")
518 9101 : CALL section_vals_get(gw_section, explicit=do_gw)
519 : END IF
520 : do_tddfpt = .FALSE.
521 9235 : do_tddfpt_unsupported_kpoints = .FALSE.
522 9235 : do_bse = .FALSE.
523 9235 : tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
524 9235 : CALL section_vals_get(tddfpt_section, explicit=do_tddfpt)
525 9235 : IF (do_tddfpt) THEN
526 688 : CALL section_vals_val_get(tddfpt_section, "KERNEL", i_val=tddfpt_kernel)
527 688 : do_tddfpt_unsupported_kpoints = tddfpt_kernel /= tddfpt_kernel_none
528 688 : IF (.NOT. do_tddfpt_unsupported_kpoints) THEN
529 58 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
530 58 : IF (use_real_wfn) THEN
531 0 : CALL cp_abort(__LOCATION__, "K-point TDDFPT requires complex wavefunctions.")
532 : END IF
533 : END IF
534 688 : CALL section_vals_val_get(tddfpt_section, "DO_BSE", l_val=do_bse)
535 688 : IF (.NOT. do_bse) THEN
536 686 : CALL section_vals_val_get(tddfpt_section, "DO_BSE_W_ONLY", l_val=do_bse)
537 : END IF
538 688 : IF (.NOT. do_bse) THEN
539 684 : CALL section_vals_val_get(tddfpt_section, "DO_BSE_GW_ONLY", l_val=do_bse)
540 : END IF
541 : END IF
542 : do_active_space = .FALSE.
543 9235 : active_space_section => section_vals_get_subs_vals(qs_env%input, "DFT%ACTIVE_SPACE")
544 9235 : CALL section_vals_get(active_space_section, explicit=do_active_space)
545 9235 : do_xtb_tblite = .FALSE.
546 9235 : IF (method_id == do_method_xtb) THEN
547 : CALL section_vals_val_get(qs_env%input, "DFT%QS%XTB%TBLITE%_SECTION_PARAMETERS_", &
548 1232 : l_val=do_xtb_tblite)
549 : END IF
550 9235 : do_dftb_scc = .FALSE.
551 9235 : IF (method_id == do_method_dftb) THEN
552 : CALL section_vals_val_get(qs_env%input, "DFT%QS%DFTB%SELF_CONSISTENT", &
553 298 : l_val=do_dftb_scc)
554 : END IF
555 9235 : do_linear_response = .FALSE.
556 9235 : IF (PRESENT(globenv)) do_linear_response = globenv%run_type_id == linear_response_run
557 4 : do_debug_fdiff = .FALSE.
558 9231 : IF (PRESENT(globenv)) do_debug_fdiff = globenv%run_type_id == debug_run
559 9235 : IF (do_debug_fdiff .AND. PRESENT(root_section)) THEN
560 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
561 892 : l_val=do_debug_forces)
562 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_STRESS_TENSOR", &
563 892 : l_val=do_debug_stress_tensor)
564 1090 : do_debug_fdiff = do_debug_forces .OR. do_debug_stress_tensor
565 : END IF
566 9235 : do_mp2 = .FALSE.
567 9235 : do_ri_mp2 = .FALSE.
568 9235 : do_ri_sos_mp2 = .FALSE.
569 9235 : do_ri_rpa = .FALSE.
570 9235 : do_wfc_low_scaling = .FALSE.
571 9235 : do_wfc_low_scaling_kpoints = .FALSE.
572 9235 : mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
573 9235 : CALL section_vals_get(mp2_section, explicit=mp2_present)
574 9235 : IF (mp2_present) THEN
575 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%MP2%_SECTION_PARAMETERS_", &
576 512 : l_val=do_mp2)
577 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_MP2%_SECTION_PARAMETERS_", &
578 512 : l_val=do_ri_mp2)
579 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_SOS_MP2%_SECTION_PARAMETERS_", &
580 512 : l_val=do_ri_sos_mp2)
581 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%_SECTION_PARAMETERS_", &
582 512 : l_val=do_ri_rpa)
583 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%_SECTION_PARAMETERS_", &
584 512 : l_val=do_wfc_low_scaling)
585 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%DO_KPOINTS", &
586 512 : l_val=do_wfc_low_scaling_kpoints)
587 512 : IF (.NOT. do_bse) THEN
588 : CALL section_vals_val_get(qs_env%input, &
589 : "DFT%XC%WF_CORRELATION%RI_RPA%GW%BSE%_SECTION_PARAMETERS_", &
590 508 : l_val=do_bse)
591 : END IF
592 : END IF
593 : CALL restrict_unsupported_atomic_kpoint_symmetry(kpoints, method_id, do_hfx, do_exx, do_gw, &
594 : do_tddfpt_unsupported_kpoints, &
595 : do_active_space, do_linear_response, &
596 : do_debug_fdiff, &
597 : do_mp2 .OR. do_ri_mp2 .OR. do_ri_sos_mp2, &
598 : do_ri_rpa .AND. .NOT. do_gw, do_bse, &
599 : do_wfc_low_scaling, do_wfc_low_scaling_kpoints, &
600 27333 : do_xtb_tblite, do_admm, .FALSE.)
601 9235 : CALL kpoint_initialize(kpoints, particle_set, my_cell)
602 : END IF
603 :
604 : CALL qs_init_subsys(qs_env, para_env, subsys, my_cell, my_cell_ref, use_ref_cell, &
605 9247 : subsys_section, silent=silent, multip=multip, charge=charge)
606 :
607 9247 : CALL get_qs_env(qs_env, dft_control=dft_control)
608 9247 : IF (owned_kpoints) THEN
609 9235 : do_dftb_scc_high_l = .FALSE.
610 9235 : IF (method_id == do_method_dftb .AND. do_dftb_scc) THEN
611 224 : do_dftb_scc_high_l = dftb_kind_set_has_high_l(qs_kind_set)
612 : END IF
613 : CALL restrict_unsupported_atomic_kpoint_symmetry(kpoints, method_id, do_hfx, do_exx, do_gw, &
614 : do_tddfpt_unsupported_kpoints, &
615 : do_active_space, do_linear_response, &
616 : do_debug_fdiff, &
617 : do_mp2 .OR. do_ri_mp2 .OR. do_ri_sos_mp2, &
618 : do_ri_rpa .AND. .NOT. do_gw, do_bse, &
619 : do_wfc_low_scaling, do_wfc_low_scaling_kpoints, &
620 : do_xtb_tblite, do_admm, do_dftb_scc_high_l, &
621 27333 : restricted=dftb_kpoint_sym_restricted)
622 9235 : final_kpoint_reinit = dftb_kpoint_sym_restricted .OR. kpoint_verbose
623 : IF (final_kpoint_reinit) THEN
624 314 : CALL kpoint_reset_initialization(kpoints)
625 314 : CALL set_kpoint_info(kpoints, verbose=kpoint_verbose)
626 314 : CALL kpoint_initialize(kpoints, particle_set, my_cell)
627 : END IF
628 9235 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
629 9235 : CALL write_kpoint_info(kpoints, dft_section=dft_section)
630 : END IF
631 9247 : IF (method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
632 48 : CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
633 48 : CALL lri_env_basis("LRI", qs_env, lri_env, qs_kind_set)
634 9199 : ELSE IF (method_id == do_method_rigpw) THEN
635 : CALL cp_warn(__LOCATION__, "Experimental code: "// &
636 2 : "RIGPW should only be used for testing.")
637 2 : CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
638 2 : CALL lri_env_basis("RI", qs_env, lri_env, qs_kind_set)
639 : END IF
640 :
641 9247 : IF (my_qmmm .AND. PRESENT(qmmm_env_qm) .AND. .NOT. dft_control%qs_control%commensurate_mgrids) THEN
642 : IF ((qmmm_env_qm%qmmm_coupl_type == do_qmmm_gauss .OR. &
643 : qmmm_env_qm%qmmm_coupl_type == do_qmmm_swave) .AND. &
644 134 : method_id /= do_method_dftb .AND. method_id /= do_method_xtb) THEN
645 : CALL cp_abort(__LOCATION__, "QM/MM with coupling GAUSS or S-WAVE requires "// &
646 0 : "keyword FORCE_EVAL/DFT/MGRID/COMMENSURATE to be enabled.")
647 : END IF
648 : END IF
649 :
650 : ! more kpoint stuff
651 9247 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, blacs_env=blacs_env)
652 9247 : IF (do_kpoints) THEN
653 750 : IF (dft_control%qs_control%do_ls_scf) THEN
654 0 : CPABORT("DFT%KPOINTS are not implemented with QS/LS_SCF; use a real-space supercell instead.")
655 : END IF
656 750 : CALL kpoint_env_initialize(kpoints, para_env, blacs_env, with_aux_fit=dft_control%do_admm)
657 750 : CALL kpoint_initialize_mos(kpoints, qs_env%mos)
658 750 : CALL get_qs_env(qs_env=qs_env, wf_history=wf_history)
659 750 : CALL wfi_create_for_kp(wf_history)
660 : END IF
661 : ! basis set symmetry rotations
662 9247 : IF (do_kpoints) THEN
663 750 : CALL qs_basis_rotation(qs_env, kpoints)
664 : END IF
665 :
666 : do_hfx = .FALSE.
667 9247 : hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
668 9247 : CALL section_vals_get(hfx_section, explicit=do_hfx)
669 9247 : CALL get_qs_env(qs_env, dft_control=dft_control, scf_control=scf_control, nelectron_total=nelectron_total)
670 9247 : IF (do_hfx) THEN
671 : ! Retrieve particle_set and atomic_kind_set (needed for both kinds of initialization)
672 5808 : nkp_grid = 1
673 1452 : IF (do_kpoints) CALL get_kpoint_info(kpoints, nkp_grid=nkp_grid)
674 1452 : IF (dft_control%do_admm) THEN
675 520 : basis_type = 'AUX_FIT'
676 : ELSE
677 932 : basis_type = 'ORB'
678 : END IF
679 : CALL hfx_create(qs_env%x_data, para_env, hfx_section, atomic_kind_set, &
680 : qs_kind_set, particle_set, dft_control, my_cell, orb_basis=basis_type, &
681 1452 : nelectron_total=nelectron_total, nkp_grid=nkp_grid)
682 : END IF
683 :
684 9247 : mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
685 9247 : CALL section_vals_get(mp2_section, explicit=mp2_present)
686 9247 : IF (mp2_present) THEN
687 512 : CPASSERT(ASSOCIATED(qs_env%mp2_env))
688 512 : CALL read_mp2_section(qs_env%input, qs_env%mp2_env)
689 : ! create the EXX section if necessary
690 : do_exx = .FALSE.
691 512 : rpa_hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
692 512 : CALL section_vals_get(rpa_hfx_section, explicit=do_exx)
693 512 : IF (do_exx) THEN
694 :
695 : ! do_exx in call of hfx_create decides whether to go without ADMM (do_exx=.TRUE.) or with
696 : ! ADMM (do_exx=.FALSE.)
697 148 : CALL section_vals_val_get(mp2_section, "RI_RPA%ADMM", l_val=do_admm_rpa)
698 :
699 : ! Reuse the HFX integrals from the qs_env if applicable
700 148 : qs_env%mp2_env%ri_rpa%reuse_hfx = .TRUE.
701 148 : IF (.NOT. do_hfx) qs_env%mp2_env%ri_rpa%reuse_hfx = .FALSE.
702 148 : CALL compare_hfx_sections(hfx_section, rpa_hfx_section, is_identical, same_except_frac)
703 148 : IF (.NOT. (is_identical .OR. same_except_frac)) qs_env%mp2_env%ri_rpa%reuse_hfx = .FALSE.
704 148 : IF (dft_control%do_admm .AND. .NOT. do_admm_rpa) qs_env%mp2_env%ri_rpa%reuse_hfx = .FALSE.
705 :
706 148 : IF (.NOT. qs_env%mp2_env%ri_rpa%reuse_hfx) THEN
707 130 : IF (do_admm_rpa) THEN
708 10 : basis_type = 'AUX_FIT'
709 : ELSE
710 120 : basis_type = 'ORB'
711 : END IF
712 : CALL hfx_create(qs_env%mp2_env%ri_rpa%x_data, para_env, rpa_hfx_section, atomic_kind_set, &
713 : qs_kind_set, particle_set, dft_control, my_cell, orb_basis=basis_type, &
714 130 : nelectron_total=nelectron_total)
715 : ELSE
716 18 : qs_env%mp2_env%ri_rpa%x_data => qs_env%x_data
717 : END IF
718 : END IF
719 : END IF
720 :
721 9247 : IF (dft_control%qs_control%do_kg) THEN
722 82 : CALL cite_reference(Iannuzzi2006)
723 82 : CALL kg_env_create(qs_env, qs_env%kg_env, qs_kind_set, qs_env%input)
724 : END IF
725 :
726 9247 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
727 : CALL section_vals_val_get(dft_section, "EXCITED_STATES%_SECTION_PARAMETERS_", &
728 9247 : l_val=qs_env%excited_state)
729 9247 : NULLIFY (exstate_env)
730 9247 : CALL exstate_create(exstate_env, qs_env%excited_state, dft_section)
731 9247 : CALL set_qs_env(qs_env, exstate_env=exstate_env)
732 :
733 : et_coupling_section => section_vals_get_subs_vals(qs_env%input, &
734 9247 : "PROPERTIES%ET_COUPLING")
735 9247 : CALL section_vals_get(et_coupling_section, explicit=do_et)
736 9247 : IF (do_et) CALL et_coupling_create(qs_env%et_coupling)
737 :
738 9247 : transport_section => section_vals_get_subs_vals(qs_env%input, "DFT%TRANSPORT")
739 9247 : CALL section_vals_get(transport_section, explicit=qs_env%do_transport)
740 9247 : IF (qs_env%do_transport) THEN
741 0 : CALL transport_env_create(qs_env)
742 : END IF
743 :
744 9247 : CALL get_qs_env(qs_env, harris_env=harris_env)
745 9247 : IF (qs_env%harris_method) THEN
746 : ! initialize the Harris input density and potential integrals
747 28 : IF (harris_env%density_source == hden_atomic) THEN
748 8 : CALL get_qs_env(qs_env, local_particles=local_particles)
749 : CALL harris_rhoin_init(harris_env%rhoin, "RHOIN", qs_kind_set, atomic_kind_set, &
750 8 : local_particles, dft_control%nspins)
751 : ELSE
752 20 : harris_env%rhoin%nspin = dft_control%nspins
753 : END IF
754 : ! Print information of the HARRIS section
755 28 : CALL harris_write_input(harris_env)
756 : END IF
757 :
758 9247 : NULLIFY (ec_env)
759 9247 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
760 : CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%_SECTION_PARAMETERS_", &
761 9247 : l_val=qs_env%energy_correction)
762 9247 : ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
763 9247 : CALL ec_env_create(qs_env, ec_env, dft_section, ec_section)
764 9247 : CALL set_qs_env(qs_env, ec_env=ec_env)
765 :
766 9247 : IF (qs_env%energy_correction) THEN
767 : ! Energy correction with Hartree-Fock exchange
768 306 : ec_hfx_section => section_vals_get_subs_vals(ec_section, "XC%HF")
769 306 : CALL section_vals_get(ec_hfx_section, explicit=do_ec_hfx)
770 :
771 306 : IF (ec_env%do_ec_hfx) THEN
772 :
773 : ! kpoints and HFX not yet compatible
774 28 : IF (ec_env%do_kpoints) THEN
775 : CALL cp_abort(__LOCATION__, &
776 : "Energy correction methods with hybrid functionals "// &
777 0 : "and kpoints is not yet available.")
778 : END IF
779 :
780 : ! Hybrid functionals require same basis
781 28 : IF (ec_env%basis_inconsistent) THEN
782 : CALL cp_abort(__LOCATION__, &
783 : "Energy correction methods with hybrid functionals: "// &
784 : "correction and ground state need to use the same basis. "// &
785 0 : "Checked by comparing basis set names only.")
786 : END IF
787 :
788 : ! Similar to RPA_HFX we can check if HFX integrals from the qs_env can be reused
789 28 : IF (ec_env%do_ec_admm .AND. .NOT. dft_control%do_admm) THEN
790 0 : CALL cp_abort(__LOCATION__, "Need an ADMM input section for ADMM EC to work")
791 : END IF
792 :
793 28 : ec_env%reuse_hfx = .TRUE.
794 28 : IF (.NOT. do_hfx) ec_env%reuse_hfx = .FALSE.
795 28 : CALL compare_hfx_sections(hfx_section, ec_hfx_section, is_identical, same_except_frac)
796 28 : IF (.NOT. (is_identical .OR. same_except_frac)) ec_env%reuse_hfx = .FALSE.
797 28 : IF (dft_control%do_admm .AND. .NOT. ec_env%do_ec_admm) ec_env%reuse_hfx = .FALSE.
798 :
799 28 : IF (.NOT. ec_env%reuse_hfx) THEN
800 12 : IF (ec_env%do_ec_admm) THEN
801 2 : basis_type = 'AUX_FIT'
802 : ELSE
803 10 : basis_type = 'ORB'
804 : END IF
805 : CALL hfx_create(ec_env%x_data, para_env, ec_hfx_section, atomic_kind_set, &
806 : qs_kind_set, particle_set, dft_control, my_cell, orb_basis=basis_type, &
807 12 : nelectron_total=nelectron_total)
808 : ELSE
809 16 : ec_env%x_data => qs_env%x_data
810 : END IF
811 : END IF
812 :
813 : ! Print information of the EC section
814 306 : CALL ec_write_input(ec_env)
815 :
816 : END IF
817 :
818 9247 : IF (dft_control%qs_control%do_almo_scf) THEN
819 72 : CALL almo_scf_env_create(qs_env)
820 : END IF
821 :
822 : ! see if we have atomic relativistic corrections
823 9247 : CALL get_qs_env(qs_env, rel_control=rel_control)
824 9247 : IF (rel_control%rel_method /= rel_none) THEN
825 18 : IF (rel_control%rel_transformation == rel_trans_atom) THEN
826 18 : nkind = SIZE(atomic_kind_set)
827 46 : DO ikind = 1, nkind
828 28 : NULLIFY (rtmat)
829 28 : CALL calculate_atomic_relkin(atomic_kind_set(ikind), qs_kind_set(ikind), rel_control, rtmat)
830 46 : IF (ASSOCIATED(rtmat)) CALL set_qs_kind(qs_kind_set(ikind), reltmat=rtmat)
831 : END DO
832 : END IF
833 : END IF
834 :
835 9247 : END SUBROUTINE qs_init
836 :
837 : ! **************************************************************************************************
838 : !> \brief Restrict atomic k-point symmetry for methods not supporting it yet
839 : !> \param kpoints ...
840 : !> \param method_id ...
841 : !> \param do_hfx ...
842 : !> \param do_exx ...
843 : !> \param do_gw ...
844 : !> \param do_tddfpt ...
845 : !> \param do_active_space ...
846 : !> \param do_linear_response ...
847 : !> \param do_debug_fdiff ...
848 : !> \param do_mp2 ...
849 : !> \param do_rpa ...
850 : !> \param do_bse ...
851 : !> \param do_wfc_low_scaling ...
852 : !> \param do_wfc_low_scaling_kpoints ...
853 : !> \param do_xtb_tblite ...
854 : !> \param do_admm ...
855 : !> \param do_dftb_scc_high_l ...
856 : !> \param restricted ...
857 : ! **************************************************************************************************
858 18470 : SUBROUTINE restrict_unsupported_atomic_kpoint_symmetry(kpoints, method_id, do_hfx, do_exx, do_gw, &
859 : do_tddfpt, do_active_space, do_linear_response, &
860 : do_debug_fdiff, &
861 : do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
862 : do_wfc_low_scaling_kpoints, do_xtb_tblite, &
863 : do_admm, do_dftb_scc_high_l, restricted)
864 : TYPE(kpoint_type), POINTER :: kpoints
865 : INTEGER, INTENT(IN) :: method_id
866 : LOGICAL, INTENT(IN) :: do_hfx, do_exx, do_gw, do_tddfpt, do_active_space, &
867 : do_linear_response, do_debug_fdiff, do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
868 : do_wfc_low_scaling_kpoints, do_xtb_tblite, do_admm, do_dftb_scc_high_l
869 : LOGICAL, INTENT(OUT), OPTIONAL :: restricted
870 :
871 : CHARACTER(LEN=default_string_length) :: kp_scheme, reason
872 : LOGICAL :: full_grid, inversion_symmetry_only, &
873 : kpoint_symmetry
874 :
875 18470 : IF (PRESENT(restricted)) restricted = .FALSE.
876 :
877 : reason = unsupported_kpoint_method_reason(method_id, do_gw, do_tddfpt, do_linear_response, &
878 18470 : do_mp2, do_bse, do_xtb_tblite)
879 18470 : IF (LEN_TRIM(reason) > 0) THEN
880 3720 : CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme)
881 3720 : IF (LEN_TRIM(kp_scheme) > 0 .AND. TRIM(kp_scheme) /= "NONE") THEN
882 0 : IF (TRIM(reason) == "GW") THEN
883 : CALL cp_abort(__LOCATION__, &
884 : "DFT%KPOINTS are not supported with GW; use "// &
885 : "WF_CORRELATION%LOW_SCALING%KPOINTS and RI_RPA%GW%KPOINTS_SELF_ENERGY "// &
886 0 : "for GW k-point sampling.")
887 : ELSE
888 : CALL cp_abort(__LOCATION__, &
889 : "DFT%KPOINTS are not supported with "//TRIM(reason)// &
890 0 : "; remove DFT%KPOINTS for these calculations.")
891 : END IF
892 : END IF
893 : END IF
894 18470 : IF (do_active_space) THEN
895 216 : CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme)
896 216 : IF (LEN_TRIM(kp_scheme) > 0 .AND. TRIM(kp_scheme) /= "NONE" .AND. &
897 : TRIM(kp_scheme) /= "GAMMA") THEN
898 : CALL cp_abort(__LOCATION__, &
899 : "Only Gamma-point DFT%KPOINTS are supported with ACTIVE_SPACE; "// &
900 0 : "use SCHEME GAMMA, SCHEME NONE, or remove DFT%KPOINTS.")
901 : END IF
902 : END IF
903 :
904 : CALL get_kpoint_info(kpoints, symmetry=kpoint_symmetry, full_grid=full_grid, &
905 18470 : inversion_symmetry_only=inversion_symmetry_only)
906 19210 : IF (.NOT. (kpoint_symmetry .AND. .NOT. full_grid .AND. .NOT. inversion_symmetry_only)) RETURN
907 :
908 : reason = unsupported_atomic_kpoint_symmetry_reason(method_id, do_hfx, do_exx, do_gw, &
909 : do_tddfpt, do_active_space, do_linear_response, &
910 : do_debug_fdiff, &
911 : do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
912 : do_wfc_low_scaling_kpoints, do_xtb_tblite, &
913 754 : do_admm, do_dftb_scc_high_l)
914 754 : IF (LEN_TRIM(reason) == 0) RETURN
915 :
916 : CALL cp_warn(__LOCATION__, &
917 : "Atomic k-point symmetry is currently not implemented for "//TRIM(reason)// &
918 14 : "; restricting to inversion/time-reversal symmetry.")
919 14 : CALL set_kpoint_info(kpoints, inversion_symmetry_only=.TRUE.)
920 14 : IF (PRESENT(restricted)) restricted = .TRUE.
921 :
922 : END SUBROUTINE restrict_unsupported_atomic_kpoint_symmetry
923 :
924 : ! **************************************************************************************************
925 : !> \brief Return the reason why k-points are not enabled for a method
926 : !> \param method_id ...
927 : !> \param do_gw ...
928 : !> \param do_tddfpt ...
929 : !> \param do_linear_response ...
930 : !> \param do_mp2 ...
931 : !> \param do_bse ...
932 : !> \param do_xtb_tblite ...
933 : !> \return reason
934 : ! **************************************************************************************************
935 18470 : FUNCTION unsupported_kpoint_method_reason(method_id, do_gw, do_tddfpt, do_linear_response, &
936 : do_mp2, do_bse, do_xtb_tblite) RESULT(reason)
937 : INTEGER, INTENT(IN) :: method_id
938 : LOGICAL, INTENT(IN) :: do_gw, do_tddfpt, do_linear_response, &
939 : do_mp2, do_bse, do_xtb_tblite
940 : CHARACTER(LEN=default_string_length) :: reason
941 :
942 : reason = ""
943 : MARK_USED(do_gw)
944 : MARK_USED(do_mp2)
945 : MARK_USED(do_xtb_tblite)
946 :
947 18470 : IF (do_bse) THEN
948 112 : reason = "BSE"
949 112 : RETURN
950 : END IF
951 18358 : IF (do_tddfpt) THEN
952 1252 : reason = "TDDFPT/TDDFT"
953 1252 : RETURN
954 : END IF
955 17106 : IF (do_linear_response) THEN
956 356 : reason = "LINEAR_RESPONSE/DFPT"
957 356 : RETURN
958 : END IF
959 16754 : SELECT CASE (method_id)
960 : CASE (do_method_rigpw)
961 4 : reason = "RIGPW"
962 : CASE (do_method_ofgpw)
963 0 : reason = "OFGPW"
964 : CASE (do_method_mndo, do_method_mndod, do_method_am1, do_method_pm3, &
965 : do_method_pm6, do_method_pm6fm, do_method_pdg, do_method_rm1, do_method_pnnl)
966 1996 : reason = "semiempirical methods"
967 : CASE DEFAULT
968 16750 : reason = ""
969 : END SELECT
970 :
971 : END FUNCTION unsupported_kpoint_method_reason
972 :
973 : ! **************************************************************************************************
974 : !> \brief Return the reason why atomic k-point symmetry is not enabled
975 : !> \param method_id ...
976 : !> \param do_hfx ...
977 : !> \param do_exx ...
978 : !> \param do_gw ...
979 : !> \param do_tddfpt ...
980 : !> \param do_active_space ...
981 : !> \param do_linear_response ...
982 : !> \param do_debug_fdiff ...
983 : !> \param do_mp2 ...
984 : !> \param do_rpa ...
985 : !> \param do_bse ...
986 : !> \param do_wfc_low_scaling ...
987 : !> \param do_wfc_low_scaling_kpoints ...
988 : !> \param do_xtb_tblite ...
989 : !> \param do_admm ...
990 : !> \param do_dftb_scc_high_l ...
991 : !> \return reason
992 : ! **************************************************************************************************
993 754 : FUNCTION unsupported_atomic_kpoint_symmetry_reason(method_id, do_hfx, do_exx, do_gw, do_tddfpt, &
994 : do_active_space, do_linear_response, do_debug_fdiff, &
995 : do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
996 : do_wfc_low_scaling_kpoints, do_xtb_tblite, &
997 : do_admm, do_dftb_scc_high_l) RESULT(reason)
998 : INTEGER, INTENT(IN) :: method_id
999 : LOGICAL, INTENT(IN) :: do_hfx, do_exx, do_gw, do_tddfpt, do_active_space, &
1000 : do_linear_response, do_debug_fdiff, do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
1001 : do_wfc_low_scaling_kpoints, do_xtb_tblite, do_admm, do_dftb_scc_high_l
1002 : CHARACTER(LEN=default_string_length) :: reason
1003 :
1004 754 : reason = ""
1005 : MARK_USED(do_debug_fdiff)
1006 : MARK_USED(do_xtb_tblite)
1007 :
1008 826 : SELECT CASE (method_id)
1009 : CASE (do_method_dftb)
1010 72 : IF (do_dftb_scc_high_l) reason = "SCC-DFTB with d orbitals"
1011 : CASE (do_method_lrigpw)
1012 2 : reason = "LRIGPW"
1013 : CASE (do_method_rigpw)
1014 0 : reason = "RIGPW"
1015 : CASE (do_method_mndo, do_method_mndod, do_method_am1, do_method_pm3, &
1016 : do_method_pm6, do_method_pm6fm, do_method_pdg, do_method_rm1, do_method_pnnl)
1017 0 : reason = "semiempirical methods"
1018 : CASE DEFAULT
1019 754 : reason = ""
1020 : END SELECT
1021 :
1022 754 : IF (LEN_TRIM(reason) > 0) RETURN
1023 748 : IF ((do_hfx .OR. do_exx) .AND. do_admm) THEN
1024 2 : reason = "HFX/HF with ADMM"
1025 746 : ELSE IF (do_bse) THEN
1026 0 : reason = "BSE"
1027 746 : ELSE IF (do_gw) THEN
1028 4 : reason = "GW"
1029 742 : ELSE IF (do_tddfpt) THEN
1030 0 : reason = "TDDFPT/TDDFT"
1031 742 : ELSE IF (do_active_space) THEN
1032 0 : reason = "ACTIVE_SPACE"
1033 742 : ELSE IF (do_linear_response) THEN
1034 0 : reason = "LINEAR_RESPONSE/DFPT"
1035 742 : ELSE IF (do_mp2) THEN
1036 0 : reason = "MP2"
1037 742 : ELSE IF (do_rpa .AND. do_wfc_low_scaling_kpoints) THEN
1038 2 : reason = "LOW_SCALING RPA"
1039 740 : ELSE IF (do_wfc_low_scaling) THEN
1040 0 : reason = "LOW_SCALING WF_CORRELATION"
1041 740 : ELSE IF (do_rpa) THEN
1042 0 : reason = "RPA"
1043 : END IF
1044 :
1045 : END FUNCTION unsupported_atomic_kpoint_symmetry_reason
1046 :
1047 : ! **************************************************************************************************
1048 : !> \brief Return whether the DFTB kind set contains d orbitals
1049 : !> \param qs_kind_set ...
1050 : !> \return has_high_l
1051 : ! **************************************************************************************************
1052 224 : FUNCTION dftb_kind_set_has_high_l(qs_kind_set) RESULT(has_high_l)
1053 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1054 : LOGICAL :: has_high_l
1055 :
1056 : INTEGER :: ikind, lmax
1057 : LOGICAL :: any_defined, defined
1058 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
1059 :
1060 224 : has_high_l = .TRUE.
1061 224 : IF (.NOT. ASSOCIATED(qs_kind_set)) RETURN
1062 :
1063 224 : any_defined = .FALSE.
1064 704 : DO ikind = 1, SIZE(qs_kind_set)
1065 484 : NULLIFY (dftb_parameter)
1066 484 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_parameter)
1067 484 : IF (.NOT. ASSOCIATED(dftb_parameter)) CYCLE
1068 : defined = .FALSE.
1069 : lmax = -1
1070 484 : CALL get_dftb_atom_param(dftb_parameter, defined=defined, lmax=lmax)
1071 484 : IF (.NOT. defined) CYCLE
1072 484 : any_defined = .TRUE.
1073 1188 : IF (lmax > 1) RETURN
1074 : END DO
1075 :
1076 220 : IF (any_defined) has_high_l = .FALSE.
1077 :
1078 : END FUNCTION dftb_kind_set_has_high_l
1079 :
1080 : ! **************************************************************************************************
1081 : !> \brief Initialize the qs environment (subsys)
1082 : !> \param qs_env ...
1083 : !> \param para_env ...
1084 : !> \param subsys ...
1085 : !> \param cell ...
1086 : !> \param cell_ref ...
1087 : !> \param use_ref_cell ...
1088 : !> \param subsys_section ...
1089 : !> \param silent ...
1090 : !> \param multip ...
1091 : !> \param charge ...
1092 : !> \author Creation (22.05.2000,MK)
1093 : ! **************************************************************************************************
1094 9247 : SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell, subsys_section, &
1095 : silent, multip, charge)
1096 :
1097 : TYPE(qs_environment_type), POINTER :: qs_env
1098 : TYPE(mp_para_env_type), POINTER :: para_env
1099 : TYPE(qs_subsys_type), POINTER :: subsys
1100 : TYPE(cell_type), POINTER :: cell, cell_ref
1101 : LOGICAL, INTENT(in) :: use_ref_cell
1102 : TYPE(section_vals_type), POINTER :: subsys_section
1103 : LOGICAL, INTENT(in), OPTIONAL :: silent
1104 : INTEGER, INTENT(IN), OPTIONAL :: multip, charge
1105 :
1106 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_init_subsys'
1107 :
1108 : CHARACTER(len=2) :: element_symbol
1109 : INTEGER :: gfn_type, handle, ikind, ispin, iw, lmax_sphere, maxl, maxlgto, maxlgto_lri, &
1110 : maxlgto_nuc, maxlppl, maxlppnl, method_id, multiplicity, my_ival, n_ao, n_mo_add, &
1111 : n_mo_auto, n_mo_target, natom, nelectron, ngauss, nkind, nlumo_dos, nlumo_molden, &
1112 : nlumo_required, output_unit, sort_basis, tnadd_method
1113 : INTEGER, DIMENSION(2) :: n_mo, nelectron_spin
1114 : INTEGER, DIMENSION(5) :: ngaussflex, occ
1115 9247 : INTEGER, DIMENSION(:), POINTER :: mo_index_range
1116 : LOGICAL :: all_potential_present, be_silent, cneo_potential_present, do_kpoints, do_ri_hfx, &
1117 : do_ri_mp2, do_ri_rpa, do_ri_sos_mp2, do_rpa_ri_exx, do_wfc_im_time, e1terms, &
1118 : has_unit_metric, lribas, mp2_present, orb_gradient, paw_atom
1119 : REAL(KIND=dp) :: alpha, ccore, ewald_rcut, fxx, maxocc, &
1120 : rc, rcut, total_zeff_corr, &
1121 : verlet_skin, zeff_correction
1122 9247 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1123 : TYPE(cp_logger_type), POINTER :: logger
1124 : TYPE(dft_control_type), POINTER :: dft_control
1125 : TYPE(dftb_control_type), POINTER :: dftb_control
1126 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
1127 : TYPE(ewald_environment_type), POINTER :: ewald_env
1128 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1129 : TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env
1130 : TYPE(gapw_control_type), POINTER :: gapw_control
1131 : TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, &
1132 : rhoin_basis, ri_aux_basis_set, &
1133 : ri_hfx_basis, ri_xas_basis, &
1134 : tmp_basis_set
1135 : TYPE(harris_type), POINTER :: harris_env
1136 : TYPE(local_rho_type), POINTER :: local_rho_set
1137 : TYPE(lri_environment_type), POINTER :: lri_env
1138 9247 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_last_converged
1139 9247 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1140 9247 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1141 : TYPE(mp2_type), POINTER :: mp2_env
1142 : TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole
1143 9247 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1144 : TYPE(pw_env_type), POINTER :: pw_env
1145 : TYPE(qs_control_type), POINTER :: qs_control
1146 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
1147 9247 : POINTER :: dftb_potential
1148 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1149 : TYPE(qs_energy_type), POINTER :: energy
1150 9247 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1151 : TYPE(qs_gcp_type), POINTER :: gcp_env
1152 9247 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1153 : TYPE(qs_kind_type), POINTER :: qs_kind
1154 : TYPE(qs_ks_env_type), POINTER :: ks_env
1155 : TYPE(qs_wf_history_type), POINTER :: wf_history
1156 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
1157 9247 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1158 : TYPE(scf_control_type), POINTER :: scf_control
1159 : TYPE(se_taper_type), POINTER :: se_taper
1160 : TYPE(section_vals_type), POINTER :: dft_section, et_coupling_section, et_ddapc_section, &
1161 : ewald_section, harris_section, lri_section, mp2_section, nl_section, poisson_section, &
1162 : pp_section, print_section, qs_section, rixs_section, se_section, tddfpt_section, &
1163 : xc_section
1164 : TYPE(semi_empirical_control_type), POINTER :: se_control
1165 : TYPE(semi_empirical_si_type), POINTER :: se_store_int_env
1166 : TYPE(xtb_control_type), POINTER :: xtb_control
1167 :
1168 9247 : CALL timeset(routineN, handle)
1169 9247 : NULLIFY (logger)
1170 9247 : logger => cp_get_default_logger()
1171 9247 : output_unit = cp_logger_get_default_io_unit(logger)
1172 :
1173 9247 : be_silent = .FALSE.
1174 9247 : IF (PRESENT(silent)) be_silent = silent
1175 :
1176 9247 : CALL cite_reference(cp2kqs2020)
1177 :
1178 : ! Initialise the Quickstep environment
1179 9247 : NULLIFY (mos, se_taper)
1180 9247 : NULLIFY (dft_control)
1181 9247 : NULLIFY (energy)
1182 9247 : NULLIFY (force)
1183 9247 : NULLIFY (local_molecules)
1184 9247 : NULLIFY (local_particles)
1185 9247 : NULLIFY (scf_control)
1186 9247 : NULLIFY (dft_section)
1187 9247 : NULLIFY (et_coupling_section)
1188 9247 : NULLIFY (ks_env)
1189 9247 : NULLIFY (mos_last_converged)
1190 9247 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
1191 9247 : qs_section => section_vals_get_subs_vals(dft_section, "QS")
1192 9247 : et_coupling_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%ET_COUPLING")
1193 : ! reimplemented TDDFPT
1194 9247 : tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
1195 9247 : rixs_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS")
1196 :
1197 : CALL qs_subsys_get(subsys, particle_set=particle_set, &
1198 : qs_kind_set=qs_kind_set, &
1199 : atomic_kind_set=atomic_kind_set, &
1200 : molecule_set=molecule_set, &
1201 9247 : molecule_kind_set=molecule_kind_set)
1202 :
1203 : ! Read the input section with the DFT control parameters
1204 9247 : CALL read_dft_control(dft_control, dft_section, cell)
1205 :
1206 : ! Set periodicity flag
1207 36988 : dft_control%qs_control%periodicity = SUM(cell%perd)
1208 :
1209 : ! Read the input section with the Quickstep control parameters
1210 9247 : CALL read_qs_section(dft_control%qs_control, qs_section, cell)
1211 :
1212 9247 : CALL check_gapw_low_spin_roks(dft_control, dft_section, subsys_section)
1213 :
1214 : ! Print the Quickstep program banner (copyright and version number)
1215 9247 : IF (.NOT. be_silent) THEN
1216 9219 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PROGRAM_BANNER", extension=".Log")
1217 9219 : CALL section_vals_val_get(qs_section, "METHOD", i_val=method_id)
1218 6695 : SELECT CASE (method_id)
1219 : CASE DEFAULT
1220 6695 : CALL qs_header(iw)
1221 : CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, &
1222 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
1223 998 : CALL se_header(iw)
1224 : CASE (do_method_dftb)
1225 298 : CALL dftb_header(iw)
1226 : CASE (do_method_xtb)
1227 9219 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
1228 198 : CALL tblite_header(iw, dft_control%qs_control%xtb_control%tblite_method)
1229 : ELSE
1230 1030 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
1231 1030 : CALL xtb_header(iw, gfn_type)
1232 : END IF
1233 : END SELECT
1234 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
1235 9219 : "PRINT%PROGRAM_BANNER")
1236 : END IF
1237 :
1238 9247 : IF (dft_control%do_sccs .AND. dft_control%qs_control%gapw) THEN
1239 0 : CPABORT("SCCS is not yet implemented with GAPW")
1240 : END IF
1241 9247 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
1242 9247 : IF (do_kpoints) THEN
1243 : IF (dft_control%nspins == 2 .AND. dft_control%qs_control%xtb .AND. &
1244 : .NOT. dft_control%qs_control%xtb_control%do_tblite .AND. &
1245 : dft_control%qs_control%xtb_control%gfn_type == gfn1xtb .AND. &
1246 750 : dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_tblite .AND. &
1247 : .NOT. dft_control%qs_control%xtb_control%tblite_mixer_damping_explicit) THEN
1248 : CALL cp_warn(__LOCATION__, &
1249 : "Reducing XTB/TBLITE_MIXER/DAMPING to 0.25 for CP2K-internal GFN1-xTB "// &
1250 : "UKS k-point calculations with SCC_MIXER TBLITE. Set XTB/TBLITE_MIXER/DAMPING "// &
1251 0 : "explicitly to override this conservative fallback.")
1252 0 : dft_control%qs_control%xtb_control%tblite_mixer_damping = 0.25_dp
1253 : END IF
1254 : ! reset some of the settings for wfn extrapolation for kpoints
1255 750 : SELECT CASE (dft_control%qs_control%wf_interpolation_method_nr)
1256 : CASE (wfi_linear_wf_method_nr, wfi_linear_ps_method_nr)
1257 : CALL cp_warn(__LOCATION__, "Linear WFN-based extrapolation methods are not "// &
1258 0 : "implemented for k-points. Switching to USE_PREV_WF.")
1259 750 : dft_control%qs_control%wf_interpolation_method_nr = wfi_use_prev_wf_method_nr
1260 : END SELECT
1261 : END IF
1262 :
1263 : ! Check if any kind of electron transfer calculation has to be performed
1264 9247 : CALL section_vals_val_get(et_coupling_section, "TYPE_OF_CONSTRAINT", i_val=my_ival)
1265 9247 : dft_control%qs_control%et_coupling_calc = .FALSE.
1266 9247 : IF (my_ival == do_et_ddapc) THEN
1267 0 : et_ddapc_section => section_vals_get_subs_vals(et_coupling_section, "DDAPC_RESTRAINT_A")
1268 0 : dft_control%qs_control%et_coupling_calc = .TRUE.
1269 0 : dft_control%qs_control%ddapc_restraint = .TRUE.
1270 0 : CALL read_ddapc_section(dft_control%qs_control, ddapc_restraint_section=et_ddapc_section)
1271 : END IF
1272 :
1273 9247 : CALL read_mgrid_section(dft_control%qs_control, dft_section)
1274 :
1275 : ! Reimplemented TDDFPT
1276 9247 : CALL read_tddfpt2_control(dft_control%tddfpt2_control, tddfpt_section, dft_control%qs_control)
1277 :
1278 : ! RIXS
1279 9247 : CALL section_vals_get(rixs_section, explicit=qs_env%do_rixs)
1280 9247 : IF (qs_env%do_rixs) THEN
1281 16 : CALL read_rixs_control(dft_control%rixs_control, rixs_section, dft_control%qs_control)
1282 : END IF
1283 :
1284 : ! Create relativistic control section
1285 : BLOCK
1286 : TYPE(rel_control_type), POINTER :: rel_control
1287 9247 : ALLOCATE (rel_control)
1288 9247 : CALL rel_c_create(rel_control)
1289 9247 : CALL rel_c_read_parameters(rel_control, dft_section)
1290 9247 : CALL set_qs_env(qs_env, rel_control=rel_control)
1291 : END BLOCK
1292 :
1293 : ! Read DFTB parameter files
1294 9247 : IF (dft_control%qs_control%method_id == do_method_dftb) THEN
1295 298 : NULLIFY (ewald_env, ewald_pw, dftb_potential)
1296 298 : dftb_control => dft_control%qs_control%dftb_control
1297 : CALL qs_dftb_param_init(atomic_kind_set, qs_kind_set, dftb_control, dftb_potential, &
1298 298 : subsys_section=subsys_section, para_env=para_env)
1299 298 : CALL set_qs_env(qs_env, dftb_potential=dftb_potential)
1300 : ! check for Ewald
1301 298 : IF (dftb_control%do_ewald) THEN
1302 2432 : ALLOCATE (ewald_env)
1303 152 : CALL ewald_env_create(ewald_env, para_env)
1304 152 : poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
1305 152 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1306 152 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1307 152 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
1308 152 : CALL get_qs_kind_set(qs_kind_set, basis_rcut=ewald_rcut)
1309 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
1310 152 : cell_periodic=cell%perd)
1311 152 : ALLOCATE (ewald_pw)
1312 152 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
1313 152 : CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
1314 : END IF
1315 8949 : ELSE IF (dft_control%qs_control%method_id == do_method_xtb) THEN
1316 : ! Read xTB parameter file
1317 1232 : xtb_control => dft_control%qs_control%xtb_control
1318 1232 : CALL get_qs_env(qs_env, nkind=nkind)
1319 1232 : IF (xtb_control%do_tblite) THEN
1320 : ! put geometry to tblite
1321 198 : CALL tb_init_geometry(qs_env, qs_env%tb_tblite)
1322 : ! select tblite method
1323 : CALL tb_set_calculator(qs_env%tb_tblite, xtb_control%tblite_method, &
1324 198 : xtb_control%tblite_accuracy, xtb_control%tblite_param_file)
1325 : !set up wave function
1326 198 : CALL tb_init_wf(qs_env%tb_tblite, dft_control)
1327 : !get basis set
1328 532 : DO ikind = 1, nkind
1329 334 : qs_kind => qs_kind_set(ikind)
1330 : ! Setup proper xTB parameters
1331 334 : CPASSERT(.NOT. ASSOCIATED(qs_kind%xtb_parameter))
1332 334 : CALL allocate_xtb_atom_param(qs_kind%xtb_parameter)
1333 : ! Set default parameters
1334 334 : CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
1335 :
1336 334 : NULLIFY (tmp_basis_set)
1337 334 : CALL tb_get_basis(qs_env%tb_tblite, tmp_basis_set, element_symbol, qs_kind%xtb_parameter, occ)
1338 334 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
1339 334 : CALL set_xtb_atom_param(qs_kind%xtb_parameter, occupation=occ)
1340 :
1341 : !setting the potential for the computation
1342 334 : zeff_correction = 0.0_dp
1343 : CALL init_potential(qs_kind%all_potential, itype="BARE", &
1344 2202 : zeff=REAL(SUM(occ), dp), zeff_correction=zeff_correction)
1345 : END DO
1346 : ELSE
1347 1034 : NULLIFY (ewald_env, ewald_pw)
1348 3314 : DO ikind = 1, nkind
1349 2280 : qs_kind => qs_kind_set(ikind)
1350 : ! Setup proper xTB parameters
1351 2280 : CPASSERT(.NOT. ASSOCIATED(qs_kind%xtb_parameter))
1352 2280 : CALL allocate_xtb_atom_param(qs_kind%xtb_parameter)
1353 : ! Set default parameters
1354 2280 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
1355 2280 : CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
1356 : CALL xtb_parameters_init(qs_kind%xtb_parameter, gfn_type, element_symbol, &
1357 : xtb_control%parameter_file_path, xtb_control%parameter_file_name, &
1358 2280 : para_env)
1359 2280 : IF (xtb_control%do_spinpol) THEN
1360 : CALL xtb_spinpol_init(qs_kind%xtb_parameter, gfn_type, element_symbol, &
1361 : xtb_control%parameter_file_path, xtb_control%spinpol_param_file_name, &
1362 58 : para_env)
1363 58 : CALL xtb_spinpol_ext(qs_kind%xtb_parameter, gfn_type, xtb_control)
1364 : END IF
1365 : ! set dependent parameters
1366 2280 : CALL xtb_parameters_set(qs_kind%xtb_parameter)
1367 : ! Generate basis set
1368 2280 : NULLIFY (tmp_basis_set)
1369 2280 : IF (qs_kind%xtb_parameter%z == 1) THEN
1370 : ! special case hydrogen
1371 516 : ngauss = xtb_control%h_sto_ng
1372 : ELSE
1373 1764 : ngauss = xtb_control%sto_ng
1374 : END IF
1375 2280 : IF (qs_kind%xtb_parameter%defined) THEN
1376 2278 : IF (xtb_control%sto_flex) THEN
1377 8 : CALL get_xtb_atom_param(qs_kind%xtb_parameter, ngauss=ngaussflex)
1378 8 : CALL init_xtb_basis(qs_kind%xtb_parameter, tmp_basis_set, ngauss, ngaussflex)
1379 : ELSE
1380 2270 : CALL init_xtb_basis(qs_kind%xtb_parameter, tmp_basis_set, ngauss)
1381 : END IF
1382 2278 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
1383 : ELSE
1384 2 : CALL set_qs_kind(qs_kind, ghost=.TRUE.)
1385 2 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
1386 2 : DEALLOCATE (qs_kind%all_potential%elec_conf)
1387 2 : DEALLOCATE (qs_kind%all_potential)
1388 : END IF
1389 : END IF
1390 : ! potential
1391 2280 : IF (qs_kind%xtb_parameter%defined) THEN
1392 2278 : zeff_correction = 0.0_dp
1393 : CALL init_potential(qs_kind%all_potential, itype="BARE", &
1394 2278 : zeff=qs_kind%xtb_parameter%zeff, zeff_correction=zeff_correction)
1395 2278 : CALL get_potential(qs_kind%all_potential, alpha_core_charge=alpha)
1396 2278 : ccore = qs_kind%xtb_parameter%zeff*SQRT((alpha/pi)**3)
1397 2278 : CALL set_potential(qs_kind%all_potential, ccore_charge=ccore)
1398 2278 : qs_kind%xtb_parameter%zeff = qs_kind%xtb_parameter%zeff - zeff_correction
1399 : END IF
1400 3314 : CALL write_xtb_atom_param(qs_kind%xtb_parameter, gfn_type, subsys_section)
1401 : END DO
1402 : !
1403 : ! set repulsive potential range
1404 : !
1405 4136 : ALLOCATE (xtb_control%rcpair(nkind, nkind))
1406 1034 : CALL xtb_pp_radius(qs_kind_set, xtb_control%rcpair, xtb_control%eps_pair, xtb_control%kf)
1407 : ! check for Ewald
1408 1034 : IF (xtb_control%do_ewald) THEN
1409 3648 : ALLOCATE (ewald_env)
1410 228 : CALL ewald_env_create(ewald_env, para_env)
1411 228 : poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
1412 228 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1413 228 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1414 228 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
1415 228 : IF (gfn_type == 0) THEN
1416 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
1417 48 : silent=silent, pset="EEQ", cell_periodic=cell%perd)
1418 : ELSE
1419 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
1420 180 : silent=silent, cell_periodic=cell%perd)
1421 : END IF
1422 228 : ALLOCATE (ewald_pw)
1423 228 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
1424 228 : CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
1425 : END IF
1426 : END IF
1427 : END IF
1428 : ! lri or ri env initialization
1429 9247 : lri_section => section_vals_get_subs_vals(qs_section, "LRIGPW")
1430 : IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. &
1431 9247 : dft_control%qs_control%lri_optbas .OR. &
1432 : dft_control%qs_control%method_id == do_method_rigpw) THEN
1433 50 : CALL lri_env_init(lri_env, lri_section)
1434 50 : CALL set_qs_env(qs_env, lri_env=lri_env)
1435 : END IF
1436 :
1437 : ! Check basis and fill in missing parts
1438 9247 : CALL check_qs_kind_set(qs_kind_set, dft_control, subsys_section=subsys_section)
1439 :
1440 : ! Check that no all-electron potential is present if GPW or GAPW_XC
1441 9247 : CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
1442 : IF ((dft_control%qs_control%method_id == do_method_gpw) .OR. &
1443 9247 : (dft_control%qs_control%method_id == do_method_gapw_xc) .OR. &
1444 : (dft_control%qs_control%method_id == do_method_ofgpw)) THEN
1445 5321 : IF (all_potential_present) THEN
1446 0 : CPABORT("All-electron calculations with GPW, GAPW_XC, and OFGPW are not implemented")
1447 : END IF
1448 : END IF
1449 :
1450 : ! Check that no cneo potential is present if not GAPW
1451 9247 : CALL get_qs_kind_set(qs_kind_set, cneo_potential_present=cneo_potential_present)
1452 9247 : IF (cneo_potential_present .AND. &
1453 : dft_control%qs_control%method_id /= do_method_gapw) THEN
1454 0 : CPABORT("CNEO calculations require GAPW method")
1455 : END IF
1456 :
1457 : ! DFT+U
1458 9247 : CALL get_qs_kind_set(qs_kind_set, dft_plus_u_atom_present=dft_control%dft_plus_u)
1459 :
1460 : ! Minimum tracking linear response U and J
1461 9247 : CALL get_qs_kind_set(qs_kind_set, do_mtlr_present=dft_control%mtlr_u_j)
1462 :
1463 9247 : IF (dft_control%do_admm) THEN
1464 : ! Check if ADMM basis is available
1465 528 : CALL get_qs_env(qs_env, nkind=nkind)
1466 1506 : DO ikind = 1, nkind
1467 978 : NULLIFY (aux_fit_basis)
1468 978 : qs_kind => qs_kind_set(ikind)
1469 978 : CALL get_qs_kind(qs_kind, basis_set=aux_fit_basis, basis_type="AUX_FIT")
1470 1506 : IF (.NOT. (ASSOCIATED(aux_fit_basis))) THEN
1471 : ! AUX_FIT basis set is not available
1472 0 : CPABORT("AUX_FIT basis set is not defined. ")
1473 : END IF
1474 : END DO
1475 : END IF
1476 :
1477 9247 : lribas = .FALSE.
1478 9247 : e1terms = .FALSE.
1479 9247 : IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN
1480 42 : lribas = .TRUE.
1481 42 : CALL get_qs_env(qs_env, lri_env=lri_env)
1482 42 : e1terms = lri_env%exact_1c_terms
1483 : END IF
1484 9247 : IF (dft_control%qs_control%do_kg) THEN
1485 82 : CALL section_vals_val_get(dft_section, "KG_METHOD%TNADD_METHOD", i_val=tnadd_method)
1486 82 : IF (tnadd_method == kg_tnadd_embed_ri) lribas = .TRUE.
1487 : END IF
1488 9241 : IF (lribas) THEN
1489 : ! Check if LRI_AUX basis is available, auto-generate if needed
1490 48 : CALL get_qs_env(qs_env, nkind=nkind)
1491 138 : DO ikind = 1, nkind
1492 90 : NULLIFY (lri_aux_basis)
1493 90 : qs_kind => qs_kind_set(ikind)
1494 90 : CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX")
1495 138 : IF (.NOT. (ASSOCIATED(lri_aux_basis))) THEN
1496 : ! LRI_AUX basis set is not yet loaded
1497 : CALL cp_warn(__LOCATION__, "Automatic Generation of LRI_AUX basis. "// &
1498 28 : "This is experimental code.")
1499 : ! Generate a default basis
1500 28 : CALL create_lri_aux_basis_set(lri_aux_basis, qs_kind, dft_control%auto_basis_lri_aux, e1terms)
1501 28 : CALL add_basis_set_to_container(qs_kind%basis_sets, lri_aux_basis, "LRI_AUX")
1502 : END IF
1503 : END DO
1504 : END IF
1505 :
1506 9247 : CALL section_vals_val_get(qs_env%input, "DFT%XC%HF%RI%_SECTION_PARAMETERS_", l_val=do_ri_hfx)
1507 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF%RI%_SECTION_PARAMETERS_", &
1508 9247 : l_val=do_rpa_ri_exx)
1509 9247 : IF (do_ri_hfx .OR. do_rpa_ri_exx) THEN
1510 114 : CALL get_qs_env(qs_env, nkind=nkind)
1511 114 : CALL section_vals_val_get(qs_env%input, "DFT%SORT_BASIS", i_val=sort_basis)
1512 306 : DO ikind = 1, nkind
1513 192 : NULLIFY (ri_hfx_basis)
1514 192 : qs_kind => qs_kind_set(ikind)
1515 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=ri_hfx_basis, &
1516 192 : basis_type="RI_HFX")
1517 9439 : IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN
1518 186 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
1519 186 : IF (dft_control%do_admm) THEN
1520 : CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hfx, &
1521 62 : basis_type="AUX_FIT", basis_sort=sort_basis)
1522 : ELSE
1523 : CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hfx, &
1524 124 : basis_sort=sort_basis)
1525 : END IF
1526 186 : CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HFX")
1527 : END IF
1528 : END DO
1529 : END IF
1530 :
1531 9247 : IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
1532 : ! Check if RI_HXC basis is available, auto-generate if needed
1533 2 : CALL get_qs_env(qs_env, nkind=nkind)
1534 4 : DO ikind = 1, nkind
1535 2 : NULLIFY (ri_hfx_basis)
1536 2 : qs_kind => qs_kind_set(ikind)
1537 2 : CALL get_qs_kind(qs_kind, basis_set=ri_hfx_basis, basis_type="RI_HXC")
1538 4 : IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN
1539 : ! Generate a default basis
1540 2 : CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hxc)
1541 2 : CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HXC")
1542 : END IF
1543 : END DO
1544 : END IF
1545 :
1546 : ! Harris method
1547 9247 : NULLIFY (harris_env)
1548 : CALL section_vals_val_get(dft_section, "HARRIS_METHOD%_SECTION_PARAMETERS_", &
1549 9247 : l_val=qs_env%harris_method)
1550 9247 : harris_section => section_vals_get_subs_vals(dft_section, "HARRIS_METHOD")
1551 9247 : CALL harris_env_create(qs_env, harris_env, harris_section)
1552 9247 : CALL set_qs_env(qs_env, harris_env=harris_env)
1553 : !
1554 9247 : IF (qs_env%harris_method .AND. qs_env%harris_env%density_source == hden_atomic) THEN
1555 8 : CALL get_qs_env(qs_env, nkind=nkind)
1556 : ! Check if RI_HXC basis is available, auto-generate if needed
1557 30 : DO ikind = 1, nkind
1558 22 : NULLIFY (tmp_basis_set)
1559 22 : qs_kind => qs_kind_set(ikind)
1560 22 : CALL get_qs_kind(qs_kind, basis_set=rhoin_basis, basis_type="RHOIN")
1561 30 : IF (.NOT. (ASSOCIATED(rhoin_basis))) THEN
1562 : ! Generate a default basis
1563 22 : CALL create_ri_aux_basis_set(tmp_basis_set, qs_kind, dft_control%auto_basis_ri_hxc)
1564 22 : CALL create_primitive_basis_set(tmp_basis_set, rhoin_basis, lmax=0)
1565 22 : CALL deallocate_gto_basis_set(tmp_basis_set)
1566 22 : CALL add_basis_set_to_container(qs_kind%basis_sets, rhoin_basis, "RHOIN")
1567 : END IF
1568 : END DO
1569 : END IF
1570 :
1571 9247 : mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
1572 9247 : CALL section_vals_get(mp2_section, explicit=mp2_present)
1573 9247 : IF (mp2_present) THEN
1574 :
1575 : ! basis should be sorted for imaginary time RPA/GW
1576 512 : CALL section_vals_val_get(qs_env%input, "DFT%SORT_BASIS", i_val=sort_basis)
1577 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%_SECTION_PARAMETERS_", &
1578 512 : l_val=do_wfc_im_time)
1579 :
1580 512 : IF (do_wfc_im_time .AND. sort_basis /= basis_sort_zet) THEN
1581 : CALL cp_warn(__LOCATION__, &
1582 10 : "Low-scaling RPA requires SORT_BASIS EXP keyword (in DFT input section) for good performance")
1583 : END IF
1584 :
1585 : ! Check if RI_AUX basis (for MP2/RPA) is given, auto-generate if not
1586 512 : CALL mp2_env_create(qs_env%mp2_env)
1587 512 : CALL get_qs_env(qs_env, mp2_env=mp2_env, nkind=nkind)
1588 512 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_MP2%_SECTION_PARAMETERS_", l_val=do_ri_mp2)
1589 512 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_SOS_MP2%_SECTION_PARAMETERS_", l_val=do_ri_sos_mp2)
1590 512 : CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%_SECTION_PARAMETERS_", l_val=do_ri_rpa)
1591 512 : IF (do_ri_mp2 .OR. do_ri_sos_mp2 .OR. do_ri_rpa) THEN
1592 1376 : DO ikind = 1, nkind
1593 902 : NULLIFY (ri_aux_basis_set)
1594 902 : qs_kind => qs_kind_set(ikind)
1595 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=ri_aux_basis_set, &
1596 902 : basis_type="RI_AUX")
1597 1414 : IF (.NOT. (ASSOCIATED(ri_aux_basis_set))) THEN
1598 : ! RI_AUX basis set is not yet loaded
1599 : ! Generate a default basis
1600 8 : CALL create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, dft_control%auto_basis_ri_aux, basis_sort=sort_basis)
1601 8 : CALL add_basis_set_to_container(qs_kind%basis_sets, ri_aux_basis_set, "RI_AUX")
1602 : ! Add a flag, which allows to check if the basis was generated
1603 : ! when applying ERI_METHOD OS to mp2, ri-rpa, gw etc
1604 8 : qs_env%mp2_env%ri_aux_auto_generated = .TRUE.
1605 : END IF
1606 : END DO
1607 : END IF
1608 :
1609 : END IF
1610 :
1611 9247 : IF (dft_control%do_xas_tdp_calculation .OR. qs_env%do_rixs) THEN
1612 : ! Check if RI_XAS basis is given, auto-generate if not
1613 68 : CALL get_qs_env(qs_env, nkind=nkind)
1614 178 : DO ikind = 1, nkind
1615 110 : NULLIFY (ri_xas_basis)
1616 110 : qs_kind => qs_kind_set(ikind)
1617 110 : CALL get_qs_kind(qs_kind, basis_Set=ri_xas_basis, basis_type="RI_XAS")
1618 9357 : IF (.NOT. ASSOCIATED(ri_xas_basis)) THEN
1619 : ! Generate a default basis
1620 106 : CALL create_ri_aux_basis_set(ri_xas_basis, qs_kind, dft_control%auto_basis_ri_xas)
1621 106 : CALL add_basis_set_to_container(qs_kind%basis_sets, ri_xas_basis, "RI_XAS")
1622 : END IF
1623 : END DO
1624 : END IF
1625 :
1626 : ! Initialize the spherical harmonics and the orbital transformation matrices
1627 9247 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, maxlppl=maxlppl, maxlppnl=maxlppnl)
1628 :
1629 : ! CNEO nuclear basis contributes to GAPW rho0
1630 9247 : IF (cneo_potential_present) THEN
1631 8 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_nuc, basis_type="NUC")
1632 8 : maxlgto = MAX(maxlgto, maxlgto_nuc)
1633 : END IF
1634 9247 : lmax_sphere = dft_control%qs_control%gapw_control%lmax_sphere
1635 9247 : IF (lmax_sphere < 0) THEN
1636 9089 : lmax_sphere = 2*maxlgto
1637 9089 : dft_control%qs_control%gapw_control%lmax_sphere = lmax_sphere
1638 : END IF
1639 9247 : IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
1640 48 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="LRI_AUX")
1641 : !take maxlgto from lri basis if larger (usually)
1642 48 : maxlgto = MAX(maxlgto, maxlgto_lri)
1643 9199 : ELSE IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
1644 2 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_HXC")
1645 2 : maxlgto = MAX(maxlgto, maxlgto_lri)
1646 : END IF
1647 9247 : IF (dft_control%do_xas_tdp_calculation .OR. qs_env%do_rixs) THEN
1648 : !done as a precaution
1649 68 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_XAS")
1650 68 : maxlgto = MAX(maxlgto, maxlgto_lri)
1651 : END IF
1652 9247 : maxl = MAX(2*maxlgto, maxlppl, maxlppnl, lmax_sphere) + 1
1653 :
1654 9247 : CALL init_orbital_pointers(maxl)
1655 :
1656 9247 : CALL init_spherical_harmonics(maxl, 0)
1657 :
1658 : ! Initialise the qs_kind_set
1659 9247 : CALL init_qs_kind_set(qs_kind_set)
1660 :
1661 : ! Initialise GAPW soft basis and projectors
1662 9247 : IF (dft_control%qs_control%method_id == do_method_gapw .OR. &
1663 : dft_control%qs_control%method_id == do_method_gapw_xc) THEN
1664 1544 : qs_control => dft_control%qs_control
1665 1544 : CALL init_gapw_basis_set(qs_kind_set, qs_control, qs_env%input)
1666 : END IF
1667 :
1668 : ! Initialise CNEO nuclear soft basis
1669 9247 : IF (cneo_potential_present) THEN
1670 8 : CALL init_cneo_basis_set(qs_kind_set, qs_control)
1671 : END IF
1672 :
1673 : ! Initialize the pretabulation for the calculation of the
1674 : ! incomplete Gamma function F_n(t) after McMurchie-Davidson
1675 9247 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
1676 9247 : maxl = MAX(3*maxlgto + 1, 0)
1677 9247 : CALL init_md_ftable(maxl)
1678 :
1679 : ! Initialize the atomic interaction radii
1680 9247 : CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
1681 : !
1682 9247 : IF (dft_control%qs_control%method_id == do_method_xtb) THEN
1683 1232 : IF (.NOT. dft_control%qs_control%xtb_control%do_tblite) THEN
1684 : ! cutoff radius
1685 1034 : CALL get_qs_env(qs_env, nkind=nkind)
1686 3314 : DO ikind = 1, nkind
1687 2280 : qs_kind => qs_kind_set(ikind)
1688 3314 : IF (qs_kind%xtb_parameter%defined) THEN
1689 2278 : CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set)
1690 2278 : rcut = xtb_control%coulomb_sr_cut
1691 2278 : fxx = 2.0_dp*xtb_control%coulomb_sr_eps*qs_kind%xtb_parameter%eta**2
1692 2278 : fxx = 0.80_dp*(1.0_dp/fxx)**0.3333_dp
1693 2278 : rcut = MIN(rcut, xtb_control%coulomb_sr_cut)
1694 2278 : qs_kind%xtb_parameter%rcut = MIN(rcut, fxx)
1695 : ELSE
1696 2 : qs_kind%xtb_parameter%rcut = 0.0_dp
1697 : END IF
1698 : END DO
1699 : END IF
1700 : END IF
1701 :
1702 9247 : IF (.NOT. be_silent) THEN
1703 9219 : CALL write_pgf_orb_radii("orb", atomic_kind_set, qs_kind_set, subsys_section)
1704 9219 : CALL write_pgf_orb_radii("aux", atomic_kind_set, qs_kind_set, subsys_section)
1705 9219 : CALL write_pgf_orb_radii("lri", atomic_kind_set, qs_kind_set, subsys_section)
1706 9219 : CALL write_pgf_orb_radii("nuc", atomic_kind_set, qs_kind_set, subsys_section)
1707 9219 : CALL write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section)
1708 9219 : CALL write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section)
1709 9219 : CALL write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section)
1710 9219 : CALL write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section)
1711 : END IF
1712 :
1713 : ! Distribute molecules and atoms using the new data structures
1714 : CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, &
1715 : particle_set=particle_set, &
1716 : local_particles=local_particles, &
1717 : molecule_kind_set=molecule_kind_set, &
1718 : molecule_set=molecule_set, &
1719 : local_molecules=local_molecules, &
1720 9247 : force_env_section=qs_env%input)
1721 :
1722 : ! SCF parameters
1723 286657 : ALLOCATE (scf_control)
1724 : ! set (non)-self consistency
1725 9247 : IF (dft_control%qs_control%dftb) THEN
1726 298 : scf_control%non_selfconsistent = .NOT. dft_control%qs_control%dftb_control%self_consistent
1727 : END IF
1728 9247 : IF (dft_control%qs_control%xtb) THEN
1729 1232 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
1730 198 : scf_control%non_selfconsistent = .FALSE.
1731 : ELSE
1732 1034 : scf_control%non_selfconsistent = (dft_control%qs_control%xtb_control%gfn_type == 0)
1733 : END IF
1734 : END IF
1735 9247 : IF (qs_env%harris_method) THEN
1736 28 : scf_control%non_selfconsistent = .TRUE.
1737 : END IF
1738 9247 : CALL scf_c_create(scf_control)
1739 9247 : CALL scf_c_read_parameters(scf_control, dft_section)
1740 9247 : IF (scf_control%gce%do_gce) THEN
1741 32 : IF (.NOT. ALL(cell%perd == 1)) THEN
1742 0 : CPABORT("Grand canonical SCF is only implemented for 3D periodic calculations.")
1743 : END IF
1744 8 : IF (.NOT. scf_control%smear%do_smear) THEN
1745 0 : CPABORT("Grand canonical SCF requires smearing.")
1746 : END IF
1747 8 : IF (scf_control%smear%method == smear_energy_window .OR. &
1748 : scf_control%smear%method == smear_list) THEN
1749 : CALL cp_abort(__LOCATION__, &
1750 : "Grand canonical SCF requires Fermi-Dirac, Gaussian, "// &
1751 0 : "Methfessel-Paxton or Marzari-Vanderbilt smearing.")
1752 : END IF
1753 8 : IF (scf_control%use_ot .OR. .NOT. scf_control%use_diag .OR. &
1754 : scf_control%diagonalization%method == diag_ot) THEN
1755 : CALL cp_abort(__LOCATION__, &
1756 : "Grand canonical SCF requires standard diagonalization. "// &
1757 0 : "It is not implemented with OT.")
1758 : END IF
1759 : END IF
1760 9247 : IF (.NOT. dft_control%qs_control%do_ls_scf) THEN
1761 9149 : SELECT CASE (dft_control%qs_control%method_id)
1762 : CASE (do_method_dftb)
1763 254 : IF (dft_control%qs_control%dftb_control%tblite_scc_mixer == tblite_scc_mixer_tblite) THEN
1764 2 : scf_control%max_scf = dft_control%qs_control%dftb_control%tblite_mixer_iterations
1765 : END IF
1766 : CASE (do_method_xtb)
1767 8895 : IF (dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_tblite) THEN
1768 34 : scf_control%max_scf = dft_control%qs_control%xtb_control%tblite_mixer_iterations
1769 : END IF
1770 : END SELECT
1771 : END IF
1772 :
1773 : ! Allocate the data structure for Quickstep energies
1774 9247 : CALL allocate_qs_energy(energy)
1775 :
1776 : ! Check for orthogonal basis
1777 9247 : has_unit_metric = .FALSE.
1778 9247 : IF (dft_control%qs_control%semi_empirical) THEN
1779 998 : IF (dft_control%qs_control%se_control%orthogonal_basis) has_unit_metric = .TRUE.
1780 : END IF
1781 9247 : IF (dft_control%qs_control%dftb) THEN
1782 298 : IF (dft_control%qs_control%dftb_control%orthogonal_basis) has_unit_metric = .TRUE.
1783 : END IF
1784 9247 : CALL set_qs_env(qs_env, has_unit_metric=has_unit_metric)
1785 :
1786 : ! Select the MTLR SCF initialization after reading SCF_GUESS.
1787 9247 : IF (dft_control%mtlr_u_j) THEN
1788 28 : SELECT CASE (scf_control%density_guess)
1789 : CASE (atomic_guess)
1790 12 : IF (dft_control%mtlr_reference_scf_explicit .AND. &
1791 4 : dft_control%mtlr_reference_scf) THEN
1792 8 : dft_control%mtlr_initialization_mode = mtlr_reference_from_atomic
1793 : ELSE
1794 4 : dft_control%mtlr_initialization_mode = mtlr_atomic_perturbations
1795 : END IF
1796 : CASE (restart_guess)
1797 4 : IF (dft_control%mtlr_reference_scf_explicit .AND. &
1798 : .NOT. dft_control%mtlr_reference_scf) THEN
1799 : CALL cp_abort(__LOCATION__, &
1800 : "MTLR_REFERENCE_SCF OFF is incompatible with SCF_GUESS RESTART. "// &
1801 0 : "Use MTLR_REFERENCE_SCF ON or SCF_GUESS ATOMIC.")
1802 : END IF
1803 4 : dft_control%mtlr_initialization_mode = mtlr_reference_from_restart
1804 : CASE DEFAULT
1805 16 : CPABORT("MTLR requires SCF_GUESS RESTART or SCF_GUESS ATOMIC.")
1806 : END SELECT
1807 :
1808 : dft_control%mtlr_reference_scf = &
1809 16 : dft_control%mtlr_initialization_mode /= mtlr_atomic_perturbations
1810 :
1811 4 : SELECT CASE (dft_control%mtlr_initialization_mode)
1812 : CASE (mtlr_atomic_perturbations)
1813 : ! Initialize every perturbation SCF from the atomic guess.
1814 4 : dft_control%qs_control%wf_interpolation_method_nr = wfi_use_guess_method_nr
1815 : CASE (mtlr_reference_from_atomic, mtlr_reference_from_restart)
1816 : ! Initialize every perturbation SCF from the restored reference MOs.
1817 12 : dft_control%qs_control%wf_interpolation_method_nr = wfi_use_prev_wf_method_nr
1818 : CASE DEFAULT
1819 16 : CPABORT("The MTLR SCF initialization mode was not resolved.")
1820 : END SELECT
1821 : END IF
1822 :
1823 : ! Activate the interpolation
1824 : CALL wfi_create(wf_history, &
1825 : interpolation_method_nr= &
1826 : dft_control%qs_control%wf_interpolation_method_nr, &
1827 : extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
1828 9247 : has_unit_metric=has_unit_metric)
1829 :
1830 : ! Set the current Quickstep environment
1831 : CALL set_qs_env(qs_env=qs_env, &
1832 : scf_control=scf_control, &
1833 9247 : wf_history=wf_history)
1834 :
1835 : CALL qs_subsys_set(subsys, &
1836 : cell_ref=cell_ref, &
1837 : use_ref_cell=use_ref_cell, &
1838 : energy=energy, &
1839 9247 : force=force)
1840 :
1841 9247 : CALL get_qs_env(qs_env, ks_env=ks_env)
1842 9247 : CALL set_ks_env(ks_env, dft_control=dft_control)
1843 :
1844 : CALL qs_subsys_set(subsys, local_molecules=local_molecules, &
1845 9247 : local_particles=local_particles, cell=cell)
1846 :
1847 9247 : CALL distribution_1d_release(local_particles)
1848 9247 : CALL distribution_1d_release(local_molecules)
1849 9247 : CALL wfi_release(wf_history)
1850 :
1851 : CALL get_qs_env(qs_env=qs_env, &
1852 : atomic_kind_set=atomic_kind_set, &
1853 : dft_control=dft_control, &
1854 9247 : scf_control=scf_control)
1855 :
1856 : ! Decide what conditions need mo_derivs
1857 : ! right now, this only appears to be OT
1858 9247 : IF (dft_control%qs_control%do_ls_scf .OR. &
1859 : dft_control%qs_control%do_almo_scf) THEN
1860 424 : CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.FALSE.)
1861 : ELSE
1862 8823 : IF (scf_control%use_ot) THEN
1863 2491 : CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.TRUE.)
1864 : ELSE
1865 6332 : CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.FALSE.)
1866 : END IF
1867 : END IF
1868 :
1869 : ! XXXXXXX this is backwards XXXXXXXX
1870 9247 : IF (dft_control%qs_control%xtb_control%do_tblite .AND. .NOT. scf_control%use_ot) THEN
1871 180 : IF (.NOT. scf_control%smear%do_smear) THEN
1872 : ! set tblite default smearing
1873 120 : scf_control%smear%do_smear = .TRUE.
1874 120 : scf_control%smear%method = smear_fermi_dirac
1875 120 : scf_control%smear%electronic_temperature = 300._dp/kelvin
1876 120 : scf_control%smear%eps_fermi_dirac = 1.E-6_dp
1877 : END IF
1878 : END IF
1879 9247 : dft_control%smear = scf_control%smear%do_smear
1880 :
1881 : ! Periodic efield needs equal occupation and orbital gradients
1882 9247 : IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)) THEN
1883 7717 : IF (dft_control%apply_period_efield) THEN
1884 30 : CALL get_qs_env(qs_env=qs_env, requires_mo_derivs=orb_gradient)
1885 30 : IF (.NOT. orb_gradient) THEN
1886 : CALL cp_abort(__LOCATION__, "Periodic Efield needs orbital gradient and direct optimization."// &
1887 0 : " Use the OT optimization method.")
1888 : END IF
1889 30 : IF (dft_control%smear) THEN
1890 : CALL cp_abort(__LOCATION__, "Periodic Efield needs equal occupation numbers."// &
1891 0 : " Smearing option is not possible.")
1892 : END IF
1893 : END IF
1894 : END IF
1895 :
1896 : ! Initialize the GAPW local densities and potentials
1897 9247 : IF (dft_control%qs_control%method_id == do_method_gapw .OR. &
1898 : dft_control%qs_control%method_id == do_method_gapw_xc) THEN
1899 : ! Allocate and initialize the set of atomic densities
1900 1544 : NULLIFY (rho_atom_set)
1901 1544 : gapw_control => dft_control%qs_control%gapw_control
1902 1544 : CALL init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1903 1544 : CALL set_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
1904 1544 : IF (dft_control%qs_control%method_id /= do_method_gapw_xc) THEN
1905 1354 : CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, natom=natom)
1906 : ! Allocate and initialize the compensation density rho0
1907 1354 : CALL init_rho0(local_rho_set, qs_env, gapw_control)
1908 : ! Allocate and Initialize the local coulomb term
1909 1354 : CALL init_coulomb_local(qs_env%hartree_local, natom)
1910 : END IF
1911 : ! NLCC
1912 1544 : CALL init_gapw_nlcc(qs_kind_set)
1913 : ! Accurate XC integration
1914 1544 : IF (gapw_control%accurate_xcint) THEN
1915 288 : CPASSERT(gapw_control%oweights >= 0)
1916 288 : CPASSERT(.NOT. ASSOCIATED(gapw_control%aw))
1917 288 : CALL get_qs_env(qs_env, nkind=nkind)
1918 864 : ALLOCATE (gapw_control%aw(nkind))
1919 288 : alpha = gapw_control%aweights
1920 816 : DO ikind = 1, nkind
1921 528 : qs_kind => qs_kind_set(ikind)
1922 528 : CALL get_qs_kind(qs_kind, hard_radius=rc, paw_atom=paw_atom)
1923 816 : IF (paw_atom) THEN
1924 510 : gapw_control%aw(ikind) = alpha*(1.2_dp/rc)**2
1925 : ELSE
1926 18 : gapw_control%aw(ikind) = 0.0_dp
1927 : END IF
1928 : END DO
1929 : END IF
1930 7703 : ELSE IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN
1931 : ! allocate local ri environment
1932 : ! nothing to do here?
1933 7661 : ELSE IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
1934 : ! allocate ri environment
1935 : ! nothing to do here?
1936 7659 : ELSE IF (dft_control%qs_control%semi_empirical) THEN
1937 998 : NULLIFY (se_store_int_env, se_nddo_mpole, se_nonbond_env)
1938 998 : natom = SIZE(particle_set)
1939 998 : se_section => section_vals_get_subs_vals(qs_section, "SE")
1940 998 : se_control => dft_control%qs_control%se_control
1941 :
1942 : ! Make the cutoff radii choice a bit smarter
1943 998 : CALL se_cutoff_compatible(se_control, se_section, cell, output_unit)
1944 :
1945 1994 : SELECT CASE (dft_control%qs_control%method_id)
1946 : CASE DEFAULT
1947 : CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pm3, &
1948 : do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
1949 : ! Neighbor lists have to be MAX(interaction range, orbital range)
1950 : ! set new kind radius
1951 998 : CALL init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, subsys_section)
1952 : END SELECT
1953 : ! Initialize to zero the max multipole to treat in the EWALD scheme..
1954 998 : se_control%max_multipole = do_multipole_none
1955 : ! check for Ewald
1956 998 : IF (se_control%do_ewald .OR. se_control%do_ewald_gks) THEN
1957 512 : ALLOCATE (ewald_env)
1958 32 : CALL ewald_env_create(ewald_env, para_env)
1959 32 : poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
1960 32 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1961 32 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1962 : print_section => section_vals_get_subs_vals(qs_env%input, &
1963 32 : "PRINT%GRID_INFORMATION")
1964 32 : CALL read_ewald_section(ewald_env, ewald_section)
1965 : ! Create ewald grids
1966 32 : ALLOCATE (ewald_pw)
1967 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, &
1968 32 : print_section=print_section)
1969 : ! Initialize ewald grids
1970 32 : CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
1971 : ! Setup the nonbond environment (real space part of Ewald)
1972 32 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
1973 : ! Setup the maximum level of multipoles to be treated in the periodic SE scheme
1974 32 : IF (se_control%do_ewald) THEN
1975 30 : CALL ewald_env_get(ewald_env, max_multipole=se_control%max_multipole)
1976 : END IF
1977 : CALL section_vals_val_get(se_section, "NEIGHBOR_LISTS%VERLET_SKIN", &
1978 32 : r_val=verlet_skin)
1979 32 : ALLOCATE (se_nonbond_env)
1980 : CALL fist_nonbond_env_create(se_nonbond_env, atomic_kind_set, do_nonbonded=.TRUE., &
1981 : do_electrostatics=.TRUE., verlet_skin=verlet_skin, ewald_rcut=ewald_rcut, &
1982 32 : ei_scale14=0.0_dp, vdw_scale14=0.0_dp, shift_cutoff=.FALSE.)
1983 : ! Create and Setup NDDO multipole environment
1984 32 : CALL nddo_mpole_setup(se_nddo_mpole, natom)
1985 : CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw, &
1986 32 : se_nonbond_env=se_nonbond_env, se_nddo_mpole=se_nddo_mpole)
1987 : ! Handle the residual integral part 1/R^3
1988 : CALL semi_empirical_expns3_setup(qs_kind_set, se_control, &
1989 32 : dft_control%qs_control%method_id)
1990 : END IF
1991 : ! Taper function
1992 : CALL se_taper_create(se_taper, se_control%integral_screening, se_control%do_ewald, &
1993 : se_control%taper_cou, se_control%range_cou, &
1994 : se_control%taper_exc, se_control%range_exc, &
1995 : se_control%taper_scr, se_control%range_scr, &
1996 998 : se_control%taper_lrc, se_control%range_lrc)
1997 998 : CALL set_qs_env(qs_env, se_taper=se_taper)
1998 : ! Store integral environment
1999 998 : CALL semi_empirical_si_create(se_store_int_env, se_section)
2000 998 : CALL set_qs_env(qs_env, se_store_int_env=se_store_int_env)
2001 : END IF
2002 :
2003 : ! Initialize possible dispersion parameters
2004 : IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
2005 : dft_control%qs_control%method_id == do_method_gapw .OR. &
2006 : dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
2007 : dft_control%qs_control%method_id == do_method_lrigpw .OR. &
2008 9247 : dft_control%qs_control%method_id == do_method_rigpw .OR. &
2009 : dft_control%qs_control%method_id == do_method_ofgpw) THEN
2010 33595 : ALLOCATE (dispersion_env)
2011 6719 : NULLIFY (xc_section)
2012 6719 : xc_section => section_vals_get_subs_vals(dft_section, "XC")
2013 6719 : CALL qs_dispersion_env_set(dispersion_env, xc_section)
2014 6719 : IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
2015 242 : NULLIFY (pp_section)
2016 242 : pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL")
2017 242 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
2018 6477 : ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
2019 62 : NULLIFY (nl_section)
2020 62 : nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
2021 62 : CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
2022 : END IF
2023 6719 : CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
2024 2528 : ELSE IF (dft_control%qs_control%method_id == do_method_dftb) THEN
2025 1490 : ALLOCATE (dispersion_env)
2026 : ! set general defaults
2027 : dispersion_env%doabc = .FALSE.
2028 : dispersion_env%c9cnst = .FALSE.
2029 : dispersion_env%lrc = .FALSE.
2030 : dispersion_env%srb = .FALSE.
2031 : dispersion_env%verbose = .FALSE.
2032 : NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
2033 : dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
2034 : dispersion_env%d3_exclude_pair)
2035 : NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel_table, &
2036 : dispersion_env%d2y_dx2, dispersion_env%dftd_section)
2037 : NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
2038 298 : IF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d3) THEN
2039 14 : dispersion_env%type = xc_vdw_fun_pairpot
2040 14 : dispersion_env%pp_type = vdw_pairpot_dftd3
2041 14 : dispersion_env%eps_cn = dftb_control%epscn
2042 14 : dispersion_env%s6 = dftb_control%sd3(1)
2043 14 : dispersion_env%sr6 = dftb_control%sd3(2)
2044 14 : dispersion_env%s8 = dftb_control%sd3(3)
2045 14 : dispersion_env%domol = .FALSE.
2046 14 : dispersion_env%kgc8 = 0._dp
2047 14 : dispersion_env%rc_disp = dftb_control%rcdisp
2048 14 : dispersion_env%exp_pre = 0._dp
2049 14 : dispersion_env%scaling = 0._dp
2050 14 : dispersion_env%nd3_exclude_pair = 0
2051 14 : dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
2052 14 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
2053 284 : ELSE IF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d3bj) THEN
2054 2 : dispersion_env%type = xc_vdw_fun_pairpot
2055 2 : dispersion_env%pp_type = vdw_pairpot_dftd3bj
2056 2 : dispersion_env%eps_cn = dftb_control%epscn
2057 2 : dispersion_env%s6 = dftb_control%sd3bj(1)
2058 2 : dispersion_env%a1 = dftb_control%sd3bj(2)
2059 2 : dispersion_env%s8 = dftb_control%sd3bj(3)
2060 2 : dispersion_env%a2 = dftb_control%sd3bj(4)
2061 2 : dispersion_env%domol = .FALSE.
2062 2 : dispersion_env%kgc8 = 0._dp
2063 2 : dispersion_env%rc_disp = dftb_control%rcdisp
2064 2 : dispersion_env%exp_pre = 0._dp
2065 2 : dispersion_env%scaling = 0._dp
2066 2 : dispersion_env%nd3_exclude_pair = 0
2067 2 : dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
2068 2 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
2069 282 : ELSE IF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d2) THEN
2070 2 : dispersion_env%type = xc_vdw_fun_pairpot
2071 2 : dispersion_env%pp_type = vdw_pairpot_dftd2
2072 2 : dispersion_env%exp_pre = dftb_control%exp_pre
2073 2 : dispersion_env%scaling = dftb_control%scaling
2074 2 : dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
2075 2 : dispersion_env%rc_disp = dftb_control%rcdisp
2076 2 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
2077 : ELSE
2078 280 : dispersion_env%type = xc_vdw_fun_none
2079 : END IF
2080 298 : CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
2081 2230 : ELSE IF (dft_control%qs_control%method_id == do_method_xtb) THEN
2082 1232 : IF (.NOT. (dft_control%qs_control%xtb_control%do_tblite)) THEN
2083 5170 : ALLOCATE (dispersion_env)
2084 : ! set general defaults
2085 : dispersion_env%doabc = .FALSE.
2086 : dispersion_env%c9cnst = .FALSE.
2087 : dispersion_env%lrc = .FALSE.
2088 : dispersion_env%srb = .FALSE.
2089 : dispersion_env%verbose = .FALSE.
2090 : NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, &
2091 : dispersion_env%r0ab, dispersion_env%rcov, &
2092 : dispersion_env%r2r4, dispersion_env%cn, &
2093 : dispersion_env%cnkind, dispersion_env%cnlist, &
2094 : dispersion_env%d3_exclude_pair)
2095 : NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel_table, &
2096 : dispersion_env%d2y_dx2, dispersion_env%dftd_section)
2097 : NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
2098 1034 : dispersion_env%type = xc_vdw_fun_pairpot
2099 1034 : dispersion_env%eps_cn = xtb_control%epscn
2100 1034 : dispersion_env%s6 = xtb_control%s6
2101 1034 : dispersion_env%s8 = xtb_control%s8
2102 1034 : dispersion_env%a1 = xtb_control%a1
2103 1034 : dispersion_env%a2 = xtb_control%a2
2104 1034 : dispersion_env%domol = .FALSE.
2105 1034 : dispersion_env%kgc8 = 0._dp
2106 1034 : dispersion_env%rc_disp = xtb_control%rcdisp
2107 1034 : dispersion_env%rc_d4 = xtb_control%rcdisp
2108 1034 : dispersion_env%exp_pre = 0._dp
2109 1034 : dispersion_env%scaling = 0._dp
2110 1034 : dispersion_env%nd3_exclude_pair = 0
2111 1034 : dispersion_env%parameter_file_name = xtb_control%dispersion_parameter_file
2112 : !
2113 1426 : SELECT CASE (xtb_control%vdw_type)
2114 : CASE (xtb_vdw_type_none, xtb_vdw_type_d3)
2115 392 : dispersion_env%pp_type = vdw_pairpot_dftd3bj
2116 392 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
2117 392 : IF (xtb_control%vdw_type == xtb_vdw_type_none) dispersion_env%type = xc_vdw_fun_none
2118 : CASE (xtb_vdw_type_d4)
2119 642 : dispersion_env%pp_type = vdw_pairpot_dftd4
2120 642 : dispersion_env%ref_functional = "none"
2121 : ! GFN0 supplies its EEQ charges directly to D4. GFN1 does not, so propagate the
2122 : ! common EEQ solver settings for the D4 charge model used in that case.
2123 642 : dispersion_env%eeq_sparam = xtb_control%eeq_sparam
2124 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, &
2125 642 : dispersion_env, para_env=para_env)
2126 642 : dispersion_env%cnfun = 2
2127 : CASE DEFAULT
2128 1034 : CPABORT("vdw type")
2129 : END SELECT
2130 1034 : CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
2131 : END IF
2132 998 : ELSE IF (dft_control%qs_control%semi_empirical) THEN
2133 4990 : ALLOCATE (dispersion_env)
2134 : ! set general defaults
2135 : dispersion_env%doabc = .FALSE.
2136 : dispersion_env%c9cnst = .FALSE.
2137 : dispersion_env%lrc = .FALSE.
2138 : dispersion_env%srb = .FALSE.
2139 : dispersion_env%verbose = .FALSE.
2140 : NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
2141 : dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
2142 : dispersion_env%d3_exclude_pair)
2143 : NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel_table, &
2144 : dispersion_env%d2y_dx2, dispersion_env%dftd_section)
2145 : NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
2146 998 : IF (se_control%dispersion) THEN
2147 6 : dispersion_env%type = xc_vdw_fun_pairpot
2148 6 : dispersion_env%pp_type = vdw_pairpot_dftd3
2149 6 : dispersion_env%eps_cn = se_control%epscn
2150 6 : dispersion_env%s6 = se_control%sd3(1)
2151 6 : dispersion_env%sr6 = se_control%sd3(2)
2152 6 : dispersion_env%s8 = se_control%sd3(3)
2153 6 : dispersion_env%domol = .FALSE.
2154 6 : dispersion_env%kgc8 = 0._dp
2155 6 : dispersion_env%rc_disp = se_control%rcdisp
2156 6 : dispersion_env%exp_pre = 0._dp
2157 6 : dispersion_env%scaling = 0._dp
2158 6 : dispersion_env%nd3_exclude_pair = 0
2159 6 : dispersion_env%parameter_file_name = se_control%dispersion_parameter_file
2160 6 : CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
2161 : ELSE
2162 992 : dispersion_env%type = xc_vdw_fun_none
2163 : END IF
2164 998 : CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
2165 : END IF
2166 :
2167 : ! Initialize possible geomertical counterpoise correction potential
2168 : IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
2169 : dft_control%qs_control%method_id == do_method_gapw .OR. &
2170 : dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
2171 : dft_control%qs_control%method_id == do_method_lrigpw .OR. &
2172 9247 : dft_control%qs_control%method_id == do_method_rigpw .OR. &
2173 : dft_control%qs_control%method_id == do_method_ofgpw) THEN
2174 6719 : ALLOCATE (gcp_env)
2175 6719 : NULLIFY (xc_section)
2176 6719 : xc_section => section_vals_get_subs_vals(dft_section, "XC")
2177 6719 : CALL qs_gcp_env_set(gcp_env, xc_section)
2178 6719 : CALL qs_gcp_init(qs_env, gcp_env)
2179 6719 : CALL set_qs_env(qs_env, gcp_env=gcp_env)
2180 : END IF
2181 :
2182 : ! Allocate the MO data types
2183 9247 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron)
2184 :
2185 : ! The total number of electrons
2186 9247 : IF (PRESENT(charge)) THEN
2187 44 : dft_control%charge = charge
2188 44 : nelectron = nelectron - dft_control%charge
2189 : ELSE
2190 9203 : nelectron = nelectron - dft_control%charge
2191 : END IF
2192 :
2193 9247 : IF (dft_control%multiplicity == 0) THEN
2194 7595 : IF (MODULO(nelectron, 2) == 0) THEN
2195 7050 : dft_control%multiplicity = 1
2196 : ELSE
2197 545 : dft_control%multiplicity = 2
2198 : END IF
2199 : END IF
2200 :
2201 9247 : multiplicity = dft_control%multiplicity
2202 :
2203 9247 : IF (PRESENT(multip)) THEN
2204 44 : multiplicity = multip
2205 : END IF
2206 :
2207 9247 : IF ((dft_control%nspins < 1) .OR. (dft_control%nspins > 2)) THEN
2208 0 : CPABORT("nspins should be 1 or 2 for the time being ...")
2209 : END IF
2210 :
2211 9247 : IF ((MODULO(nelectron, 2) /= 0) .AND. (dft_control%nspins == 1)) THEN
2212 50 : IF (.NOT. dft_control%qs_control%ofgpw .AND. .NOT. dft_control%smear) THEN
2213 0 : CPABORT("Use the LSD option for an odd number of electrons")
2214 : END IF
2215 : END IF
2216 :
2217 : ! The transition potential method to calculate XAS needs LSD
2218 9247 : IF (dft_control%do_xas_calculation) THEN
2219 42 : IF (dft_control%nspins == 1) THEN
2220 0 : CPABORT("Use the LSD option for XAS with transition potential")
2221 : END IF
2222 : END IF
2223 :
2224 : ! assigning the number of states per spin initial version, not yet very
2225 : ! general. Should work for an even number of electrons and a single
2226 : ! additional electron this set of options that requires full matrices,
2227 : ! however, makes things a bit ugly right now.... we try to make a
2228 : ! distinction between the number of electrons per spin and the number of
2229 : ! MOs per spin this should allow the use of fractional occupations later on
2230 9247 : IF (dft_control%qs_control%ofgpw) THEN
2231 :
2232 0 : IF (dft_control%nspins == 1) THEN
2233 0 : maxocc = nelectron
2234 0 : nelectron_spin(1) = nelectron
2235 0 : nelectron_spin(2) = 0
2236 0 : n_mo(1) = 1
2237 0 : n_mo(2) = 0
2238 : ELSE
2239 0 : nelectron_spin(1) = (nelectron + multiplicity - 1)/2
2240 0 : nelectron_spin(2) = (nelectron - multiplicity + 1)/2
2241 0 : IF (nelectron_spin(1) < 0) THEN
2242 0 : CPABORT("LSD: too few electrons for this multiplicity")
2243 : END IF
2244 0 : maxocc = MAXVAL(nelectron_spin)
2245 0 : n_mo(1) = MIN(nelectron_spin(1), 1)
2246 0 : n_mo(2) = MIN(nelectron_spin(2), 1)
2247 : END IF
2248 :
2249 : ELSE
2250 :
2251 9247 : IF (dft_control%nspins == 1) THEN
2252 7285 : maxocc = 2.0_dp
2253 7285 : nelectron_spin(1) = nelectron
2254 7285 : nelectron_spin(2) = 0
2255 7285 : IF (MODULO(nelectron, 2) == 0) THEN
2256 7235 : n_mo(1) = nelectron/2
2257 : ELSE
2258 50 : n_mo(1) = INT(nelectron/2._dp) + 1
2259 : END IF
2260 7285 : n_mo(2) = 0
2261 : ELSE
2262 1962 : maxocc = 1.0_dp
2263 :
2264 : ! The simplist spin distribution is written here. Special cases will
2265 : ! need additional user input
2266 1962 : IF (MODULO(nelectron + multiplicity - 1, 2) /= 0) THEN
2267 0 : CPABORT("LSD: try to use a different multiplicity")
2268 : END IF
2269 :
2270 1962 : nelectron_spin(1) = (nelectron + multiplicity - 1)/2
2271 1962 : nelectron_spin(2) = (nelectron - multiplicity + 1)/2
2272 :
2273 1962 : IF (nelectron_spin(2) < 0) THEN
2274 0 : CPABORT("LSD: too few electrons for this multiplicity")
2275 : END IF
2276 :
2277 1962 : n_mo(1) = nelectron_spin(1)
2278 1962 : n_mo(2) = nelectron_spin(2)
2279 :
2280 : END IF
2281 :
2282 : END IF
2283 :
2284 : ! Read the total_zeff_corr here [SGh]
2285 9247 : CALL get_qs_kind_set(qs_kind_set, total_zeff_corr=total_zeff_corr)
2286 : ! store it in qs_env
2287 9247 : qs_env%total_zeff_corr = total_zeff_corr
2288 :
2289 : ! Store the number of electrons once and for all
2290 : CALL qs_subsys_set(subsys, &
2291 : nelectron_total=nelectron, &
2292 9247 : nelectron_spin=nelectron_spin)
2293 :
2294 27661 : IF (ANY(scf_control%added_mos_auto) .AND. .NOT. scf_control%smear%do_smear) THEN
2295 0 : CPABORT("ADDED_MOS AUTO requires smearing")
2296 : END IF
2297 27661 : IF (ANY(scf_control%added_mos_auto) .AND. .NOT. do_kpoints) THEN
2298 0 : CPABORT("ADDED_MOS AUTO currently requires a K-point calculation")
2299 : END IF
2300 :
2301 27661 : IF (scf_control%smear%do_smear .AND. ANY(scf_control%added_mos_auto)) THEN
2302 88 : n_mo_auto = MAX(4, CEILING(0.2_dp*REAL(MAXVAL(n_mo(1:dft_control%nspins)), KIND=dp)))
2303 120 : IF (dft_control%nspins == 2 .AND. ALL(scf_control%added_mos_auto(1:2))) THEN
2304 24 : n_mo_target = MIN(n_ao, MAXVAL(n_mo(1:2)) + n_mo_auto)
2305 24 : DO ispin = 1, 2
2306 16 : scf_control%added_mos(ispin) = MAX(0, n_mo_target - n_mo(ispin))
2307 24 : scf_control%added_mos_auto(ispin) = .TRUE.
2308 : END DO
2309 : ELSE
2310 64 : DO ispin = 1, dft_control%nspins
2311 64 : IF (scf_control%added_mos_auto(ispin)) THEN
2312 32 : scf_control%added_mos(ispin) = MIN(n_mo_auto, MAX(0, n_ao - n_mo(ispin)))
2313 32 : scf_control%added_mos_auto(ispin) = .TRUE.
2314 : END IF
2315 : END DO
2316 : END IF
2317 40 : IF (output_unit > 0) THEN
2318 20 : IF (dft_control%nspins == 2) THEN
2319 : WRITE (UNIT=output_unit, FMT="(T2,A,2I5)") &
2320 4 : "SCF smearing: automatically selected ADDED_MOS:", &
2321 16 : scf_control%added_mos(1:dft_control%nspins)
2322 : ELSE
2323 : WRITE (UNIT=output_unit, FMT="(T2,A,I0)") &
2324 16 : "SCF smearing: automatically selected ADDED_MOS: ", scf_control%added_mos(1)
2325 : END IF
2326 : END IF
2327 : END IF
2328 :
2329 : ! Ensure that all orbitals requested for printout are added even
2330 : ! if the keyword ADDED_MOS was not specified or set properly
2331 9247 : mo_index_range => section_get_ivals(dft_section, "PRINT%MO%MO_INDEX_RANGE")
2332 9247 : CPASSERT(ASSOCIATED(mo_index_range))
2333 9279 : IF (ALL(mo_index_range > 0)) THEN
2334 16 : IF (mo_index_range(1) > mo_index_range(2)) THEN
2335 : CALL cp_abort(__LOCATION__, &
2336 : "The upper orbital index ("// &
2337 : TRIM(ADJUSTL(cp_to_string(mo_index_range(2))))// &
2338 : ") of the MO_INDEX_RANGE should be equal or larger "// &
2339 : "than the lower orbital index ("// &
2340 : TRIM(ADJUSTL(cp_to_string(mo_index_range(1))))// &
2341 0 : ") for printout.")
2342 : END IF
2343 : ! Adapt ADDED_MOS automatically if needed for printout
2344 16 : IF (.NOT. scf_control%use_ot) THEN
2345 : scf_control%added_mos(1) = MIN(MAX(scf_control%added_mos(1), &
2346 : mo_index_range(2) - n_mo(1)), &
2347 10 : n_ao - n_mo(1))
2348 10 : IF (dft_control%nspins == 2) THEN
2349 : scf_control%added_mos(2) = MIN(MAX(scf_control%added_mos(2), &
2350 : mo_index_range(2) - n_mo(2)), &
2351 8 : n_ao - n_mo(2))
2352 : END IF
2353 : END IF
2354 9231 : ELSE IF (mo_index_range(2) < 0) THEN
2355 0 : IF (.NOT. scf_control%use_ot) THEN
2356 : ! Add all available orbitals
2357 0 : scf_control%added_mos(1) = n_ao - n_mo(1)
2358 0 : IF (dft_control%nspins == 2) THEN
2359 : ! Ensure the same number for the spin-down (beta) orbitals
2360 0 : scf_control%added_mos(2) = n_ao - n_mo(2)
2361 : END IF
2362 : END IF
2363 : END IF
2364 :
2365 9247 : nlumo_dos = section_get_ival(dft_section, "PRINT%DOS%NLUMO")
2366 9247 : nlumo_molden = section_get_ival(dft_section, "PRINT%MO_MOLDEN%NLUMO")
2367 9247 : nlumo_required = MAX(nlumo_dos, nlumo_molden)
2368 9247 : IF (nlumo_dos == -1 .OR. nlumo_molden == -1) nlumo_required = -1
2369 9247 : IF (.NOT. scf_control%use_ot .AND. nlumo_required /= 0) THEN
2370 10 : IF (nlumo_required == -1) THEN
2371 4 : IF (scf_control%added_mos(1) /= -1 .OR. &
2372 : (dft_control%nspins == 2 .AND. scf_control%added_mos(2) /= -1)) THEN
2373 : CALL cp_warn(__LOCATION__, &
2374 : "NLUMO requested by DOS/PDOS/Molden exceeds SCF%ADDED_MOS. "// &
2375 : "For diagonalization calculations, SCF%ADDED_MOS is "// &
2376 2 : "increased to provide the requested unoccupied orbitals.")
2377 : END IF
2378 4 : scf_control%added_mos(1) = -1
2379 4 : IF (dft_control%nspins == 2) scf_control%added_mos(2) = -1
2380 : ELSE
2381 6 : IF (scf_control%added_mos(1) >= 0 .AND. &
2382 : nlumo_required > scf_control%added_mos(1)) THEN
2383 : CALL cp_warn(__LOCATION__, &
2384 : "NLUMO requested by DOS/PDOS/Molden exceeds SCF%ADDED_MOS. "// &
2385 : "For diagonalization calculations, SCF%ADDED_MOS is "// &
2386 6 : "increased to provide the requested unoccupied orbitals.")
2387 6 : scf_control%added_mos(1) = nlumo_required
2388 : END IF
2389 6 : IF (dft_control%nspins == 2 .AND. scf_control%added_mos(2) > 0 .AND. &
2390 : nlumo_required > scf_control%added_mos(2)) THEN
2391 0 : scf_control%added_mos(2) = nlumo_required
2392 : END IF
2393 : END IF
2394 : END IF
2395 :
2396 9247 : IF (dft_control%nspins == 2) THEN
2397 : ! Check and set number of added (unoccupied) orbitals for beta spin
2398 1962 : IF (scf_control%added_mos(2) < 0) THEN
2399 154 : n_mo_add = n_ao - n_mo(2) ! use all available MOs
2400 1808 : ELSE IF (scf_control%added_mos(2) > 0) THEN
2401 : n_mo_add = scf_control%added_mos(2)
2402 : ELSE
2403 1632 : n_mo_add = scf_control%added_mos(1)
2404 : END IF
2405 1962 : IF (n_mo_add > n_ao - n_mo(2)) THEN
2406 22 : CPWARN("More ADDED_MOs requested for beta spin than available.")
2407 : END IF
2408 1962 : scf_control%added_mos(2) = MIN(n_mo_add, n_ao - n_mo(2))
2409 1962 : n_mo(2) = n_mo(2) + scf_control%added_mos(2)
2410 : END IF
2411 :
2412 : ! proceed alpha orbitals after the beta orbitals; this is essential to avoid
2413 : ! reduction in the number of available unoccupied molecular orbitals.
2414 : ! E.g. n_ao = 10, nelectrons = 10, multiplicity = 3 implies n_mo(1) = 6, n_mo(2) = 4;
2415 : ! added_mos(1:2) = (6,undef) should increase the number of molecular orbitals as
2416 : ! n_mo(1) = min(n_ao, n_mo(1) + added_mos(1)) = 10, n_mo(2) = 10.
2417 : ! However, if we try to proceed alpha orbitals first, this leads us n_mo(1:2) = (10,8)
2418 : ! due to the following assignment instruction above:
2419 : ! IF (scf_control%added_mos(2) > 0) THEN ... ELSE; n_mo_add = scf_control%added_mos(1); END IF
2420 9247 : IF (dft_control%qs_control%xtb_control%do_tblite .AND. .NOT. scf_control%use_ot) THEN
2421 180 : scf_control%added_mos(1) = n_ao - n_mo(1) ! tblite needs all MO's
2422 9067 : ELSE IF (scf_control%added_mos(1) < 0) THEN
2423 800 : scf_control%added_mos(1) = n_ao - n_mo(1) ! use all available MOs
2424 8267 : ELSE IF (scf_control%added_mos(1) > n_ao - n_mo(1)) THEN
2425 : CALL cp_warn(__LOCATION__, &
2426 : "More added MOs requested than available. "// &
2427 : "The full set of unoccupied MOs will be used. "// &
2428 : "Use 'ADDED_MOS -1' to always use all available MOs "// &
2429 132 : "and to get rid of this warning.")
2430 : END IF
2431 9247 : scf_control%added_mos(1) = MIN(scf_control%added_mos(1), n_ao - n_mo(1))
2432 9247 : n_mo(1) = n_mo(1) + scf_control%added_mos(1)
2433 :
2434 9247 : IF (dft_control%nspins == 2) THEN
2435 1962 : IF (n_mo(2) > n_mo(1)) THEN
2436 : CALL cp_warn(__LOCATION__, &
2437 : "More beta than alpha MOs requested. "// &
2438 0 : "The number of beta MOs will be reduced to the number alpha MOs.")
2439 : END IF
2440 1962 : n_mo(2) = MIN(n_mo(1), n_mo(2))
2441 1962 : CPASSERT(n_mo(1) >= nelectron_spin(1))
2442 1962 : CPASSERT(n_mo(2) >= nelectron_spin(2))
2443 : END IF
2444 :
2445 : ! kpoints
2446 9247 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
2447 9247 : IF (do_kpoints .AND. dft_control%nspins == 2) THEN
2448 : ! we need equal number of calculated states
2449 96 : IF (n_mo(2) /= n_mo(1)) THEN
2450 : CALL cp_warn(__LOCATION__, &
2451 : "Kpoints: Different number of MOs requested. "// &
2452 52 : "The number of beta MOs will be set to the number alpha MOs.")
2453 : END IF
2454 96 : n_mo(2) = n_mo(1)
2455 96 : CPASSERT(n_mo(1) >= nelectron_spin(1))
2456 96 : CPASSERT(n_mo(2) >= nelectron_spin(2))
2457 : END IF
2458 :
2459 : ! Compatibility checks for smearing
2460 9247 : IF (scf_control%smear%do_smear) THEN
2461 1212 : IF (scf_control%added_mos(1) == 0) THEN
2462 0 : CPABORT("Extra MOs (ADDED_MOS) are required for smearing")
2463 : END IF
2464 : END IF
2465 :
2466 : ! Some options require that all MOs are computed ...
2467 : IF ((scf_control%level_shift /= 0.0_dp) .OR. &
2468 9247 : (scf_control%diagonalization%eps_jacobi /= 0.0_dp) .OR. &
2469 : (dft_control%roks .AND. (.NOT. scf_control%use_ot))) THEN
2470 9383 : n_mo(:) = n_ao
2471 : END IF
2472 :
2473 : ! Compatibility checks for ROKS
2474 9247 : IF (dft_control%roks .AND. (.NOT. scf_control%use_ot)) THEN
2475 44 : IF (scf_control%roks_scheme == general_roks) THEN
2476 0 : CPWARN("General ROKS scheme is not yet tested!")
2477 : END IF
2478 44 : IF (scf_control%smear%do_smear) THEN
2479 : CALL cp_abort(__LOCATION__, &
2480 : "The options ROKS and SMEAR are not compatible. "// &
2481 0 : "Try UKS instead of ROKS")
2482 : END IF
2483 : END IF
2484 9247 : IF (dft_control%low_spin_roks) THEN
2485 40 : SELECT CASE (dft_control%qs_control%method_id)
2486 : CASE DEFAULT
2487 : CASE (do_method_xtb, do_method_dftb)
2488 : CALL cp_abort(__LOCATION__, &
2489 0 : "xTB/DFTB methods are not compatible with low spin ROKS.")
2490 : CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pm3, &
2491 : do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
2492 : CALL cp_abort(__LOCATION__, &
2493 40 : "SE methods are not compatible with low spin ROKS.")
2494 : END SELECT
2495 : END IF
2496 :
2497 : ! in principle the restricted calculation could be performed
2498 : ! using just one set of MOs and special casing most of the code
2499 : ! right now we'll just take care of what is effectively an additional constraint
2500 : ! at as few places as possible, just duplicating the beta orbitals
2501 9247 : IF (dft_control%restricted .AND. (output_unit > 0)) THEN
2502 : ! it is really not yet tested till the end ! Joost
2503 33 : WRITE (output_unit, *) ""
2504 33 : WRITE (output_unit, *) " **************************************"
2505 33 : WRITE (output_unit, *) " restricted calculation cutting corners"
2506 33 : WRITE (output_unit, *) " experimental feature, check code "
2507 33 : WRITE (output_unit, *) " **************************************"
2508 : END IF
2509 :
2510 : ! no point in allocating these things here ?
2511 9247 : IF (dft_control%qs_control%do_ls_scf) THEN
2512 352 : NULLIFY (mos)
2513 : ELSE
2514 37520 : ALLOCATE (mos(dft_control%nspins))
2515 19730 : DO ispin = 1, dft_control%nspins
2516 : CALL allocate_mo_set(mo_set=mos(ispin), &
2517 : nao=n_ao, &
2518 : nmo=n_mo(ispin), &
2519 : nelectron=nelectron_spin(ispin), &
2520 : n_el_f=REAL(nelectron_spin(ispin), dp), &
2521 : maxocc=maxocc, &
2522 19730 : flexible_electron_count=dft_control%relax_multiplicity)
2523 : END DO
2524 : END IF
2525 :
2526 9247 : CALL set_qs_env(qs_env, mos=mos)
2527 :
2528 : ! allocate mos when switch_surf_dip is triggered [SGh]
2529 9247 : IF (dft_control%switch_surf_dip) THEN
2530 8 : ALLOCATE (mos_last_converged(dft_control%nspins))
2531 4 : DO ispin = 1, dft_control%nspins
2532 : CALL allocate_mo_set(mo_set=mos_last_converged(ispin), &
2533 : nao=n_ao, &
2534 : nmo=n_mo(ispin), &
2535 : nelectron=nelectron_spin(ispin), &
2536 : n_el_f=REAL(nelectron_spin(ispin), dp), &
2537 : maxocc=maxocc, &
2538 4 : flexible_electron_count=dft_control%relax_multiplicity)
2539 : END DO
2540 2 : CALL set_qs_env(qs_env, mos_last_converged=mos_last_converged)
2541 : END IF
2542 :
2543 9247 : IF (.NOT. be_silent) THEN
2544 : ! Print the DFT control parameters
2545 9219 : IF (PRESENT(multip)) THEN
2546 44 : dft_control%multiplicity = multiplicity
2547 : END IF
2548 9219 : CALL write_dft_control(dft_control, dft_section)
2549 :
2550 : ! Print the vdW control parameters
2551 : IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
2552 : dft_control%qs_control%method_id == do_method_gapw .OR. &
2553 : dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
2554 : dft_control%qs_control%method_id == do_method_lrigpw .OR. &
2555 : dft_control%qs_control%method_id == do_method_rigpw .OR. &
2556 : dft_control%qs_control%method_id == do_method_dftb .OR. &
2557 : (dft_control%qs_control%method_id == do_method_xtb .AND. &
2558 9219 : (.NOT. dft_control%qs_control%xtb_control%do_tblite)) .OR. &
2559 : dft_control%qs_control%method_id == do_method_ofgpw) THEN
2560 8023 : CALL get_qs_env(qs_env, dispersion_env=dispersion_env)
2561 8023 : CALL qs_write_dispersion(qs_env, dispersion_env)
2562 : END IF
2563 :
2564 : ! Print the Quickstep control parameters
2565 9219 : CALL write_qs_control(dft_control%qs_control, dft_section)
2566 :
2567 : ! Print the ADMM control parameters
2568 9219 : IF (dft_control%do_admm) THEN
2569 528 : CALL write_admm_control(dft_control%admm_control, dft_section)
2570 : END IF
2571 :
2572 : ! Print XES/XAS control parameters
2573 9219 : IF (dft_control%do_xas_calculation) THEN
2574 42 : CALL cite_reference(Iannuzzi2007)
2575 : !CALL write_xas_control(dft_control%xas_control,dft_section)
2576 : END IF
2577 :
2578 : ! Print the unnormalized basis set information (input data)
2579 9219 : CALL write_gto_basis_sets(qs_kind_set, subsys_section)
2580 :
2581 : ! Print the atomic kind set
2582 9219 : CALL write_qs_kind_set(qs_kind_set, subsys_section)
2583 :
2584 : ! Print the molecule kind set
2585 9219 : CALL write_molecule_kind_set(molecule_kind_set, subsys_section)
2586 :
2587 : ! Print the total number of kinds, atoms, basis functions etc.
2588 9219 : CALL write_total_numbers(qs_kind_set, particle_set, qs_env%input)
2589 :
2590 : ! Print the atomic coordinates
2591 9219 : CALL write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label="QUICKSTEP")
2592 :
2593 : ! Print the interatomic distances
2594 9219 : CALL write_particle_distances(particle_set, cell, subsys_section)
2595 :
2596 : ! Print the requested structure data
2597 9219 : CALL write_structure_data(particle_set, cell, subsys_section)
2598 :
2599 : ! Print symmetry information
2600 9219 : CALL write_symmetry(particle_set, cell, subsys_section)
2601 :
2602 : ! Print the SCF parameters
2603 9219 : IF ((.NOT. dft_control%qs_control%do_ls_scf) .AND. &
2604 : (.NOT. dft_control%qs_control%do_almo_scf)) THEN
2605 8795 : CALL scf_c_write_parameters(scf_control, dft_section)
2606 : END IF
2607 : END IF
2608 :
2609 : ! Sets up pw_env, qs_charges, mpools ...
2610 9247 : CALL qs_env_setup(qs_env)
2611 :
2612 : ! Allocate and initialise rho0 soft on the global grid
2613 9247 : IF (dft_control%qs_control%method_id == do_method_gapw) THEN
2614 1354 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho0_mpole=rho0_mpole)
2615 1354 : CALL rho0_s_grid_create(pw_env, rho0_mpole)
2616 : END IF
2617 :
2618 9247 : IF (output_unit > 0) CALL m_flush(output_unit)
2619 9247 : CALL timestop(handle)
2620 :
2621 101717 : END SUBROUTINE qs_init_subsys
2622 :
2623 : ! **************************************************************************************************
2624 : !> \brief Write the total number of kinds, atoms, etc. to the logical unit
2625 : !> number lunit.
2626 : !> \param qs_kind_set ...
2627 : !> \param particle_set ...
2628 : !> \param force_env_section ...
2629 : !> \author Creation (06.10.2000)
2630 : ! **************************************************************************************************
2631 9219 : SUBROUTINE write_total_numbers(qs_kind_set, particle_set, force_env_section)
2632 :
2633 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2634 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2635 : TYPE(section_vals_type), POINTER :: force_env_section
2636 :
2637 : INTEGER :: maxlgto, maxlppl, maxlppnl, natom, &
2638 : natom_q, ncgf, nkind, nkind_q, npgf, &
2639 : nset, nsgf, nshell, output_unit
2640 : TYPE(cp_logger_type), POINTER :: logger
2641 :
2642 9219 : NULLIFY (logger)
2643 9219 : logger => cp_get_default_logger()
2644 : output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%TOTAL_NUMBERS", &
2645 9219 : extension=".Log")
2646 :
2647 9219 : IF (output_unit > 0) THEN
2648 4615 : natom = SIZE(particle_set)
2649 4615 : nkind = SIZE(qs_kind_set)
2650 :
2651 : CALL get_qs_kind_set(qs_kind_set, &
2652 : maxlgto=maxlgto, &
2653 : ncgf=ncgf, &
2654 : npgf=npgf, &
2655 : nset=nset, &
2656 : nsgf=nsgf, &
2657 : nshell=nshell, &
2658 : maxlppl=maxlppl, &
2659 4615 : maxlppnl=maxlppnl)
2660 :
2661 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
2662 4615 : "TOTAL NUMBERS AND MAXIMUM NUMBERS"
2663 :
2664 4615 : IF (nset + npgf + ncgf > 0) THEN
2665 : WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") &
2666 4615 : "Total number of", &
2667 4615 : "- Atomic kinds: ", nkind, &
2668 4615 : "- Atoms: ", natom, &
2669 4615 : "- Shell sets: ", nset, &
2670 4615 : "- Shells: ", nshell, &
2671 4615 : "- Primitive Cartesian functions: ", npgf, &
2672 4615 : "- Cartesian basis functions: ", ncgf, &
2673 9230 : "- Spherical basis functions: ", nsgf
2674 0 : ELSE IF (nshell + nsgf > 0) THEN
2675 : WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") &
2676 0 : "Total number of", &
2677 0 : "- Atomic kinds: ", nkind, &
2678 0 : "- Atoms: ", natom, &
2679 0 : "- Shells: ", nshell, &
2680 0 : "- Spherical basis functions: ", nsgf
2681 : ELSE
2682 : WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") &
2683 0 : "Total number of", &
2684 0 : "- Atomic kinds: ", nkind, &
2685 0 : "- Atoms: ", natom
2686 : END IF
2687 :
2688 4615 : IF ((maxlppl > -1) .AND. (maxlppnl > -1)) THEN
2689 : WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T75,I6))") &
2690 2269 : "Maximum angular momentum of the", &
2691 2269 : "- Orbital basis functions: ", maxlgto, &
2692 2269 : "- Local part of the GTH pseudopotential: ", maxlppl, &
2693 4538 : "- Non-local part of the GTH pseudopotential: ", maxlppnl
2694 2346 : ELSE IF (maxlppl > -1) THEN
2695 : WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T75,I6))") &
2696 661 : "Maximum angular momentum of the", &
2697 661 : "- Orbital basis functions: ", maxlgto, &
2698 1322 : "- Local part of the GTH pseudopotential: ", maxlppl
2699 : ELSE
2700 : WRITE (UNIT=output_unit, FMT="(/,T3,A,T75,I6)") &
2701 1685 : "Maximum angular momentum of the orbital basis functions: ", maxlgto
2702 : END IF
2703 :
2704 : ! LRI_AUX BASIS
2705 : CALL get_qs_kind_set(qs_kind_set, &
2706 : maxlgto=maxlgto, &
2707 : ncgf=ncgf, &
2708 : npgf=npgf, &
2709 : nset=nset, &
2710 : nsgf=nsgf, &
2711 : nshell=nshell, &
2712 4615 : basis_type="LRI_AUX")
2713 4615 : IF (nset + npgf + ncgf > 0) THEN
2714 : WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2715 176 : "LRI_AUX Basis: ", &
2716 176 : "Total number of", &
2717 176 : "- Shell sets: ", nset, &
2718 176 : "- Shells: ", nshell, &
2719 176 : "- Primitive Cartesian functions: ", npgf, &
2720 176 : "- Cartesian basis functions: ", ncgf, &
2721 352 : "- Spherical basis functions: ", nsgf
2722 : WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
2723 176 : " Maximum angular momentum ", maxlgto
2724 : END IF
2725 :
2726 : ! RI_HXC BASIS
2727 : CALL get_qs_kind_set(qs_kind_set, &
2728 : maxlgto=maxlgto, &
2729 : ncgf=ncgf, &
2730 : npgf=npgf, &
2731 : nset=nset, &
2732 : nsgf=nsgf, &
2733 : nshell=nshell, &
2734 4615 : basis_type="RI_HXC")
2735 4615 : IF (nset + npgf + ncgf > 0) THEN
2736 : WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2737 150 : "RI_HXC Basis: ", &
2738 150 : "Total number of", &
2739 150 : "- Shell sets: ", nset, &
2740 150 : "- Shells: ", nshell, &
2741 150 : "- Primitive Cartesian functions: ", npgf, &
2742 150 : "- Cartesian basis functions: ", ncgf, &
2743 300 : "- Spherical basis functions: ", nsgf
2744 : WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
2745 150 : " Maximum angular momentum ", maxlgto
2746 : END IF
2747 :
2748 : ! AUX_FIT BASIS
2749 : CALL get_qs_kind_set(qs_kind_set, &
2750 : maxlgto=maxlgto, &
2751 : ncgf=ncgf, &
2752 : npgf=npgf, &
2753 : nset=nset, &
2754 : nsgf=nsgf, &
2755 : nshell=nshell, &
2756 4615 : basis_type="AUX_FIT")
2757 4615 : IF (nset + npgf + ncgf > 0) THEN
2758 : WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2759 416 : "AUX_FIT ADMM-Basis: ", &
2760 416 : "Total number of", &
2761 416 : "- Shell sets: ", nset, &
2762 416 : "- Shells: ", nshell, &
2763 416 : "- Primitive Cartesian functions: ", npgf, &
2764 416 : "- Cartesian basis functions: ", ncgf, &
2765 832 : "- Spherical basis functions: ", nsgf
2766 : WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
2767 416 : " Maximum angular momentum ", maxlgto
2768 : END IF
2769 :
2770 : ! NUCLEAR BASIS
2771 : CALL get_qs_kind_set(qs_kind_set, &
2772 : nkind_q=nkind_q, &
2773 : natom_q=natom_q, &
2774 : maxlgto=maxlgto, &
2775 : ncgf=ncgf, &
2776 : npgf=npgf, &
2777 : nset=nset, &
2778 : nsgf=nsgf, &
2779 : nshell=nshell, &
2780 4615 : basis_type="NUC")
2781 4615 : IF (nset + npgf + ncgf > 0) THEN
2782 : WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2783 153 : "Nuclear Basis: ", &
2784 153 : "Total number of", &
2785 153 : "- Quantum atomic kinds: ", nkind_q, &
2786 153 : "- Quantum atoms: ", natom_q, &
2787 153 : "- Shell sets: ", nset, &
2788 153 : "- Shells: ", nshell, &
2789 153 : "- Primitive Cartesian functions: ", npgf, &
2790 153 : "- Cartesian basis functions: ", ncgf, &
2791 306 : "- Spherical basis functions: ", nsgf
2792 : WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
2793 153 : " Maximum angular momentum ", maxlgto
2794 : END IF
2795 :
2796 : END IF
2797 : CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
2798 9219 : "PRINT%TOTAL_NUMBERS")
2799 :
2800 9219 : END SUBROUTINE write_total_numbers
2801 :
2802 : END MODULE qs_environment
|
