mix.F90 13.5 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 39 40 41 42 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 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
    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
    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

    !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

    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

    ! 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
Gregor Michalicek's avatar
Gregor Michalicek committed
113
      !In an LDA+U caclulation, also the density matrix is included in the
114
      !supervectors (sm,fsm) if no linear mixing is performed on it.
115 116 117 118 119 120
      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
    else
121
      atoms%n_u = 0
122 123 124 125 126 127 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 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 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
    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'')' )
    end if

    !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 ) 

    !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 )

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

    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 )

    call resDen%Residual( outDen, inDen )

    ! KERKER PRECONDITIONER
    if( input%preconditioning_param /= 0 ) then
      !resDen = inDen ! just for the general setup - values are next: 
      call resDen%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 1001 )
      call resDen%Residual( outDen, inDen )
      !call brysh2( input, stars, atoms, sphhar, noco, vacuum, sym, fsm, oneD, resDen )   
      if( input%jspins == 2 ) call resDen%SpinsToChargeAndMagnetisation( )
      call vYukawa%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 4 )
      call vgen_coulomb( 1, mpi, dimension, oneD, input, field, vacuum, sym, stars, cell, &
                         sphhar, atoms, resDen, vYukawa, results )
      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 preconditioning section

    !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)
204
!
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 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 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
!       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) 

    !calculate the distance of charge densities...

    !induce metric in fsm use sm as an output array: |sm> = w |fsm>
!    CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                mmap,nmaph,mapmt,mapvac2,fsm, sm)

    !calculate the charge density distance for each spin
    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)))
    DEALLOCATE (sm,fsm,smMet,fmMet)
    CALL closeXMLElement('densityConvergence')

    !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)

    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)

  end subroutine mix

end module m_mix