LCOV - code coverage report
Current view: top level - src - sirius_interface.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:3158929) Lines: 330 385 85.7 %
Date: 2025-07-22 22:41:15 Functions: 7 7 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : !***************************************************************************************************
       9             : !> \brief Interface to the SIRIUS Library
      10             : !> \par History
      11             : !>      07.2018 initial create
      12             : !> \author JHU
      13             : !***************************************************************************************************
      14             : #if defined(__SIRIUS)
      15             : MODULE sirius_interface
      16             :    USE ISO_C_BINDING,                   ONLY: C_DOUBLE,&
      17             :                                               C_INT,&
      18             :                                               C_LOC
      19             :    USE atom_kind_orbitals,              ONLY: calculate_atomic_orbitals,&
      20             :                                               gth_potential_conversion
      21             :    USE atom_types,                      ONLY: atom_gthpot_type
      22             :    USE atom_upf,                        ONLY: atom_upfpot_type
      23             :    USE atom_utils,                      ONLY: atom_local_potential
      24             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      25             :                                               get_atomic_kind
      26             :    USE cell_types,                      ONLY: cell_type,&
      27             :                                               real_to_scaled
      28             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      29             :                                               cp_logger_get_default_io_unit,&
      30             :                                               cp_logger_type
      31             :    USE cp_output_handling,              ONLY: cp_p_file,&
      32             :                                               cp_print_key_finished_output,&
      33             :                                               cp_print_key_should_output,&
      34             :                                               cp_print_key_unit_nr
      35             :    USE external_potential_types,        ONLY: gth_potential_type
      36             :    USE input_constants,                 ONLY: do_gapw_log
      37             :    USE input_cp2k_pwdft,                ONLY: SIRIUS_FUNC_VDWDF,&
      38             :                                               SIRIUS_FUNC_VDWDF2,&
      39             :                                               SIRIUS_FUNC_VDWDFCX
      40             :    USE input_section_types,             ONLY: section_vals_get,&
      41             :                                               section_vals_get_subs_vals,&
      42             :                                               section_vals_get_subs_vals2,&
      43             :                                               section_vals_type,&
      44             :                                               section_vals_val_get
      45             :    USE kinds,                           ONLY: default_string_length,&
      46             :                                               dp
      47             :    USE machine,                         ONLY: m_flush
      48             :    USE mathconstants,                   ONLY: fourpi,&
      49             :                                               gamma1
      50             :    USE message_passing,                 ONLY: mp_para_env_type
      51             :    USE particle_types,                  ONLY: particle_type
      52             :    USE physcon,                         ONLY: massunit
      53             :    USE pwdft_environment_types,         ONLY: pwdft_energy_type,&
      54             :                                               pwdft_env_get,&
      55             :                                               pwdft_env_set,&
      56             :                                               pwdft_environment_type
      57             :    USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
      58             :                                               create_grid_atom,&
      59             :                                               deallocate_grid_atom,&
      60             :                                               grid_atom_type
      61             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      62             :                                               qs_kind_type
      63             :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
      64             :                                               qs_subsys_type
      65             :    USE sirius,                          ONLY: &
      66             :         SIRIUS_INTEGER_ARRAY_TYPE, SIRIUS_INTEGER_TYPE, SIRIUS_LOGICAL_ARRAY_TYPE, &
      67             :         SIRIUS_LOGICAL_TYPE, SIRIUS_NUMBER_ARRAY_TYPE, SIRIUS_NUMBER_TYPE, &
      68             :         SIRIUS_STRING_ARRAY_TYPE, SIRIUS_STRING_TYPE, sirius_add_atom, sirius_add_atom_type, &
      69             :         sirius_add_atom_type_radial_function, sirius_add_xc_functional, sirius_context_handler, &
      70             :         sirius_create_context, sirius_create_ground_state, sirius_create_kset_from_grid, &
      71             :         sirius_finalize, sirius_find_ground_state, sirius_get_band_energies, &
      72             :         sirius_get_band_occupancies, sirius_get_energy, sirius_get_forces, &
      73             :         sirius_get_kpoint_properties, sirius_get_num_kpoints, sirius_get_parameters, &
      74             :         sirius_get_stress_tensor, sirius_ground_state_handler, sirius_import_parameters, &
      75             :         sirius_initialize, sirius_initialize_context, sirius_kpoint_set_handler, &
      76             :         sirius_option_get_info, sirius_option_get_section_length, sirius_option_set, &
      77             :         sirius_set_atom_position, sirius_set_atom_type_dion, sirius_set_atom_type_hubbard, &
      78             :         sirius_set_atom_type_radial_grid, sirius_set_lattice_vectors, sirius_set_mpi_grid_dims, &
      79             :         sirius_update_ground_state
      80             : #include "./base/base_uses.f90"
      81             : 
      82             :    IMPLICIT NONE
      83             : 
      84             :    PRIVATE
      85             : 
      86             : !     *** Global parameters ***
      87             : 
      88             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'sirius_interface'
      89             : 
      90             : !     *** Public subroutines ***
      91             : 
      92             :    PUBLIC :: cp_sirius_init, cp_sirius_finalize
      93             :    PUBLIC :: cp_sirius_create_env, cp_sirius_energy_force, cp_sirius_update_context
      94             : 
      95             : CONTAINS
      96             : 
      97             : !***************************************************************************************************
      98             : !> \brief ...
      99             : !> \param
     100             : !> \par History
     101             : !>      07.2018 start the Sirius library
     102             : !> \author JHU
     103             : ! **************************************************************************************************
     104        9238 :    SUBROUTINE cp_sirius_init()
     105        9238 :       CALL sirius_initialize(.FALSE.)
     106        9238 :    END SUBROUTINE cp_sirius_init
     107             : 
     108             : !***************************************************************************************************
     109             : !> \brief ...
     110             : !> \param
     111             : !> \par History
     112             : !>      07.2018 stop the Sirius library
     113             : !> \author JHU
     114             : ! **************************************************************************************************
     115        9238 :    SUBROUTINE cp_sirius_finalize()
     116        9238 :       CALL sirius_finalize(.FALSE., .FALSE., .FALSE.)
     117        9238 :    END SUBROUTINE cp_sirius_finalize
     118             : 
     119             : !***************************************************************************************************
     120             : !> \brief ...
     121             : !> \param pwdft_env ...
     122             : !> \param
     123             : !> \par History
     124             : !>      07.2018 Create the Sirius environment
     125             : !> \author JHU
     126             : ! **************************************************************************************************
     127          16 :    SUBROUTINE cp_sirius_create_env(pwdft_env)
     128             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     129             : #if defined(__SIRIUS)
     130             : 
     131             :       CHARACTER(len=2)                                   :: element_symbol
     132             :       CHARACTER(len=default_string_length)               :: label
     133             :       INTEGER                                            :: i, iatom, ibeta, ifun, ikind, iwf, j, l, &
     134             :                                                             n, natom, nbeta, nkind, nmesh, &
     135             :                                                             num_mag_dims, sirius_mpi_comm, vdw_func, nu, lu, output_unit
     136          16 :       INTEGER, DIMENSION(:), POINTER                     :: mpi_grid_dims
     137             :       INTEGER(KIND=C_INT), DIMENSION(3)                  :: k_grid, k_shift
     138          16 :       INTEGER, DIMENSION(:), POINTER                     :: kk
     139             :       LOGICAL                                            :: up, use_ref_cell
     140             :       LOGICAL(4)                                         :: use_so, use_symmetry, dft_plus_u_atom
     141          16 :       REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:)     :: fun
     142          16 :       REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :)  :: dion
     143             :       REAL(KIND=C_DOUBLE), DIMENSION(3)                  :: a1, a2, a3, v1, v2
     144             :       REAL(KIND=dp)                                      :: al, angle1, angle2, cval, focc, &
     145             :                                                             magnetization, mass, pf, rl, zeff, alpha_u, beta_u, &
     146             :                                                             J0_u, J_u, U_u, occ_u, u_minus_J
     147          16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: beta, corden, ef, fe, locpot, rc, rp
     148             :       REAL(KIND=dp), DIMENSION(3)                        :: vr, vs, j_t
     149          16 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: density
     150          16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: wavefunction, wfninfo
     151             :       TYPE(atom_gthpot_type), POINTER                    :: gth_atompot
     152             :       TYPE(atom_upfpot_type), POINTER                    :: upf_pot
     153          16 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     154             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     155             :       TYPE(cell_type), POINTER                           :: my_cell
     156             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     157             :       TYPE(grid_atom_type), POINTER                      :: atom_grid
     158             :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     159          16 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     160          16 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     161             :       TYPE(qs_subsys_type), POINTER                      :: qs_subsys
     162             :       TYPE(section_vals_type), POINTER                   :: pwdft_section, pwdft_sub_section, &
     163             :                                                             xc_fun, xc_section
     164             :       TYPE(sirius_context_handler)                       :: sctx
     165             :       TYPE(sirius_ground_state_handler)                  :: gs_handler
     166             :       TYPE(sirius_kpoint_set_handler)                    :: ks_handler
     167             : 
     168           0 :       CPASSERT(ASSOCIATED(pwdft_env))
     169             : 
     170          16 :       output_unit = cp_logger_get_default_io_unit()
     171             :       ! create context of simulation
     172          16 :       CALL pwdft_env_get(pwdft_env, para_env=para_env)
     173          16 :       sirius_mpi_comm = para_env%get_handle()
     174          16 :       CALL sirius_create_context(sirius_mpi_comm, sctx)
     175             : 
     176             : !     the "fun" starts.
     177             : 
     178          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_section, xc_input=xc_section)
     179             : 
     180             :       CALL section_vals_val_get(pwdft_section, "ignore_convergence_failure", &
     181          16 :                                 l_val=pwdft_env%ignore_convergence_failure)
     182             :       ! cp2k should *have* a function that return all xc_functionals. Doing
     183             :       ! manually is prone to errors
     184             : 
     185          16 :       IF (ASSOCIATED(xc_section)) THEN
     186          16 :          ifun = 0
     187          30 :          DO
     188          46 :             ifun = ifun + 1
     189          46 :             xc_fun => section_vals_get_subs_vals2(xc_section, i_section=ifun)
     190          46 :             IF (.NOT. ASSOCIATED(xc_fun)) EXIT
     191             :             ! Here, we do not have to check whether the functional name starts with XC_
     192             :             ! because we only allow the shorter form w/o XC_
     193          46 :             CALL sirius_add_xc_functional(sctx, "XC_"//TRIM(xc_fun%section%name))
     194             :          END DO
     195             :       END IF
     196             : 
     197             :       !     import control section
     198          16 :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "control")
     199          16 :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     200          16 :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "control")
     201          16 :          CALL section_vals_val_get(pwdft_sub_section, "mpi_grid_dims", i_vals=mpi_grid_dims)
     202             :       END IF
     203             : 
     204             : !     import parameters section
     205          16 :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "parameters")
     206             : 
     207          16 :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     208          16 :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "parameters")
     209          16 :          CALL section_vals_val_get(pwdft_sub_section, "ngridk", i_vals=kk)
     210          16 :          k_grid(1) = kk(1)
     211          16 :          k_grid(2) = kk(2)
     212          16 :          k_grid(3) = kk(3)
     213             : 
     214          16 :          CALL section_vals_val_get(pwdft_sub_section, "shiftk", i_vals=kk)
     215          16 :          k_shift(1) = kk(1)
     216          16 :          k_shift(2) = kk(2)
     217          16 :          k_shift(3) = kk(3)
     218          16 :          CALL section_vals_val_get(pwdft_sub_section, "num_mag_dims", i_val=num_mag_dims)
     219          16 :          CALL section_vals_val_get(pwdft_sub_section, "use_symmetry", l_val=use_symmetry)
     220          16 :          CALL section_vals_val_get(pwdft_sub_section, "so_correction", l_val=use_so)
     221             : 
     222             : ! now check if van der walls corrections are needed
     223             :          vdw_func = -1
     224             : #ifdef __LIBVDWXC
     225          16 :          CALL section_vals_val_get(pwdft_sub_section, "vdw_functional", i_val=vdw_func)
     226           0 :          SELECT CASE (vdw_func)
     227             :          CASE (SIRIUS_FUNC_VDWDF)
     228           0 :             CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF")
     229             :          CASE (SIRIUS_FUNC_VDWDF2)
     230           0 :             CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF2")
     231             :          CASE (SIRIUS_FUNC_VDWDFCX)
     232          16 :             CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF2")
     233             :          CASE default
     234             :          END SELECT
     235             : #endif
     236             : 
     237             :       END IF
     238             : 
     239             : !     import mixer section
     240          16 :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "mixer")
     241          16 :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     242          16 :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "mixer")
     243             :       END IF
     244             : 
     245             : !     import settings section
     246          16 :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "settings")
     247          16 :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     248          16 :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "settings")
     249             :       END IF
     250             : 
     251             :       !     import solver section
     252          16 :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "iterative_solver")
     253          16 :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     254          16 :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "iterative_solver")
     255             :       END IF
     256             : 
     257             : #if defined(__SIRIUS_DFTD4)
     258             :       ! import dftd3 and dftd4 section
     259             :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd4")
     260             :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     261             :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "dftd4")
     262             :       END IF
     263             : 
     264             :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd3")
     265             :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     266             :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "dftd3")
     267             :       END IF
     268             : #endif
     269             : 
     270             :       !
     271             :       ! uncomment these lines when nlcg is officially supported
     272             :       !
     273             : #if defined(__SIRIUS_NLCG)
     274             :       !     import nlcg section
     275             :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "nlcg")
     276             :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     277             :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "nlcg")
     278             :       END IF
     279             : #endif
     280             : 
     281             : #if defined(__SIRIUS_VCSQNM)
     282             :       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "vcsqnm")
     283             :       IF (ASSOCIATED(pwdft_sub_section)) THEN
     284             :          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "vcsqnm")
     285             :       END IF
     286             : #endif
     287             : 
     288             :       !CALL sirius_dump_runtime_setup(sctx, "runtime.json")
     289          16 :       CALL sirius_import_parameters(sctx, '{}')
     290             : 
     291             : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
     292          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
     293          16 :       CALL qs_subsys_get(qs_subsys, cell=my_cell, use_ref_cell=use_ref_cell)
     294          64 :       a1(:) = my_cell%hmat(:, 1)
     295          64 :       a2(:) = my_cell%hmat(:, 2)
     296          64 :       a3(:) = my_cell%hmat(:, 3)
     297          16 :       CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
     298             : 
     299          16 :       IF (use_ref_cell) THEN
     300           0 :          CPWARN("SIRIUS| The specified CELL_REF will be ignored for PW_DFT calculations")
     301             :       END IF
     302             : 
     303             : ! set up the atomic type definitions
     304             :       CALL qs_subsys_get(qs_subsys, &
     305             :                          atomic_kind_set=atomic_kind_set, &
     306             :                          qs_kind_set=qs_kind_set, &
     307          16 :                          particle_set=particle_set)
     308          16 :       nkind = SIZE(atomic_kind_set)
     309          38 :       DO ikind = 1, nkind
     310             :          CALL get_atomic_kind(atomic_kind_set(ikind), &
     311          22 :                               name=label, element_symbol=element_symbol, mass=mass)
     312          22 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
     313          22 :          NULLIFY (upf_pot, gth_potential)
     314          22 :          CALL get_qs_kind(qs_kind_set(ikind), upf_potential=upf_pot, gth_potential=gth_potential)
     315             : 
     316          22 :          IF (ASSOCIATED(upf_pot)) THEN
     317             :             CALL sirius_add_atom_type(sctx, label, fname=upf_pot%filename, &
     318             :                                       symbol=element_symbol, &
     319          20 :                                       mass=REAL(mass/massunit, KIND=C_DOUBLE))
     320             : 
     321           2 :          ELSEIF (ASSOCIATED(gth_potential)) THEN
     322             : !
     323           2 :             NULLIFY (atom_grid)
     324           2 :             CALL allocate_grid_atom(atom_grid)
     325           2 :             nmesh = 929
     326           2 :             atom_grid%nr = nmesh
     327           2 :             CALL create_grid_atom(atom_grid, nmesh, 1, 1, 0, do_gapw_log)
     328           8 :             ALLOCATE (rp(nmesh), fun(nmesh))
     329           2 :             IF (atom_grid%rad(1) < atom_grid%rad(nmesh)) THEN
     330             :                up = .TRUE.
     331             :             ELSE
     332             :                up = .FALSE.
     333             :             END IF
     334             :             IF (up) THEN
     335           0 :                rp(1:nmesh) = atom_grid%rad(1:nmesh)
     336             :             ELSE
     337        1860 :                DO i = 1, nmesh
     338        1860 :                   rp(i) = atom_grid%rad(nmesh - i + 1)
     339             :                END DO
     340             :             END IF
     341             : ! add new atom type
     342             :             CALL sirius_add_atom_type(sctx, label, &
     343             :                                       zn=NINT(zeff + 0.001d0), &
     344             :                                       symbol=element_symbol, &
     345             :                                       mass=REAL(mass/massunit, KIND=C_DOUBLE), &
     346           2 :                                       spin_orbit=.FALSE.)
     347             : !
     348         720 :             ALLOCATE (gth_atompot)
     349           2 :             CALL gth_potential_conversion(gth_potential, gth_atompot)
     350             : ! set radial grid
     351        1860 :             fun(1:nmesh) = rp(1:nmesh)
     352           2 :             CALL sirius_set_atom_type_radial_grid(sctx, label, nmesh, fun(1))
     353             : ! set beta-projectors
     354           8 :             ALLOCATE (ef(nmesh), beta(nmesh))
     355           2 :             ibeta = 0
     356          10 :             DO l = 0, 3
     357           8 :                IF (gth_atompot%nl(l) == 0) CYCLE
     358           2 :                rl = gth_atompot%rcnl(l)
     359             : ! we need to multiply by r so that data transferred to sirius are r \beta(r) not beta(r)
     360        1860 :                ef(1:nmesh) = EXP(-0.5_dp*rp(1:nmesh)*rp(1:nmesh)/(rl*rl))
     361           6 :                DO i = 1, gth_atompot%nl(l)
     362           2 :                   pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
     363           2 :                   j = l + 2*i - 1
     364           2 :                   pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
     365        1860 :                   beta(:) = pf*rp**(l + 2*i - 2)*ef
     366           2 :                   ibeta = ibeta + 1
     367        1860 :                   fun(1:nmesh) = beta(1:nmesh)*rp(1:nmesh)
     368             :                   CALL sirius_add_atom_type_radial_function(sctx, label, &
     369          10 :                                                             "beta", fun(1), nmesh, l=l)
     370             :                END DO
     371             :             END DO
     372           2 :             DEALLOCATE (ef, beta)
     373           2 :             nbeta = ibeta
     374             : 
     375             : ! nonlocal PP matrix elements
     376           8 :             ALLOCATE (dion(nbeta, nbeta))
     377           6 :             dion = 0.0_dp
     378          10 :             DO l = 0, 3
     379           8 :                IF (gth_atompot%nl(l) == 0) CYCLE
     380           2 :                ibeta = SUM(gth_atompot%nl(0:l - 1)) + 1
     381           2 :                i = ibeta + gth_atompot%nl(l) - 1
     382           8 :                dion(ibeta:i, ibeta:i) = gth_atompot%hnl(1:gth_atompot%nl(l), 1:gth_atompot%nl(l), l)
     383             :             END DO
     384           2 :             CALL sirius_set_atom_type_dion(sctx, label, nbeta, dion(1, 1))
     385           2 :             DEALLOCATE (dion)
     386             : 
     387             : ! set non-linear core correction
     388           2 :             IF (gth_atompot%nlcc) THEN
     389           0 :                ALLOCATE (corden(nmesh), fe(nmesh), rc(nmesh))
     390           0 :                corden(:) = 0.0_dp
     391           0 :                n = gth_atompot%nexp_nlcc
     392           0 :                DO i = 1, n
     393           0 :                   al = gth_atompot%alpha_nlcc(i)
     394           0 :                   rc(:) = rp(:)/al
     395           0 :                   fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
     396           0 :                   DO j = 1, gth_atompot%nct_nlcc(i)
     397           0 :                      cval = gth_atompot%cval_nlcc(j, i)
     398           0 :                      corden(:) = corden(:) + fe(:)*rc(:)**(2*j - 2)*cval
     399             :                   END DO
     400             :                END DO
     401           0 :                fun(1:nmesh) = corden(1:nmesh)*rp(1:nmesh)
     402             :                CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_core", &
     403           0 :                                                          fun(1), nmesh)
     404           0 :                DEALLOCATE (corden, fe, rc)
     405             :             END IF
     406             : 
     407             : ! local potential
     408           6 :             ALLOCATE (locpot(nmesh))
     409        1860 :             locpot(:) = 0.0_dp
     410           2 :             CALL atom_local_potential(locpot, gth_atompot, rp)
     411        1860 :             fun(1:nmesh) = locpot(1:nmesh)
     412             :             CALL sirius_add_atom_type_radial_function(sctx, label, "vloc", &
     413           2 :                                                       fun(1), nmesh)
     414           2 :             DEALLOCATE (locpot)
     415             : !
     416           2 :             NULLIFY (density, wavefunction, wfninfo)
     417             :             CALL calculate_atomic_orbitals(atomic_kind_set(ikind), qs_kind_set(ikind), &
     418             :                                            density=density, wavefunction=wavefunction, &
     419           2 :                                            wfninfo=wfninfo, agrid=atom_grid)
     420             : 
     421             : ! set the atomic radial functions
     422           6 :             DO iwf = 1, SIZE(wavefunction, 2)
     423           4 :                focc = wfninfo(1, iwf)
     424           4 :                l = NINT(wfninfo(2, iwf))
     425             : ! we can not easily get the principal quantum number
     426           4 :                nu = -1
     427           4 :                IF (up) THEN
     428           0 :                   fun(1:nmesh) = wavefunction(1:nmesh, iwf)*rp(i)
     429             :                ELSE
     430        3720 :                   DO i = 1, nmesh
     431        3720 :                      fun(i) = wavefunction(nmesh - i + 1, iwf)*rp(i)
     432             :                   END DO
     433             :                END IF
     434             :                CALL sirius_add_atom_type_radial_function(sctx, &
     435             :                                                          label, "ps_atomic_wf", &
     436           6 :                                                          fun(1), nmesh, l=l, occ=REAL(focc, KIND=C_DOUBLE), n=nu)
     437             :             END DO
     438             : 
     439             : ! set total charge density of a free atom (to compute initial rho(r))
     440           2 :             IF (up) THEN
     441           0 :                fun(1:nmesh) = fourpi*density(1:nmesh)*atom_grid%rad(1:nmesh)**2
     442             :             ELSE
     443        1860 :                DO i = 1, nmesh
     444        1860 :                   fun(i) = fourpi*density(nmesh - i + 1)*atom_grid%rad(nmesh - i + 1)**2
     445             :                END DO
     446             :             END IF
     447             :             CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_total", &
     448           2 :                                                       fun(1), nmesh)
     449             : 
     450           2 :             IF (ASSOCIATED(density)) DEALLOCATE (density)
     451           2 :             IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
     452           2 :             IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
     453             : 
     454           2 :             CALL deallocate_grid_atom(atom_grid)
     455           2 :             DEALLOCATE (rp, fun)
     456           2 :             DEALLOCATE (gth_atompot)
     457             : !
     458             :          ELSE
     459             :             CALL cp_abort(__LOCATION__, &
     460           0 :                           "CP2K/SIRIUS: atomic kind needs UPF or GTH potential definition")
     461             :          END IF
     462             : 
     463             :          CALL get_qs_kind(qs_kind_set(ikind), &
     464             :                           dft_plus_u_atom=dft_plus_u_atom, &
     465             :                           l_of_dft_plus_u=lu, &
     466             :                           n_of_dft_plus_u=nu, &
     467             :                           u_minus_j_target=u_minus_j, &
     468             :                           U_of_dft_plus_u=U_u, &
     469             :                           J_of_dft_plus_u=J_u, &
     470             :                           alpha_of_dft_plus_u=alpha_u, &
     471             :                           beta_of_dft_plus_u=beta_u, &
     472             :                           J0_of_dft_plus_u=J0_u, &
     473          22 :                           occupation_of_dft_plus_u=occ_u)
     474             : 
     475          60 :          IF (dft_plus_u_atom) THEN
     476           0 :             IF (nu < 1) THEN
     477           0 :                CPABORT("CP2K/SIRIUS (hubbard): principal quantum number not specified")
     478             :             END IF
     479             : 
     480           0 :             IF (lu < 0) THEN
     481           0 :                CPABORT("CP2K/SIRIUS (hubbard): l can not be negative.")
     482             :             END IF
     483             : 
     484           0 :             IF (occ_u < 0.0) THEN
     485           0 :                CPABORT("CP2K/SIRIUS (hubbard): the occupation number can not be negative.")
     486             :             END IF
     487             : 
     488           0 :             j_t(:) = 0.0
     489           0 :             IF (ABS(u_minus_j) < 1e-8) THEN
     490           0 :                j_t(1) = J_u
     491             :                CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
     492           0 :                                                  occ_u, U_u, j_t, alpha_u, beta_u, J0_u)
     493             :             ELSE
     494             :                CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
     495           0 :                                                  occ_u, u_minus_j, j_t, alpha_u, beta_u, J0_u)
     496             :             END IF
     497             :          END IF
     498             : 
     499             :       END DO
     500             : 
     501             : ! add atoms to the unit cell
     502             : ! WARNING: sirius accepts only fractional coordinates;
     503          16 :       natom = SIZE(particle_set)
     504          50 :       DO iatom = 1, natom
     505         136 :          vr(1:3) = particle_set(iatom)%r(1:3)
     506          34 :          CALL real_to_scaled(vs, vr, my_cell)
     507          34 :          atomic_kind => particle_set(iatom)%atomic_kind
     508          34 :          ikind = atomic_kind%kind_number
     509          34 :          CALL get_atomic_kind(atomic_kind, name=label)
     510          34 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, magnetization=magnetization)
     511             : ! angle of magnetization might come from input Atom x y z mx my mz
     512             : ! or as an angle?
     513             : ! Answer : SIRIUS only accept the magnetization as mx, my, mz
     514          34 :          IF (num_mag_dims .EQ. 3) THEN
     515           2 :             angle1 = 0.0_dp
     516           2 :             angle2 = 0.0_dp
     517           2 :             v1(1) = magnetization*SIN(angle1)*COS(angle2)
     518           2 :             v1(2) = magnetization*SIN(angle1)*SIN(angle2)
     519           2 :             v1(3) = magnetization*COS(angle1)
     520             :          ELSE
     521          32 :             v1 = 0._dp
     522          32 :             v1(3) = magnetization
     523             :          END IF
     524          34 :          v2(1:3) = vs(1:3)
     525          84 :          CALL sirius_add_atom(sctx, label, v2(1), v1(1))
     526             :       END DO
     527             : 
     528          16 :       CALL sirius_set_mpi_grid_dims(sctx, 2, mpi_grid_dims)
     529             : 
     530             : ! initialize global variables/indices/arrays/etc. of the simulation
     531          16 :       CALL sirius_initialize_context(sctx)
     532             : 
     533             :       ! strictly speaking the parameter use_symmetry is initialized at the
     534             :       ! beginning but it does no harm to do it that way
     535          16 :       IF (use_symmetry) THEN
     536          14 :          CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.TRUE., kset_handler=ks_handler)
     537             :       ELSE
     538           2 :          CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.FALSE., kset_handler=ks_handler)
     539             :       END IF
     540             : ! create ground-state class
     541          16 :       CALL sirius_create_ground_state(ks_handler, gs_handler)
     542             : 
     543          16 :       CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler, ks_handler=ks_handler)
     544             : #endif
     545          32 :    END SUBROUTINE cp_sirius_create_env
     546             : 
     547             : !***************************************************************************************************
     548             : !> \brief ...
     549             : !> \param pwdft_env ...
     550             : !> \param
     551             : !> \par History
     552             : !>      07.2018 Update the Sirius environment
     553             : !> \author JHU
     554             : ! **************************************************************************************************
     555          16 :    SUBROUTINE cp_sirius_update_context(pwdft_env)
     556             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     557             : 
     558             :       INTEGER                                            :: iatom, natom
     559             :       REAL(KIND=C_DOUBLE), DIMENSION(3)                  :: a1, a2, a3, v2
     560             :       REAL(KIND=dp), DIMENSION(3)                        :: vr, vs
     561             :       TYPE(cell_type), POINTER                           :: my_cell
     562          16 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     563             :       TYPE(qs_subsys_type), POINTER                      :: qs_subsys
     564             :       TYPE(sirius_context_handler)                       :: sctx
     565             :       TYPE(sirius_ground_state_handler)                  :: gs_handler
     566             : 
     567           0 :       CPASSERT(ASSOCIATED(pwdft_env))
     568          16 :       CALL pwdft_env_get(pwdft_env, sctx=sctx, gs_handler=gs_handler)
     569             : 
     570             : ! get current positions and lattice vectors
     571          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
     572             : 
     573             : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
     574          16 :       CALL qs_subsys_get(qs_subsys, cell=my_cell)
     575          64 :       a1(:) = my_cell%hmat(:, 1)
     576          64 :       a2(:) = my_cell%hmat(:, 2)
     577          64 :       a3(:) = my_cell%hmat(:, 3)
     578          16 :       CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
     579             : 
     580             : ! new atomic positions
     581          16 :       CALL qs_subsys_get(qs_subsys, particle_set=particle_set)
     582          16 :       natom = SIZE(particle_set)
     583          50 :       DO iatom = 1, natom
     584         136 :          vr(1:3) = particle_set(iatom)%r(1:3)
     585          34 :          CALL real_to_scaled(vs, vr, my_cell)
     586          34 :          v2(1:3) = vs(1:3)
     587          50 :          CALL sirius_set_atom_position(sctx, iatom, v2(1))
     588             :       END DO
     589             : 
     590             : ! update ground-state class
     591          16 :       CALL sirius_update_ground_state(gs_handler)
     592             : 
     593          16 :       CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler)
     594             : 
     595          16 :    END SUBROUTINE cp_sirius_update_context
     596             : 
     597             : ! **************************************************************************************************
     598             : !> \brief ...
     599             : !> \param sctx ...
     600             : !> \param section ...
     601             : !> \param section_name ...
     602             : ! **************************************************************************************************
     603          80 :    SUBROUTINE cp_sirius_fill_in_section(sctx, section, section_name)
     604             :       TYPE(sirius_context_handler), INTENT(INOUT)        :: sctx
     605             :       TYPE(section_vals_type), POINTER                   :: section
     606             :       CHARACTER(*), INTENT(in)                           :: section_name
     607             : 
     608             :       CHARACTER(len=256), TARGET                         :: option_name
     609             :       CHARACTER(len=4096)                                :: description, usage
     610          80 :       CHARACTER(len=80), DIMENSION(:), POINTER           :: tmp
     611             :       CHARACTER(len=80), TARGET                          :: str
     612             :       INTEGER                                            :: ctype, elem, ic, j
     613          80 :       INTEGER, DIMENSION(:), POINTER                     :: ivals
     614             :       INTEGER, TARGET                                    :: enum_length, ival, length, &
     615             :                                                             num_possible_values, number_of_options
     616             :       LOGICAL                                            :: explicit
     617          80 :       LOGICAL, DIMENSION(:), POINTER                     :: lvals
     618             :       LOGICAL, TARGET                                    :: found, lval
     619          80 :       REAL(kind=dp), DIMENSION(:), POINTER               :: rvals
     620             :       REAL(kind=dp), TARGET                              :: rval
     621             : 
     622          80 :       NULLIFY (rvals)
     623          80 :       NULLIFY (ivals)
     624          80 :       CALL sirius_option_get_section_length(section_name, number_of_options)
     625             : 
     626        1712 :       DO elem = 1, number_of_options
     627        1632 :          option_name = ''
     628             :          CALL sirius_option_get_info(section_name, &
     629             :                                      elem, &
     630             :                                      option_name, &
     631             :                                      256, &
     632             :                                      ctype, &
     633             :                                      num_possible_values, &
     634             :                                      enum_length, &
     635             :                                      description, &
     636             :                                      4096, &
     637             :                                      usage, &
     638        1632 :                                      4096)
     639        1712 :          IF ((option_name /= 'memory_usage') .AND. (option_name /= 'xc_functionals') .AND. (option_name /= 'vk')) THEN
     640        1600 :             CALL section_vals_val_get(section, option_name, explicit=found)
     641        1600 :             IF (found) THEN
     642         128 :                SELECT CASE (ctype)
     643             :                CASE (SIRIUS_INTEGER_TYPE)
     644         128 :                   CALL section_vals_val_get(section, option_name, i_val=ival)
     645         128 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ival))
     646             :                CASE (SIRIUS_NUMBER_TYPE)
     647         104 :                   CALL section_vals_val_get(section, option_name, r_val=rval)
     648         104 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rval))
     649             :                CASE (SIRIUS_LOGICAL_TYPE)
     650          24 :                   CALL section_vals_val_get(section, option_name, l_val=lval)
     651          24 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lval))
     652             :                CASE (SIRIUS_STRING_TYPE)      ! string nightmare
     653         100 :                   str = ''
     654         100 :                   CALL section_vals_val_get(section, option_name, explicit=explicit, c_val=str)
     655         100 :                   str = TRIM(ADJUSTL(str))
     656        8100 :                   DO j = 1, LEN(str)
     657        8000 :                      ic = ICHAR(str(j:j))
     658        8100 :                      IF (ic >= 65 .AND. ic < 90) str(j:j) = CHAR(ic + 32)
     659             :                   END DO
     660             : 
     661         100 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), max_length=LEN_TRIM(str))
     662             :                CASE (SIRIUS_INTEGER_ARRAY_TYPE)
     663          16 :                   CALL section_vals_val_get(section, option_name, i_vals=ivals)
     664             :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ivals(1)), &
     665          16 :                                          max_length=num_possible_values)
     666             :                CASE (SIRIUS_NUMBER_ARRAY_TYPE)
     667           0 :                   CALL section_vals_val_get(section, option_name, r_vals=rvals)
     668             :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rvals(1)), &
     669           0 :                                          max_length=num_possible_values)
     670             :                CASE (SIRIUS_LOGICAL_ARRAY_TYPE)
     671           0 :                   CALL section_vals_val_get(section, option_name, l_vals=lvals)
     672             :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lvals(1)), &
     673           0 :                                          max_length=num_possible_values)
     674             :                CASE (SIRIUS_STRING_ARRAY_TYPE)
     675           0 :                   CALL section_vals_val_get(section, option_name, explicit=explicit, n_rep_val=length)
     676         372 :                   DO j = 1, length
     677           0 :                      str = ''
     678           0 :                      CALL section_vals_val_get(section, option_name, i_rep_val=j, explicit=explicit, c_vals=tmp)
     679           0 :                      str = TRIM(ADJUSTL(tmp(j)))
     680             :                      CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), &
     681           0 :                                             max_length=LEN_TRIM(str), append=.TRUE.)
     682             :                   END DO
     683             :                CASE DEFAULT
     684             :                END SELECT
     685             :             END IF
     686             :          END IF
     687             :       END DO
     688          80 :    END SUBROUTINE cp_sirius_fill_in_section
     689             : 
     690             : !***************************************************************************************************
     691             : !> \brief ...
     692             : !> \param pwdft_env ...
     693             : !> \param calculate_forces ...
     694             : !> \param calculate_stress_tensor ...
     695             : !> \param
     696             : !> \par History
     697             : !>      07.2018 start the Sirius library
     698             : !> \author JHU
     699             : ! **************************************************************************************************
     700          16 :    SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_tensor)
     701             :       TYPE(pwdft_environment_type), INTENT(INOUT), &
     702             :          POINTER                                         :: pwdft_env
     703             :       LOGICAL, INTENT(IN)                                :: calculate_forces, calculate_stress_tensor
     704             : 
     705             :       INTEGER                                            :: iw, n1, n2
     706             :       LOGICAL                                            :: do_print, gs_converged
     707             :       REAL(KIND=C_DOUBLE)                                :: etotal
     708          16 :       REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :)  :: cforces
     709             :       REAL(KIND=C_DOUBLE), DIMENSION(3, 3)               :: cstress
     710             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: stress
     711          16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: forces
     712             :       TYPE(cp_logger_type), POINTER                      :: logger
     713             :       TYPE(pwdft_energy_type), POINTER                   :: energy
     714             :       TYPE(section_vals_type), POINTER                   :: print_section, pwdft_input
     715             :       TYPE(sirius_ground_state_handler)                  :: gs_handler
     716             : 
     717           0 :       CPASSERT(ASSOCIATED(pwdft_env))
     718             : 
     719          16 :       NULLIFY (logger)
     720          16 :       logger => cp_get_default_logger()
     721          16 :       iw = cp_logger_get_default_io_unit(logger)
     722             : 
     723          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, gs_handler=gs_handler)
     724          16 :       CALL sirius_find_ground_state(gs_handler, converged=gs_converged)
     725             : 
     726          16 :       IF (gs_converged) THEN
     727          16 :          IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS: ground state is converged"
     728             :       ELSE
     729           0 :          IF (pwdft_env%ignore_convergence_failure) THEN
     730           0 :             IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS Warning: ground state is not converged"
     731             :          ELSE
     732           0 :             CPABORT("CP2K/SIRIUS (ground state): SIRIUS did not converge.")
     733             :          END IF
     734             :       END IF
     735          16 :       IF (iw > 0) CALL m_flush(iw)
     736             : 
     737          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, energy=energy)
     738             :       etotal = 0.0_C_DOUBLE
     739             : 
     740          16 :       CALL sirius_get_energy(gs_handler, 'band-gap', etotal)
     741          16 :       energy%band_gap = etotal
     742             : 
     743             :       etotal = 0.0_C_DOUBLE
     744          16 :       CALL sirius_get_energy(gs_handler, 'total', etotal)
     745          16 :       energy%etotal = etotal
     746             : 
     747             :       ! extract entropy (TS returned by sirius is always negative, sign
     748             :       ! convention in QE)
     749             :       etotal = 0.0_C_DOUBLE
     750          16 :       CALL sirius_get_energy(gs_handler, 'demet', etotal)
     751          16 :       energy%entropy = -etotal
     752             : 
     753          16 :       IF (calculate_forces) THEN
     754          12 :          CALL pwdft_env_get(pwdft_env=pwdft_env, forces=forces)
     755          12 :          n1 = SIZE(forces, 1)
     756          12 :          n2 = SIZE(forces, 2)
     757             : 
     758          48 :          ALLOCATE (cforces(n2, n1))
     759         132 :          cforces = 0.0_C_DOUBLE
     760          12 :          CALL sirius_get_forces(gs_handler, 'total', cforces)
     761             :          ! Sirius computes the forces but cp2k use the gradient everywhere
     762             :          ! so a minus sign is needed.
     763             :          ! note also that sirius and cp2k store the forces transpose to each other
     764             :          ! sirius : forces(coordinates, atoms)
     765             :          ! cp2k : forces(atoms, coordinates)
     766         138 :          forces = -TRANSPOSE(cforces(:, :))
     767          12 :          DEALLOCATE (cforces)
     768             :       END IF
     769             : 
     770          16 :       IF (calculate_stress_tensor) THEN
     771           0 :          cstress = 0.0_C_DOUBLE
     772           0 :          CALL sirius_get_stress_tensor(gs_handler, 'total', cstress)
     773           0 :          stress(1:3, 1:3) = cstress(1:3, 1:3)
     774           0 :          CALL pwdft_env_set(pwdft_env=pwdft_env, stress=stress)
     775             :       END IF
     776             : 
     777          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_input)
     778          16 :       print_section => section_vals_get_subs_vals(pwdft_input, "PRINT")
     779          16 :       CALL section_vals_get(print_section, explicit=do_print)
     780          16 :       IF (do_print) THEN
     781           2 :          CALL cp_sirius_print_results(pwdft_env, print_section)
     782             :       END IF
     783          32 :    END SUBROUTINE cp_sirius_energy_force
     784             : 
     785             : !***************************************************************************************************
     786             : !> \brief ...
     787             : !> \param pwdft_env ...
     788             : !> \param print_section ...
     789             : !> \param
     790             : !> \par History
     791             : !>      12.2019 init
     792             : !> \author JHU
     793             : ! **************************************************************************************************
     794           2 :    SUBROUTINE cp_sirius_print_results(pwdft_env, print_section)
     795             :       TYPE(pwdft_environment_type), INTENT(INOUT), &
     796             :          POINTER                                         :: pwdft_env
     797             :       TYPE(section_vals_type), POINTER                   :: print_section
     798             : 
     799             :       CHARACTER(LEN=default_string_length)               :: my_act, my_pos
     800             :       INTEGER                                            :: i, ik, iounit, ispn, iterstep, iv, iw, &
     801             :                                                             nbands, nhist, nkpts, nspins
     802             :       INTEGER(KIND=C_INT)                                :: cint
     803             :       LOGICAL                                            :: append, dos, ionode
     804             :       REAL(KIND=C_DOUBLE)                                :: creal
     805           2 :       REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:)     :: slist
     806             :       REAL(KIND=dp)                                      :: de, e_fermi(2), emax, emin, eval
     807           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: wkpt
     808           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ehist, hist, occval
     809           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: energies, occupations
     810             :       TYPE(cp_logger_type), POINTER                      :: logger
     811             :       TYPE(sirius_context_handler)                       :: sctx
     812             :       TYPE(sirius_ground_state_handler)                  :: gs_handler
     813             :       TYPE(sirius_kpoint_set_handler)                    :: ks_handler
     814             : 
     815           2 :       NULLIFY (logger)
     816           4 :       logger => cp_get_default_logger()
     817             :       ionode = logger%para_env%is_source()
     818           2 :       iounit = cp_logger_get_default_io_unit(logger)
     819             : 
     820             :       ! Density of States
     821           2 :       dos = BTEST(cp_print_key_should_output(logger%iter_info, print_section, "DOS"), cp_p_file)
     822           2 :       IF (dos) THEN
     823           2 :          CALL pwdft_env_get(pwdft_env, ks_handler=ks_handler)
     824           2 :          CALL pwdft_env_get(pwdft_env, gs_handler=gs_handler)
     825           2 :          CALL pwdft_env_get(pwdft_env, sctx=sctx)
     826             : 
     827           2 :          CALL section_vals_val_get(print_section, "DOS%DELTA_E", r_val=de)
     828           2 :          CALL section_vals_val_get(print_section, "DOS%APPEND", l_val=append)
     829             : 
     830           2 :          CALL sirius_get_num_kpoints(ks_handler, cint)
     831           2 :          nkpts = cint
     832           2 :          CALL sirius_get_parameters(sctx, num_bands=cint)
     833           2 :          nbands = cint
     834           2 :          CALL sirius_get_parameters(sctx, num_spins=cint)
     835           2 :          nspins = cint
     836           2 :          e_fermi(:) = 0.0_dp
     837          10 :          ALLOCATE (energies(nbands, nspins, nkpts))
     838         442 :          energies = 0.0_dp
     839           8 :          ALLOCATE (occupations(nbands, nspins, nkpts))
     840         442 :          occupations = 0.0_dp
     841           6 :          ALLOCATE (wkpt(nkpts))
     842           6 :          ALLOCATE (slist(nbands))
     843          10 :          DO ik = 1, nkpts
     844           8 :             CALL sirius_get_kpoint_properties(ks_handler, ik, creal)
     845          10 :             wkpt(ik) = creal
     846             :          END DO
     847          10 :          DO ik = 1, nkpts
     848          26 :             DO ispn = 1, nspins
     849          16 :                CALL sirius_get_band_energies(ks_handler, ik, ispn, slist)
     850         432 :                energies(1:nbands, ispn, ik) = slist(1:nbands)
     851          16 :                CALL sirius_get_band_occupancies(ks_handler, ik, ispn, slist)
     852         440 :                occupations(1:nbands, ispn, ik) = slist(1:nbands)
     853             :             END DO
     854             :          END DO
     855         442 :          emin = MINVAL(energies)
     856         442 :          emax = MAXVAL(energies)
     857           2 :          nhist = NINT((emax - emin)/de) + 1
     858          16 :          ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
     859       16110 :          hist = 0.0_dp
     860       16110 :          occval = 0.0_dp
     861       16110 :          ehist = 0.0_dp
     862             : 
     863          10 :          DO ik = 1, nkpts
     864          26 :             DO ispn = 1, nspins
     865         440 :                DO i = 1, nbands
     866         416 :                   eval = energies(i, ispn, ik) - emin
     867         416 :                   iv = NINT(eval/de) + 1
     868         416 :                   CPASSERT((iv > 0) .AND. (iv <= nhist))
     869         416 :                   hist(iv, ispn) = hist(iv, ispn) + wkpt(ik)
     870         432 :                   occval(iv, ispn) = occval(iv, ispn) + wkpt(ik)*occupations(i, ispn, ik)
     871             :                END DO
     872             :             END DO
     873             :          END DO
     874       16110 :          hist = hist/REAL(nbands, KIND=dp)
     875        8054 :          DO i = 1, nhist
     876       24158 :             ehist(i, 1:nspins) = emin + (i - 1)*de
     877             :          END DO
     878             : 
     879           2 :          iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
     880           2 :          my_act = "WRITE"
     881           2 :          IF (append .AND. iterstep > 1) THEN
     882           0 :             my_pos = "APPEND"
     883             :          ELSE
     884           2 :             my_pos = "REWIND"
     885             :          END IF
     886             : 
     887             :          iw = cp_print_key_unit_nr(logger, print_section, "DOS", &
     888             :                                    extension=".dos", file_position=my_pos, file_action=my_act, &
     889           2 :                                    file_form="FORMATTED")
     890           2 :          IF (iw > 0) THEN
     891           1 :             IF (nspins == 2) THEN
     892             :                WRITE (UNIT=iw, FMT="(T2,A,I0,A,2F12.6)") &
     893           1 :                   "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
     894           1 :                WRITE (UNIT=iw, FMT="(T2,A, A)") "   Energy[a.u.]  Alpha_Density     Occupation", &
     895           2 :                   "   Beta_Density      Occupation"
     896             :             ELSE
     897             :                WRITE (UNIT=iw, FMT="(T2,A,I0,A,F12.6)") &
     898           0 :                   "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
     899           0 :                WRITE (UNIT=iw, FMT="(T2,A)") "   Energy[a.u.]       Density     Occupation"
     900             :             END IF
     901        4027 :             DO i = 1, nhist
     902        4026 :                eval = emin + (i - 1)*de
     903        4027 :                IF (nspins == 2) THEN
     904        4026 :                   WRITE (UNIT=iw, FMT="(F15.8,4F15.4)") eval, hist(i, 1), occval(i, 1), &
     905        8052 :                      hist(i, 2), occval(i, 2)
     906             :                ELSE
     907           0 :                   WRITE (UNIT=iw, FMT="(F15.8,2F15.4)") eval, hist(i, 1), occval(i, 1)
     908             :                END IF
     909             :             END DO
     910             :          END IF
     911           2 :          CALL cp_print_key_finished_output(iw, logger, print_section, "DOS")
     912             : 
     913           2 :          DEALLOCATE (energies, occupations, wkpt, slist)
     914           4 :          DEALLOCATE (hist, occval, ehist)
     915             :       END IF
     916           4 :    END SUBROUTINE cp_sirius_print_results
     917             : 
     918             : END MODULE sirius_interface
     919             : 
     920             : #else
     921             : 
     922             : !***************************************************************************************************
     923             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     924             : !***************************************************************************************************
     925             : MODULE sirius_interface
     926             :    USE pwdft_environment_types, ONLY: pwdft_environment_type
     927             : #include "./base/base_uses.f90"
     928             : 
     929             :    IMPLICIT NONE
     930             :    PRIVATE
     931             : 
     932             :    PUBLIC :: cp_sirius_init, cp_sirius_finalize
     933             :    PUBLIC :: cp_sirius_create_env, cp_sirius_energy_force, cp_sirius_update_context
     934             : 
     935             : CONTAINS
     936             : 
     937             : ! **************************************************************************************************
     938             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     939             : ! **************************************************************************************************
     940             :    SUBROUTINE cp_sirius_init()
     941             :    END SUBROUTINE cp_sirius_init
     942             : 
     943             : ! **************************************************************************************************
     944             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     945             : ! **************************************************************************************************
     946             :    SUBROUTINE cp_sirius_finalize()
     947             :    END SUBROUTINE cp_sirius_finalize
     948             : 
     949             : ! **************************************************************************************************
     950             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     951             : !> \param pwdft_env ...
     952             : ! **************************************************************************************************
     953             :    SUBROUTINE cp_sirius_create_env(pwdft_env)
     954             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     955             : 
     956             :       MARK_USED(pwdft_env)
     957             :       CPABORT("Sirius library is missing")
     958             :    END SUBROUTINE cp_sirius_create_env
     959             : 
     960             : ! **************************************************************************************************
     961             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     962             : !> \param pwdft_env ...
     963             : !> \param calculate_forces ...
     964             : !> \param calculate_stress ...
     965             : ! **************************************************************************************************
     966             :    SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
     967             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     968             :       LOGICAL                                            :: calculate_forces, calculate_stress
     969             : 
     970             :       MARK_USED(pwdft_env)
     971             :       MARK_USED(calculate_forces)
     972             :       MARK_USED(calculate_stress)
     973             :       CPABORT("Sirius library is missing")
     974             :    END SUBROUTINE cp_sirius_energy_force
     975             : 
     976             : ! **************************************************************************************************
     977             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     978             : !> \param pwdft_env ...
     979             : ! **************************************************************************************************
     980             :    SUBROUTINE cp_sirius_update_context(pwdft_env)
     981             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     982             : 
     983             :       MARK_USED(pwdft_env)
     984             :       CPABORT("Sirius library is missing")
     985             :    END SUBROUTINE cp_sirius_update_context
     986             : 
     987             : END MODULE sirius_interface
     988             : 
     989             : #endif

Generated by: LCOV version 1.15