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