LCOV - code coverage report
Current view: top level - force - fix_by_gaussian.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 32 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 3 0.0 %

          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,stars,mpi,sym,vacuum,sphhar,input,oned,cell,noco,den)
      12             :     !The idea of this fix is to add an Gaussian to the INT which make the charge flat at the MT-boundary and to 
      13             :     !shift this Gaussian with the displacement.
      14             :     USE m_qfix
      15             :     USE m_spgrot
      16             :     USE m_constants
      17             :     USE m_types
      18             :     IMPLICIT NONE
      19             :     REAL,INTENT(in)             :: shift(:,:)
      20             :     TYPE(t_mpi),INTENT(IN)      :: mpi
      21             :     TYPE(t_atoms),INTENT(IN)    :: atoms
      22             :     TYPE(t_sym),INTENT(IN)      :: sym
      23             :     TYPE(t_vacuum),INTENT(IN)   :: vacuum
      24             :     TYPE(t_sphhar),INTENT(IN)   :: sphhar
      25             :     TYPE(t_input),INTENT(IN)    :: input
      26             :     TYPE(t_oneD),INTENT(IN)     :: oneD
      27             :     TYPE(t_cell),INTENT(IN)     :: cell
      28             :     TYPE(t_noco),INTENT(IN)     :: noco
      29             :     TYPE(t_stars),INTENT(IN)    :: stars
      30             :     TYPE(t_potden),INTENT(INOUT):: den
      31             : 
      32             :     REAL    :: slope1,slope2,dr,alpha
      33             :     REAL    :: sigma,a_fac,gauss,x,fix
      34           0 :     INTEGER :: kr(3,sym%nop)
      35           0 :     COMPLEX :: sf,phas(sym%nop)
      36             :     INTEGER :: js,n,l,k,nat,j
      37             : 
      38           0 :     DO js=1,input%jspins
      39           0 :        DO n=1,atoms%ntype
      40           0 :           DO l=0,0 !Currently only l=0
      41             : !             alpha = LOG( den%mt(atoms%jri(n)-1,l,n,js) / den%mt(atoms%jri(n),l,n,js) )
      42             : !             alpha = SQRT(alpha / ( atoms%rmt(n)*atoms%rmt(n)*( 1.0-EXP( -2.0*atoms%dx(n) ) ) ))
      43             : !             A_fac= den%mt(atoms%jri(n),l,n,js)/gaussian_r(atoms%rmt(n),alpha)
      44           0 :              dr=atoms%rmsh(atoms%jri(n)-1,n)-atoms%rmsh(atoms%jri(n),n)
      45           0 :              slope1=(den%mt(atoms%jri(n)-1,l,n,js)/atoms%rmsh(atoms%jri(n)-1,n)**2-den%mt(atoms%jri(n),l,n,js)/atoms%rmsh(atoms%jri(n),n)**2)/dr
      46           0 :              slope1=(den%mt(atoms%jri(n)-1,l,n,js)-den%mt(atoms%jri(n),l,n,js))/dr/atoms%rmt(n)**2
      47           0 :              sigma=atoms%rmt(n)/2.0
      48           0 :              alpha=1/sigma
      49           0 :              slope2=(gaussian_r(atoms%rmsh(atoms%jri(n)-1,n),alpha)-gaussian_r(atoms%rmsh(atoms%jri(n),n),alpha))/dr
      50           0 :              A_fac=slope1/slope2
      51           0 :              PRINT *, a_fac, 1/alpha
      52           0 :              DO k=2,stars%ng3
      53           0 :                 gauss=A_fac*gaussian_g(stars%sk3(k),alpha)
      54           0 :                 CALL spgrot(sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,stars%kv3(:,k),kr,phas)
      55           0 :                 DO nat=SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
      56             :                    sf=0.0
      57           0 :                    DO  j = 1,sym%nop
      58           0 :                       x=-tpi_const*DOT_PRODUCT(1.*kr(:,j),atoms%taual(:,nat))
      59           0 :                       sf = sf + CMPLX(COS(x),SIN(x))*CONJG(phas(j))
      60           0 :                       x=-tpi_const*DOT_PRODUCT(1.*kr(:,j),atoms%taual(:,nat)+shift(:,nat))
      61           0 :                       sf = sf - CMPLX(COS(x),SIN(x))*CONJG(phas(j))
      62             :                    ENDDO
      63             :                 ENDDO
      64           0 :                 den%pw(k,js)=den%pw(k,js)+gauss*sf/sym%nop/cell%omtil
      65             :              ENDDO
      66             :           ENDDO
      67             :        END DO
      68             :     END DO
      69             :     CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
      70           0 :          den,noco%l_noco,mpi%isize==1,force_fix=.TRUE.,fix=fix,fix_pw_only=.true.)
      71           0 :   END SUBROUTINE fix_by_gaussian
      72             : 
      73           0 :   FUNCTION gaussian_r(r,alpha)
      74             :    
      75             :     REAL,INTENT(IN) :: r,alpha
      76             :     real            :: gaussian_r
      77           0 :     gaussian_r=EXP(-r**2*alpha**2)
      78           0 :   END FUNCTION gaussian_r
      79             : 
      80           0 :   FUNCTION gaussian_g(g,alpha)
      81             :     USE m_constants
      82             :     REAL,INTENT(IN) :: g,alpha
      83             :     REAL            :: gaussian_g
      84           0 :     gaussian_g=SQRT(pi_const**3)/alpha**3*EXP(-0.25*g**2/alpha**2)
      85           0 :   END FUNCTION gaussian_g
      86             : 
      87             : END MODULE m_fix_by_gaussian

Generated by: LCOV version 1.13