mix.F90 14.4 KB
Newer Older
1 2 3 4 5
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
6
module m_mix
Gregor Michalicek's avatar
Gregor Michalicek committed
7

8
  !------------------------------------------------------------------------
9
  !  mixing of charge densities or potentials:
10
  !    IMIX = 0 : linear mixing                                     
11 12 13 14
  !    IMIX = 3 : Broyden's First method                            
  !    IMIX = 5 : Broyden's Second method                           
  !    IMIX = 7 : Generalized Anderson method                       
  !------------------------------------------------------------------------
Gregor Michalicek's avatar
Gregor Michalicek committed
15

16
contains
Gregor Michalicek's avatar
Gregor Michalicek committed
17

18
  subroutine mix( field, xcpot, dimension, obsolete, sliceplot, mpi, &
19 20
                stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
                oneD, hybrid, archiveType, inDen, outDen, results )
Gregor Michalicek's avatar
Gregor Michalicek committed
21

22
#include"cpp_double.h"
Gregor Michalicek's avatar
Gregor Michalicek committed
23

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
    use m_juDFT
    use m_constants
    use m_cdn_io
    use m_broyd_io
    use m_brysh1
    use m_stmix
    use m_broyden
    use m_broyden2
    use m_brysh2
    use m_metric
    use m_qfix
    use m_types
    use m_xmlOutput
    use m_umix
    use m_vgen_coulomb
39
    use m_VYukawaFilm
40 41 42
#ifdef CPP_MPI
    use m_mpi_bc_potden
#endif
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
    implicit none

    type(t_oneD),      intent(in)    :: oneD
    type(t_hybrid),    intent(in)    :: hybrid 
    type(t_input),     intent(in)    :: input
    type(t_vacuum),    intent(in)    :: vacuum
    type(t_noco),      intent(in)    :: noco
    type(t_sym),       intent(in)    :: sym
    type(t_stars),     intent(in)    :: stars
    type(t_cell),      intent(in)    :: cell
    type(t_sphhar),    intent(in)    :: sphhar
    type(t_field),     intent(inout) :: field
    class(t_xcpot),    intent(in)    :: xcpot
    type(t_dimension), intent(in)    :: dimension
    type(t_obsolete),  intent(in)    :: obsolete
    type(t_sliceplot), intent(in)    :: sliceplot
    type(t_mpi),       intent(in)    :: mpi
    type(t_atoms),     intent(inout) :: atoms !n_u is modified temporarily
    type(t_potden),    intent(inout) :: outDen
    type(t_results),   intent(inout) :: results
    type(t_potden),    intent(inout) :: inDen
    integer,           intent(in)    :: archiveType

    real                             :: fix, intfac, vacfac
    integer                          :: i, imap, js, n, lh
    integer                          :: mmap, mmaph, nmaph, nmap, mapmt, mapvac, mapvac2
    integer                          :: iofl, n_u_keep
    logical                          :: l_exist, l_ldaU, l_densityMatrixPresent, l_pot
    real                             :: dist(6)
    real, allocatable                :: sm(:), fsm(:), fmMet(:), smMet(:)
    character(len=20)                :: attributes(2)
    complex                          :: n_mmpTemp(-3:3,-3:3,max(1,atoms%n_u),input%jspins)
    type(t_potden)                   :: resDen, vYukawa
76
    integer                          :: ierr(2)
77 78 79 80 81 82 83 84 85 86 87 88

    !External functions
    real :: CPP_BLAS_sdot
    external :: CPP_BLAS_sdot

    ! YM: I have exported 'vol' from outside, be aware
    !     IF (film) THEN
    !        vol = 2.0 * z1 * area
    !     ELSE
    !        vol = omtil
    !     ENDIF

89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
    MPI0_a: if( mpi%irank == 0 ) then

      l_densityMatrixPresent = any( inDen%mmpMat(:,:,:,:) /= 0.0 )

      !In systems without inversions symmetry the interstitial star-
      !coefficients are complex. Thus twice as many numbers have to be
      !stored.
      intfac = 2.0
      if ( sym%invs ) intfac = 1.0

      !The corresponding is true for the coeff. of the warping vacuum
      !density depending on the two dimensional inversion.
      vacfac = 2.0
      if ( sym%invs2 ) vacfac = 1.0

      mmaph = intfac * stars%ng3 + atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd + &
              vacfac * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + vacuum%nmzd * vacuum%nvac
      mmap  = mmaph * input%jspins
      !in a non-collinear calculations extra space is needed for the
      !off-diag. part of the density matrix. these coeff. are generally
      !complex independ of invs and invs2.
      if ( noco%l_noco ) then
         mmap = mmap + 2 * stars%ng3 + 2 * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + &
                2 * vacuum%nmzd * vacuum%nvac
      end if
114

115 116 117 118 119 120 121 122 123 124 125 126
      ! LDA+U (start)
      n_mmpTemp = inDen%mmpMat
      n_u_keep = atoms%n_u
      if ( atoms%n_u > 0 ) call u_mix( input, atoms, inDen%mmpMat, outDen%mmpMat )
      if ( l_densityMatrixPresent ) then
        !In an LDA+U caclulation, also the density matrix is included in the
        !supervectors (sm,fsm) if no linear mixing is performed on it.
        if ( input%ldauLinMix ) then
          atoms%n_u = 0
        else
          mmap = mmap + 7 * 7 * 2 * atoms%n_u * input%jspins ! add 7*7 complex numbers per atoms%n_u and spin
        end if
127
      else
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
        atoms%n_u = 0
      end if
      ! LDA+U (end)

      allocate( sm(mmap), fsm(mmap) )
      allocate( smMet(mmap), fmMet(mmap) )
      dist(:) = 0.0

      !determine type of mixing:
      !imix=0:straight, imix=o broyden first, imix=5:broyden second
      !imix=:generalozed anderson mixing
      select case( input%imix )
        case( 0 )
          write( 16, fmt='(a,2f10.5)' ) 'STRAIGHT MIXING',input%alpha
        case( 3 )
          write( 16, fmt='(a,f10.5)' ) 'BROYDEN FIRST MIXING',input%alpha
        case( 5 )
          write( 16, fmt='(a,f10.5)' ) 'BROYDEN SECOND MIXING',input%alpha
        case( 7 )
          write( 16, fmt='(a,f10.5)' ) 'ANDERSON GENERALIZED',input%alpha
        case default
          call juDFT_error( "mix: input%imix =/= 0,3,5,7 ", calledby ="mix" )
      end select

      if ( input%jspins == 2 .and. input%imix /= 0 ) then
        write( 6, '(''WARNING : for QUASI-NEWTON METHODS SPINF=1'')' )
154 155
      end if

156 157 158 159
      !put input charge density into array sm
      !(in the spin polarized case the arrays sm and fsm consist of spin up and spin down densities)
      call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
                   intfac, vacfac, inDen, nmap, nmaph, mapmt, mapvac, mapvac2, sm ) 
160

161 162 163
      !put output charge density into array fsm
      call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
                   intfac, vacfac, outDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
164

165 166
      !store the difference fsm - sm in fsm
      fsm(:nmap) = fsm(:nmap) - sm(:nmap)
167

168 169 170 171
      l_pot = .false.
      ! Apply metric w to fsm and store in fmMet:  w |fsm>
      call metric( cell, atoms, vacuum, sphhar, input, noco, stars, sym, oneD, &
                   mmap, nmaph, mapmt, mapvac2, fsm, fmMet, l_pot )
172

173
      !calculate the distance of charge densities for each spin
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
      IF(hybrid%l_calhf) THEN
         CALL openXMLElement('densityConvergence',(/'units  ','comment'/),(/'me/bohr^3','HF       '/))
      ELSE
         CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
      END IF

      DO js = 1,input%jspins
         dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, fmMet(nmaph*(js-1)+1),1)

         attributes = ''
         WRITE(attributes(1),'(i0)') js
         WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/cell%vol))
         CALL writeXMLElementForm('chargeDensity',(/'spin    ','distance'/),attributes,reshape((/4,8,1,20/),(/2,2/)))
         IF( hybrid%l_calhf ) THEN
            WRITE (16,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
            WRITE ( 6,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
         ELSE
            WRITE (16,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
            WRITE ( 6,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
         END IF
      END DO
      IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,fmMet(nmaph*2+1),1)
      IF (noco%l_noco) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/cell%vol))

      !calculate the distance of total charge and spin density
      !|rho/m(o) - rho/m(i)| = |rh1(o) -rh1(i)|+ |rh2(o) -rh2(i)| +/_
      !                        +/_2<rh2(o) -rh2(i)|rh1(o) -rh1(i)>
      IF (input%jspins.EQ.2) THEN
         dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,fmMet(nmaph+1),1)
         dist(4) = dist(1) + dist(2) + 2.0e0*dist(3)
         dist(5) = dist(1) + dist(2) - 2.0e0*dist(3)
         CALL writeXMLElementFormPoly('overallChargeDensity',(/'distance'/),&
                                      (/1000*SQRT(ABS(dist(4)/cell%vol))/),reshape((/10,20/),(/1,2/)))
         CALL writeXMLElementFormPoly('spinDensity',(/'distance'/),&
                                      (/1000*SQRT(ABS(dist(5)/cell%vol))/),reshape((/19,20/),(/1,2/)))
         IF( hybrid%l_calhf ) THEN
            WRITE (16,FMT=8001) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
            WRITE (16,FMT=8011) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
            WRITE ( 6,FMT=8001) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
            WRITE ( 6,FMT=8011) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
         ELSE
            WRITE (16,FMT=8000) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
            WRITE (16,FMT=8010) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
            WRITE ( 6,FMT=8000) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
            WRITE ( 6,FMT=8010) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
         END IF

         !dist/vol should always be >= 0 ,
         !but for dist=0 numerically you might obtain dist/vol < 0
         !(e.g. when calculating non-magnetic systems with jspins=2).
      END IF
      results%last_distance=maxval(1000*SQRT(ABS(dist/cell%vol)))
226
      DEALLOCATE (smMet,fmMet)
227 228
      CALL closeXMLElement('densityConvergence')

229 230 231 232 233 234 235 236 237 238 239 240 241
    end if MPI0_a

    ! KERKER PRECONDITIONER
    if( input%preconditioning_param /= 0 ) then
      call resDen%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 1001 )
      call vYukawa%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 4 )
      MPI0_b: if( mpi%irank == 0 ) then 
        call resDen%Residual( outDen, inDen )
        if( input%jspins == 2 ) call resDen%SpinsToChargeAndMagnetisation()
      end if MPI0_b
#ifdef CPP_MPI
      call mpi_bc_potden( mpi, stars, sphhar, atoms, input, vacuum, oneD, noco, resDen )
#endif
242
      if ( .not. input%film ) then
243 244
        call vgen_coulomb( 1, mpi, dimension, oneD, input, field, vacuum, sym, stars, cell, &
                           sphhar, atoms, resDen, vYukawa )
245 246 247 248 249
      else
        vYukawa%iter = resDen%iter
        call VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, resDen, &
                         vYukawa )
      end if
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290
    end if
    MPI0_c: if( mpi%irank == 0 ) then
      if( input%preconditioning_param /= 0 ) then
        resDen%pw(1:stars%ng3,1) = resDen%pw(1:stars%ng3,1) - input%preconditioning_param ** 2 / fpi_const * vYukawa%pw(1:stars%ng3,1)
        do n = 1, atoms%ntype
          do lh = 0, sphhar%nlhd
            resDen%mt(1:atoms%jri(n),lh,n,1) = resDen%mt(1:atoms%jri(n),lh,n,1) &
                    - input%preconditioning_param ** 2 / fpi_const &
                    * vYukawa%mt(1:atoms%jri(n),lh,n,1) * atoms%rmsh(1:atoms%jri(n),n) ** 2
          end do
        end do
        if( input%jspins == 2 ) call resDen%ChargeAndMagnetisationToSpins()
        call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
                     intfac, vacfac, resDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
      end if
    ! end of preconditioner


    !mixing of the densities
      IF (input%imix.EQ.0) THEN
         CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
      ELSE
         CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
                      hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)

!        Replace the broyden call above by the commented metric and broyden2 calls
!        below to switch on the continuous restart of the Broyden method.
         ! Apply metric w to sm and store in smMet:  w |sm>
!         CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                     mmap,nmaph,mapmt,mapvac2,sm,smMet,l_pot)
!  
!         CALL broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
!                       hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm,fmMet,smMet)
      END IF

      !initiatlize mixed density and extract it with brysh2 call
      inDen%mmpMat = CMPLX(0.0,0.0)

      CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen) 
      DEALLOCATE (sm,fsm)

291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
      !fix charge of the new density
      CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false., fix)

      IF(atoms%n_u.NE.n_u_keep) THEN
         inDen%mmpMat = n_mmpTemp
      END IF

      atoms%n_u=n_u_keep

      IF(vacuum%nvac.EQ.1) THEN
         inDen%vacz(:,2,:) = inDen%vacz(:,1,:)
         IF (sym%invs) THEN
            inDen%vacxy(:,:,2,:) = CONJG(inDen%vacxy(:,:,1,:))
         ELSE
            inDen%vacxy(:,:,2,:) = inDen%vacxy(:,:,1,:)
         END IF
      END IF

      IF (atoms%n_u > 0) THEN
         IF (.NOT.l_densityMatrixPresent) THEN
            inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
            CALL resetBroydenHistory()
         END IF
      ENDIF

      !write out mixed density
      CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
                        1,results%last_distance,results%ef,.TRUE.,inDen)

320 321 322 323 324 325 326 327 328
#ifdef CPP_HDF
      IF (judft_was_argument("-last_extra")) THEN
         CALL system("rm cdn_last.hdf")
         CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
                           1,results%last_distance,results%ef,.TRUE.,inDen,'cdn_last')

      END IF
#endif

329 330 331 332 333 334 335 336 337 338 339
      inDen%iter = inDen%iter + 1

      7900 FORMAT (/,'---->    distance of charge densities for spin ',i2,'                 it=',i5,':',f13.6,' me/bohr**3')
      7901 FORMAT (/,'----> HF distance of charge densities for spin ',i2,'                 it=',i5,':',f13.6,' me/bohr**3')
      8000 FORMAT (/,'---->    distance of charge densities for it=',i5,':', f13.6,' me/bohr**3')
      8001 FORMAT (/,'----> HF distance of charge densities for it=',i5,':', f13.6,' me/bohr**3')
      8010 FORMAT (/,'---->    distance of spin densities for it=',i5,':', f13.6,' me/bohr**3')
      8011 FORMAT (/,'----> HF distance of spin densities for it=',i5,':', f13.6,' me/bohr**3')
      8020 FORMAT (4d25.14)
      8030 FORMAT (10i10)

340
    end if MPI0_c
341 342 343 344

  end subroutine mix

end module m_mix