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