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