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