(git:c28f603)
Loading...
Searching...
No Matches
qs_environment.F
Go to the documentation of this file.
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! **************************************************************************************************
30 USE bibliography, ONLY: iannuzzi2006,&
32 cite_reference,&
34 USE cell_types, ONLY: cell_type
44 USE cp_control_utils, ONLY: &
53 USE cp_output_handling, ONLY: cp_p_file,&
63 USE ec_environment, ONLY: ec_env_create,&
82 USE gamma, ONLY: init_md_ftable
85 USE header, ONLY: dftb_header,&
86 qs_header,&
87 se_header,&
92 USE input_constants, ONLY: &
110 USE kinds, ONLY: default_string_length,&
111 dp
115 USE kpoint_types, ONLY: get_kpoint_info,&
125 USE machine, ONLY: m_flush
126 USE mathconstants, ONLY: pi
131 USE mp2_setup, ONLY: read_mp2_section
132 USE mp2_types, ONLY: mp2_env_create,&
141 USE physcon, ONLY: kelvin
142 USE pw_env_types, ONLY: pw_env_type
161 USE qs_gcp_types, ONLY: qs_gcp_type
162 USE qs_gcp_utils, ONLY: qs_gcp_env_set,&
175 USE qs_kind_types, ONLY: &
179 USE qs_ks_types, ONLY: qs_ks_env_create,&
183 USE qs_mo_types, ONLY: allocate_mo_set,&
186 USE qs_rho0_methods, ONLY: init_rho0
191 USE qs_subsys_types, ONLY: qs_subsys_get,&
214 USE tblite_interface, ONLY: tb_get_basis,&
216 tb_init_wf,&
219 USE xtb_parameters, ONLY: init_xtb_basis,&
225#include "./base/base_uses.f90"
226
227 IMPLICIT NONE
228
229 PRIVATE
230
231 ! *** Global parameters ***
232 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_environment'
233
234 ! *** Public subroutines ***
235 PUBLIC :: qs_init
236
237CONTAINS
238
239! **************************************************************************************************
240!> \brief Read the input and the database files for the setup of the
241!> QUICKSTEP environment.
242!> \param qs_env ...
243!> \param para_env ...
244!> \param root_section ...
245!> \param globenv ...
246!> \param cp_subsys ...
247!> \param kpoint_env ...
248!> \param qmmm ...
249!> \param qmmm_env_qm ...
250!> \param force_env_section ...
251!> \param subsys_section ...
252!> \param use_motion_section ...
253!> \param silent ...
254!> \param multip ...
255!> \param charge ...
256!> \author Creation (22.05.2000,MK)
257! **************************************************************************************************
258 SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, &
259 qmmm, qmmm_env_qm, force_env_section, subsys_section, &
260 use_motion_section, silent, multip, charge)
261
262 TYPE(qs_environment_type), POINTER :: qs_env
263 TYPE(mp_para_env_type), POINTER :: para_env
264 TYPE(section_vals_type), OPTIONAL, POINTER :: root_section
265 TYPE(global_environment_type), OPTIONAL, POINTER :: globenv
266 TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys
267 TYPE(kpoint_type), OPTIONAL, POINTER :: kpoint_env
268 LOGICAL, INTENT(IN), OPTIONAL :: qmmm
269 TYPE(qmmm_env_qm_type), OPTIONAL, POINTER :: qmmm_env_qm
270 TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
271 LOGICAL, INTENT(IN) :: use_motion_section
272 LOGICAL, INTENT(IN), OPTIONAL :: silent
273 INTEGER, INTENT(IN), OPTIONAL :: multip, charge
274
275 CHARACTER(LEN=default_string_length) :: basis_type
276 INTEGER :: ikind, method_id, nelectron_total, &
277 nkind, nkp_grid(3), tddfpt_kernel
278 LOGICAL :: dftb_kpoint_sym_restricted, do_active_space, do_admm, do_admm_rpa, do_bse, &
279 do_debug_fdiff, do_debug_forces, do_debug_stress_tensor, do_dftb_scc, do_dftb_scc_high_l, &
280 do_ec_hfx, do_et, do_exx, do_gw, do_hfx, do_kpoints, do_linear_response, do_mp2, &
281 do_ri_mp2, do_ri_rpa, do_ri_sos_mp2, do_tddfpt, do_tddfpt_unsupported_kpoints, &
282 do_wfc_low_scaling, do_wfc_low_scaling_kpoints, do_xtb_tblite, final_kpoint_reinit, &
283 is_identical, is_semi, kpoint_verbose, mp2_present, my_qmmm, owned_kpoints, qmmm_decoupl, &
284 same_except_frac, use_real_wfn, use_ref_cell
285 REAL(kind=dp), DIMENSION(:, :), POINTER :: rtmat
286 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
287 TYPE(cell_type), POINTER :: my_cell, my_cell_ref
288 TYPE(cp_blacs_env_type), POINTER :: blacs_env
289 TYPE(dft_control_type), POINTER :: dft_control
290 TYPE(distribution_1d_type), POINTER :: local_particles
291 TYPE(energy_correction_type), POINTER :: ec_env
292 TYPE(excited_energy_type), POINTER :: exstate_env
293 TYPE(harris_type), POINTER :: harris_env
294 TYPE(kpoint_type), POINTER :: kpoints
295 TYPE(lri_environment_type), POINTER :: lri_env
296 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
297 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
298 TYPE(qs_ks_env_type), POINTER :: ks_env
299 TYPE(qs_subsys_type), POINTER :: subsys
300 TYPE(qs_wf_history_type), POINTER :: wf_history
301 TYPE(rel_control_type), POINTER :: rel_control
302 TYPE(scf_control_type), POINTER :: scf_control
303 TYPE(section_vals_type), POINTER :: active_space_section, admm_section, dft_section, &
304 ec_hfx_section, ec_section, et_coupling_section, gw_section, hfx_section, kpoint_section, &
305 mp2_section, rpa_hfx_section, tddfpt_section, transport_section
306
307 NULLIFY (my_cell, my_cell_ref, atomic_kind_set, particle_set, &
308 qs_kind_set, kpoint_section, dft_section, ec_section, &
309 subsys, ks_env, dft_control, blacs_env)
310
311 CALL set_qs_env(qs_env, input=force_env_section)
312 IF (.NOT. ASSOCIATED(subsys_section)) THEN
313 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
314 END IF
315
316 ! QMMM
317 my_qmmm = .false.
318 IF (PRESENT(qmmm)) my_qmmm = qmmm
319 qmmm_decoupl = .false.
320 IF (PRESENT(qmmm_env_qm)) THEN
321 IF (qmmm_env_qm%qmmm_coupl_type == do_qmmm_gauss .OR. &
322 qmmm_env_qm%qmmm_coupl_type == do_qmmm_swave) THEN
323 ! For GAUSS/SWAVE methods there could be a DDAPC decoupling requested
324 qmmm_decoupl = my_qmmm .AND. qmmm_env_qm%periodic .AND. qmmm_env_qm%multipole
325 END IF
326 qs_env%qmmm_env_qm => qmmm_env_qm
327 END IF
328 CALL set_qs_env(qs_env=qs_env, qmmm=my_qmmm)
329
330 ! Possibly initialize arrays for SE
331 CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=method_id)
332 SELECT CASE (method_id)
335 CALL init_se_intd_array()
336 is_semi = .true.
338 is_semi = .true.
339 CASE DEFAULT
340 is_semi = .false.
341 END SELECT
342
343 ALLOCATE (subsys)
344 CALL qs_subsys_create(subsys, para_env, &
345 force_env_section=force_env_section, &
346 subsys_section=subsys_section, &
347 use_motion_section=use_motion_section, &
348 root_section=root_section, &
349 cp_subsys=cp_subsys, &
350 elkind=is_semi, silent=silent)
351
352 ALLOCATE (ks_env)
353 CALL qs_ks_env_create(ks_env)
354 CALL set_ks_env(ks_env, subsys=subsys)
355 CALL set_qs_env(qs_env, ks_env=ks_env)
356
357 CALL qs_subsys_get(subsys, &
358 cell=my_cell, &
359 cell_ref=my_cell_ref, &
360 use_ref_cell=use_ref_cell, &
361 atomic_kind_set=atomic_kind_set, &
362 qs_kind_set=qs_kind_set, &
363 particle_set=particle_set)
364
365 CALL set_ks_env(ks_env, para_env=para_env)
366 IF (PRESENT(globenv)) THEN
367 CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
368 globenv%blacs_repeatable)
369 ELSE
370 CALL cp_blacs_env_create(blacs_env, para_env)
371 END IF
372 CALL set_ks_env(ks_env, blacs_env=blacs_env)
373 CALL cp_blacs_env_release(blacs_env)
374
375 ! *** Setup the grids for the G-space Interpolation if any
376 CALL cp_ddapc_ewald_create(qs_env%cp_ddapc_ewald, qmmm_decoupl, my_cell, &
377 force_env_section, subsys_section, para_env)
378
379 ! kpoints
380 IF (PRESENT(kpoint_env)) THEN
381 owned_kpoints = .false.
382 kpoints => kpoint_env
383 CALL set_qs_env(qs_env=qs_env, kpoints=kpoints)
384 CALL kpoint_initialize(kpoints, particle_set, my_cell)
385 ELSE
386 owned_kpoints = .true.
387 NULLIFY (kpoints)
388 CALL kpoint_create(kpoints)
389 CALL set_qs_env(qs_env=qs_env, kpoints=kpoints)
390 kpoint_section => section_vals_get_subs_vals(qs_env%input, "DFT%KPOINTS")
391 CALL read_kpoint_section(kpoints, kpoint_section, my_cell%hmat, my_cell)
392 CALL get_kpoint_info(kpoints, verbose=kpoint_verbose)
393 IF (kpoint_verbose) CALL set_kpoint_info(kpoints, verbose=.false.)
394 do_hfx = .false.
395 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
396 CALL section_vals_get(hfx_section, explicit=do_hfx)
397 do_exx = .false.
398 rpa_hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
399 CALL section_vals_get(rpa_hfx_section, explicit=do_exx)
400 do_admm = .false.
401 admm_section => section_vals_get_subs_vals(qs_env%input, "DFT%AUXILIARY_DENSITY_MATRIX_METHOD")
402 CALL section_vals_get(admm_section, explicit=do_admm)
403 do_gw = .false.
404 gw_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%GW")
405 CALL section_vals_get(gw_section, explicit=do_gw)
406 IF (.NOT. do_gw) THEN
407 gw_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%BANDSTRUCTURE%GW")
408 CALL section_vals_get(gw_section, explicit=do_gw)
409 END IF
410 do_tddfpt = .false.
411 do_tddfpt_unsupported_kpoints = .false.
412 do_bse = .false.
413 tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
414 CALL section_vals_get(tddfpt_section, explicit=do_tddfpt)
415 IF (do_tddfpt) THEN
416 CALL section_vals_val_get(tddfpt_section, "KERNEL", i_val=tddfpt_kernel)
417 do_tddfpt_unsupported_kpoints = tddfpt_kernel /= tddfpt_kernel_none
418 IF (.NOT. do_tddfpt_unsupported_kpoints) THEN
419 CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
420 IF (use_real_wfn) &
421 CALL cp_abort(__location__, "K-point TDDFPT requires complex wavefunctions.")
422 END IF
423 CALL section_vals_val_get(tddfpt_section, "DO_BSE", l_val=do_bse)
424 IF (.NOT. do_bse) &
425 CALL section_vals_val_get(tddfpt_section, "DO_BSE_W_ONLY", l_val=do_bse)
426 IF (.NOT. do_bse) &
427 CALL section_vals_val_get(tddfpt_section, "DO_BSE_GW_ONLY", l_val=do_bse)
428 END IF
429 do_active_space = .false.
430 active_space_section => section_vals_get_subs_vals(qs_env%input, "DFT%ACTIVE_SPACE")
431 CALL section_vals_get(active_space_section, explicit=do_active_space)
432 do_xtb_tblite = .false.
433 IF (method_id == do_method_xtb) THEN
434 CALL section_vals_val_get(qs_env%input, "DFT%QS%XTB%TBLITE%_SECTION_PARAMETERS_", &
435 l_val=do_xtb_tblite)
436 END IF
437 do_dftb_scc = .false.
438 IF (method_id == do_method_dftb) THEN
439 CALL section_vals_val_get(qs_env%input, "DFT%QS%DFTB%SELF_CONSISTENT", &
440 l_val=do_dftb_scc)
441 END IF
442 do_linear_response = .false.
443 IF (PRESENT(globenv)) do_linear_response = globenv%run_type_id == linear_response_run
444 do_debug_fdiff = .false.
445 IF (PRESENT(globenv)) do_debug_fdiff = globenv%run_type_id == debug_run
446 IF (do_debug_fdiff .AND. PRESENT(root_section)) THEN
447 CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
448 l_val=do_debug_forces)
449 CALL section_vals_val_get(root_section, "DEBUG%DEBUG_STRESS_TENSOR", &
450 l_val=do_debug_stress_tensor)
451 do_debug_fdiff = do_debug_forces .OR. do_debug_stress_tensor
452 END IF
453 do_mp2 = .false.
454 do_ri_mp2 = .false.
455 do_ri_sos_mp2 = .false.
456 do_ri_rpa = .false.
457 do_wfc_low_scaling = .false.
458 do_wfc_low_scaling_kpoints = .false.
459 mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
460 CALL section_vals_get(mp2_section, explicit=mp2_present)
461 IF (mp2_present) THEN
462 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%MP2%_SECTION_PARAMETERS_", &
463 l_val=do_mp2)
464 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_MP2%_SECTION_PARAMETERS_", &
465 l_val=do_ri_mp2)
466 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_SOS_MP2%_SECTION_PARAMETERS_", &
467 l_val=do_ri_sos_mp2)
468 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%_SECTION_PARAMETERS_", &
469 l_val=do_ri_rpa)
470 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%_SECTION_PARAMETERS_", &
471 l_val=do_wfc_low_scaling)
472 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%DO_KPOINTS", &
473 l_val=do_wfc_low_scaling_kpoints)
474 IF (.NOT. do_bse) &
475 CALL section_vals_val_get(qs_env%input, &
476 "DFT%XC%WF_CORRELATION%RI_RPA%GW%BSE%_SECTION_PARAMETERS_", &
477 l_val=do_bse)
478 END IF
479 CALL restrict_unsupported_atomic_kpoint_symmetry(kpoints, method_id, do_hfx, do_exx, do_gw, &
480 do_tddfpt_unsupported_kpoints, &
481 do_active_space, do_linear_response, &
482 do_debug_fdiff, &
483 do_mp2 .OR. do_ri_mp2 .OR. do_ri_sos_mp2, &
484 do_ri_rpa .AND. .NOT. do_gw, do_bse, &
485 do_wfc_low_scaling, do_wfc_low_scaling_kpoints, &
486 do_xtb_tblite, do_admm, .false.)
487 CALL kpoint_initialize(kpoints, particle_set, my_cell)
488 END IF
489
490 CALL qs_init_subsys(qs_env, para_env, subsys, my_cell, my_cell_ref, use_ref_cell, &
491 subsys_section, silent=silent, multip=multip, charge=charge)
492
493 CALL get_qs_env(qs_env, dft_control=dft_control)
494 IF (owned_kpoints) THEN
495 do_dftb_scc_high_l = .false.
496 IF (method_id == do_method_dftb .AND. do_dftb_scc) THEN
497 do_dftb_scc_high_l = dftb_kind_set_has_high_l(qs_kind_set)
498 END IF
499 CALL restrict_unsupported_atomic_kpoint_symmetry(kpoints, method_id, do_hfx, do_exx, do_gw, &
500 do_tddfpt_unsupported_kpoints, &
501 do_active_space, do_linear_response, &
502 do_debug_fdiff, &
503 do_mp2 .OR. do_ri_mp2 .OR. do_ri_sos_mp2, &
504 do_ri_rpa .AND. .NOT. do_gw, do_bse, &
505 do_wfc_low_scaling, do_wfc_low_scaling_kpoints, &
506 do_xtb_tblite, do_admm, do_dftb_scc_high_l, &
507 restricted=dftb_kpoint_sym_restricted)
508 final_kpoint_reinit = dftb_kpoint_sym_restricted .OR. kpoint_verbose
509 IF (final_kpoint_reinit) THEN
510 CALL kpoint_reset_initialization(kpoints)
511 CALL set_kpoint_info(kpoints, verbose=kpoint_verbose)
512 CALL kpoint_initialize(kpoints, particle_set, my_cell)
513 END IF
514 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
515 CALL write_kpoint_info(kpoints, dft_section=dft_section)
516 END IF
517 IF (method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
518 CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
519 CALL lri_env_basis("LRI", qs_env, lri_env, qs_kind_set)
520 ELSE IF (method_id == do_method_rigpw) THEN
521 CALL cp_warn(__location__, "Experimental code: "// &
522 "RIGPW should only be used for testing.")
523 CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
524 CALL lri_env_basis("RI", qs_env, lri_env, qs_kind_set)
525 END IF
526
527 IF (my_qmmm .AND. PRESENT(qmmm_env_qm) .AND. .NOT. dft_control%qs_control%commensurate_mgrids) THEN
528 IF (qmmm_env_qm%qmmm_coupl_type == do_qmmm_gauss .OR. qmmm_env_qm%qmmm_coupl_type == do_qmmm_swave) THEN
529 CALL cp_abort(__location__, "QM/MM with coupling GAUSS or S-WAVE requires "// &
530 "keyword FORCE_EVAL/DFT/MGRID/COMMENSURATE to be enabled.")
531 END IF
532 END IF
533
534 ! more kpoint stuff
535 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints, blacs_env=blacs_env)
536 IF (do_kpoints) THEN
537 IF (dft_control%qs_control%do_ls_scf) &
538 cpabort("DFT%KPOINTS are not implemented with QS/LS_SCF; use a real-space supercell instead.")
539 CALL kpoint_env_initialize(kpoints, para_env, blacs_env, with_aux_fit=dft_control%do_admm)
540 CALL kpoint_initialize_mos(kpoints, qs_env%mos)
541 CALL get_qs_env(qs_env=qs_env, wf_history=wf_history)
542 CALL wfi_create_for_kp(wf_history)
543 END IF
544 ! basis set symmetry rotations
545 IF (do_kpoints) THEN
546 CALL qs_basis_rotation(qs_env, kpoints)
547 END IF
548
549 do_hfx = .false.
550 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
551 CALL section_vals_get(hfx_section, explicit=do_hfx)
552 CALL get_qs_env(qs_env, dft_control=dft_control, scf_control=scf_control, nelectron_total=nelectron_total)
553 IF (do_hfx) THEN
554 ! Retrieve particle_set and atomic_kind_set (needed for both kinds of initialization)
555 nkp_grid = 1
556 IF (do_kpoints) CALL get_kpoint_info(kpoints, nkp_grid=nkp_grid)
557 IF (dft_control%do_admm) THEN
558 basis_type = 'AUX_FIT'
559 ELSE
560 basis_type = 'ORB'
561 END IF
562 CALL hfx_create(qs_env%x_data, para_env, hfx_section, atomic_kind_set, &
563 qs_kind_set, particle_set, dft_control, my_cell, orb_basis=basis_type, &
564 nelectron_total=nelectron_total, nkp_grid=nkp_grid)
565 END IF
566
567 mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
568 CALL section_vals_get(mp2_section, explicit=mp2_present)
569 IF (mp2_present) THEN
570 cpassert(ASSOCIATED(qs_env%mp2_env))
571 CALL read_mp2_section(qs_env%input, qs_env%mp2_env)
572 ! create the EXX section if necessary
573 do_exx = .false.
574 rpa_hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
575 CALL section_vals_get(rpa_hfx_section, explicit=do_exx)
576 IF (do_exx) THEN
577
578 ! do_exx in call of hfx_create decides whether to go without ADMM (do_exx=.TRUE.) or with
579 ! ADMM (do_exx=.FALSE.)
580 CALL section_vals_val_get(mp2_section, "RI_RPA%ADMM", l_val=do_admm_rpa)
581
582 ! Reuse the HFX integrals from the qs_env if applicable
583 qs_env%mp2_env%ri_rpa%reuse_hfx = .true.
584 IF (.NOT. do_hfx) qs_env%mp2_env%ri_rpa%reuse_hfx = .false.
585 CALL compare_hfx_sections(hfx_section, rpa_hfx_section, is_identical, same_except_frac)
586 IF (.NOT. (is_identical .OR. same_except_frac)) qs_env%mp2_env%ri_rpa%reuse_hfx = .false.
587 IF (dft_control%do_admm .AND. .NOT. do_admm_rpa) qs_env%mp2_env%ri_rpa%reuse_hfx = .false.
588
589 IF (.NOT. qs_env%mp2_env%ri_rpa%reuse_hfx) THEN
590 IF (do_admm_rpa) THEN
591 basis_type = 'AUX_FIT'
592 ELSE
593 basis_type = 'ORB'
594 END IF
595 CALL hfx_create(qs_env%mp2_env%ri_rpa%x_data, para_env, rpa_hfx_section, atomic_kind_set, &
596 qs_kind_set, particle_set, dft_control, my_cell, orb_basis=basis_type, &
597 nelectron_total=nelectron_total)
598 ELSE
599 qs_env%mp2_env%ri_rpa%x_data => qs_env%x_data
600 END IF
601 END IF
602 END IF
603
604 IF (dft_control%qs_control%do_kg) THEN
605 CALL cite_reference(iannuzzi2006)
606 CALL kg_env_create(qs_env, qs_env%kg_env, qs_kind_set, qs_env%input)
607 END IF
608
609 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
610 CALL section_vals_val_get(dft_section, "EXCITED_STATES%_SECTION_PARAMETERS_", &
611 l_val=qs_env%excited_state)
612 NULLIFY (exstate_env)
613 CALL exstate_create(exstate_env, qs_env%excited_state, dft_section)
614 CALL set_qs_env(qs_env, exstate_env=exstate_env)
615
616 et_coupling_section => section_vals_get_subs_vals(qs_env%input, &
617 "PROPERTIES%ET_COUPLING")
618 CALL section_vals_get(et_coupling_section, explicit=do_et)
619 IF (do_et) CALL et_coupling_create(qs_env%et_coupling)
620
621 transport_section => section_vals_get_subs_vals(qs_env%input, "DFT%TRANSPORT")
622 CALL section_vals_get(transport_section, explicit=qs_env%do_transport)
623 IF (qs_env%do_transport) THEN
624 CALL transport_env_create(qs_env)
625 END IF
626
627 CALL get_qs_env(qs_env, harris_env=harris_env)
628 IF (qs_env%harris_method) THEN
629 ! initialize the Harris input density and potential integrals
630 CALL get_qs_env(qs_env, local_particles=local_particles)
631 CALL harris_rhoin_init(harris_env%rhoin, "RHOIN", qs_kind_set, atomic_kind_set, &
632 local_particles, dft_control%nspins)
633 ! Print information of the HARRIS section
634 CALL harris_write_input(harris_env)
635 END IF
636
637 NULLIFY (ec_env)
638 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
639 CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%_SECTION_PARAMETERS_", &
640 l_val=qs_env%energy_correction)
641 ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
642 CALL ec_env_create(qs_env, ec_env, dft_section, ec_section)
643 CALL set_qs_env(qs_env, ec_env=ec_env)
644
645 IF (qs_env%energy_correction) THEN
646 ! Energy correction with Hartree-Fock exchange
647 ec_hfx_section => section_vals_get_subs_vals(ec_section, "XC%HF")
648 CALL section_vals_get(ec_hfx_section, explicit=do_ec_hfx)
649
650 IF (ec_env%do_ec_hfx) THEN
651
652 ! kpoints and HFX not yet compatible
653 IF (ec_env%do_kpoints) THEN
654 CALL cp_abort(__location__, &
655 "Energy correction methods with hybrid functionals "// &
656 "and kpoints is not yet available.")
657 END IF
658
659 ! Hybrid functionals require same basis
660 IF (ec_env%basis_inconsistent) THEN
661 CALL cp_abort(__location__, &
662 "Energy correction methods with hybrid functionals: "// &
663 "correction and ground state need to use the same basis. "// &
664 "Checked by comparing basis set names only.")
665 END IF
666
667 ! Similar to RPA_HFX we can check if HFX integrals from the qs_env can be reused
668 IF (ec_env%do_ec_admm .AND. .NOT. dft_control%do_admm) THEN
669 CALL cp_abort(__location__, "Need an ADMM input section for ADMM EC to work")
670 END IF
671
672 ec_env%reuse_hfx = .true.
673 IF (.NOT. do_hfx) ec_env%reuse_hfx = .false.
674 CALL compare_hfx_sections(hfx_section, ec_hfx_section, is_identical, same_except_frac)
675 IF (.NOT. (is_identical .OR. same_except_frac)) ec_env%reuse_hfx = .false.
676 IF (dft_control%do_admm .AND. .NOT. ec_env%do_ec_admm) ec_env%reuse_hfx = .false.
677
678 IF (.NOT. ec_env%reuse_hfx) THEN
679 IF (ec_env%do_ec_admm) THEN
680 basis_type = 'AUX_FIT'
681 ELSE
682 basis_type = 'ORB'
683 END IF
684 CALL hfx_create(ec_env%x_data, para_env, ec_hfx_section, atomic_kind_set, &
685 qs_kind_set, particle_set, dft_control, my_cell, orb_basis=basis_type, &
686 nelectron_total=nelectron_total)
687 ELSE
688 ec_env%x_data => qs_env%x_data
689 END IF
690 END IF
691
692 ! Print information of the EC section
693 CALL ec_write_input(ec_env)
694
695 END IF
696
697 IF (dft_control%qs_control%do_almo_scf) THEN
698 CALL almo_scf_env_create(qs_env)
699 END IF
700
701 ! see if we have atomic relativistic corrections
702 CALL get_qs_env(qs_env, rel_control=rel_control)
703 IF (rel_control%rel_method /= rel_none) THEN
704 IF (rel_control%rel_transformation == rel_trans_atom) THEN
705 nkind = SIZE(atomic_kind_set)
706 DO ikind = 1, nkind
707 NULLIFY (rtmat)
708 CALL calculate_atomic_relkin(atomic_kind_set(ikind), qs_kind_set(ikind), rel_control, rtmat)
709 IF (ASSOCIATED(rtmat)) CALL set_qs_kind(qs_kind_set(ikind), reltmat=rtmat)
710 END DO
711 END IF
712 END IF
713
714 END SUBROUTINE qs_init
715
716! **************************************************************************************************
717!> \brief Restrict atomic k-point symmetry for methods not supporting it yet
718!> \param kpoints ...
719!> \param method_id ...
720!> \param do_hfx ...
721!> \param do_exx ...
722!> \param do_gw ...
723!> \param do_tddfpt ...
724!> \param do_active_space ...
725!> \param do_linear_response ...
726!> \param do_debug_fdiff ...
727!> \param do_mp2 ...
728!> \param do_rpa ...
729!> \param do_bse ...
730!> \param do_wfc_low_scaling ...
731!> \param do_wfc_low_scaling_kpoints ...
732!> \param do_xtb_tblite ...
733!> \param do_admm ...
734!> \param do_dftb_scc_high_l ...
735!> \param restricted ...
736! **************************************************************************************************
737 SUBROUTINE restrict_unsupported_atomic_kpoint_symmetry(kpoints, method_id, do_hfx, do_exx, do_gw, &
738 do_tddfpt, do_active_space, do_linear_response, &
739 do_debug_fdiff, &
740 do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
741 do_wfc_low_scaling_kpoints, do_xtb_tblite, &
742 do_admm, do_dftb_scc_high_l, restricted)
743 TYPE(kpoint_type), POINTER :: kpoints
744 INTEGER, INTENT(IN) :: method_id
745 LOGICAL, INTENT(IN) :: do_hfx, do_exx, do_gw, do_tddfpt, do_active_space, &
746 do_linear_response, do_debug_fdiff, do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
747 do_wfc_low_scaling_kpoints, do_xtb_tblite, do_admm, do_dftb_scc_high_l
748 LOGICAL, INTENT(OUT), OPTIONAL :: restricted
749
750 CHARACTER(LEN=default_string_length) :: kp_scheme, reason
751 LOGICAL :: full_grid, inversion_symmetry_only, &
752 kpoint_symmetry
753
754 IF (PRESENT(restricted)) restricted = .false.
755
756 reason = unsupported_kpoint_method_reason(method_id, do_gw, do_tddfpt, do_linear_response, &
757 do_mp2, do_bse, do_xtb_tblite)
758 IF (len_trim(reason) > 0) THEN
759 CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme)
760 IF (len_trim(kp_scheme) > 0 .AND. trim(kp_scheme) /= "NONE") THEN
761 IF (trim(reason) == "GW") THEN
762 CALL cp_abort(__location__, &
763 "DFT%KPOINTS are not supported with GW; use "// &
764 "WF_CORRELATION%LOW_SCALING%KPOINTS and RI_RPA%GW%KPOINTS_SELF_ENERGY "// &
765 "for GW k-point sampling.")
766 ELSE
767 CALL cp_abort(__location__, &
768 "DFT%KPOINTS are not supported with "//trim(reason)// &
769 "; remove DFT%KPOINTS for these calculations.")
770 END IF
771 END IF
772 END IF
773 IF (do_active_space) THEN
774 CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme)
775 IF (len_trim(kp_scheme) > 0 .AND. trim(kp_scheme) /= "NONE" .AND. &
776 trim(kp_scheme) /= "GAMMA") THEN
777 CALL cp_abort(__location__, &
778 "Only Gamma-point DFT%KPOINTS are supported with ACTIVE_SPACE; "// &
779 "use SCHEME GAMMA, SCHEME NONE, or remove DFT%KPOINTS.")
780 END IF
781 END IF
782
783 CALL get_kpoint_info(kpoints, symmetry=kpoint_symmetry, full_grid=full_grid, &
784 inversion_symmetry_only=inversion_symmetry_only)
785 IF (.NOT. (kpoint_symmetry .AND. .NOT. full_grid .AND. .NOT. inversion_symmetry_only)) RETURN
786
787 reason = unsupported_atomic_kpoint_symmetry_reason(method_id, do_hfx, do_exx, do_gw, &
788 do_tddfpt, do_active_space, do_linear_response, &
789 do_debug_fdiff, &
790 do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
791 do_wfc_low_scaling_kpoints, do_xtb_tblite, &
792 do_admm, do_dftb_scc_high_l)
793 IF (len_trim(reason) == 0) RETURN
794
795 CALL cp_warn(__location__, &
796 "Atomic k-point symmetry is currently not implemented for "//trim(reason)// &
797 "; restricting to inversion/time-reversal symmetry.")
798 CALL set_kpoint_info(kpoints, inversion_symmetry_only=.true.)
799 IF (PRESENT(restricted)) restricted = .true.
800
801 END SUBROUTINE restrict_unsupported_atomic_kpoint_symmetry
802
803! **************************************************************************************************
804!> \brief Return the reason why k-points are not enabled for a method
805!> \param method_id ...
806!> \param do_gw ...
807!> \param do_tddfpt ...
808!> \param do_linear_response ...
809!> \param do_mp2 ...
810!> \param do_bse ...
811!> \param do_xtb_tblite ...
812!> \return reason
813! **************************************************************************************************
814 FUNCTION unsupported_kpoint_method_reason(method_id, do_gw, do_tddfpt, do_linear_response, &
815 do_mp2, do_bse, do_xtb_tblite) RESULT(reason)
816 INTEGER, INTENT(IN) :: method_id
817 LOGICAL, INTENT(IN) :: do_gw, do_tddfpt, do_linear_response, &
818 do_mp2, do_bse, do_xtb_tblite
819 CHARACTER(LEN=default_string_length) :: reason
820
821 reason = ""
822 mark_used(do_gw)
823 mark_used(do_mp2)
824 mark_used(do_xtb_tblite)
825
826 IF (do_bse) THEN
827 reason = "BSE"
828 RETURN
829 END IF
830 IF (do_tddfpt) THEN
831 reason = "TDDFPT/TDDFT"
832 RETURN
833 END IF
834 IF (do_linear_response) THEN
835 reason = "LINEAR_RESPONSE/DFPT"
836 RETURN
837 END IF
838 SELECT CASE (method_id)
839 CASE (do_method_rigpw)
840 reason = "RIGPW"
841 CASE (do_method_ofgpw)
842 reason = "OFGPW"
845 reason = "semiempirical methods"
846 CASE DEFAULT
847 reason = ""
848 END SELECT
849
850 END FUNCTION unsupported_kpoint_method_reason
851
852! **************************************************************************************************
853!> \brief Return the reason why atomic k-point symmetry is not enabled
854!> \param method_id ...
855!> \param do_hfx ...
856!> \param do_exx ...
857!> \param do_gw ...
858!> \param do_tddfpt ...
859!> \param do_active_space ...
860!> \param do_linear_response ...
861!> \param do_debug_fdiff ...
862!> \param do_mp2 ...
863!> \param do_rpa ...
864!> \param do_bse ...
865!> \param do_wfc_low_scaling ...
866!> \param do_wfc_low_scaling_kpoints ...
867!> \param do_xtb_tblite ...
868!> \param do_admm ...
869!> \param do_dftb_scc_high_l ...
870!> \return reason
871! **************************************************************************************************
872 FUNCTION unsupported_atomic_kpoint_symmetry_reason(method_id, do_hfx, do_exx, do_gw, do_tddfpt, &
873 do_active_space, do_linear_response, do_debug_fdiff, &
874 do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
875 do_wfc_low_scaling_kpoints, do_xtb_tblite, &
876 do_admm, do_dftb_scc_high_l) RESULT(reason)
877 INTEGER, INTENT(IN) :: method_id
878 LOGICAL, INTENT(IN) :: do_hfx, do_exx, do_gw, do_tddfpt, do_active_space, &
879 do_linear_response, do_debug_fdiff, do_mp2, do_rpa, do_bse, do_wfc_low_scaling, &
880 do_wfc_low_scaling_kpoints, do_xtb_tblite, do_admm, do_dftb_scc_high_l
881 CHARACTER(LEN=default_string_length) :: reason
882
883 reason = ""
884 mark_used(do_debug_fdiff)
885 mark_used(do_xtb_tblite)
886
887 SELECT CASE (method_id)
888 CASE (do_method_dftb)
889 IF (do_dftb_scc_high_l) reason = "SCC-DFTB with d orbitals"
890 CASE (do_method_lrigpw)
891 reason = "LRIGPW"
892 CASE (do_method_rigpw)
893 reason = "RIGPW"
896 reason = "semiempirical methods"
897 CASE DEFAULT
898 reason = ""
899 END SELECT
900
901 IF (len_trim(reason) > 0) RETURN
902 IF ((do_hfx .OR. do_exx) .AND. do_admm) THEN
903 reason = "HFX/HF with ADMM"
904 ELSE IF (do_bse) THEN
905 reason = "BSE"
906 ELSE IF (do_gw) THEN
907 reason = "GW"
908 ELSE IF (do_tddfpt) THEN
909 reason = "TDDFPT/TDDFT"
910 ELSE IF (do_active_space) THEN
911 reason = "ACTIVE_SPACE"
912 ELSE IF (do_linear_response) THEN
913 reason = "LINEAR_RESPONSE/DFPT"
914 ELSE IF (do_mp2) THEN
915 reason = "MP2"
916 ELSE IF (do_rpa .AND. do_wfc_low_scaling_kpoints) THEN
917 reason = "LOW_SCALING RPA"
918 ELSE IF (do_wfc_low_scaling) THEN
919 reason = "LOW_SCALING WF_CORRELATION"
920 ELSE IF (do_rpa) THEN
921 reason = "RPA"
922 END IF
923
924 END FUNCTION unsupported_atomic_kpoint_symmetry_reason
925
926! **************************************************************************************************
927!> \brief Return whether the DFTB kind set contains d orbitals
928!> \param qs_kind_set ...
929!> \return has_high_l
930! **************************************************************************************************
931 FUNCTION dftb_kind_set_has_high_l(qs_kind_set) RESULT(has_high_l)
932 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
933 LOGICAL :: has_high_l
934
935 INTEGER :: ikind, lmax
936 LOGICAL :: any_defined, defined
937 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
938
939 has_high_l = .true.
940 IF (.NOT. ASSOCIATED(qs_kind_set)) RETURN
941
942 any_defined = .false.
943 DO ikind = 1, SIZE(qs_kind_set)
944 NULLIFY (dftb_parameter)
945 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_parameter)
946 IF (.NOT. ASSOCIATED(dftb_parameter)) cycle
947 defined = .false.
948 lmax = -1
949 CALL get_dftb_atom_param(dftb_parameter, defined=defined, lmax=lmax)
950 IF (.NOT. defined) cycle
951 any_defined = .true.
952 IF (lmax > 1) RETURN
953 END DO
954
955 IF (any_defined) has_high_l = .false.
956
957 END FUNCTION dftb_kind_set_has_high_l
958
959! **************************************************************************************************
960!> \brief Initialize the qs environment (subsys)
961!> \param qs_env ...
962!> \param para_env ...
963!> \param subsys ...
964!> \param cell ...
965!> \param cell_ref ...
966!> \param use_ref_cell ...
967!> \param subsys_section ...
968!> \param silent ...
969!> \param multip ...
970!> \param charge ...
971!> \author Creation (22.05.2000,MK)
972! **************************************************************************************************
973 SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell, subsys_section, &
974 silent, multip, charge)
975
976 TYPE(qs_environment_type), POINTER :: qs_env
977 TYPE(mp_para_env_type), POINTER :: para_env
978 TYPE(qs_subsys_type), POINTER :: subsys
979 TYPE(cell_type), POINTER :: cell, cell_ref
980 LOGICAL, INTENT(in) :: use_ref_cell
981 TYPE(section_vals_type), POINTER :: subsys_section
982 LOGICAL, INTENT(in), OPTIONAL :: silent
983 INTEGER, INTENT(IN), OPTIONAL :: multip, charge
984
985 CHARACTER(len=*), PARAMETER :: routinen = 'qs_init_subsys'
986
987 CHARACTER(len=2) :: element_symbol
988 INTEGER :: gfn_type, handle, ikind, ispin, iw, lmax_sphere, maxl, maxlgto, maxlgto_lri, &
989 maxlgto_nuc, maxlppl, maxlppnl, method_id, multiplicity, my_ival, n_ao, n_mo_add, natom, &
990 nelectron, ngauss, nkind, nlumo_dos, nlumo_molden, nlumo_required, output_unit, &
991 sort_basis, tnadd_method
992 INTEGER, DIMENSION(2) :: n_mo, nelectron_spin
993 INTEGER, DIMENSION(5) :: occ
994 INTEGER, DIMENSION(:), POINTER :: mo_index_range
995 LOGICAL :: all_potential_present, be_silent, cneo_potential_present, do_kpoints, do_ri_hfx, &
996 do_ri_mp2, do_ri_rpa, do_ri_sos_mp2, do_rpa_ri_exx, do_wfc_im_time, e1terms, &
997 has_unit_metric, lribas, mp2_present, orb_gradient, paw_atom
998 REAL(kind=dp) :: alpha, ccore, ewald_rcut, fxx, maxocc, &
999 rc, rcut, total_zeff_corr, &
1000 verlet_skin, zeff_correction
1001 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1002 TYPE(cp_logger_type), POINTER :: logger
1003 TYPE(dft_control_type), POINTER :: dft_control
1004 TYPE(dftb_control_type), POINTER :: dftb_control
1005 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
1006 TYPE(ewald_environment_type), POINTER :: ewald_env
1007 TYPE(ewald_pw_type), POINTER :: ewald_pw
1008 TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env
1009 TYPE(gapw_control_type), POINTER :: gapw_control
1010 TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, &
1011 rhoin_basis, ri_aux_basis_set, &
1012 ri_hfx_basis, ri_xas_basis, &
1013 tmp_basis_set
1014 TYPE(harris_type), POINTER :: harris_env
1015 TYPE(local_rho_type), POINTER :: local_rho_set
1016 TYPE(lri_environment_type), POINTER :: lri_env
1017 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_last_converged
1018 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1019 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1020 TYPE(mp2_type), POINTER :: mp2_env
1021 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole
1022 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1023 TYPE(pw_env_type), POINTER :: pw_env
1024 TYPE(qs_control_type), POINTER :: qs_control
1025 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
1026 POINTER :: dftb_potential
1027 TYPE(qs_dispersion_type), POINTER :: dispersion_env
1028 TYPE(qs_energy_type), POINTER :: energy
1029 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1030 TYPE(qs_gcp_type), POINTER :: gcp_env
1031 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1032 TYPE(qs_kind_type), POINTER :: qs_kind
1033 TYPE(qs_ks_env_type), POINTER :: ks_env
1034 TYPE(qs_wf_history_type), POINTER :: wf_history
1035 TYPE(rho0_mpole_type), POINTER :: rho0_mpole
1036 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
1037 TYPE(scf_control_type), POINTER :: scf_control
1038 TYPE(se_taper_type), POINTER :: se_taper
1039 TYPE(section_vals_type), POINTER :: dft_section, et_coupling_section, et_ddapc_section, &
1040 ewald_section, harris_section, lri_section, mp2_section, nl_section, poisson_section, &
1041 pp_section, print_section, qs_section, rixs_section, se_section, tddfpt_section, &
1042 xc_section
1043 TYPE(semi_empirical_control_type), POINTER :: se_control
1044 TYPE(semi_empirical_si_type), POINTER :: se_store_int_env
1045 TYPE(xtb_control_type), POINTER :: xtb_control
1046
1047 CALL timeset(routinen, handle)
1048 NULLIFY (logger)
1049 logger => cp_get_default_logger()
1050 output_unit = cp_logger_get_default_io_unit(logger)
1051
1052 be_silent = .false.
1053 IF (PRESENT(silent)) be_silent = silent
1054
1055 CALL cite_reference(cp2kqs2020)
1056
1057 ! Initialise the Quickstep environment
1058 NULLIFY (mos, se_taper)
1059 NULLIFY (dft_control)
1060 NULLIFY (energy)
1061 NULLIFY (force)
1062 NULLIFY (local_molecules)
1063 NULLIFY (local_particles)
1064 NULLIFY (scf_control)
1065 NULLIFY (dft_section)
1066 NULLIFY (et_coupling_section)
1067 NULLIFY (ks_env)
1068 NULLIFY (mos_last_converged)
1069 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
1070 qs_section => section_vals_get_subs_vals(dft_section, "QS")
1071 et_coupling_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%ET_COUPLING")
1072 ! reimplemented TDDFPT
1073 tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
1074 rixs_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS")
1075
1076 CALL qs_subsys_get(subsys, particle_set=particle_set, &
1077 qs_kind_set=qs_kind_set, &
1078 atomic_kind_set=atomic_kind_set, &
1079 molecule_set=molecule_set, &
1080 molecule_kind_set=molecule_kind_set)
1081
1082 ! Read the input section with the DFT control parameters
1083 CALL read_dft_control(dft_control, dft_section, cell)
1084
1085 ! Set periodicity flag
1086 dft_control%qs_control%periodicity = sum(cell%perd)
1087
1088 ! Read the input section with the Quickstep control parameters
1089 CALL read_qs_section(dft_control%qs_control, qs_section, cell)
1090
1091 ! Print the Quickstep program banner (copyright and version number)
1092 IF (.NOT. be_silent) THEN
1093 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PROGRAM_BANNER", extension=".Log")
1094 CALL section_vals_val_get(qs_section, "METHOD", i_val=method_id)
1095 SELECT CASE (method_id)
1096 CASE DEFAULT
1097 CALL qs_header(iw)
1100 CALL se_header(iw)
1101 CASE (do_method_dftb)
1102 CALL dftb_header(iw)
1103 CASE (do_method_xtb)
1104 IF (dft_control%qs_control%xtb_control%do_tblite) THEN
1105 CALL tblite_header(iw, dft_control%qs_control%xtb_control%tblite_method)
1106 ELSE
1107 gfn_type = dft_control%qs_control%xtb_control%gfn_type
1108 CALL xtb_header(iw, gfn_type)
1109 END IF
1110 END SELECT
1111 CALL cp_print_key_finished_output(iw, logger, dft_section, &
1112 "PRINT%PROGRAM_BANNER")
1113 END IF
1114
1115 IF (dft_control%do_sccs .AND. dft_control%qs_control%gapw) THEN
1116 cpabort("SCCS is not yet implemented with GAPW")
1117 END IF
1118 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
1119 IF (do_kpoints) THEN
1120 IF (dft_control%nspins == 2 .AND. dft_control%qs_control%xtb .AND. &
1121 .NOT. dft_control%qs_control%xtb_control%do_tblite .AND. &
1122 dft_control%qs_control%xtb_control%gfn_type == gfn1xtb .AND. &
1123 dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_tblite .AND. &
1124 .NOT. dft_control%qs_control%xtb_control%tblite_mixer_damping_explicit) THEN
1125 CALL cp_warn(__location__, &
1126 "Reducing XTB/TBLITE_MIXER/DAMPING to 0.25 for CP2K-internal GFN1-xTB "// &
1127 "UKS k-point calculations with SCC_MIXER TBLITE. Set XTB/TBLITE_MIXER/DAMPING "// &
1128 "explicitly to override this conservative fallback.")
1129 dft_control%qs_control%xtb_control%tblite_mixer_damping = 0.25_dp
1130 END IF
1131 ! reset some of the settings for wfn extrapolation for kpoints
1132 SELECT CASE (dft_control%qs_control%wf_interpolation_method_nr)
1134 CALL cp_warn(__location__, "Linear WFN-based extrapolation methods are not "// &
1135 "implemented for k-points. Switching to USE_PREV_WF.")
1136 dft_control%qs_control%wf_interpolation_method_nr = wfi_use_prev_wf_method_nr
1137 END SELECT
1138 END IF
1139
1140 ! Check if any kind of electron transfer calculation has to be performed
1141 CALL section_vals_val_get(et_coupling_section, "TYPE_OF_CONSTRAINT", i_val=my_ival)
1142 dft_control%qs_control%et_coupling_calc = .false.
1143 IF (my_ival == do_et_ddapc) THEN
1144 et_ddapc_section => section_vals_get_subs_vals(et_coupling_section, "DDAPC_RESTRAINT_A")
1145 dft_control%qs_control%et_coupling_calc = .true.
1146 dft_control%qs_control%ddapc_restraint = .true.
1147 CALL read_ddapc_section(dft_control%qs_control, ddapc_restraint_section=et_ddapc_section)
1148 END IF
1149
1150 CALL read_mgrid_section(dft_control%qs_control, dft_section)
1151
1152 ! Reimplemented TDDFPT
1153 CALL read_tddfpt2_control(dft_control%tddfpt2_control, tddfpt_section, dft_control%qs_control)
1154
1155 ! RIXS
1156 CALL section_vals_get(rixs_section, explicit=qs_env%do_rixs)
1157 IF (qs_env%do_rixs) THEN
1158 CALL read_rixs_control(dft_control%rixs_control, rixs_section, dft_control%qs_control)
1159 END IF
1160
1161 ! Create relativistic control section
1162 block
1163 TYPE(rel_control_type), POINTER :: rel_control
1164 ALLOCATE (rel_control)
1165 CALL rel_c_create(rel_control)
1166 CALL rel_c_read_parameters(rel_control, dft_section)
1167 CALL set_qs_env(qs_env, rel_control=rel_control)
1168 END block
1169
1170 ! Read DFTB parameter files
1171 IF (dft_control%qs_control%method_id == do_method_dftb) THEN
1172 NULLIFY (ewald_env, ewald_pw, dftb_potential)
1173 dftb_control => dft_control%qs_control%dftb_control
1174 CALL qs_dftb_param_init(atomic_kind_set, qs_kind_set, dftb_control, dftb_potential, &
1175 subsys_section=subsys_section, para_env=para_env)
1176 CALL set_qs_env(qs_env, dftb_potential=dftb_potential)
1177 ! check for Ewald
1178 IF (dftb_control%do_ewald) THEN
1179 ALLOCATE (ewald_env)
1180 CALL ewald_env_create(ewald_env, para_env)
1181 poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
1182 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1183 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1184 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
1185 CALL get_qs_kind_set(qs_kind_set, basis_rcut=ewald_rcut)
1186 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
1187 cell_periodic=cell%perd)
1188 ALLOCATE (ewald_pw)
1189 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
1190 CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
1191 END IF
1192 ELSEIF (dft_control%qs_control%method_id == do_method_xtb) THEN
1193 ! Read xTB parameter file
1194 xtb_control => dft_control%qs_control%xtb_control
1195 CALL get_qs_env(qs_env, nkind=nkind)
1196 IF (xtb_control%do_tblite) THEN
1197 ! put geometry to tblite
1198 CALL tb_init_geometry(qs_env, qs_env%tb_tblite)
1199 ! select tblite method
1200 CALL tb_set_calculator(qs_env%tb_tblite, xtb_control%tblite_method, &
1201 xtb_control%tblite_accuracy, xtb_control%tblite_param_file)
1202 !set up wave function
1203 CALL tb_init_wf(qs_env%tb_tblite, dft_control)
1204 !get basis set
1205 DO ikind = 1, nkind
1206 qs_kind => qs_kind_set(ikind)
1207 ! Setup proper xTB parameters
1208 cpassert(.NOT. ASSOCIATED(qs_kind%xtb_parameter))
1209 CALL allocate_xtb_atom_param(qs_kind%xtb_parameter)
1210 ! Set default parameters
1211 CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
1212
1213 NULLIFY (tmp_basis_set)
1214 CALL tb_get_basis(qs_env%tb_tblite, tmp_basis_set, element_symbol, qs_kind%xtb_parameter, occ)
1215 CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
1216 CALL set_xtb_atom_param(qs_kind%xtb_parameter, occupation=occ)
1217
1218 !setting the potential for the computation
1219 zeff_correction = 0.0_dp
1220 CALL init_potential(qs_kind%all_potential, itype="BARE", &
1221 zeff=real(sum(occ), dp), zeff_correction=zeff_correction)
1222 END DO
1223 ELSE
1224 NULLIFY (ewald_env, ewald_pw)
1225 DO ikind = 1, nkind
1226 qs_kind => qs_kind_set(ikind)
1227 ! Setup proper xTB parameters
1228 cpassert(.NOT. ASSOCIATED(qs_kind%xtb_parameter))
1229 CALL allocate_xtb_atom_param(qs_kind%xtb_parameter)
1230 ! Set default parameters
1231 gfn_type = dft_control%qs_control%xtb_control%gfn_type
1232 CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
1233 CALL xtb_parameters_init(qs_kind%xtb_parameter, gfn_type, element_symbol, &
1234 xtb_control%parameter_file_path, xtb_control%parameter_file_name, &
1235 para_env)
1236 ! set dependent parameters
1237 CALL xtb_parameters_set(qs_kind%xtb_parameter)
1238 ! Generate basis set
1239 NULLIFY (tmp_basis_set)
1240 IF (qs_kind%xtb_parameter%z == 1) THEN
1241 ! special case hydrogen
1242 ngauss = xtb_control%h_sto_ng
1243 ELSE
1244 ngauss = xtb_control%sto_ng
1245 END IF
1246 IF (qs_kind%xtb_parameter%defined) THEN
1247 CALL init_xtb_basis(qs_kind%xtb_parameter, tmp_basis_set, ngauss)
1248 CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
1249 ELSE
1250 CALL set_qs_kind(qs_kind, ghost=.true.)
1251 IF (ASSOCIATED(qs_kind%all_potential)) THEN
1252 DEALLOCATE (qs_kind%all_potential%elec_conf)
1253 DEALLOCATE (qs_kind%all_potential)
1254 END IF
1255 END IF
1256 ! potential
1257 IF (qs_kind%xtb_parameter%defined) THEN
1258 zeff_correction = 0.0_dp
1259 CALL init_potential(qs_kind%all_potential, itype="BARE", &
1260 zeff=qs_kind%xtb_parameter%zeff, zeff_correction=zeff_correction)
1261 CALL get_potential(qs_kind%all_potential, alpha_core_charge=alpha)
1262 ccore = qs_kind%xtb_parameter%zeff*sqrt((alpha/pi)**3)
1263 CALL set_potential(qs_kind%all_potential, ccore_charge=ccore)
1264 qs_kind%xtb_parameter%zeff = qs_kind%xtb_parameter%zeff - zeff_correction
1265 END IF
1266 END DO
1267 !
1268 ! set repulsive potential range
1269 !
1270 ALLOCATE (xtb_control%rcpair(nkind, nkind))
1271 CALL xtb_pp_radius(qs_kind_set, xtb_control%rcpair, xtb_control%eps_pair, xtb_control%kf)
1272 ! check for Ewald
1273 IF (xtb_control%do_ewald) THEN
1274 ALLOCATE (ewald_env)
1275 CALL ewald_env_create(ewald_env, para_env)
1276 poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
1277 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1278 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1279 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
1280 IF (gfn_type == 0) THEN
1281 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
1282 silent=silent, pset="EEQ", cell_periodic=cell%perd)
1283 ELSE
1284 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
1285 silent=silent, cell_periodic=cell%perd)
1286 END IF
1287 ALLOCATE (ewald_pw)
1288 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
1289 CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
1290 END IF
1291 END IF
1292 END IF
1293 ! lri or ri env initialization
1294 lri_section => section_vals_get_subs_vals(qs_section, "LRIGPW")
1295 IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. &
1296 dft_control%qs_control%lri_optbas .OR. &
1297 dft_control%qs_control%method_id == do_method_rigpw) THEN
1298 CALL lri_env_init(lri_env, lri_section)
1299 CALL set_qs_env(qs_env, lri_env=lri_env)
1300 END IF
1301
1302 ! Check basis and fill in missing parts
1303 CALL check_qs_kind_set(qs_kind_set, dft_control, subsys_section=subsys_section)
1304
1305 ! Check that no all-electron potential is present if GPW or GAPW_XC
1306 CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
1307 IF ((dft_control%qs_control%method_id == do_method_gpw) .OR. &
1308 (dft_control%qs_control%method_id == do_method_gapw_xc) .OR. &
1309 (dft_control%qs_control%method_id == do_method_ofgpw)) THEN
1310 IF (all_potential_present) THEN
1311 cpabort("All-electron calculations with GPW, GAPW_XC, and OFGPW are not implemented")
1312 END IF
1313 END IF
1314
1315 ! Check that no cneo potential is present if not GAPW
1316 CALL get_qs_kind_set(qs_kind_set, cneo_potential_present=cneo_potential_present)
1317 IF (cneo_potential_present .AND. &
1318 dft_control%qs_control%method_id /= do_method_gapw) THEN
1319 cpabort("CNEO calculations require GAPW method")
1320 END IF
1321
1322 ! DFT+U
1323 CALL get_qs_kind_set(qs_kind_set, dft_plus_u_atom_present=dft_control%dft_plus_u)
1324
1325 IF (dft_control%do_admm) THEN
1326 ! Check if ADMM basis is available
1327 CALL get_qs_env(qs_env, nkind=nkind)
1328 DO ikind = 1, nkind
1329 NULLIFY (aux_fit_basis)
1330 qs_kind => qs_kind_set(ikind)
1331 CALL get_qs_kind(qs_kind, basis_set=aux_fit_basis, basis_type="AUX_FIT")
1332 IF (.NOT. (ASSOCIATED(aux_fit_basis))) THEN
1333 ! AUX_FIT basis set is not available
1334 cpabort("AUX_FIT basis set is not defined. ")
1335 END IF
1336 END DO
1337 END IF
1338
1339 lribas = .false.
1340 e1terms = .false.
1341 IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN
1342 lribas = .true.
1343 CALL get_qs_env(qs_env, lri_env=lri_env)
1344 e1terms = lri_env%exact_1c_terms
1345 END IF
1346 IF (dft_control%qs_control%do_kg) THEN
1347 CALL section_vals_val_get(dft_section, "KG_METHOD%TNADD_METHOD", i_val=tnadd_method)
1348 IF (tnadd_method == kg_tnadd_embed_ri) lribas = .true.
1349 END IF
1350 IF (lribas) THEN
1351 ! Check if LRI_AUX basis is available, auto-generate if needed
1352 CALL get_qs_env(qs_env, nkind=nkind)
1353 DO ikind = 1, nkind
1354 NULLIFY (lri_aux_basis)
1355 qs_kind => qs_kind_set(ikind)
1356 CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX")
1357 IF (.NOT. (ASSOCIATED(lri_aux_basis))) THEN
1358 ! LRI_AUX basis set is not yet loaded
1359 CALL cp_warn(__location__, "Automatic Generation of LRI_AUX basis. "// &
1360 "This is experimental code.")
1361 ! Generate a default basis
1362 CALL create_lri_aux_basis_set(lri_aux_basis, qs_kind, dft_control%auto_basis_lri_aux, e1terms)
1363 CALL add_basis_set_to_container(qs_kind%basis_sets, lri_aux_basis, "LRI_AUX")
1364 END IF
1365 END DO
1366 END IF
1367
1368 CALL section_vals_val_get(qs_env%input, "DFT%XC%HF%RI%_SECTION_PARAMETERS_", l_val=do_ri_hfx)
1369 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF%RI%_SECTION_PARAMETERS_", &
1370 l_val=do_rpa_ri_exx)
1371 IF (do_ri_hfx .OR. do_rpa_ri_exx) THEN
1372 CALL get_qs_env(qs_env, nkind=nkind)
1373 CALL section_vals_val_get(qs_env%input, "DFT%SORT_BASIS", i_val=sort_basis)
1374 DO ikind = 1, nkind
1375 NULLIFY (ri_hfx_basis)
1376 qs_kind => qs_kind_set(ikind)
1377 CALL get_qs_kind(qs_kind=qs_kind, basis_set=ri_hfx_basis, &
1378 basis_type="RI_HFX")
1379 IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN
1380 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
1381 IF (dft_control%do_admm) THEN
1382 CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hfx, &
1383 basis_type="AUX_FIT", basis_sort=sort_basis)
1384 ELSE
1385 CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hfx, &
1386 basis_sort=sort_basis)
1387 END IF
1388 CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HFX")
1389 END IF
1390 END DO
1391 END IF
1392
1393 IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
1394 ! Check if RI_HXC basis is available, auto-generate if needed
1395 CALL get_qs_env(qs_env, nkind=nkind)
1396 DO ikind = 1, nkind
1397 NULLIFY (ri_hfx_basis)
1398 qs_kind => qs_kind_set(ikind)
1399 CALL get_qs_kind(qs_kind, basis_set=ri_hfx_basis, basis_type="RI_HXC")
1400 IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN
1401 ! Generate a default basis
1402 CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hxc)
1403 CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HXC")
1404 END IF
1405 END DO
1406 END IF
1407
1408 ! Harris method
1409 NULLIFY (harris_env)
1410 CALL section_vals_val_get(dft_section, "HARRIS_METHOD%_SECTION_PARAMETERS_", &
1411 l_val=qs_env%harris_method)
1412 harris_section => section_vals_get_subs_vals(dft_section, "HARRIS_METHOD")
1413 CALL harris_env_create(qs_env, harris_env, harris_section)
1414 CALL set_qs_env(qs_env, harris_env=harris_env)
1415 !
1416 IF (qs_env%harris_method) THEN
1417 CALL get_qs_env(qs_env, nkind=nkind)
1418 ! Check if RI_HXC basis is available, auto-generate if needed
1419 DO ikind = 1, nkind
1420 NULLIFY (tmp_basis_set)
1421 qs_kind => qs_kind_set(ikind)
1422 CALL get_qs_kind(qs_kind, basis_set=rhoin_basis, basis_type="RHOIN")
1423 IF (.NOT. (ASSOCIATED(rhoin_basis))) THEN
1424 ! Generate a default basis
1425 CALL create_ri_aux_basis_set(tmp_basis_set, qs_kind, dft_control%auto_basis_ri_hxc)
1426 IF (qs_env%harris_env%density_source == hden_atomic) THEN
1427 CALL create_primitive_basis_set(tmp_basis_set, rhoin_basis, lmax=0)
1428 CALL deallocate_gto_basis_set(tmp_basis_set)
1429 ELSE
1430 rhoin_basis => tmp_basis_set
1431 END IF
1432 CALL add_basis_set_to_container(qs_kind%basis_sets, rhoin_basis, "RHOIN")
1433 END IF
1434 END DO
1435 END IF
1436
1437 mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
1438 CALL section_vals_get(mp2_section, explicit=mp2_present)
1439 IF (mp2_present) THEN
1440
1441 ! basis should be sorted for imaginary time RPA/GW
1442 CALL section_vals_val_get(qs_env%input, "DFT%SORT_BASIS", i_val=sort_basis)
1443 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%_SECTION_PARAMETERS_", &
1444 l_val=do_wfc_im_time)
1445
1446 IF (do_wfc_im_time .AND. sort_basis /= basis_sort_zet) THEN
1447 CALL cp_warn(__location__, &
1448 "Low-scaling RPA requires SORT_BASIS EXP keyword (in DFT input section) for good performance")
1449 END IF
1450
1451 ! Check if RI_AUX basis (for MP2/RPA) is given, auto-generate if not
1452 CALL mp2_env_create(qs_env%mp2_env)
1453 CALL get_qs_env(qs_env, mp2_env=mp2_env, nkind=nkind)
1454 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_MP2%_SECTION_PARAMETERS_", l_val=do_ri_mp2)
1455 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_SOS_MP2%_SECTION_PARAMETERS_", l_val=do_ri_sos_mp2)
1456 CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%_SECTION_PARAMETERS_", l_val=do_ri_rpa)
1457 IF (do_ri_mp2 .OR. do_ri_sos_mp2 .OR. do_ri_rpa) THEN
1458 DO ikind = 1, nkind
1459 NULLIFY (ri_aux_basis_set)
1460 qs_kind => qs_kind_set(ikind)
1461 CALL get_qs_kind(qs_kind=qs_kind, basis_set=ri_aux_basis_set, &
1462 basis_type="RI_AUX")
1463 IF (.NOT. (ASSOCIATED(ri_aux_basis_set))) THEN
1464 ! RI_AUX basis set is not yet loaded
1465 ! Generate a default basis
1466 CALL create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, dft_control%auto_basis_ri_aux, basis_sort=sort_basis)
1467 CALL add_basis_set_to_container(qs_kind%basis_sets, ri_aux_basis_set, "RI_AUX")
1468 ! Add a flag, which allows to check if the basis was generated
1469 ! when applying ERI_METHOD OS to mp2, ri-rpa, gw etc
1470 qs_env%mp2_env%ri_aux_auto_generated = .true.
1471 END IF
1472 END DO
1473 END IF
1474
1475 END IF
1476
1477 IF (dft_control%do_xas_tdp_calculation .OR. qs_env%do_rixs) THEN
1478 ! Check if RI_XAS basis is given, auto-generate if not
1479 CALL get_qs_env(qs_env, nkind=nkind)
1480 DO ikind = 1, nkind
1481 NULLIFY (ri_xas_basis)
1482 qs_kind => qs_kind_set(ikind)
1483 CALL get_qs_kind(qs_kind, basis_set=ri_xas_basis, basis_type="RI_XAS")
1484 IF (.NOT. ASSOCIATED(ri_xas_basis)) THEN
1485 ! Generate a default basis
1486 CALL create_ri_aux_basis_set(ri_xas_basis, qs_kind, dft_control%auto_basis_ri_xas)
1487 CALL add_basis_set_to_container(qs_kind%basis_sets, ri_xas_basis, "RI_XAS")
1488 END IF
1489 END DO
1490 END IF
1491
1492 ! Initialize the spherical harmonics and the orbital transformation matrices
1493 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, maxlppl=maxlppl, maxlppnl=maxlppnl)
1494
1495 ! CNEO nuclear basis contributes to GAPW rho0
1496 IF (cneo_potential_present) THEN
1497 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_nuc, basis_type="NUC")
1498 maxlgto = max(maxlgto, maxlgto_nuc)
1499 END IF
1500 lmax_sphere = dft_control%qs_control%gapw_control%lmax_sphere
1501 IF (lmax_sphere < 0) THEN
1502 lmax_sphere = 2*maxlgto
1503 dft_control%qs_control%gapw_control%lmax_sphere = lmax_sphere
1504 END IF
1505 IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
1506 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="LRI_AUX")
1507 !take maxlgto from lri basis if larger (usually)
1508 maxlgto = max(maxlgto, maxlgto_lri)
1509 ELSE IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
1510 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_HXC")
1511 maxlgto = max(maxlgto, maxlgto_lri)
1512 END IF
1513 IF (dft_control%do_xas_tdp_calculation .OR. qs_env%do_rixs) THEN
1514 !done as a precaution
1515 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_XAS")
1516 maxlgto = max(maxlgto, maxlgto_lri)
1517 END IF
1518 maxl = max(2*maxlgto, maxlppl, maxlppnl, lmax_sphere) + 1
1519
1520 CALL init_orbital_pointers(maxl)
1521
1522 CALL init_spherical_harmonics(maxl, 0)
1523
1524 ! Initialise the qs_kind_set
1525 CALL init_qs_kind_set(qs_kind_set)
1526
1527 ! Initialise GAPW soft basis and projectors
1528 IF (dft_control%qs_control%method_id == do_method_gapw .OR. &
1529 dft_control%qs_control%method_id == do_method_gapw_xc) THEN
1530 qs_control => dft_control%qs_control
1531 CALL init_gapw_basis_set(qs_kind_set, qs_control, qs_env%input)
1532 END IF
1533
1534 ! Initialise CNEO nuclear soft basis
1535 IF (cneo_potential_present) THEN
1536 CALL init_cneo_basis_set(qs_kind_set, qs_control)
1537 END IF
1538
1539 ! Initialize the pretabulation for the calculation of the
1540 ! incomplete Gamma function F_n(t) after McMurchie-Davidson
1541 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
1542 maxl = max(3*maxlgto + 1, 0)
1543 CALL init_md_ftable(maxl)
1544
1545 ! Initialize the atomic interaction radii
1546 CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
1547 !
1548 IF (dft_control%qs_control%method_id == do_method_xtb) THEN
1549 IF (.NOT. dft_control%qs_control%xtb_control%do_tblite) THEN
1550 ! cutoff radius
1551 CALL get_qs_env(qs_env, nkind=nkind)
1552 DO ikind = 1, nkind
1553 qs_kind => qs_kind_set(ikind)
1554 IF (qs_kind%xtb_parameter%defined) THEN
1555 CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set)
1556 rcut = xtb_control%coulomb_sr_cut
1557 fxx = 2.0_dp*xtb_control%coulomb_sr_eps*qs_kind%xtb_parameter%eta**2
1558 fxx = 0.80_dp*(1.0_dp/fxx)**0.3333_dp
1559 rcut = min(rcut, xtb_control%coulomb_sr_cut)
1560 qs_kind%xtb_parameter%rcut = min(rcut, fxx)
1561 ELSE
1562 qs_kind%xtb_parameter%rcut = 0.0_dp
1563 END IF
1564 END DO
1565 END IF
1566 END IF
1567
1568 IF (.NOT. be_silent) THEN
1569 CALL write_pgf_orb_radii("orb", atomic_kind_set, qs_kind_set, subsys_section)
1570 CALL write_pgf_orb_radii("aux", atomic_kind_set, qs_kind_set, subsys_section)
1571 CALL write_pgf_orb_radii("lri", atomic_kind_set, qs_kind_set, subsys_section)
1572 CALL write_pgf_orb_radii("nuc", atomic_kind_set, qs_kind_set, subsys_section)
1573 CALL write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section)
1574 CALL write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section)
1575 CALL write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section)
1576 CALL write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section)
1577 END IF
1578
1579 ! Distribute molecules and atoms using the new data structures
1580 CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, &
1581 particle_set=particle_set, &
1582 local_particles=local_particles, &
1583 molecule_kind_set=molecule_kind_set, &
1584 molecule_set=molecule_set, &
1585 local_molecules=local_molecules, &
1586 force_env_section=qs_env%input)
1587
1588 ! SCF parameters
1589 ALLOCATE (scf_control)
1590 ! set (non)-self consistency
1591 IF (dft_control%qs_control%dftb) THEN
1592 scf_control%non_selfconsistent = .NOT. dft_control%qs_control%dftb_control%self_consistent
1593 END IF
1594 IF (dft_control%qs_control%xtb) THEN
1595 IF (dft_control%qs_control%xtb_control%do_tblite) THEN
1596 scf_control%non_selfconsistent = .false.
1597 ELSE
1598 scf_control%non_selfconsistent = (dft_control%qs_control%xtb_control%gfn_type == 0)
1599 END IF
1600 END IF
1601 IF (qs_env%harris_method) THEN
1602 scf_control%non_selfconsistent = .true.
1603 END IF
1604 CALL scf_c_create(scf_control)
1605 CALL scf_c_read_parameters(scf_control, dft_section)
1606 IF (.NOT. dft_control%qs_control%do_ls_scf) THEN
1607 SELECT CASE (dft_control%qs_control%method_id)
1608 CASE (do_method_dftb)
1609 IF (dft_control%qs_control%dftb_control%tblite_scc_mixer == tblite_scc_mixer_tblite) &
1610 scf_control%max_scf = dft_control%qs_control%dftb_control%tblite_mixer_iterations
1611 CASE (do_method_xtb)
1612 IF (dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_tblite) &
1613 scf_control%max_scf = dft_control%qs_control%xtb_control%tblite_mixer_iterations
1614 END SELECT
1615 END IF
1616
1617 ! Allocate the data structure for Quickstep energies
1618 CALL allocate_qs_energy(energy)
1619
1620 ! Check for orthogonal basis
1621 has_unit_metric = .false.
1622 IF (dft_control%qs_control%semi_empirical) THEN
1623 IF (dft_control%qs_control%se_control%orthogonal_basis) has_unit_metric = .true.
1624 END IF
1625 IF (dft_control%qs_control%dftb) THEN
1626 IF (dft_control%qs_control%dftb_control%orthogonal_basis) has_unit_metric = .true.
1627 END IF
1628 CALL set_qs_env(qs_env, has_unit_metric=has_unit_metric)
1629
1630 ! Activate the interpolation
1631 CALL wfi_create(wf_history, &
1632 interpolation_method_nr= &
1633 dft_control%qs_control%wf_interpolation_method_nr, &
1634 extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
1635 has_unit_metric=has_unit_metric)
1636
1637 ! Set the current Quickstep environment
1638 CALL set_qs_env(qs_env=qs_env, &
1639 scf_control=scf_control, &
1640 wf_history=wf_history)
1641
1642 CALL qs_subsys_set(subsys, &
1643 cell_ref=cell_ref, &
1644 use_ref_cell=use_ref_cell, &
1645 energy=energy, &
1646 force=force)
1647
1648 CALL get_qs_env(qs_env, ks_env=ks_env)
1649 CALL set_ks_env(ks_env, dft_control=dft_control)
1650
1651 CALL qs_subsys_set(subsys, local_molecules=local_molecules, &
1652 local_particles=local_particles, cell=cell)
1653
1654 CALL distribution_1d_release(local_particles)
1655 CALL distribution_1d_release(local_molecules)
1656 CALL wfi_release(wf_history)
1657
1658 CALL get_qs_env(qs_env=qs_env, &
1659 atomic_kind_set=atomic_kind_set, &
1660 dft_control=dft_control, &
1661 scf_control=scf_control)
1662
1663 ! Decide what conditions need mo_derivs
1664 ! right now, this only appears to be OT
1665 IF (dft_control%qs_control%do_ls_scf .OR. &
1666 dft_control%qs_control%do_almo_scf) THEN
1667 CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.false.)
1668 ELSE
1669 IF (scf_control%use_ot) THEN
1670 CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.true.)
1671 ELSE
1672 CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.false.)
1673 END IF
1674 END IF
1675
1676 ! XXXXXXX this is backwards XXXXXXXX
1677 IF (dft_control%qs_control%xtb_control%do_tblite .AND. .NOT. scf_control%use_ot) THEN
1678 IF (.NOT. scf_control%smear%do_smear) THEN
1679 ! set tblite default smearing
1680 scf_control%smear%do_smear = .true.
1681 scf_control%smear%method = smear_fermi_dirac
1682 scf_control%smear%electronic_temperature = 300._dp/kelvin
1683 scf_control%smear%eps_fermi_dirac = 1.e-6_dp
1684 END IF
1685 END IF
1686 dft_control%smear = scf_control%smear%do_smear
1687
1688 ! Periodic efield needs equal occupation and orbital gradients
1689 IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)) THEN
1690 IF (dft_control%apply_period_efield) THEN
1691 CALL get_qs_env(qs_env=qs_env, requires_mo_derivs=orb_gradient)
1692 IF (.NOT. orb_gradient) THEN
1693 CALL cp_abort(__location__, "Periodic Efield needs orbital gradient and direct optimization."// &
1694 " Use the OT optimization method.")
1695 END IF
1696 IF (dft_control%smear) THEN
1697 CALL cp_abort(__location__, "Periodic Efield needs equal occupation numbers."// &
1698 " Smearing option is not possible.")
1699 END IF
1700 END IF
1701 END IF
1702
1703 ! Initialize the GAPW local densities and potentials
1704 IF (dft_control%qs_control%method_id == do_method_gapw .OR. &
1705 dft_control%qs_control%method_id == do_method_gapw_xc) THEN
1706 ! Allocate and initialize the set of atomic densities
1707 NULLIFY (rho_atom_set)
1708 gapw_control => dft_control%qs_control%gapw_control
1709 CALL init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
1710 CALL set_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
1711 IF (dft_control%qs_control%method_id /= do_method_gapw_xc) THEN
1712 CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, natom=natom)
1713 ! Allocate and initialize the compensation density rho0
1714 CALL init_rho0(local_rho_set, qs_env, gapw_control)
1715 ! Allocate and Initialize the local coulomb term
1716 CALL init_coulomb_local(qs_env%hartree_local, natom)
1717 END IF
1718 ! NLCC
1719 CALL init_gapw_nlcc(qs_kind_set)
1720 ! Accurate XC integration
1721 IF (gapw_control%accurate_xcint) THEN
1722 cpassert(.NOT. ASSOCIATED(gapw_control%aw))
1723 CALL get_qs_env(qs_env, nkind=nkind)
1724 ALLOCATE (gapw_control%aw(nkind))
1725 alpha = gapw_control%aweights
1726 DO ikind = 1, nkind
1727 qs_kind => qs_kind_set(ikind)
1728 CALL get_qs_kind(qs_kind, hard_radius=rc, paw_atom=paw_atom)
1729 IF (paw_atom) THEN
1730 gapw_control%aw(ikind) = alpha*(1.2_dp/rc)**2
1731 ELSE
1732 gapw_control%aw(ikind) = 0.0_dp
1733 END IF
1734 END DO
1735 END IF
1736 ELSE IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN
1737 ! allocate local ri environment
1738 ! nothing to do here?
1739 ELSE IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
1740 ! allocate ri environment
1741 ! nothing to do here?
1742 ELSE IF (dft_control%qs_control%semi_empirical) THEN
1743 NULLIFY (se_store_int_env, se_nddo_mpole, se_nonbond_env)
1744 natom = SIZE(particle_set)
1745 se_section => section_vals_get_subs_vals(qs_section, "SE")
1746 se_control => dft_control%qs_control%se_control
1747
1748 ! Make the cutoff radii choice a bit smarter
1749 CALL se_cutoff_compatible(se_control, se_section, cell, output_unit)
1750
1751 SELECT CASE (dft_control%qs_control%method_id)
1752 CASE DEFAULT
1755 ! Neighbor lists have to be MAX(interaction range, orbital range)
1756 ! set new kind radius
1757 CALL init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, subsys_section)
1758 END SELECT
1759 ! Initialize to zero the max multipole to treat in the EWALD scheme..
1760 se_control%max_multipole = do_multipole_none
1761 ! check for Ewald
1762 IF (se_control%do_ewald .OR. se_control%do_ewald_gks) THEN
1763 ALLOCATE (ewald_env)
1764 CALL ewald_env_create(ewald_env, para_env)
1765 poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
1766 CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
1767 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
1768 print_section => section_vals_get_subs_vals(qs_env%input, &
1769 "PRINT%GRID_INFORMATION")
1770 CALL read_ewald_section(ewald_env, ewald_section)
1771 ! Create ewald grids
1772 ALLOCATE (ewald_pw)
1773 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, &
1774 print_section=print_section)
1775 ! Initialize ewald grids
1776 CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
1777 ! Setup the nonbond environment (real space part of Ewald)
1778 CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
1779 ! Setup the maximum level of multipoles to be treated in the periodic SE scheme
1780 IF (se_control%do_ewald) THEN
1781 CALL ewald_env_get(ewald_env, max_multipole=se_control%max_multipole)
1782 END IF
1783 CALL section_vals_val_get(se_section, "NEIGHBOR_LISTS%VERLET_SKIN", &
1784 r_val=verlet_skin)
1785 ALLOCATE (se_nonbond_env)
1786 CALL fist_nonbond_env_create(se_nonbond_env, atomic_kind_set, do_nonbonded=.true., &
1787 do_electrostatics=.true., verlet_skin=verlet_skin, ewald_rcut=ewald_rcut, &
1788 ei_scale14=0.0_dp, vdw_scale14=0.0_dp, shift_cutoff=.false.)
1789 ! Create and Setup NDDO multipole environment
1790 CALL nddo_mpole_setup(se_nddo_mpole, natom)
1791 CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw, &
1792 se_nonbond_env=se_nonbond_env, se_nddo_mpole=se_nddo_mpole)
1793 ! Handle the residual integral part 1/R^3
1794 CALL semi_empirical_expns3_setup(qs_kind_set, se_control, &
1795 dft_control%qs_control%method_id)
1796 END IF
1797 ! Taper function
1798 CALL se_taper_create(se_taper, se_control%integral_screening, se_control%do_ewald, &
1799 se_control%taper_cou, se_control%range_cou, &
1800 se_control%taper_exc, se_control%range_exc, &
1801 se_control%taper_scr, se_control%range_scr, &
1802 se_control%taper_lrc, se_control%range_lrc)
1803 CALL set_qs_env(qs_env, se_taper=se_taper)
1804 ! Store integral environment
1805 CALL semi_empirical_si_create(se_store_int_env, se_section)
1806 CALL set_qs_env(qs_env, se_store_int_env=se_store_int_env)
1807 END IF
1808
1809 ! Initialize possible dispersion parameters
1810 IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
1811 dft_control%qs_control%method_id == do_method_gapw .OR. &
1812 dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
1813 dft_control%qs_control%method_id == do_method_lrigpw .OR. &
1814 dft_control%qs_control%method_id == do_method_rigpw .OR. &
1815 dft_control%qs_control%method_id == do_method_ofgpw) THEN
1816 ALLOCATE (dispersion_env)
1817 NULLIFY (xc_section)
1818 xc_section => section_vals_get_subs_vals(dft_section, "XC")
1819 CALL qs_dispersion_env_set(dispersion_env, xc_section)
1820 IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
1821 NULLIFY (pp_section)
1822 pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL")
1823 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
1824 ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
1825 NULLIFY (nl_section)
1826 nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
1827 CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
1828 END IF
1829 CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
1830 ELSE IF (dft_control%qs_control%method_id == do_method_dftb) THEN
1831 ALLOCATE (dispersion_env)
1832 ! set general defaults
1833 dispersion_env%doabc = .false.
1834 dispersion_env%c9cnst = .false.
1835 dispersion_env%lrc = .false.
1836 dispersion_env%srb = .false.
1837 dispersion_env%verbose = .false.
1838 NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
1839 dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
1840 dispersion_env%d3_exclude_pair)
1841 NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
1842 dispersion_env%d2y_dx2, dispersion_env%dftd_section)
1843 NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
1844 IF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d3) THEN
1845 dispersion_env%type = xc_vdw_fun_pairpot
1846 dispersion_env%pp_type = vdw_pairpot_dftd3
1847 dispersion_env%eps_cn = dftb_control%epscn
1848 dispersion_env%s6 = dftb_control%sd3(1)
1849 dispersion_env%sr6 = dftb_control%sd3(2)
1850 dispersion_env%s8 = dftb_control%sd3(3)
1851 dispersion_env%domol = .false.
1852 dispersion_env%kgc8 = 0._dp
1853 dispersion_env%rc_disp = dftb_control%rcdisp
1854 dispersion_env%exp_pre = 0._dp
1855 dispersion_env%scaling = 0._dp
1856 dispersion_env%nd3_exclude_pair = 0
1857 dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
1858 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
1859 ELSEIF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d3bj) THEN
1860 dispersion_env%type = xc_vdw_fun_pairpot
1861 dispersion_env%pp_type = vdw_pairpot_dftd3bj
1862 dispersion_env%eps_cn = dftb_control%epscn
1863 dispersion_env%s6 = dftb_control%sd3bj(1)
1864 dispersion_env%a1 = dftb_control%sd3bj(2)
1865 dispersion_env%s8 = dftb_control%sd3bj(3)
1866 dispersion_env%a2 = dftb_control%sd3bj(4)
1867 dispersion_env%domol = .false.
1868 dispersion_env%kgc8 = 0._dp
1869 dispersion_env%rc_disp = dftb_control%rcdisp
1870 dispersion_env%exp_pre = 0._dp
1871 dispersion_env%scaling = 0._dp
1872 dispersion_env%nd3_exclude_pair = 0
1873 dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
1874 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
1875 ELSEIF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d2) THEN
1876 dispersion_env%type = xc_vdw_fun_pairpot
1877 dispersion_env%pp_type = vdw_pairpot_dftd2
1878 dispersion_env%exp_pre = dftb_control%exp_pre
1879 dispersion_env%scaling = dftb_control%scaling
1880 dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
1881 dispersion_env%rc_disp = dftb_control%rcdisp
1882 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
1883 ELSE
1884 dispersion_env%type = xc_vdw_fun_none
1885 END IF
1886 CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
1887 ELSE IF (dft_control%qs_control%method_id == do_method_xtb) THEN
1888 IF (.NOT. (dft_control%qs_control%xtb_control%do_tblite)) THEN
1889 ALLOCATE (dispersion_env)
1890 ! set general defaults
1891 dispersion_env%doabc = .false.
1892 dispersion_env%c9cnst = .false.
1893 dispersion_env%lrc = .false.
1894 dispersion_env%srb = .false.
1895 dispersion_env%verbose = .false.
1896 NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, &
1897 dispersion_env%r0ab, dispersion_env%rcov, &
1898 dispersion_env%r2r4, dispersion_env%cn, &
1899 dispersion_env%cnkind, dispersion_env%cnlist, &
1900 dispersion_env%d3_exclude_pair)
1901 NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
1902 dispersion_env%d2y_dx2, dispersion_env%dftd_section)
1903 NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
1904 dispersion_env%type = xc_vdw_fun_pairpot
1905 dispersion_env%eps_cn = xtb_control%epscn
1906 dispersion_env%s6 = xtb_control%s6
1907 dispersion_env%s8 = xtb_control%s8
1908 dispersion_env%a1 = xtb_control%a1
1909 dispersion_env%a2 = xtb_control%a2
1910 dispersion_env%domol = .false.
1911 dispersion_env%kgc8 = 0._dp
1912 dispersion_env%rc_disp = xtb_control%rcdisp
1913 dispersion_env%rc_d4 = xtb_control%rcdisp
1914 dispersion_env%exp_pre = 0._dp
1915 dispersion_env%scaling = 0._dp
1916 dispersion_env%nd3_exclude_pair = 0
1917 dispersion_env%parameter_file_name = xtb_control%dispersion_parameter_file
1918 !
1919 SELECT CASE (xtb_control%vdw_type)
1921 dispersion_env%pp_type = vdw_pairpot_dftd3bj
1922 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
1923 IF (xtb_control%vdw_type == xtb_vdw_type_none) dispersion_env%type = xc_vdw_fun_none
1924 CASE (xtb_vdw_type_d4)
1925 dispersion_env%pp_type = vdw_pairpot_dftd4
1926 dispersion_env%ref_functional = "none"
1927 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, &
1928 dispersion_env, para_env=para_env)
1929 dispersion_env%cnfun = 2
1930 CASE DEFAULT
1931 cpabort("vdw type")
1932 END SELECT
1933 CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
1934 END IF
1935 ELSE IF (dft_control%qs_control%semi_empirical) THEN
1936 ALLOCATE (dispersion_env)
1937 ! set general defaults
1938 dispersion_env%doabc = .false.
1939 dispersion_env%c9cnst = .false.
1940 dispersion_env%lrc = .false.
1941 dispersion_env%srb = .false.
1942 dispersion_env%verbose = .false.
1943 NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
1944 dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
1945 dispersion_env%d3_exclude_pair)
1946 NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
1947 dispersion_env%d2y_dx2, dispersion_env%dftd_section)
1948 NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
1949 IF (se_control%dispersion) THEN
1950 dispersion_env%type = xc_vdw_fun_pairpot
1951 dispersion_env%pp_type = vdw_pairpot_dftd3
1952 dispersion_env%eps_cn = se_control%epscn
1953 dispersion_env%s6 = se_control%sd3(1)
1954 dispersion_env%sr6 = se_control%sd3(2)
1955 dispersion_env%s8 = se_control%sd3(3)
1956 dispersion_env%domol = .false.
1957 dispersion_env%kgc8 = 0._dp
1958 dispersion_env%rc_disp = se_control%rcdisp
1959 dispersion_env%exp_pre = 0._dp
1960 dispersion_env%scaling = 0._dp
1961 dispersion_env%nd3_exclude_pair = 0
1962 dispersion_env%parameter_file_name = se_control%dispersion_parameter_file
1963 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
1964 ELSE
1965 dispersion_env%type = xc_vdw_fun_none
1966 END IF
1967 CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
1968 END IF
1969
1970 ! Initialize possible geomertical counterpoise correction potential
1971 IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
1972 dft_control%qs_control%method_id == do_method_gapw .OR. &
1973 dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
1974 dft_control%qs_control%method_id == do_method_lrigpw .OR. &
1975 dft_control%qs_control%method_id == do_method_rigpw .OR. &
1976 dft_control%qs_control%method_id == do_method_ofgpw) THEN
1977 ALLOCATE (gcp_env)
1978 NULLIFY (xc_section)
1979 xc_section => section_vals_get_subs_vals(dft_section, "XC")
1980 CALL qs_gcp_env_set(gcp_env, xc_section)
1981 CALL qs_gcp_init(qs_env, gcp_env)
1982 CALL set_qs_env(qs_env, gcp_env=gcp_env)
1983 END IF
1984
1985 ! Allocate the MO data types
1986 CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron)
1987
1988 ! The total number of electrons
1989 IF (PRESENT(charge)) THEN
1990 dft_control%charge = charge
1991 nelectron = nelectron - dft_control%charge
1992 ELSE
1993 nelectron = nelectron - dft_control%charge
1994 END IF
1995
1996 IF (dft_control%multiplicity == 0) THEN
1997 IF (modulo(nelectron, 2) == 0) THEN
1998 dft_control%multiplicity = 1
1999 ELSE
2000 dft_control%multiplicity = 2
2001 END IF
2002 END IF
2003
2004 multiplicity = dft_control%multiplicity
2005
2006 IF (PRESENT(multip)) THEN
2007 multiplicity = multip
2008 END IF
2009
2010 IF ((dft_control%nspins < 1) .OR. (dft_control%nspins > 2)) THEN
2011 cpabort("nspins should be 1 or 2 for the time being ...")
2012 END IF
2013
2014 IF ((modulo(nelectron, 2) /= 0) .AND. (dft_control%nspins == 1)) THEN
2015 IF (.NOT. dft_control%qs_control%ofgpw .AND. .NOT. dft_control%smear) THEN
2016 cpabort("Use the LSD option for an odd number of electrons")
2017 END IF
2018 END IF
2019
2020 ! The transition potential method to calculate XAS needs LSD
2021 IF (dft_control%do_xas_calculation) THEN
2022 IF (dft_control%nspins == 1) THEN
2023 cpabort("Use the LSD option for XAS with transition potential")
2024 END IF
2025 END IF
2026
2027 ! assigning the number of states per spin initial version, not yet very
2028 ! general. Should work for an even number of electrons and a single
2029 ! additional electron this set of options that requires full matrices,
2030 ! however, makes things a bit ugly right now.... we try to make a
2031 ! distinction between the number of electrons per spin and the number of
2032 ! MOs per spin this should allow the use of fractional occupations later on
2033 IF (dft_control%qs_control%ofgpw) THEN
2034
2035 IF (dft_control%nspins == 1) THEN
2036 maxocc = nelectron
2037 nelectron_spin(1) = nelectron
2038 nelectron_spin(2) = 0
2039 n_mo(1) = 1
2040 n_mo(2) = 0
2041 ELSE
2042 nelectron_spin(1) = (nelectron + multiplicity - 1)/2
2043 nelectron_spin(2) = (nelectron - multiplicity + 1)/2
2044 IF (nelectron_spin(1) < 0) THEN
2045 cpabort("LSD: too few electrons for this multiplicity")
2046 END IF
2047 maxocc = maxval(nelectron_spin)
2048 n_mo(1) = min(nelectron_spin(1), 1)
2049 n_mo(2) = min(nelectron_spin(2), 1)
2050 END IF
2051
2052 ELSE
2053
2054 IF (dft_control%nspins == 1) THEN
2055 maxocc = 2.0_dp
2056 nelectron_spin(1) = nelectron
2057 nelectron_spin(2) = 0
2058 IF (modulo(nelectron, 2) == 0) THEN
2059 n_mo(1) = nelectron/2
2060 ELSE
2061 n_mo(1) = int(nelectron/2._dp) + 1
2062 END IF
2063 n_mo(2) = 0
2064 ELSE
2065 maxocc = 1.0_dp
2066
2067 ! The simplist spin distribution is written here. Special cases will
2068 ! need additional user input
2069 IF (modulo(nelectron + multiplicity - 1, 2) /= 0) THEN
2070 cpabort("LSD: try to use a different multiplicity")
2071 END IF
2072
2073 nelectron_spin(1) = (nelectron + multiplicity - 1)/2
2074 nelectron_spin(2) = (nelectron - multiplicity + 1)/2
2075
2076 IF (nelectron_spin(2) < 0) THEN
2077 cpabort("LSD: too few electrons for this multiplicity")
2078 END IF
2079
2080 n_mo(1) = nelectron_spin(1)
2081 n_mo(2) = nelectron_spin(2)
2082
2083 END IF
2084
2085 END IF
2086
2087 ! Read the total_zeff_corr here [SGh]
2088 CALL get_qs_kind_set(qs_kind_set, total_zeff_corr=total_zeff_corr)
2089 ! store it in qs_env
2090 qs_env%total_zeff_corr = total_zeff_corr
2091
2092 ! Store the number of electrons once and for all
2093 CALL qs_subsys_set(subsys, &
2094 nelectron_total=nelectron, &
2095 nelectron_spin=nelectron_spin)
2096
2097 ! Ensure that all orbitals requested for printout are added even
2098 ! if the keyword ADDED_MOS was not specified or set properly
2099 mo_index_range => section_get_ivals(dft_section, "PRINT%MO%MO_INDEX_RANGE")
2100 cpassert(ASSOCIATED(mo_index_range))
2101 IF (all(mo_index_range > 0)) THEN
2102 IF (mo_index_range(1) > mo_index_range(2)) THEN
2103 CALL cp_abort(__location__, &
2104 "The upper orbital index ("// &
2105 trim(adjustl(cp_to_string(mo_index_range(2))))// &
2106 ") of the MO_INDEX_RANGE should be equal or larger "// &
2107 "than the lower orbital index ("// &
2108 trim(adjustl(cp_to_string(mo_index_range(1))))// &
2109 ") for printout.")
2110 END IF
2111 ! Adapt ADDED_MOS automatically if needed for printout
2112 IF (.NOT. scf_control%use_ot) THEN
2113 scf_control%added_mos(1) = min(max(scf_control%added_mos(1), &
2114 mo_index_range(2) - n_mo(1)), &
2115 n_ao - n_mo(1))
2116 IF (dft_control%nspins == 2) THEN
2117 scf_control%added_mos(2) = min(max(scf_control%added_mos(2), &
2118 mo_index_range(2) - n_mo(2)), &
2119 n_ao - n_mo(2))
2120 END IF
2121 END IF
2122 ELSE IF (mo_index_range(2) < 0) THEN
2123 IF (.NOT. scf_control%use_ot) THEN
2124 ! Add all available orbitals
2125 scf_control%added_mos(1) = n_ao - n_mo(1)
2126 IF (dft_control%nspins == 2) THEN
2127 ! Ensure the same number for the spin-down (beta) orbitals
2128 scf_control%added_mos(2) = n_ao - n_mo(2)
2129 END IF
2130 END IF
2131 END IF
2132
2133 nlumo_dos = section_get_ival(dft_section, "PRINT%DOS%NLUMO")
2134 nlumo_molden = section_get_ival(dft_section, "PRINT%MO_MOLDEN%NLUMO")
2135 nlumo_required = max(nlumo_dos, nlumo_molden)
2136 IF (nlumo_dos == -1 .OR. nlumo_molden == -1) nlumo_required = -1
2137 IF (.NOT. scf_control%use_ot .AND. nlumo_required /= 0) THEN
2138 IF (nlumo_required == -1) THEN
2139 IF (scf_control%added_mos(1) /= -1 .OR. &
2140 (dft_control%nspins == 2 .AND. scf_control%added_mos(2) /= -1)) THEN
2141 CALL cp_warn(__location__, &
2142 "NLUMO requested by DOS/PDOS/Molden exceeds SCF%ADDED_MOS. "// &
2143 "For diagonalization calculations, SCF%ADDED_MOS is "// &
2144 "increased to provide the requested unoccupied orbitals.")
2145 END IF
2146 scf_control%added_mos(1) = -1
2147 IF (dft_control%nspins == 2) scf_control%added_mos(2) = -1
2148 ELSE
2149 IF (scf_control%added_mos(1) >= 0 .AND. &
2150 nlumo_required > scf_control%added_mos(1)) THEN
2151 CALL cp_warn(__location__, &
2152 "NLUMO requested by DOS/PDOS/Molden exceeds SCF%ADDED_MOS. "// &
2153 "For diagonalization calculations, SCF%ADDED_MOS is "// &
2154 "increased to provide the requested unoccupied orbitals.")
2155 scf_control%added_mos(1) = nlumo_required
2156 END IF
2157 IF (dft_control%nspins == 2 .AND. scf_control%added_mos(2) > 0 .AND. &
2158 nlumo_required > scf_control%added_mos(2)) THEN
2159 scf_control%added_mos(2) = nlumo_required
2160 END IF
2161 END IF
2162 END IF
2163
2164 IF (dft_control%nspins == 2) THEN
2165 ! Check and set number of added (unoccupied) orbitals for beta spin
2166 IF (scf_control%added_mos(2) < 0) THEN
2167 n_mo_add = n_ao - n_mo(2) ! use all available MOs
2168 ELSE IF (scf_control%added_mos(2) > 0) THEN
2169 n_mo_add = scf_control%added_mos(2)
2170 ELSE
2171 n_mo_add = scf_control%added_mos(1)
2172 END IF
2173 IF (n_mo_add > n_ao - n_mo(2)) THEN
2174 cpwarn("More ADDED_MOs requested for beta spin than available.")
2175 END IF
2176 scf_control%added_mos(2) = min(n_mo_add, n_ao - n_mo(2))
2177 n_mo(2) = n_mo(2) + scf_control%added_mos(2)
2178 END IF
2179
2180 ! proceed alpha orbitals after the beta orbitals; this is essential to avoid
2181 ! reduction in the number of available unoccupied molecular orbitals.
2182 ! E.g. n_ao = 10, nelectrons = 10, multiplicity = 3 implies n_mo(1) = 6, n_mo(2) = 4;
2183 ! added_mos(1:2) = (6,undef) should increase the number of molecular orbitals as
2184 ! n_mo(1) = min(n_ao, n_mo(1) + added_mos(1)) = 10, n_mo(2) = 10.
2185 ! However, if we try to proceed alpha orbitals first, this leads us n_mo(1:2) = (10,8)
2186 ! due to the following assignment instruction above:
2187 ! IF (scf_control%added_mos(2) > 0) THEN ... ELSE; n_mo_add = scf_control%added_mos(1); END IF
2188 IF (dft_control%qs_control%xtb_control%do_tblite .AND. .NOT. scf_control%use_ot) THEN
2189 scf_control%added_mos(1) = n_ao - n_mo(1) ! tblite needs all MO's
2190 ELSE IF (scf_control%added_mos(1) < 0) THEN
2191 scf_control%added_mos(1) = n_ao - n_mo(1) ! use all available MOs
2192 ELSE IF (scf_control%added_mos(1) > n_ao - n_mo(1)) THEN
2193 CALL cp_warn(__location__, &
2194 "More added MOs requested than available. "// &
2195 "The full set of unoccupied MOs will be used. "// &
2196 "Use 'ADDED_MOS -1' to always use all available MOs "// &
2197 "and to get rid of this warning.")
2198 END IF
2199 scf_control%added_mos(1) = min(scf_control%added_mos(1), n_ao - n_mo(1))
2200 n_mo(1) = n_mo(1) + scf_control%added_mos(1)
2201
2202 IF (dft_control%nspins == 2) THEN
2203 IF (n_mo(2) > n_mo(1)) &
2204 CALL cp_warn(__location__, &
2205 "More beta than alpha MOs requested. "// &
2206 "The number of beta MOs will be reduced to the number alpha MOs.")
2207 n_mo(2) = min(n_mo(1), n_mo(2))
2208 cpassert(n_mo(1) >= nelectron_spin(1))
2209 cpassert(n_mo(2) >= nelectron_spin(2))
2210 END IF
2211
2212 ! kpoints
2213 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
2214 IF (do_kpoints .AND. dft_control%nspins == 2) THEN
2215 ! we need equal number of calculated states
2216 IF (n_mo(2) /= n_mo(1)) &
2217 CALL cp_warn(__location__, &
2218 "Kpoints: Different number of MOs requested. "// &
2219 "The number of beta MOs will be set to the number alpha MOs.")
2220 n_mo(2) = n_mo(1)
2221 cpassert(n_mo(1) >= nelectron_spin(1))
2222 cpassert(n_mo(2) >= nelectron_spin(2))
2223 END IF
2224
2225 ! Compatibility checks for smearing
2226 IF (scf_control%smear%do_smear) THEN
2227 IF (dft_control%qs_control%xtb_control%do_tblite .AND. scf_control%use_ot) THEN
2228 cpabort("CP2K/tblite with OT does not support smearing.")
2229 END IF
2230 IF (scf_control%added_mos(1) == 0) THEN
2231 cpabort("Extra MOs (ADDED_MOS) are required for smearing")
2232 END IF
2233 END IF
2234
2235 ! Some options require that all MOs are computed ...
2236 IF (btest(cp_print_key_should_output(logger%iter_info, dft_section, &
2237 "PRINT%MO/CARTESIAN"), &
2238 cp_p_file) .OR. &
2239 (scf_control%level_shift /= 0.0_dp) .OR. &
2240 (scf_control%diagonalization%eps_jacobi /= 0.0_dp) .OR. &
2241 (dft_control%roks .AND. (.NOT. scf_control%use_ot))) THEN
2242 n_mo(:) = n_ao
2243 END IF
2244
2245 ! Compatibility checks for ROKS
2246 IF (dft_control%roks .AND. (.NOT. scf_control%use_ot)) THEN
2247 IF (scf_control%roks_scheme == general_roks) THEN
2248 cpwarn("General ROKS scheme is not yet tested!")
2249 END IF
2250 IF (scf_control%smear%do_smear) THEN
2251 CALL cp_abort(__location__, &
2252 "The options ROKS and SMEAR are not compatible. "// &
2253 "Try UKS instead of ROKS")
2254 END IF
2255 END IF
2256 IF (dft_control%low_spin_roks) THEN
2257 SELECT CASE (dft_control%qs_control%method_id)
2258 CASE DEFAULT
2260 CALL cp_abort(__location__, &
2261 "xTB/DFTB methods are not compatible with low spin ROKS.")
2264 CALL cp_abort(__location__, &
2265 "SE methods are not compatible with low spin ROKS.")
2266 END SELECT
2267 END IF
2268
2269 ! in principle the restricted calculation could be performed
2270 ! using just one set of MOs and special casing most of the code
2271 ! right now we'll just take care of what is effectively an additional constraint
2272 ! at as few places as possible, just duplicating the beta orbitals
2273 IF (dft_control%restricted .AND. (output_unit > 0)) THEN
2274 ! it is really not yet tested till the end ! Joost
2275 WRITE (output_unit, *) ""
2276 WRITE (output_unit, *) " **************************************"
2277 WRITE (output_unit, *) " restricted calculation cutting corners"
2278 WRITE (output_unit, *) " experimental feature, check code "
2279 WRITE (output_unit, *) " **************************************"
2280 END IF
2281
2282 ! no point in allocating these things here ?
2283 IF (dft_control%qs_control%do_ls_scf) THEN
2284 NULLIFY (mos)
2285 ELSE
2286 ALLOCATE (mos(dft_control%nspins))
2287 DO ispin = 1, dft_control%nspins
2288 CALL allocate_mo_set(mo_set=mos(ispin), &
2289 nao=n_ao, &
2290 nmo=n_mo(ispin), &
2291 nelectron=nelectron_spin(ispin), &
2292 n_el_f=real(nelectron_spin(ispin), dp), &
2293 maxocc=maxocc, &
2294 flexible_electron_count=dft_control%relax_multiplicity)
2295 END DO
2296 END IF
2297
2298 CALL set_qs_env(qs_env, mos=mos)
2299
2300 ! allocate mos when switch_surf_dip is triggered [SGh]
2301 IF (dft_control%switch_surf_dip) THEN
2302 ALLOCATE (mos_last_converged(dft_control%nspins))
2303 DO ispin = 1, dft_control%nspins
2304 CALL allocate_mo_set(mo_set=mos_last_converged(ispin), &
2305 nao=n_ao, &
2306 nmo=n_mo(ispin), &
2307 nelectron=nelectron_spin(ispin), &
2308 n_el_f=real(nelectron_spin(ispin), dp), &
2309 maxocc=maxocc, &
2310 flexible_electron_count=dft_control%relax_multiplicity)
2311 END DO
2312 CALL set_qs_env(qs_env, mos_last_converged=mos_last_converged)
2313 END IF
2314
2315 IF (.NOT. be_silent) THEN
2316 ! Print the DFT control parameters
2317 IF (PRESENT(multip)) THEN
2318 dft_control%multiplicity = multiplicity
2319 END IF
2320 CALL write_dft_control(dft_control, dft_section)
2321
2322 ! Print the vdW control parameters
2323 IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
2324 dft_control%qs_control%method_id == do_method_gapw .OR. &
2325 dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
2326 dft_control%qs_control%method_id == do_method_lrigpw .OR. &
2327 dft_control%qs_control%method_id == do_method_rigpw .OR. &
2328 dft_control%qs_control%method_id == do_method_dftb .OR. &
2329 (dft_control%qs_control%method_id == do_method_xtb .AND. &
2330 (.NOT. dft_control%qs_control%xtb_control%do_tblite)) .OR. &
2331 dft_control%qs_control%method_id == do_method_ofgpw) THEN
2332 CALL get_qs_env(qs_env, dispersion_env=dispersion_env)
2333 CALL qs_write_dispersion(qs_env, dispersion_env)
2334 END IF
2335
2336 ! Print the Quickstep control parameters
2337 CALL write_qs_control(dft_control%qs_control, dft_section)
2338
2339 ! Print the ADMM control parameters
2340 IF (dft_control%do_admm) THEN
2341 CALL write_admm_control(dft_control%admm_control, dft_section)
2342 END IF
2343
2344 ! Print XES/XAS control parameters
2345 IF (dft_control%do_xas_calculation) THEN
2346 CALL cite_reference(iannuzzi2007)
2347 !CALL write_xas_control(dft_control%xas_control,dft_section)
2348 END IF
2349
2350 ! Print the unnormalized basis set information (input data)
2351 CALL write_gto_basis_sets(qs_kind_set, subsys_section)
2352
2353 ! Print the atomic kind set
2354 CALL write_qs_kind_set(qs_kind_set, subsys_section)
2355
2356 ! Print the molecule kind set
2357 CALL write_molecule_kind_set(molecule_kind_set, subsys_section)
2358
2359 ! Print the total number of kinds, atoms, basis functions etc.
2360 CALL write_total_numbers(qs_kind_set, particle_set, qs_env%input)
2361
2362 ! Print the atomic coordinates
2363 CALL write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label="QUICKSTEP")
2364
2365 ! Print the interatomic distances
2366 CALL write_particle_distances(particle_set, cell, subsys_section)
2367
2368 ! Print the requested structure data
2369 CALL write_structure_data(particle_set, cell, subsys_section)
2370
2371 ! Print symmetry information
2372 CALL write_symmetry(particle_set, cell, subsys_section)
2373
2374 ! Print the SCF parameters
2375 IF ((.NOT. dft_control%qs_control%do_ls_scf) .AND. &
2376 (.NOT. dft_control%qs_control%do_almo_scf)) THEN
2377 CALL scf_c_write_parameters(scf_control, dft_section)
2378 END IF
2379 END IF
2380
2381 ! Sets up pw_env, qs_charges, mpools ...
2382 CALL qs_env_setup(qs_env)
2383
2384 ! Allocate and initialise rho0 soft on the global grid
2385 IF (dft_control%qs_control%method_id == do_method_gapw) THEN
2386 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho0_mpole=rho0_mpole)
2387 CALL rho0_s_grid_create(pw_env, rho0_mpole)
2388 END IF
2389
2390 IF (output_unit > 0) CALL m_flush(output_unit)
2391 CALL timestop(handle)
2392
2393 END SUBROUTINE qs_init_subsys
2394
2395! **************************************************************************************************
2396!> \brief Write the total number of kinds, atoms, etc. to the logical unit
2397!> number lunit.
2398!> \param qs_kind_set ...
2399!> \param particle_set ...
2400!> \param force_env_section ...
2401!> \author Creation (06.10.2000)
2402! **************************************************************************************************
2403 SUBROUTINE write_total_numbers(qs_kind_set, particle_set, force_env_section)
2404
2405 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2406 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2407 TYPE(section_vals_type), POINTER :: force_env_section
2408
2409 INTEGER :: maxlgto, maxlppl, maxlppnl, natom, &
2410 natom_q, ncgf, nkind, nkind_q, npgf, &
2411 nset, nsgf, nshell, output_unit
2412 TYPE(cp_logger_type), POINTER :: logger
2413
2414 NULLIFY (logger)
2415 logger => cp_get_default_logger()
2416 output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%TOTAL_NUMBERS", &
2417 extension=".Log")
2418
2419 IF (output_unit > 0) THEN
2420 natom = SIZE(particle_set)
2421 nkind = SIZE(qs_kind_set)
2422
2423 CALL get_qs_kind_set(qs_kind_set, &
2424 maxlgto=maxlgto, &
2425 ncgf=ncgf, &
2426 npgf=npgf, &
2427 nset=nset, &
2428 nsgf=nsgf, &
2429 nshell=nshell, &
2430 maxlppl=maxlppl, &
2431 maxlppnl=maxlppnl)
2432
2433 WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
2434 "TOTAL NUMBERS AND MAXIMUM NUMBERS"
2435
2436 IF (nset + npgf + ncgf > 0) THEN
2437 WRITE (unit=output_unit, fmt="(/,T3,A,(T30,A,T71,I10))") &
2438 "Total number of", &
2439 "- Atomic kinds: ", nkind, &
2440 "- Atoms: ", natom, &
2441 "- Shell sets: ", nset, &
2442 "- Shells: ", nshell, &
2443 "- Primitive Cartesian functions: ", npgf, &
2444 "- Cartesian basis functions: ", ncgf, &
2445 "- Spherical basis functions: ", nsgf
2446 ELSE IF (nshell + nsgf > 0) THEN
2447 WRITE (unit=output_unit, fmt="(/,T3,A,(T30,A,T71,I10))") &
2448 "Total number of", &
2449 "- Atomic kinds: ", nkind, &
2450 "- Atoms: ", natom, &
2451 "- Shells: ", nshell, &
2452 "- Spherical basis functions: ", nsgf
2453 ELSE
2454 WRITE (unit=output_unit, fmt="(/,T3,A,(T30,A,T71,I10))") &
2455 "Total number of", &
2456 "- Atomic kinds: ", nkind, &
2457 "- Atoms: ", natom
2458 END IF
2459
2460 IF ((maxlppl > -1) .AND. (maxlppnl > -1)) THEN
2461 WRITE (unit=output_unit, fmt="(/,T3,A,(T30,A,T75,I6))") &
2462 "Maximum angular momentum of the", &
2463 "- Orbital basis functions: ", maxlgto, &
2464 "- Local part of the GTH pseudopotential: ", maxlppl, &
2465 "- Non-local part of the GTH pseudopotential: ", maxlppnl
2466 ELSEIF (maxlppl > -1) THEN
2467 WRITE (unit=output_unit, fmt="(/,T3,A,(T30,A,T75,I6))") &
2468 "Maximum angular momentum of the", &
2469 "- Orbital basis functions: ", maxlgto, &
2470 "- Local part of the GTH pseudopotential: ", maxlppl
2471 ELSE
2472 WRITE (unit=output_unit, fmt="(/,T3,A,T75,I6)") &
2473 "Maximum angular momentum of the orbital basis functions: ", maxlgto
2474 END IF
2475
2476 ! LRI_AUX BASIS
2477 CALL get_qs_kind_set(qs_kind_set, &
2478 maxlgto=maxlgto, &
2479 ncgf=ncgf, &
2480 npgf=npgf, &
2481 nset=nset, &
2482 nsgf=nsgf, &
2483 nshell=nshell, &
2484 basis_type="LRI_AUX")
2485 IF (nset + npgf + ncgf > 0) THEN
2486 WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2487 "LRI_AUX Basis: ", &
2488 "Total number of", &
2489 "- Shell sets: ", nset, &
2490 "- Shells: ", nshell, &
2491 "- Primitive Cartesian functions: ", npgf, &
2492 "- Cartesian basis functions: ", ncgf, &
2493 "- Spherical basis functions: ", nsgf
2494 WRITE (unit=output_unit, fmt="(T30,A,T75,I6)") &
2495 " Maximum angular momentum ", maxlgto
2496 END IF
2497
2498 ! RI_HXC BASIS
2499 CALL get_qs_kind_set(qs_kind_set, &
2500 maxlgto=maxlgto, &
2501 ncgf=ncgf, &
2502 npgf=npgf, &
2503 nset=nset, &
2504 nsgf=nsgf, &
2505 nshell=nshell, &
2506 basis_type="RI_HXC")
2507 IF (nset + npgf + ncgf > 0) THEN
2508 WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2509 "RI_HXC Basis: ", &
2510 "Total number of", &
2511 "- Shell sets: ", nset, &
2512 "- Shells: ", nshell, &
2513 "- Primitive Cartesian functions: ", npgf, &
2514 "- Cartesian basis functions: ", ncgf, &
2515 "- Spherical basis functions: ", nsgf
2516 WRITE (unit=output_unit, fmt="(T30,A,T75,I6)") &
2517 " Maximum angular momentum ", maxlgto
2518 END IF
2519
2520 ! AUX_FIT BASIS
2521 CALL get_qs_kind_set(qs_kind_set, &
2522 maxlgto=maxlgto, &
2523 ncgf=ncgf, &
2524 npgf=npgf, &
2525 nset=nset, &
2526 nsgf=nsgf, &
2527 nshell=nshell, &
2528 basis_type="AUX_FIT")
2529 IF (nset + npgf + ncgf > 0) THEN
2530 WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2531 "AUX_FIT ADMM-Basis: ", &
2532 "Total number of", &
2533 "- Shell sets: ", nset, &
2534 "- Shells: ", nshell, &
2535 "- Primitive Cartesian functions: ", npgf, &
2536 "- Cartesian basis functions: ", ncgf, &
2537 "- Spherical basis functions: ", nsgf
2538 WRITE (unit=output_unit, fmt="(T30,A,T75,I6)") &
2539 " Maximum angular momentum ", maxlgto
2540 END IF
2541
2542 ! NUCLEAR BASIS
2543 CALL get_qs_kind_set(qs_kind_set, &
2544 nkind_q=nkind_q, &
2545 natom_q=natom_q, &
2546 maxlgto=maxlgto, &
2547 ncgf=ncgf, &
2548 npgf=npgf, &
2549 nset=nset, &
2550 nsgf=nsgf, &
2551 nshell=nshell, &
2552 basis_type="NUC")
2553 IF (nset + npgf + ncgf > 0) THEN
2554 WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
2555 "Nuclear Basis: ", &
2556 "Total number of", &
2557 "- Quantum atomic kinds: ", nkind_q, &
2558 "- Quantum atoms: ", natom_q, &
2559 "- Shell sets: ", nset, &
2560 "- Shells: ", nshell, &
2561 "- Primitive Cartesian functions: ", npgf, &
2562 "- Cartesian basis functions: ", ncgf, &
2563 "- Spherical basis functions: ", nsgf
2564 WRITE (unit=output_unit, fmt="(T30,A,T75,I6)") &
2565 " Maximum angular momentum ", maxlgto
2566 END IF
2567
2568 END IF
2569 CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
2570 "PRINT%TOTAL_NUMBERS")
2571
2572 END SUBROUTINE write_total_numbers
2573
2574END MODULE qs_environment
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
almo_scf_env methods
subroutine, public almo_scf_env_create(qs_env)
Creation and basic initialization of the almo environment.
calculate the orbitals for a given atomic kind type
subroutine, public calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
...
Define the atomic kind types and their sub types.
Automatic generation of auxiliary basis sets of different kind.
Definition auto_basis.F:15
subroutine, public create_lri_aux_basis_set(lri_aux_basis_set, qs_kind, basis_cntrl, exact_1c_terms, tda_kernel)
Create a LRI_AUX basis set using some heuristics.
Definition auto_basis.F:233
subroutine, public create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis_type, basis_sort)
Create a RI_AUX basis set using some heuristics.
Definition auto_basis.F:57
subroutine, public add_basis_set_to_container(container, basis_set, basis_set_type)
...
integer, parameter, public basis_sort_zet
subroutine, public deallocate_gto_basis_set(gto_basis_set)
...
subroutine, public create_primitive_basis_set(basis_set, pbasis, lmax)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public cp2kqs2020
integer, save, public iannuzzi2007
integer, save, public iannuzzi2006
Handles all functions related to the CELL.
Definition cell_types.F:15
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Utilities to set up the control types.
subroutine, public write_qs_control(qs_control, dft_section)
Purpose: Write the QS control parameters to the output unit.
subroutine, public read_rixs_control(rixs_control, rixs_section, qs_control)
Reads the input and stores in the rixs_control_type.
subroutine, public read_qs_section(qs_control, qs_section, cell)
...
subroutine, public read_tddfpt2_control(t_control, t_section, qs_control)
Read TDDFPT-related input parameters.
subroutine, public read_dft_control(dft_control, dft_section, cell)
...
subroutine, public write_admm_control(admm_control, dft_section)
Write the ADMM control parameters to the output unit.
subroutine, public write_dft_control(dft_control, dft_section)
Write the DFT control parameters to the output unit.
subroutine, public read_ddapc_section(qs_control, qs_section, ddapc_restraint_section)
reads the input parameters needed for ddapc.
subroutine, public read_mgrid_section(qs_control, dft_section)
...
contains information regarding the decoupling/recoupling method of Bloechl
subroutine, public cp_ddapc_ewald_create(cp_ddapc_ewald, qmmm_decoupl, qm_cell, force_env_section, subsys_section, para_env)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
types that represent a subsys, i.e. a part of the system
Work with symmetry.
Definition cp_symmetry.F:13
subroutine, public write_symmetry(particle_set, cell, input_section)
Write symmetry information to output.
Definition cp_symmetry.F:61
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public distribution_1d_release(distribution_1d)
releases the given distribution_1d
Distribution methods for atoms, particles, or molecules.
subroutine, public distribute_molecules_1d(atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, force_env_section, prev_molecule_kind_set, prev_local_molecules)
Distribute molecules and particles.
Types needed for a for a Energy Correction.
Energy correction environment setup and handling.
subroutine, public ec_write_input(ec_env)
Print out the energy correction input section.
subroutine, public ec_env_create(qs_env, ec_env, dft_section, ec_section)
Allocates and intitializes ec_env.
Definition and initialisation of the et_coupling data type.
subroutine, public et_coupling_create(et_coupling)
...
subroutine, public ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, gmax, ns_max, precs, o_spline, para_env, poisson_section, interaction_cutoffs, cell_hmat)
Purpose: Set the EWALD environment.
subroutine, public ewald_env_create(ewald_env, para_env)
allocates and intitializes a ewald_env
subroutine, public read_ewald_section(ewald_env, ewald_section)
Purpose: read the EWALD section.
subroutine, public read_ewald_section_tb(ewald_env, ewald_section, hmat, silent, pset, cell_periodic)
Purpose: read the EWALD section for TB methods.
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
subroutine, public ewald_pw_grid_update(ewald_pw, ewald_env, cell_hmat)
Rescales pw_grids for given box, if necessary.
subroutine, public ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section)
creates the structure ewald_pw_type
Types for excited states potential energies.
subroutine, public exstate_create(ex_env, excited_state, dft_section)
Allocates and intitializes exstate_env.
Definition of the atomic potential types.
subroutine, public fist_nonbond_env_create(fist_nonbond_env, atomic_kind_set, potparm14, potparm, do_nonbonded, do_electrostatics, verlet_skin, ewald_rcut, ei_scale14, vdw_scale14, shift_cutoff)
allocates and intitializes a fist_nonbond_env
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
subroutine, public init_md_ftable(nmax)
Initialize a table of F_n(t) values in the range 0 <= t <= 12 with a stepsize of 0....
Definition gamma.F:540
Define type storing the global information of a run. Keep the amount of stored data small....
subroutine, public init_coulomb_local(hartree_local, natom)
...
subroutine, public dftb_header(iw)
...
Definition header.F:193
subroutine, public qs_header(iw)
...
Definition header.F:310
subroutine, public se_header(iw)
...
Definition header.F:286
subroutine, public xtb_header(iw, gfn_type)
...
Definition header.F:217
subroutine, public tblite_header(iw, tb_type)
...
Definition header.F:252
Types and set/get functions for HFX.
Definition hfx_types.F:15
subroutine, public hfx_create(x_data, para_env, hfx_section, atomic_kind_set, qs_kind_set, particle_set, dft_control, cell, orb_basis, ri_basis, nelectron_total, nkp_grid)
This routine allocates and initializes all types in hfx_data
Definition hfx_types.F:604
subroutine, public compare_hfx_sections(hfx_section1, hfx_section2, is_identical, same_except_frac)
Compares the non-technical parts of two HFX input section and check whether they are the same Ignore ...
Definition hfx_types.F:3025
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public smear_fermi_dirac
integer, parameter, public do_method_ofgpw
integer, parameter, public xc_vdw_fun_none
integer, parameter, public dispersion_d3
integer, parameter, public gfn1xtb
integer, parameter, public do_method_rigpw
integer, parameter, public do_method_gpw
integer, parameter, public vdw_pairpot_dftd3
integer, parameter, public do_method_pdg
integer, parameter, public wfi_linear_wf_method_nr
integer, parameter, public tddfpt_kernel_none
integer, parameter, public linear_response_run
integer, parameter, public wfi_linear_ps_method_nr
integer, parameter, public do_method_pnnl
integer, parameter, public xc_vdw_fun_nonloc
integer, parameter, public vdw_pairpot_dftd4
integer, parameter, public xtb_vdw_type_d3
integer, parameter, public kg_tnadd_embed_ri
integer, parameter, public dispersion_d3bj
integer, parameter, public debug_run
integer, parameter, public tblite_scc_mixer_tblite
integer, parameter, public do_method_rm1
integer, parameter, public do_qmmm_swave
integer, parameter, public do_method_pm3
integer, parameter, public do_method_mndo
integer, parameter, public xtb_vdw_type_d4
integer, parameter, public do_method_gapw
integer, parameter, public do_method_mndod
integer, parameter, public rel_trans_atom
integer, parameter, public xtb_vdw_type_none
integer, parameter, public do_method_am1
integer, parameter, public do_method_dftb
integer, parameter, public wfi_use_prev_wf_method_nr
integer, parameter, public do_qmmm_gauss
integer, parameter, public rel_none
integer, parameter, public general_roks
integer, parameter, public vdw_pairpot_dftd2
integer, parameter, public do_method_lrigpw
integer, parameter, public dispersion_d2
integer, parameter, public do_method_xtb
integer, parameter, public do_method_pm6fm
integer, parameter, public xc_vdw_fun_pairpot
integer, parameter, public do_et_ddapc
integer, parameter, public vdw_pairpot_dftd3bj
integer, parameter, public do_method_gapw_xc
integer, parameter, public do_method_pm6
integer, parameter, public hden_atomic
objects that represent the structure of input sections and the data contained in an input section
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
integer function, public section_get_ival(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Routines for a Kim-Gordon-like partitioning into molecular subunits.
subroutine, public kg_env_create(qs_env, kg_env, qs_kind_set, input)
Allocates and intitializes kg_env.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Routines needed for kpoint calculation.
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group)
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
subroutine, public kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
Initialize the kpoint environment.
Types and basic routines needed for a kpoint calculation.
subroutine, public set_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Set information in a kpoint environment.
subroutine, public kpoint_reset_initialization(kpoint)
Reset all data derived from a concrete k-point initialization. Input options such as the scheme,...
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
subroutine, public write_kpoint_info(kpoint, iounit, dft_section)
Write information on the kpoints to output.
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
subroutine, public read_kpoint_section(kpoint, kpoint_section, a_vec, cell)
Read the kpoint input section.
initializes the environment for lri lri : local resolution of the identity
subroutine, public lri_env_init(lri_env, lri_section)
initializes the lri env
subroutine, public lri_env_basis(ri_type, qs_env, lri_env, qs_kind_set)
initializes the lri env
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
Define the molecule kind structure types and the corresponding functionality.
subroutine, public write_molecule_kind_set(molecule_kind_set, subsys_section)
Write a moleculeatomic kind set data set to the output unit.
Define the data structure for the molecule information.
Types needed for MP2 calculations.
Definition mp2_setup.F:14
subroutine, public read_mp2_section(input, mp2_env)
...
Definition mp2_setup.F:62
Types needed for MP2 calculations.
Definition mp2_types.F:14
subroutine, public mp2_env_create(mp2_env)
...
Definition mp2_types.F:471
Multipole structure: for multipole (fixed and induced) in FF based MD.
integer, parameter, public do_multipole_none
Provides Cartesian and spherical orbital pointers and indices.
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
subroutine, public init_spherical_harmonics(maxl, output_unit)
Initialize or update the orbital transformation matrices.
Define methods related to particle_type.
subroutine, public write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label)
Write the atomic coordinates to the output unit.
subroutine, public write_structure_data(particle_set, cell, input_section)
Write structure data requested by a separate structure data input section to the output unit....
subroutine, public write_particle_distances(particle_set, cell, subsys_section)
Write the matrix of the particle distances to the output unit.
Define the data structure for the particle information.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public kelvin
Definition physcon.F:165
container for various plainwaves related things
subroutine, public qs_basis_rotation(qs_env, kpoints, basis_type)
Construct basis set rotation matrices.
subroutine, public qs_dftb_param_init(atomic_kind_set, qs_kind_set, dftb_control, dftb_potential, subsys_section, para_env)
...
Definition of the DFTB parameter types.
Working with the DFTB parameter types.
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Calculation of non local dispersion functionals Some routines adapted from: Copyright (C) 2001-2009 Q...
subroutine, public qs_dispersion_nonloc_init(dispersion_env, para_env)
...
Calculation of dispersion using pair potentials.
subroutine, public qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
...
Definition of disperson types for DFT calculations.
Set disperson types for DFT calculations.
subroutine, public qs_dispersion_env_set(dispersion_env, xc_section)
...
subroutine, public qs_write_dispersion(qs_env, dispersion_env, ounit)
...
subroutine, public allocate_qs_energy(qs_energy)
Allocate and/or initialise a Quickstep energy data structure.
qs_environment methods that use many other modules
subroutine, public qs_env_setup(qs_env)
initializes various components of the qs_env, that need only atomic_kind_set, cell,...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
subroutine, public qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, qmmm, qmmm_env_qm, force_env_section, subsys_section, use_motion_section, silent, multip, charge)
Read the input and the database files for the setup of the QUICKSTEP environment.
Definition of gCP types for DFT calculations.
Set disperson types for DFT calculations.
subroutine, public qs_gcp_env_set(gcp_env, xc_section)
...
subroutine, public qs_gcp_init(qs_env, gcp_env)
...
Types needed for a for a Harris model calculation.
subroutine, public harris_rhoin_init(rhoin, basis_type, qs_kind_set, atomic_kind_set, local_particles, nspin)
...
Harris method environment setup and handling.
subroutine, public harris_write_input(harris_env)
Print out the Harris method input section.
subroutine, public harris_env_create(qs_env, harris_env, harris_section)
Allocates and intitializes harris_env.
Calculate the interaction radii for the operator matrix calculation.
subroutine, public write_pgf_orb_radii(basis, atomic_kind_set, qs_kind_set, subsys_section)
Write the orbital basis function radii to the output unit.
subroutine, public write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section)
Write the radii of the one center projector.
subroutine, public write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section)
Write the radii of the projector functions of the Goedecker pseudopotential (GTH, non-local part) to ...
subroutine, public init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, subsys_section)
...
subroutine, public write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section)
Write the radii of the exponential functions of the Goedecker pseudopotential (GTH,...
subroutine, public init_interaction_radii(qs_control, qs_kind_set)
Initialize all the atomic kind radii for a given threshold value.
subroutine, public write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section)
Write the radii of the core charge distributions to the output unit.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
subroutine, public init_gapw_nlcc(qs_kind_set)
...
subroutine, public init_cneo_basis_set(qs_kind_set, qs_control)
...
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
subroutine, public write_qs_kind_set(qs_kind_set, subsys_section)
Write an atomic kind set data set to the output unit.
subroutine, public init_gapw_basis_set(qs_kind_set, qs_control, force_env_section, modify_qs_control)
...
subroutine, public init_qs_kind_set(qs_kind_set)
Initialise an atomic kind set data set.
subroutine, public check_qs_kind_set(qs_kind_set, dft_control, subsys_section)
...
subroutine, public write_gto_basis_sets(qs_kind_set, subsys_section)
Write all the GTO basis sets of an atomic kind set to the output unit (for the printing of the unnorm...
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
subroutine, public qs_ks_env_create(ks_env)
Allocates a new instance of ks_env.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public rho0_s_grid_create(pw_env, rho0_mpole)
...
subroutine, public init_rho0(local_rho_set, qs_env, gapw_control, zcore)
...
subroutine, public init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
Routines that work on qs_subsys_type.
subroutine, public qs_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, use_motion_section, cp_subsys, elkind, silent)
Creates a qs_subsys. Optionally an existsing cp_subsys is used.
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
subroutine, public qs_subsys_set(subsys, cp_subsys, local_particles, local_molecules, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, nelectron_total, nelectron_spin)
...
Storage of past states of the qs_env. Methods to interpolate (or actually normally extrapolate) the n...
subroutine, public wfi_create_for_kp(wf_history)
Adapts wf_history storage flags for k-point calculations. For ASPC, switches from Gamma WFN storage t...
subroutine, public wfi_create(wf_history, interpolation_method_nr, extrapolation_order, has_unit_metric)
...
interpolate the wavefunctions to speed up the convergence when doing MD
subroutine, public wfi_release(wf_history)
releases a wf_history of a wavefunction (see doc/ReferenceCounting.html)
parameters that control a relativistic calculation
subroutine, public rel_c_create(rel_control)
allocates and initializes an rel control object with the default values
subroutine, public rel_c_read_parameters(rel_control, dft_section)
reads the parameters of the relativistic section into the given rel_control
parameters that control an scf iteration
subroutine, public scf_c_read_parameters(scf_control, inp_section)
reads the parameters of the scf section into the given scf_control
subroutine, public scf_c_write_parameters(scf_control, dft_section)
writes out the scf parameters
subroutine, public scf_c_create(scf_control)
allocates and initializes an scf control object with the default values
Methods for handling the 1/R^3 residual integral part.
subroutine, public semi_empirical_expns3_setup(qs_kind_set, se_control, method_id)
Setup the quantity necessary to handle the slowly convergent residual integral term 1/R^3.
Arrays of parameters used in the semi-empirical calculations \References Everywhere in this module TC...
subroutine, public init_se_intd_array()
Initialize all arrays used for the evaluation of the integrals.
Setup and Methods for semi-empirical multipole types.
subroutine, public nddo_mpole_setup(nddo_mpole, natom)
Setup NDDO multipole type.
Definition of the semi empirical multipole integral expansions types.
Type to store integrals for semi-empirical calculations.
subroutine, public semi_empirical_si_create(store_int_env, se_section, compression)
Allocate semi-empirical store integrals type.
Definition of the semi empirical parameter types.
subroutine, public se_taper_create(se_taper, integral_screening, do_ewald, taper_cou, range_cou, taper_exc, range_exc, taper_scr, range_scr, taper_lrc, range_lrc)
Creates the taper type used in SE calculations.
Working with the semi empirical parameter types.
subroutine, public se_cutoff_compatible(se_control, se_section, cell, output_unit)
Reset cutoffs trying to be somehow a bit smarter.
interface to tblite
subroutine, public tb_init_wf(tb, dft_control)
initialize wavefunction ...
subroutine, public tb_init_geometry(qs_env, tb)
intialize geometry objects ...
subroutine, public tb_set_calculator(tb, typ, accuracy, param_file)
...
subroutine, public tb_get_basis(tb, gto_basis_set, element_symbol, param, occ)
...
routines for DFT+NEGF calculations (coupling with the quantum transport code OMEN)
Definition transport.F:19
subroutine, public transport_env_create(qs_env)
creates the transport environment
Definition transport.F:122
Read xTB parameters.
subroutine, public xtb_parameters_set(param)
Read atom parameters for xTB Hamiltonian from input file.
subroutine, public xtb_parameters_init(param, gfn_type, element_symbol, parameter_file_path, parameter_file_name, para_env)
...
subroutine, public init_xtb_basis(param, gto_basis_set, ngauss)
...
xTB (repulsive) pair potentials Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov JCTC 1...
subroutine, public xtb_pp_radius(qs_kind_set, ppradius, eps_pair, kfparam)
...
Definition of the xTB parameter types.
Definition xtb_types.F:20
subroutine, public allocate_xtb_atom_param(xtb_parameter)
...
Definition xtb_types.F:99
subroutine, public set_xtb_atom_param(xtb_parameter, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, xi, kappa0, alpg, electronegativity, occupation, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:299
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represents a system: atoms, molecules, their pos,vel,...
structure to store local (to a processor) ordered lists of integers.
Contains information on the energy correction functional for KG.
Contains information on the excited states energy.
contains the initially parsed file and the initial parallel environment
Contains information about kpoints.
stores all the informations relevant to an mpi environment
contained for different pw related things
Contains information on the Harris method.
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps track of the previous wavefunctions and can extrapolate them for the next step of md
contains the parameters needed by a relativistic calculation
Global Multipolar NDDO information type.
Taper type use in semi-empirical calculations.