Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2019 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_fix_by_gaussian
       8             :    USE m_judft
       9             :    IMPLICIT NONE
      10             : CONTAINS
      11           0 :    SUBROUTINE fix_by_gaussian(shift,atoms,nococonv,stars,fmpi,sym,vacuum,sphhar,input ,cell,noco,den)
      12             :       ! The idea of this fix is to add an Gaussian to the INT which make the charge flat at the
      13             :       ! MT-boundary and to shift this Gaussian with the displacement.
      14             :       USE m_qfix
      15             :       USE m_spgrot
      16             :       USE m_constants
      17             :       USE m_types
      18             : 
      19             :       REAL,           INTENT(IN)    :: shift(:,:)
      20             :       TYPE(t_mpi),    INTENT(IN)    :: fmpi
      21             :       TYPE(t_atoms),  INTENT(IN)    :: atoms
      22             :       type(t_nococonv),INTENT(IN)   :: nococonv
      23             :       TYPE(t_sym),    INTENT(IN)    :: sym
      24             :       TYPE(t_vacuum), INTENT(IN)    :: vacuum
      25             :       TYPE(t_sphhar), INTENT(IN)    :: sphhar
      26             :       TYPE(t_input),  INTENT(IN)    :: input
      27             :        
      28             :       TYPE(t_cell),   INTENT(IN)    :: cell
      29             :       TYPE(t_noco),   INTENT(IN)    :: noco
      30             :       TYPE(t_stars),  INTENT(IN)    :: stars
      31             :       TYPE(t_potden), INTENT(INOUT) :: den
      32             : 
      33             :       REAL    :: slope1, slope2, dr, alpha
      34             :       REAL    :: sigma, a_fac, gauss, x, fix
      35           0 :       INTEGER :: kr(3,sym%nop)
      36           0 :       COMPLEX :: sf, phas(sym%nop)
      37             :       INTEGER :: js, n, l, k, nat, j
      38             : 
      39           0 :       DO js=1, input%jspins
      40           0 :          DO n=1, atoms%ntype
      41           0 :             DO l=0, 0 ! Currently only l=0
      42             :                ! alpha = LOG( den%mt(atoms%jri(n)-1,l,n,js) / den%mt(atoms%jri(n),l,n,js) )
      43             :                ! alpha = SQRT(alpha / ( atoms%rmt(n)*atoms%rmt(n)*( 1.0-EXP( -2.0*atoms%dx(n) ) ) ))
      44             :                ! A_fac= den%mt(atoms%jri(n),l,n,js)/gaussian_r(atoms%rmt(n),alpha)
      45           0 :                dr=atoms%rmsh(atoms%jri(n)-1,n)-atoms%rmsh(atoms%jri(n),n)
      46             :                slope1=(den%mt(atoms%jri(n)-1,l,n,js)/atoms%rmsh(atoms%jri(n)-1,n)**2 &
      47           0 :                       -den%mt(atoms%jri(n),l,n,js)/atoms%rmsh(atoms%jri(n),n)**2)/dr
      48             :                slope1=(den%mt(atoms%jri(n)-1,l,n,js) &
      49           0 :                       -den%mt(atoms%jri(n),l,n,js))/dr/atoms%rmt(n)**2 
      50             :                ! TODO: Only one of those slopes can be right
      51           0 :                sigma=atoms%rmt(n)/2.0
      52           0 :                alpha=1/sigma
      53           0 :                slope2=(gaussian_r(atoms%rmsh(atoms%jri(n)-1,n),alpha)-gaussian_r(atoms%rmsh(atoms%jri(n),n),alpha))/dr
      54           0 :                A_fac=slope1/slope2
      55           0 :                PRINT *, a_fac, 1/alpha
      56           0 :                DO k=2,stars%ng3
      57           0 :                   gauss=A_fac*gaussian_g(stars%sk3(k),alpha)
      58           0 :                   CALL spgrot(sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,stars%kv3(:,k),kr,phas)
      59           0 :                   DO nat = atoms%firstAtom(n), atoms%firstAtom(n) + atoms%neq(n) - 1
      60             :                      sf=0.0
      61           0 :                      DO  j = 1,sym%nop
      62           0 :                         x=-tpi_const*DOT_PRODUCT(1.*kr(:,j),atoms%taual(:,nat))
      63           0 :                         sf = sf + CMPLX(COS(x),SIN(x))*CONJG(phas(j))
      64           0 :                         x=-tpi_const*DOT_PRODUCT(1.*kr(:,j),atoms%taual(:,nat)+shift(:,nat))
      65           0 :                         sf = sf - CMPLX(COS(x),SIN(x))*CONJG(phas(j))
      66             :                      END DO
      67             :                   END DO
      68           0 :                   den%pw(k,js)=den%pw(k,js)+gauss*sf/sym%nop/cell%omtil
      69             :                END DO
      70             :             END DO
      71             :          END DO
      72             :       END DO
      73           0 :       CALL qfix(fmpi,stars,nococonv,atoms,sym,vacuum,sphhar,input,cell ,den,noco%l_noco,.FALSE.,l_par=.FALSE.,force_fix=.TRUE.,fix=fix,fix_pw_only=.true.)
      74           0 :    END SUBROUTINE fix_by_gaussian
      75             : 
      76           0 :    FUNCTION gaussian_r(r,alpha)   
      77             :       REAL, INTENT(IN) :: r,alpha
      78             :       REAL             :: gaussian_r
      79           0 :       gaussian_r=EXP(-r**2*alpha**2)
      80           0 :    END FUNCTION gaussian_r
      81             : 
      82           0 :    FUNCTION gaussian_g(g,alpha)
      83             :       USE m_constants
      84             :       REAL,INTENT(IN) :: g,alpha
      85             :       REAL            :: gaussian_g
      86             : 
      87           0 :       gaussian_g=SQRT(pi_const**3)/alpha**3*EXP(-0.25*g**2/alpha**2)
      88           0 :    END FUNCTION gaussian_g
      89             : 
      90             : END MODULE m_fix_by_gaussian