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