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 Does all kind of post scf calculations for GPW/GAPW
      10             : !> \par History
      11             : !>      Taken out from qs_scf_post_gpw
      12             : !> \author JGH
      13             : ! **************************************************************************************************
      14             : MODULE qs_elf_methods
      15             :    USE dbcsr_api,                       ONLY: dbcsr_p_type
      16             :    USE kinds,                           ONLY: dp
      17             :    USE mathconstants,                   ONLY: pi
      18             :    USE pw_env_types,                    ONLY: pw_env_get,&
      19             :                                               pw_env_type
      20             :    USE pw_methods,                      ONLY: pw_derive,&
      21             :                                               pw_transfer,&
      22             :                                               pw_zero
      23             :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      24             :                                               pw_pool_type
      25             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      26             :                                               pw_r3d_rs_type
      27             :    USE qs_collocate_density,            ONLY: calculate_rho_elec
      28             :    USE qs_environment_types,            ONLY: get_qs_env,&
      29             :                                               qs_environment_type
      30             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      31             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      32             :                                               qs_rho_type
      33             : #include "./base/base_uses.f90"
      34             : 
      35             :    IMPLICIT NONE
      36             :    PRIVATE
      37             : 
      38             :    ! Global parameters
      39             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_elf_methods'
      40             : 
      41             :    PUBLIC :: qs_elf_calc
      42             : 
      43             : ! **************************************************************************************************
      44             : 
      45             : CONTAINS
      46             : 
      47             : ! **************************************************************************************************
      48             : !> \brief ...
      49             : !> \param qs_env ...
      50             : !> \param elf_r ...
      51             : !> \param rho_cutoff ...
      52             : ! **************************************************************************************************
      53          82 :    SUBROUTINE qs_elf_calc(qs_env, elf_r, rho_cutoff)
      54             : 
      55             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      56             :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN)     :: elf_r
      57             :       REAL(kind=dp), INTENT(IN)                          :: rho_cutoff
      58             : 
      59             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_elf_calc'
      60             :       INTEGER, DIMENSION(3, 3), PARAMETER :: nd = RESHAPE((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/))
      61             :       REAL(KIND=dp), PARAMETER                           :: ELFCUT = 0.0001_dp, &
      62             :                                                             f18 = (1.0_dp/8.0_dp), &
      63             :                                                             f23 = (2.0_dp/3.0_dp), &
      64             :                                                             f53 = (5.0_dp/3.0_dp)
      65             : 
      66             :       INTEGER                                            :: handle, i, idir, ispin, j, k, nspin
      67             :       INTEGER, DIMENSION(2, 3)                           :: bo
      68             :       LOGICAL                                            :: deriv_pw, drho_r_valid, tau_r_valid
      69             :       REAL(kind=dp)                                      :: cfermi, elf_kernel, norm_drho, rho_53, &
      70             :                                                             udvol
      71          82 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
      72          82 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_struct_ao
      73             :       TYPE(pw_c1d_gs_type)                               :: tmp_g
      74             :       TYPE(pw_env_type), POINTER                         :: pw_env
      75          82 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
      76             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      77         328 :       TYPE(pw_r3d_rs_type), DIMENSION(3)                 :: drho_r
      78          82 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_struct_r, tau_struct_r
      79          82 :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER     :: drho_struct_r
      80             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_r, tau_r
      81             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
      82             :       TYPE(qs_rho_type), POINTER                         :: rho_struct
      83             : 
      84          82 :       CALL timeset(routineN, handle)
      85             : 
      86          82 :       NULLIFY (rho_struct, rho_r, tau_r, pw_env, auxbas_pw_pool, pw_pools, ks_env)
      87          82 :       NULLIFY (rho_struct_ao, rho_struct_r, tau_struct_r, drho_struct_r)
      88             : 
      89          82 :       CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env, rho=rho_struct)
      90             : 
      91             :       CALL qs_rho_get(rho_struct, &
      92             :                       rho_ao_kp=rho_struct_ao, &
      93             :                       rho_r=rho_struct_r, &
      94             :                       tau_r=tau_struct_r, &
      95             :                       drho_r=drho_struct_r, &
      96             :                       tau_r_valid=tau_r_valid, &
      97          82 :                       drho_r_valid=drho_r_valid)
      98             : 
      99             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     100          82 :                       pw_pools=pw_pools)
     101          82 :       nspin = SIZE(rho_struct_r)
     102         820 :       bo = rho_struct_r(1)%pw_grid%bounds_local
     103          82 :       cfermi = (3.0_dp/10.0_dp)*(pi*pi*3.0_dp)**f23
     104             : 
     105             :       ! In this case, we need a work matrix containing tau in g space
     106             :       ! We will not have further use for it, so we will need only one
     107          82 :       IF (.NOT. tau_r_valid) THEN
     108          82 :          ALLOCATE (tau_r)
     109          82 :          CALL auxbas_pw_pool%create_pw(tau_r)
     110             :       END IF
     111          82 :       IF (.NOT. tau_r_valid .OR. .NOT. drho_r_valid) THEN
     112          82 :          CALL auxbas_pw_pool%create_pw(tmp_g)
     113             :       END IF
     114          82 :       IF (.NOT. drho_r_valid) THEN
     115         328 :          DO idir = 1, 3
     116         328 :             CALL auxbas_pw_pool%create_pw(drho_r(idir))
     117             :          END DO
     118             :       END IF
     119             : 
     120         168 :       DO ispin = 1, nspin
     121          86 :          rho_r => rho_struct_r(ispin)
     122          86 :          IF (tau_r_valid) THEN
     123           0 :             tau_r => tau_struct_r(ispin)
     124             :          ELSE
     125          86 :             rho_ao => rho_struct_ao(ispin, :)
     126          86 :             CALL pw_zero(tau_r)
     127             :             CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     128             :                                     rho=tau_r, &
     129             :                                     rho_gspace=tmp_g, &
     130             :                                     ks_env=ks_env, soft_valid=.FALSE., &
     131          86 :                                     compute_tau=.TRUE.)
     132             :          END IF
     133             : 
     134          86 :          IF (drho_r_valid) THEN
     135           0 :             drho_r(:) = drho_struct_r(:, ispin)
     136             :          ELSE
     137          86 :             deriv_pw = .FALSE.
     138             :             IF (deriv_pw) THEN
     139             :                udvol = 1.0_dp/rho_r%pw_grid%dvol
     140             :                DO idir = 1, 3
     141             :                   CALL pw_transfer(rho_r, tmp_g)
     142             :                   CALL pw_derive(tmp_g, nd(:, idir))
     143             :                   CALL pw_transfer(tmp_g, drho_r(idir))
     144             :                END DO
     145             : 
     146             :             ELSE
     147         344 :                DO idir = 1, 3
     148         258 :                   rho_ao => rho_struct_ao(ispin, :)
     149             :                   CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     150             :                                           rho=drho_r(idir), &
     151             :                                           rho_gspace=tmp_g, &
     152             :                                           ks_env=ks_env, soft_valid=.FALSE., &
     153         344 :                                           compute_tau=.FALSE., compute_grad=.TRUE., idir=idir)
     154             : 
     155             :                END DO
     156             :             END IF
     157             :          END IF
     158             : 
     159             :          ! Calculate elf_r
     160             : !$OMP        PARALLEL DO DEFAULT(NONE) SHARED(bo,elf_r, ispin, drho_r,rho_r, tau_r, cfermi, rho_cutoff)&
     161         168 : !$OMP                    PRIVATE(k,j,i, norm_drho, rho_53, elf_kernel)
     162             :          DO k = bo(1, 3), bo(2, 3)
     163             :             DO j = bo(1, 2), bo(2, 2)
     164             :                DO i = bo(1, 1), bo(2, 1)
     165             :                   norm_drho = drho_r(1)%array(i, j, k)**2 + &
     166             :                               drho_r(2)%array(i, j, k)**2 + &
     167             :                               drho_r(3)%array(i, j, k)**2
     168             :                   norm_drho = norm_drho/MAX(rho_r%array(i, j, k), rho_cutoff)
     169             :                   rho_53 = cfermi*MAX(rho_r%array(i, j, k), rho_cutoff)**f53
     170             :                   elf_kernel = (tau_r%array(i, j, k) - f18*norm_drho) + 2.87E-5_dp
     171             :                   elf_kernel = (elf_kernel/rho_53)**2
     172             :                   elf_r(ispin)%array(i, j, k) = 1.0_dp/(1.0_dp + elf_kernel)
     173             :                   IF (elf_r(ispin)%array(i, j, k) < ELFCUT) elf_r(ispin)%array(i, j, k) = 0.0_dp
     174             :                END DO
     175             :             END DO
     176             :          END DO
     177             :       END DO
     178             : 
     179          82 :       IF (.NOT. drho_r_valid) THEN
     180         328 :          DO idir = 1, 3
     181         328 :             CALL auxbas_pw_pool%give_back_pw(drho_r(idir))
     182             :          END DO
     183             :       END IF
     184          82 :       IF (.NOT. tau_r_valid) THEN
     185          82 :          CALL auxbas_pw_pool%give_back_pw(tau_r)
     186          82 :          DEALLOCATE (tau_r)
     187             :       END IF
     188          82 :       IF (.NOT. tau_r_valid .OR. .NOT. drho_r_valid) THEN
     189          82 :          CALL auxbas_pw_pool%give_back_pw(tmp_g)
     190             :       END IF
     191             : 
     192          82 :       CALL timestop(handle)
     193             : 
     194          82 :    END SUBROUTINE qs_elf_calc
     195             : 
     196             : END MODULE qs_elf_methods