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