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