LCOV - code coverage report
Current view: top level - src - sirius_interface.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:4dc10b3) Lines: 330 383 86.2 %
Date: 2024-11-21 06:45:46 Functions: 7 7 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 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        8530 :    SUBROUTINE cp_sirius_init()
     105        8530 :       CALL sirius_initialize(.FALSE.)
     106        8530 :    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        8530 :    SUBROUTINE cp_sirius_finalize()
     116        8530 :       CALL sirius_finalize(.FALSE., .FALSE., .FALSE.)
     117        8530 :    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
     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             :       !
     258             :       ! uncomment these lines when nlcg is officially supported
     259             :       !
     260             : 
     261             :       !     import nlcg section
     262             :       !      pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "nlcg")
     263             :       !      IF (ASSOCIATED(pwdft_sub_section)) THEN
     264             :       !         CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "nlcg")
     265             :       !      ENDIF
     266             : 
     267             :       !CALL sirius_dump_runtime_setup(sctx, "runtime.json")
     268          16 :       CALL sirius_import_parameters(sctx, '{}')
     269             : 
     270             : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
     271          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
     272          16 :       CALL qs_subsys_get(qs_subsys, cell=my_cell, use_ref_cell=use_ref_cell)
     273          64 :       a1(:) = my_cell%hmat(:, 1)
     274          64 :       a2(:) = my_cell%hmat(:, 2)
     275          64 :       a3(:) = my_cell%hmat(:, 3)
     276          16 :       CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
     277             : 
     278          16 :       IF (use_ref_cell) THEN
     279           0 :          CPWARN("SIRIUS| The specified CELL_REF will be ignored for PW_DFT calculations")
     280             :       END IF
     281             : 
     282             : ! set up the atomic type definitions
     283             :       CALL qs_subsys_get(qs_subsys, &
     284             :                          atomic_kind_set=atomic_kind_set, &
     285             :                          qs_kind_set=qs_kind_set, &
     286          16 :                          particle_set=particle_set)
     287          16 :       nkind = SIZE(atomic_kind_set)
     288          38 :       DO ikind = 1, nkind
     289             :          CALL get_atomic_kind(atomic_kind_set(ikind), &
     290          22 :                               name=label, element_symbol=element_symbol, mass=mass)
     291          22 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
     292          22 :          NULLIFY (upf_pot, gth_potential)
     293          22 :          CALL get_qs_kind(qs_kind_set(ikind), upf_potential=upf_pot, gth_potential=gth_potential)
     294             : 
     295          22 :          IF (ASSOCIATED(upf_pot)) THEN
     296             :             CALL sirius_add_atom_type(sctx, label, fname=upf_pot%filename, &
     297             :                                       symbol=element_symbol, &
     298          20 :                                       mass=REAL(mass/massunit, KIND=C_DOUBLE))
     299             : 
     300           2 :          ELSEIF (ASSOCIATED(gth_potential)) THEN
     301             : !
     302           2 :             NULLIFY (atom_grid)
     303           2 :             CALL allocate_grid_atom(atom_grid)
     304           2 :             nmesh = 929
     305           2 :             atom_grid%nr = nmesh
     306           2 :             CALL create_grid_atom(atom_grid, nmesh, 1, 1, 0, do_gapw_log)
     307           8 :             ALLOCATE (rp(nmesh), fun(nmesh))
     308           2 :             IF (atom_grid%rad(1) < atom_grid%rad(nmesh)) THEN
     309             :                up = .TRUE.
     310             :             ELSE
     311             :                up = .FALSE.
     312             :             END IF
     313             :             IF (up) THEN
     314           0 :                rp(1:nmesh) = atom_grid%rad(1:nmesh)
     315             :             ELSE
     316        1860 :                DO i = 1, nmesh
     317        1860 :                   rp(i) = atom_grid%rad(nmesh - i + 1)
     318             :                END DO
     319             :             END IF
     320             : ! add new atom type
     321             :             CALL sirius_add_atom_type(sctx, label, &
     322             :                                       zn=NINT(zeff + 0.001d0), &
     323             :                                       symbol=element_symbol, &
     324             :                                       mass=REAL(mass/massunit, KIND=C_DOUBLE), &
     325           2 :                                       spin_orbit=.FALSE.)
     326             : !
     327         720 :             ALLOCATE (gth_atompot)
     328           2 :             CALL gth_potential_conversion(gth_potential, gth_atompot)
     329             : ! set radial grid
     330        1860 :             fun(1:nmesh) = rp(1:nmesh)
     331           2 :             CALL sirius_set_atom_type_radial_grid(sctx, label, nmesh, fun(1))
     332             : ! set beta-projectors
     333           8 :             ALLOCATE (ef(nmesh), beta(nmesh))
     334           2 :             ibeta = 0
     335          10 :             DO l = 0, 3
     336           8 :                IF (gth_atompot%nl(l) == 0) CYCLE
     337           2 :                rl = gth_atompot%rcnl(l)
     338             : ! we need to multiply by r so that data transferred to sirius are r \beta(r) not beta(r)
     339        1860 :                ef(1:nmesh) = EXP(-0.5_dp*rp(1:nmesh)*rp(1:nmesh)/(rl*rl))
     340           6 :                DO i = 1, gth_atompot%nl(l)
     341           2 :                   pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
     342           2 :                   j = l + 2*i - 1
     343           2 :                   pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
     344        1860 :                   beta(:) = pf*rp**(l + 2*i - 2)*ef
     345           2 :                   ibeta = ibeta + 1
     346        1860 :                   fun(1:nmesh) = beta(1:nmesh)*rp(1:nmesh)
     347             :                   CALL sirius_add_atom_type_radial_function(sctx, label, &
     348          10 :                                                             "beta", fun(1), nmesh, l=l)
     349             :                END DO
     350             :             END DO
     351           2 :             DEALLOCATE (ef, beta)
     352           2 :             nbeta = ibeta
     353             : 
     354             : ! nonlocal PP matrix elements
     355           8 :             ALLOCATE (dion(nbeta, nbeta))
     356           6 :             dion = 0.0_dp
     357          10 :             DO l = 0, 3
     358           8 :                IF (gth_atompot%nl(l) == 0) CYCLE
     359           2 :                ibeta = SUM(gth_atompot%nl(0:l - 1)) + 1
     360           2 :                i = ibeta + gth_atompot%nl(l) - 1
     361           8 :                dion(ibeta:i, ibeta:i) = gth_atompot%hnl(1:gth_atompot%nl(l), 1:gth_atompot%nl(l), l)
     362             :             END DO
     363           2 :             CALL sirius_set_atom_type_dion(sctx, label, nbeta, dion(1, 1))
     364           2 :             DEALLOCATE (dion)
     365             : 
     366             : ! set non-linear core correction
     367           2 :             IF (gth_atompot%nlcc) THEN
     368           0 :                ALLOCATE (corden(nmesh), fe(nmesh), rc(nmesh))
     369           0 :                corden(:) = 0.0_dp
     370           0 :                n = gth_atompot%nexp_nlcc
     371           0 :                DO i = 1, n
     372           0 :                   al = gth_atompot%alpha_nlcc(i)
     373           0 :                   rc(:) = rp(:)/al
     374           0 :                   fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
     375           0 :                   DO j = 1, gth_atompot%nct_nlcc(i)
     376           0 :                      cval = gth_atompot%cval_nlcc(j, i)
     377           0 :                      corden(:) = corden(:) + fe(:)*rc(:)**(2*j - 2)*cval
     378             :                   END DO
     379             :                END DO
     380           0 :                fun(1:nmesh) = corden(1:nmesh)*rp(1:nmesh)
     381             :                CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_core", &
     382           0 :                                                          fun(1), nmesh)
     383           0 :                DEALLOCATE (corden, fe, rc)
     384             :             END IF
     385             : 
     386             : ! local potential
     387           6 :             ALLOCATE (locpot(nmesh))
     388        1860 :             locpot(:) = 0.0_dp
     389           2 :             CALL atom_local_potential(locpot, gth_atompot, rp)
     390        1860 :             fun(1:nmesh) = locpot(1:nmesh)
     391             :             CALL sirius_add_atom_type_radial_function(sctx, label, "vloc", &
     392           2 :                                                       fun(1), nmesh)
     393           2 :             DEALLOCATE (locpot)
     394             : !
     395           2 :             NULLIFY (density, wavefunction, wfninfo)
     396             :             CALL calculate_atomic_orbitals(atomic_kind_set(ikind), qs_kind_set(ikind), &
     397             :                                            density=density, wavefunction=wavefunction, &
     398           2 :                                            wfninfo=wfninfo, agrid=atom_grid)
     399             : 
     400             : ! set the atomic radial functions
     401           6 :             DO iwf = 1, SIZE(wavefunction, 2)
     402           4 :                focc = wfninfo(1, iwf)
     403           4 :                l = NINT(wfninfo(2, iwf))
     404             : ! we can not easily get the principal quantum number
     405           4 :                nu = -1
     406           4 :                IF (up) THEN
     407           0 :                   fun(1:nmesh) = wavefunction(1:nmesh, iwf)*rp(i)
     408             :                ELSE
     409        3720 :                   DO i = 1, nmesh
     410        3720 :                      fun(i) = wavefunction(nmesh - i + 1, iwf)*rp(i)
     411             :                   END DO
     412             :                END IF
     413             :                CALL sirius_add_atom_type_radial_function(sctx, &
     414             :                                                          label, "ps_atomic_wf", &
     415           6 :                                                          fun(1), nmesh, l=l, occ=REAL(focc, KIND=C_DOUBLE), n=nu)
     416             :             END DO
     417             : 
     418             : ! set total charge density of a free atom (to compute initial rho(r))
     419           2 :             IF (up) THEN
     420           0 :                fun(1:nmesh) = fourpi*density(1:nmesh)*atom_grid%rad(1:nmesh)**2
     421             :             ELSE
     422        1860 :                DO i = 1, nmesh
     423        1860 :                   fun(i) = fourpi*density(nmesh - i + 1)*atom_grid%rad(nmesh - i + 1)**2
     424             :                END DO
     425             :             END IF
     426             :             CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_total", &
     427           2 :                                                       fun(1), nmesh)
     428             : 
     429           2 :             IF (ASSOCIATED(density)) DEALLOCATE (density)
     430           2 :             IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
     431           2 :             IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
     432             : 
     433           2 :             CALL deallocate_grid_atom(atom_grid)
     434           2 :             DEALLOCATE (rp, fun)
     435           2 :             DEALLOCATE (gth_atompot)
     436             : !
     437             :          ELSE
     438             :             CALL cp_abort(__LOCATION__, &
     439           0 :                           "CP2K/SIRIUS: atomic kind needs UPF or GTH potential definition")
     440             :          END IF
     441             : 
     442             :          CALL get_qs_kind(qs_kind_set(ikind), &
     443             :                           dft_plus_u_atom=dft_plus_u_atom, &
     444             :                           l_of_dft_plus_u=lu, &
     445             :                           n_of_dft_plus_u=nu, &
     446             :                           u_minus_j_target=u_minus_j, &
     447             :                           U_of_dft_plus_u=U_u, &
     448             :                           J_of_dft_plus_u=J_u, &
     449             :                           alpha_of_dft_plus_u=alpha_u, &
     450             :                           beta_of_dft_plus_u=beta_u, &
     451             :                           J0_of_dft_plus_u=J0_u, &
     452          22 :                           occupation_of_dft_plus_u=occ_u)
     453             : 
     454          60 :          IF (dft_plus_u_atom) THEN
     455           0 :             IF (nu < 1) THEN
     456           0 :                CPABORT("CP2K/SIRIUS (hubbard): principal quantum number not specified")
     457             :             END IF
     458             : 
     459           0 :             IF (lu < 0) THEN
     460           0 :                CPABORT("CP2K/SIRIUS (hubbard): l can not be negative.")
     461             :             END IF
     462             : 
     463           0 :             IF (occ_u < 0.0) THEN
     464           0 :                CPABORT("CP2K/SIRIUS (hubbard): the occupation number can not be negative.")
     465             :             END IF
     466             : 
     467           0 :             IF (ABS(u_minus_j) < 1e-8) THEN
     468             :                CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
     469           0 :                                                  occ_u, U_u, J_u, alpha_u, beta_u, J0_u)
     470             :             ELSE
     471             :                CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
     472           0 :                                                  occ_u, u_minus_j, 0.0_dp, alpha_u, beta_u, J0_u)
     473             :             END IF
     474             :          END IF
     475             : 
     476             :       END DO
     477             : 
     478             : ! add atoms to the unit cell
     479             : ! WARNING: sirius accepts only fractional coordinates;
     480          16 :       natom = SIZE(particle_set)
     481          50 :       DO iatom = 1, natom
     482         136 :          vr(1:3) = particle_set(iatom)%r(1:3)
     483          34 :          CALL real_to_scaled(vs, vr, my_cell)
     484          34 :          atomic_kind => particle_set(iatom)%atomic_kind
     485          34 :          ikind = atomic_kind%kind_number
     486          34 :          CALL get_atomic_kind(atomic_kind, name=label)
     487          34 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, magnetization=magnetization)
     488             : ! angle of magnetization might come from input Atom x y z mx my mz
     489             : ! or as an angle?
     490             : ! Answer : SIRIUS only accept the magnetization as mx, my, mz
     491          34 :          IF (num_mag_dims .EQ. 3) THEN
     492           2 :             angle1 = 0.0_dp
     493           2 :             angle2 = 0.0_dp
     494           2 :             v1(1) = magnetization*SIN(angle1)*COS(angle2)
     495           2 :             v1(2) = magnetization*SIN(angle1)*SIN(angle2)
     496           2 :             v1(3) = magnetization*COS(angle1)
     497             :          ELSE
     498          32 :             v1 = 0._dp
     499          32 :             v1(3) = magnetization
     500             :          END IF
     501          34 :          v2(1:3) = vs(1:3)
     502          84 :          CALL sirius_add_atom(sctx, label, v2(1), v1(1))
     503             :       END DO
     504             : 
     505          16 :       CALL sirius_set_mpi_grid_dims(sctx, 2, mpi_grid_dims)
     506             : 
     507             : ! initialize global variables/indices/arrays/etc. of the simulation
     508          16 :       CALL sirius_initialize_context(sctx)
     509             : 
     510             :       ! strictly speaking the parameter use_symmetry is initialized at the
     511             :       ! beginning but it does no harm to do it that way
     512          16 :       IF (use_symmetry) THEN
     513          14 :          CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.TRUE., kset_handler=ks_handler)
     514             :       ELSE
     515           2 :          CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.FALSE., kset_handler=ks_handler)
     516             :       END IF
     517             : ! create ground-state class
     518          16 :       CALL sirius_create_ground_state(ks_handler, gs_handler)
     519             : 
     520          16 :       CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler, ks_handler=ks_handler)
     521             : #endif
     522          32 :    END SUBROUTINE cp_sirius_create_env
     523             : 
     524             : !***************************************************************************************************
     525             : !> \brief ...
     526             : !> \param pwdft_env ...
     527             : !> \param
     528             : !> \par History
     529             : !>      07.2018 Update the Sirius environment
     530             : !> \author JHU
     531             : ! **************************************************************************************************
     532          16 :    SUBROUTINE cp_sirius_update_context(pwdft_env)
     533             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     534             : 
     535             :       INTEGER                                            :: iatom, natom
     536             :       REAL(KIND=C_DOUBLE), DIMENSION(3)                  :: a1, a2, a3, v2
     537             :       REAL(KIND=dp), DIMENSION(3)                        :: vr, vs
     538             :       TYPE(cell_type), POINTER                           :: my_cell
     539          16 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     540             :       TYPE(qs_subsys_type), POINTER                      :: qs_subsys
     541             :       TYPE(sirius_context_handler)                       :: sctx
     542             :       TYPE(sirius_ground_state_handler)                  :: gs_handler
     543             : 
     544           0 :       CPASSERT(ASSOCIATED(pwdft_env))
     545          16 :       CALL pwdft_env_get(pwdft_env, sctx=sctx, gs_handler=gs_handler)
     546             : 
     547             : ! get current positions and lattice vectors
     548          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
     549             : 
     550             : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
     551          16 :       CALL qs_subsys_get(qs_subsys, cell=my_cell)
     552          64 :       a1(:) = my_cell%hmat(:, 1)
     553          64 :       a2(:) = my_cell%hmat(:, 2)
     554          64 :       a3(:) = my_cell%hmat(:, 3)
     555          16 :       CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
     556             : 
     557             : ! new atomic positions
     558          16 :       CALL qs_subsys_get(qs_subsys, particle_set=particle_set)
     559          16 :       natom = SIZE(particle_set)
     560          50 :       DO iatom = 1, natom
     561         136 :          vr(1:3) = particle_set(iatom)%r(1:3)
     562          34 :          CALL real_to_scaled(vs, vr, my_cell)
     563          34 :          v2(1:3) = vs(1:3)
     564          50 :          CALL sirius_set_atom_position(sctx, iatom, v2(1))
     565             :       END DO
     566             : 
     567             : ! update ground-state class
     568          16 :       CALL sirius_update_ground_state(gs_handler)
     569             : 
     570          16 :       CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler)
     571             : 
     572          16 :    END SUBROUTINE cp_sirius_update_context
     573             : 
     574             : ! **************************************************************************************************
     575             : !> \brief ...
     576             : !> \param sctx ...
     577             : !> \param section ...
     578             : !> \param section_name ...
     579             : ! **************************************************************************************************
     580          80 :    SUBROUTINE cp_sirius_fill_in_section(sctx, section, section_name)
     581             :       TYPE(sirius_context_handler), INTENT(INOUT)        :: sctx
     582             :       TYPE(section_vals_type), POINTER                   :: section
     583             :       CHARACTER(*), INTENT(in)                           :: section_name
     584             : 
     585             :       CHARACTER(len=256), TARGET                         :: option_name
     586             :       CHARACTER(len=4096)                                :: description, usage
     587          80 :       CHARACTER(len=80), DIMENSION(:), POINTER           :: tmp
     588             :       CHARACTER(len=80), TARGET                          :: str
     589             :       INTEGER                                            :: ctype, elem, ic, j
     590          80 :       INTEGER, DIMENSION(:), POINTER                     :: ivals
     591             :       INTEGER, TARGET                                    :: enum_length, ival, length, &
     592             :                                                             num_possible_values, number_of_options
     593             :       LOGICAL                                            :: explicit
     594          80 :       LOGICAL, DIMENSION(:), POINTER                     :: lvals
     595             :       LOGICAL, TARGET                                    :: found, lval
     596          80 :       REAL(kind=dp), DIMENSION(:), POINTER               :: rvals
     597             :       REAL(kind=dp), TARGET                              :: rval
     598             : 
     599          80 :       NULLIFY (rvals)
     600          80 :       NULLIFY (ivals)
     601          80 :       CALL sirius_option_get_section_length(section_name, number_of_options)
     602             : 
     603        1696 :       DO elem = 1, number_of_options
     604        1616 :          option_name = ''
     605             :          CALL sirius_option_get_info(section_name, &
     606             :                                      elem, &
     607             :                                      option_name, &
     608             :                                      256, &
     609             :                                      ctype, &
     610             :                                      num_possible_values, &
     611             :                                      enum_length, &
     612             :                                      description, &
     613             :                                      4096, &
     614             :                                      usage, &
     615        1616 :                                      4096)
     616        1696 :          IF ((option_name /= 'memory_usage') .AND. (option_name /= 'xc_functionals') .AND. (option_name /= 'vk')) THEN
     617        1584 :             CALL section_vals_val_get(section, option_name, explicit=found)
     618        1584 :             IF (found) THEN
     619         128 :                SELECT CASE (ctype)
     620             :                CASE (SIRIUS_INTEGER_TYPE)
     621         128 :                   CALL section_vals_val_get(section, option_name, i_val=ival)
     622         128 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ival))
     623             :                CASE (SIRIUS_NUMBER_TYPE)
     624         104 :                   CALL section_vals_val_get(section, option_name, r_val=rval)
     625         104 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rval))
     626             :                CASE (SIRIUS_LOGICAL_TYPE)
     627          24 :                   CALL section_vals_val_get(section, option_name, l_val=lval)
     628          24 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lval))
     629             :                CASE (SIRIUS_STRING_TYPE)      ! string nightmare
     630         100 :                   str = ''
     631         100 :                   CALL section_vals_val_get(section, option_name, explicit=explicit, c_val=str)
     632         100 :                   str = TRIM(ADJUSTL(str))
     633        8100 :                   DO j = 1, LEN(str)
     634        8000 :                      ic = ICHAR(str(j:j))
     635        8100 :                      IF (ic >= 65 .AND. ic < 90) str(j:j) = CHAR(ic + 32)
     636             :                   END DO
     637             : 
     638         100 :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), max_length=LEN_TRIM(str))
     639             :                CASE (SIRIUS_INTEGER_ARRAY_TYPE)
     640          16 :                   CALL section_vals_val_get(section, option_name, i_vals=ivals)
     641             :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ivals(1)), &
     642          16 :                                          max_length=num_possible_values)
     643             :                CASE (SIRIUS_NUMBER_ARRAY_TYPE)
     644           0 :                   CALL section_vals_val_get(section, option_name, r_vals=rvals)
     645             :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rvals(1)), &
     646           0 :                                          max_length=num_possible_values)
     647             :                CASE (SIRIUS_LOGICAL_ARRAY_TYPE)
     648           0 :                   CALL section_vals_val_get(section, option_name, l_vals=lvals)
     649             :                   CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lvals(1)), &
     650           0 :                                          max_length=num_possible_values)
     651             :                CASE (SIRIUS_STRING_ARRAY_TYPE)
     652           0 :                   CALL section_vals_val_get(section, option_name, explicit=explicit, n_rep_val=length)
     653         372 :                   DO j = 1, length
     654           0 :                      str = ''
     655           0 :                      CALL section_vals_val_get(section, option_name, i_rep_val=j, explicit=explicit, c_vals=tmp)
     656           0 :                      str = TRIM(ADJUSTL(tmp(j)))
     657             :                      CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), &
     658           0 :                                             max_length=LEN_TRIM(str), append=.TRUE.)
     659             :                   END DO
     660             :                CASE DEFAULT
     661             :                END SELECT
     662             :             END IF
     663             :          END IF
     664             :       END DO
     665          80 :    END SUBROUTINE cp_sirius_fill_in_section
     666             : 
     667             : !***************************************************************************************************
     668             : !> \brief ...
     669             : !> \param pwdft_env ...
     670             : !> \param calculate_forces ...
     671             : !> \param calculate_stress_tensor ...
     672             : !> \param
     673             : !> \par History
     674             : !>      07.2018 start the Sirius library
     675             : !> \author JHU
     676             : ! **************************************************************************************************
     677          16 :    SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_tensor)
     678             :       TYPE(pwdft_environment_type), INTENT(INOUT), &
     679             :          POINTER                                         :: pwdft_env
     680             :       LOGICAL, INTENT(IN)                                :: calculate_forces, calculate_stress_tensor
     681             : 
     682             :       INTEGER                                            :: iw, n1, n2
     683             :       LOGICAL                                            :: do_print, gs_converged
     684             :       REAL(KIND=C_DOUBLE)                                :: etotal
     685          16 :       REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :)  :: cforces
     686             :       REAL(KIND=C_DOUBLE), DIMENSION(3, 3)               :: cstress
     687             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: stress
     688          16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: forces
     689             :       TYPE(cp_logger_type), POINTER                      :: logger
     690             :       TYPE(pwdft_energy_type), POINTER                   :: energy
     691             :       TYPE(section_vals_type), POINTER                   :: print_section, pwdft_input
     692             :       TYPE(sirius_ground_state_handler)                  :: gs_handler
     693             : 
     694           0 :       CPASSERT(ASSOCIATED(pwdft_env))
     695             : 
     696          16 :       NULLIFY (logger)
     697          16 :       logger => cp_get_default_logger()
     698          16 :       iw = cp_logger_get_default_io_unit(logger)
     699             : 
     700          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, gs_handler=gs_handler)
     701          16 :       CALL sirius_find_ground_state(gs_handler, converged=gs_converged)
     702             : 
     703          16 :       IF (gs_converged) THEN
     704          16 :          IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS: ground state is converged"
     705             :       ELSE
     706           0 :          IF (pwdft_env%ignore_convergence_failure) THEN
     707           0 :             IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS Warning: ground state is not converged"
     708             :          ELSE
     709           0 :             CPABORT("CP2K/SIRIUS (ground state): SIRIUS did not converge.")
     710             :          END IF
     711             :       END IF
     712          16 :       IF (iw > 0) CALL m_flush(iw)
     713             : 
     714          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, energy=energy)
     715             :       etotal = 0.0_C_DOUBLE
     716             : 
     717          16 :       CALL sirius_get_energy(gs_handler, 'band-gap', etotal)
     718          16 :       energy%band_gap = etotal
     719             : 
     720             :       etotal = 0.0_C_DOUBLE
     721          16 :       CALL sirius_get_energy(gs_handler, 'total', etotal)
     722          16 :       energy%etotal = etotal
     723             : 
     724             :       ! extract entropy (TS returned by sirius is always negative, sign
     725             :       ! convention in QE)
     726             :       etotal = 0.0_C_DOUBLE
     727          16 :       CALL sirius_get_energy(gs_handler, 'demet', etotal)
     728          16 :       energy%entropy = -etotal
     729             : 
     730          16 :       IF (calculate_forces) THEN
     731          12 :          CALL pwdft_env_get(pwdft_env=pwdft_env, forces=forces)
     732          12 :          n1 = SIZE(forces, 1)
     733          12 :          n2 = SIZE(forces, 2)
     734             : 
     735          48 :          ALLOCATE (cforces(n2, n1))
     736         132 :          cforces = 0.0_C_DOUBLE
     737          12 :          CALL sirius_get_forces(gs_handler, 'total', cforces)
     738             :          ! Sirius computes the forces but cp2k use the gradient everywhere
     739             :          ! so a minus sign is needed.
     740             :          ! note also that sirius and cp2k store the forces transpose to each other
     741             :          ! sirius : forces(coordinates, atoms)
     742             :          ! cp2k : forces(atoms, coordinates)
     743         138 :          forces = -TRANSPOSE(cforces(:, :))
     744          12 :          DEALLOCATE (cforces)
     745             :       END IF
     746             : 
     747          16 :       IF (calculate_stress_tensor) THEN
     748           0 :          cstress = 0.0_C_DOUBLE
     749           0 :          CALL sirius_get_stress_tensor(gs_handler, 'total', cstress)
     750           0 :          stress(1:3, 1:3) = cstress(1:3, 1:3)
     751           0 :          CALL pwdft_env_set(pwdft_env=pwdft_env, stress=stress)
     752             :       END IF
     753             : 
     754          16 :       CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_input)
     755          16 :       print_section => section_vals_get_subs_vals(pwdft_input, "PRINT")
     756          16 :       CALL section_vals_get(print_section, explicit=do_print)
     757          16 :       IF (do_print) THEN
     758           2 :          CALL cp_sirius_print_results(pwdft_env, print_section)
     759             :       END IF
     760          32 :    END SUBROUTINE cp_sirius_energy_force
     761             : 
     762             : !***************************************************************************************************
     763             : !> \brief ...
     764             : !> \param pwdft_env ...
     765             : !> \param print_section ...
     766             : !> \param
     767             : !> \par History
     768             : !>      12.2019 init
     769             : !> \author JHU
     770             : ! **************************************************************************************************
     771           2 :    SUBROUTINE cp_sirius_print_results(pwdft_env, print_section)
     772             :       TYPE(pwdft_environment_type), INTENT(INOUT), &
     773             :          POINTER                                         :: pwdft_env
     774             :       TYPE(section_vals_type), POINTER                   :: print_section
     775             : 
     776             :       CHARACTER(LEN=default_string_length)               :: my_act, my_pos
     777             :       INTEGER                                            :: i, ik, iounit, ispn, iterstep, iv, iw, &
     778             :                                                             nbands, nhist, nkpts, nspins
     779             :       INTEGER(KIND=C_INT)                                :: cint
     780             :       LOGICAL                                            :: append, dos, ionode
     781             :       REAL(KIND=C_DOUBLE)                                :: creal
     782           2 :       REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:)     :: slist
     783             :       REAL(KIND=dp)                                      :: de, e_fermi(2), emax, emin, eval
     784           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: wkpt
     785           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ehist, hist, occval
     786           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: energies, occupations
     787             :       TYPE(cp_logger_type), POINTER                      :: logger
     788             :       TYPE(sirius_context_handler)                       :: sctx
     789             :       TYPE(sirius_ground_state_handler)                  :: gs_handler
     790             :       TYPE(sirius_kpoint_set_handler)                    :: ks_handler
     791             : 
     792           2 :       NULLIFY (logger)
     793           4 :       logger => cp_get_default_logger()
     794             :       ionode = logger%para_env%is_source()
     795           2 :       iounit = cp_logger_get_default_io_unit(logger)
     796             : 
     797             :       ! Density of States
     798           2 :       dos = BTEST(cp_print_key_should_output(logger%iter_info, print_section, "DOS"), cp_p_file)
     799           2 :       IF (dos) THEN
     800           2 :          CALL pwdft_env_get(pwdft_env, ks_handler=ks_handler)
     801           2 :          CALL pwdft_env_get(pwdft_env, gs_handler=gs_handler)
     802           2 :          CALL pwdft_env_get(pwdft_env, sctx=sctx)
     803             : 
     804           2 :          CALL section_vals_val_get(print_section, "DOS%DELTA_E", r_val=de)
     805           2 :          CALL section_vals_val_get(print_section, "DOS%APPEND", l_val=append)
     806             : 
     807           2 :          CALL sirius_get_num_kpoints(ks_handler, cint)
     808           2 :          nkpts = cint
     809           2 :          CALL sirius_get_parameters(sctx, num_bands=cint)
     810           2 :          nbands = cint
     811           2 :          CALL sirius_get_parameters(sctx, num_spins=cint)
     812           2 :          nspins = cint
     813           2 :          e_fermi(:) = 0.0_dp
     814          10 :          ALLOCATE (energies(nbands, nspins, nkpts))
     815         442 :          energies = 0.0_dp
     816           8 :          ALLOCATE (occupations(nbands, nspins, nkpts))
     817         442 :          occupations = 0.0_dp
     818           6 :          ALLOCATE (wkpt(nkpts))
     819           6 :          ALLOCATE (slist(nbands))
     820          10 :          DO ik = 1, nkpts
     821           8 :             CALL sirius_get_kpoint_properties(ks_handler, ik, creal)
     822          10 :             wkpt(ik) = creal
     823             :          END DO
     824          10 :          DO ik = 1, nkpts
     825          26 :             DO ispn = 1, nspins
     826          16 :                CALL sirius_get_band_energies(ks_handler, ik, ispn, slist)
     827         432 :                energies(1:nbands, ispn, ik) = slist(1:nbands)
     828          16 :                CALL sirius_get_band_occupancies(ks_handler, ik, ispn, slist)
     829         440 :                occupations(1:nbands, ispn, ik) = slist(1:nbands)
     830             :             END DO
     831             :          END DO
     832         442 :          emin = MINVAL(energies)
     833         442 :          emax = MAXVAL(energies)
     834           2 :          nhist = NINT((emax - emin)/de) + 1
     835          16 :          ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
     836       16110 :          hist = 0.0_dp
     837       16110 :          occval = 0.0_dp
     838       16110 :          ehist = 0.0_dp
     839             : 
     840          10 :          DO ik = 1, nkpts
     841          26 :             DO ispn = 1, nspins
     842         440 :                DO i = 1, nbands
     843         416 :                   eval = energies(i, ispn, ik) - emin
     844         416 :                   iv = NINT(eval/de) + 1
     845         416 :                   CPASSERT((iv > 0) .AND. (iv <= nhist))
     846         416 :                   hist(iv, ispn) = hist(iv, ispn) + wkpt(ik)
     847         432 :                   occval(iv, ispn) = occval(iv, ispn) + wkpt(ik)*occupations(i, ispn, ik)
     848             :                END DO
     849             :             END DO
     850             :          END DO
     851       16110 :          hist = hist/REAL(nbands, KIND=dp)
     852        8054 :          DO i = 1, nhist
     853       24158 :             ehist(i, 1:nspins) = emin + (i - 1)*de
     854             :          END DO
     855             : 
     856           2 :          iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
     857           2 :          my_act = "WRITE"
     858           2 :          IF (append .AND. iterstep > 1) THEN
     859           0 :             my_pos = "APPEND"
     860             :          ELSE
     861           2 :             my_pos = "REWIND"
     862             :          END IF
     863             : 
     864             :          iw = cp_print_key_unit_nr(logger, print_section, "DOS", &
     865             :                                    extension=".dos", file_position=my_pos, file_action=my_act, &
     866           2 :                                    file_form="FORMATTED")
     867           2 :          IF (iw > 0) THEN
     868           1 :             IF (nspins == 2) THEN
     869             :                WRITE (UNIT=iw, FMT="(T2,A,I0,A,2F12.6)") &
     870           1 :                   "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
     871           1 :                WRITE (UNIT=iw, FMT="(T2,A, A)") "   Energy[a.u.]  Alpha_Density     Occupation", &
     872           2 :                   "   Beta_Density      Occupation"
     873             :             ELSE
     874             :                WRITE (UNIT=iw, FMT="(T2,A,I0,A,F12.6)") &
     875           0 :                   "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
     876           0 :                WRITE (UNIT=iw, FMT="(T2,A)") "   Energy[a.u.]       Density     Occupation"
     877             :             END IF
     878        4027 :             DO i = 1, nhist
     879        4026 :                eval = emin + (i - 1)*de
     880        4027 :                IF (nspins == 2) THEN
     881        4026 :                   WRITE (UNIT=iw, FMT="(F15.8,4F15.4)") eval, hist(i, 1), occval(i, 1), &
     882        8052 :                      hist(i, 2), occval(i, 2)
     883             :                ELSE
     884           0 :                   WRITE (UNIT=iw, FMT="(F15.8,2F15.4)") eval, hist(i, 1), occval(i, 1)
     885             :                END IF
     886             :             END DO
     887             :          END IF
     888           2 :          CALL cp_print_key_finished_output(iw, logger, print_section, "DOS")
     889             : 
     890           2 :          DEALLOCATE (energies, occupations, wkpt, slist)
     891           4 :          DEALLOCATE (hist, occval, ehist)
     892             :       END IF
     893           4 :    END SUBROUTINE cp_sirius_print_results
     894             : 
     895             : END MODULE sirius_interface
     896             : 
     897             : #else
     898             : 
     899             : !***************************************************************************************************
     900             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     901             : !***************************************************************************************************
     902             : MODULE sirius_interface
     903             :    USE pwdft_environment_types, ONLY: pwdft_environment_type
     904             : #include "./base/base_uses.f90"
     905             : 
     906             :    IMPLICIT NONE
     907             :    PRIVATE
     908             : 
     909             :    PUBLIC :: cp_sirius_init, cp_sirius_finalize
     910             :    PUBLIC :: cp_sirius_create_env, cp_sirius_energy_force, cp_sirius_update_context
     911             : 
     912             : CONTAINS
     913             : 
     914             : ! **************************************************************************************************
     915             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     916             : ! **************************************************************************************************
     917             :    SUBROUTINE cp_sirius_init()
     918             :    END SUBROUTINE cp_sirius_init
     919             : 
     920             : ! **************************************************************************************************
     921             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     922             : ! **************************************************************************************************
     923             :    SUBROUTINE cp_sirius_finalize()
     924             :    END SUBROUTINE cp_sirius_finalize
     925             : 
     926             : ! **************************************************************************************************
     927             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     928             : !> \param pwdft_env ...
     929             : ! **************************************************************************************************
     930             :    SUBROUTINE cp_sirius_create_env(pwdft_env)
     931             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     932             : 
     933             :       MARK_USED(pwdft_env)
     934             :       CPABORT("Sirius library is missing")
     935             :    END SUBROUTINE cp_sirius_create_env
     936             : 
     937             : ! **************************************************************************************************
     938             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     939             : !> \param pwdft_env ...
     940             : !> \param calculate_forces ...
     941             : !> \param calculate_stress ...
     942             : ! **************************************************************************************************
     943             :    SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
     944             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     945             :       LOGICAL                                            :: calculate_forces, calculate_stress
     946             : 
     947             :       MARK_USED(pwdft_env)
     948             :       MARK_USED(calculate_forces)
     949             :       MARK_USED(calculate_stress)
     950             :       CPABORT("Sirius library is missing")
     951             :    END SUBROUTINE cp_sirius_energy_force
     952             : 
     953             : ! **************************************************************************************************
     954             : !> \brief Empty implementation in case SIRIUS is not compiled in.
     955             : !> \param pwdft_env ...
     956             : ! **************************************************************************************************
     957             :    SUBROUTINE cp_sirius_update_context(pwdft_env)
     958             :       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
     959             : 
     960             :       MARK_USED(pwdft_env)
     961             :       CPABORT("Sirius library is missing")
     962             :    END SUBROUTINE cp_sirius_update_context
     963             : 
     964             : END MODULE sirius_interface
     965             : 
     966             : #endif

Generated by: LCOV version 1.15