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