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