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