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