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