Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 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             : MODULE m_kerker
       7             : 
       8             : CONTAINS
       9             : 
      10           6 :   SUBROUTINE kerker( field,  fmpi, &
      11             :        stars, atoms, sphhar, vacuum, input, sym, juphon, cell, noco, nococonv,&
      12             :          inDen, outDen, precon_v  )
      13             : 
      14             :     !Implementation of the Kerker preconditioner by M.Hinzen
      15             : 
      16             :     USE m_vgen_coulomb
      17             :     USE m_VYukawaFilm
      18             :     USE m_juDFT
      19             :     USE m_qfix
      20             :     USE m_types
      21             :     USE m_types_mixvector
      22             :     USE m_constants
      23             : 
      24             :     IMPLICIT NONE
      25             : 
      26             :      
      27             :     TYPE(t_input),     INTENT(in)    :: input
      28             :     TYPE(t_vacuum),    INTENT(in)    :: vacuum
      29             :     TYPE(t_noco),      INTENT(in)    :: noco
      30             :     TYPE(t_nococonv),      INTENT(in)    :: nococonv
      31             :     TYPE(t_sym),       INTENT(in)    :: sym
      32             :     TYPE(t_juphon),    INTENT(in)    :: juphon
      33             :     TYPE(t_stars),     INTENT(in)    :: stars
      34             :     TYPE(t_cell),      INTENT(in)    :: cell
      35             :     TYPE(t_sphhar),    INTENT(in)    :: sphhar
      36             :     TYPE(t_field),     INTENT(inout) :: field
      37             : 
      38             :     TYPE(t_mpi),       INTENT(in)    :: fmpi
      39             :     TYPE(t_atoms),     INTENT(in)    :: atoms
      40             :     TYPE(t_potden),    INTENT(inout) :: outDen
      41             :     TYPE(t_potden),    INTENT(in)    :: inDen
      42             :     TYPE(t_mixvector), INTENT(INOUT) :: precon_v
      43             : 
      44           6 :     type(t_potden)                   :: resDen, vYukawa, resDenMod
      45             :     real                             :: fix
      46             :     integer                          :: lh,n
      47             :     complex                           :: sigma_loc(2)
      48             : 
      49           6 :     if (sym%invs) then
      50             :       !This is for easier debugging of the preconditioner. The imaginary part
      51             :       !of the output density in the interstitial is never constrained to be
      52             :       !0 in the case of inversion symmetric systems. This leads to all numerical
      53             :       !noise leaking through making comparisons especially with different parallelizations
      54             :       !more difficult. The input density is implicitly constrained since the
      55             :       !mixvector%from_density subroutine throws away the imaginary part if sym%invs is .true.
      56       15366 :       outDen%pw(:,:input%jspins) = real(outDen%pw(:,:input%jspins))
      57             :     endif
      58             : 
      59           6 :     CALL resDen%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
      60           6 :     CALL vYukawa%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_POTYUK )
      61           6 :     MPI0_b: IF( fmpi%irank == 0 ) THEN
      62           3 :        CALL resDen%subPotDen( outDen, inDen )
      63           3 :        IF( input%jspins == 2 ) CALL resDen%SpinsToChargeAndMagnetisation()
      64             :     END IF MPI0_b
      65           6 :     CALL resDen%distribute(fmpi%mpi_comm)
      66           6 :     IF ( .NOT. input%film ) THEN
      67           6 :        sigma_loc = cmplx(0.0,0.0)
      68             :        CALL vgen_coulomb( 1, fmpi,    input, field, vacuum, sym, juphon, stars, cell, &
      69           6 :             sphhar, atoms, .FALSE., resDen, vYukawa, sigma_loc )
      70             :     ELSE
      71           0 :        call resDenMod%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
      72           0 :        if( fmpi%irank == 0 ) then
      73           0 :           call resDenMod%copyPotDen( resDen )
      74             :        end if
      75           0 :        CALL resDenMod%distribute(fmpi%mpi_comm)
      76           0 :        vYukawa%iter = resDen%iter
      77             :        CALL VYukawaFilm( stars, vacuum, cell, sym, juphon, input, fmpi, atoms, sphhar,   noco, nococonv,resDenMod, &
      78           0 :             vYukawa )
      79             :     END IF
      80             : 
      81           6 :     MPI0_c: IF( fmpi%irank == 0 ) THEN
      82        3840 :        resDen%pw(1:stars%ng3,1) = resDen%pw(1:stars%ng3,1) - input%preconditioning_param ** 2 / fpi_const * vYukawa%pw(1:stars%ng3,1)
      83          18 :        DO n = 1, atoms%ntype
      84         225 :           DO lh = 0, sphhar%nlh(sym%ntypsy(atoms%firstAtom(n)))
      85             :              resDen%mt(:atoms%jri(n),lh,n,1) = resDen%mt(:atoms%jri(n),lh,n,1) &
      86             :                   - input%preconditioning_param ** 2 / fpi_const &
      87      156921 :                   * vYukawa%mt(:atoms%jri(n),lh,n,1) * atoms%rmsh(:atoms%jri(n),n) ** 2
      88             :           END DO
      89             :        END DO
      90           3 :        IF (input%film) THEN
      91           0 :           resDen%vac(:,1,:,:) = resDen%vac(:,1,:,:) - input%preconditioning_param ** 2 / fpi_const * REAL(vYukawa%vac(:,1,:,:)) ! TODO: AN TB; REAL to COMPLEX OK?
      92           0 :           resDen%vac(:vacuum%nmzxyd,2:,:,:) = resDen%vac(:vacuum%nmzxyd,2:,:,:) - input%preconditioning_param ** 2 / fpi_const * vYukawa%vac(:vacuum%nmzxyd,2:,:,:)
      93             :        END IF
      94           3 :        IF( input%jspins == 2 ) CALL resDen%ChargeAndMagnetisationToSpins()
      95             :        ! fix the preconditioned density
      96           3 :        CALL outDen%addPotDen( resDen, inDen )
      97           3 :        CALL qfix(fmpi,stars,nococonv, atoms, sym, vacuum, sphhar, input, cell,   outDen, noco%l_noco, .FALSE., l_par=.FALSE., force_fix=.TRUE., fix=fix )
      98           3 :        CALL resDen%subPotDen( outDen, inDen )
      99         354 :        resDen%mmpMat = outDen%mmpMat - inDen%mmpMat
     100             :     END IF MPI0_c
     101           6 :     CALL precon_v%from_density(resden,vacuum%nmzxyd)
     102             : 
     103           6 :   END SUBROUTINE kerker
     104             : 
     105             : END MODULE m_kerker