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

Generated by: LCOV version 1.15