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