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