mix.F90 14.1 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 40 41
#ifdef CPP_MPI
    use m_mpi_bc_potden
#endif
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
    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
75
    integer                          :: ierr(2)
76 77 78 79 80 81 82 83 84 85 86 87

    !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

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
    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
113

114 115 116 117 118 119 120 121 122 123 124 125
      ! 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
126
      else
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
        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'')' )
153 154
      end if

155 156 157 158
      !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 ) 
159

160 161 162
      !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 )
163

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

167 168 169 170
      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 )
171

172
    end if MPI0_a
173 174 175 176 177

    ! 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 )
178
      MPI0_b: if( mpi%irank == 0 ) then 
179 180
        call resDen%Residual( outDen, inDen )
        if( input%jspins == 2 ) call resDen%SpinsToChargeAndMagnetisation()
181 182 183 184
      end if MPI0_b
#ifdef CPP_MPI
      call mpi_bc_potden( mpi, stars, sphhar, atoms, input, vacuum, oneD, noco, resDen )
#endif
185
      call vgen_coulomb( 1, mpi, dimension, oneD, input, field, vacuum, sym, stars, cell, &
186
                         sphhar, atoms, resDen, vYukawa )
187
    end if
188
    MPI0_c: if( mpi%irank == 0 ) then
189
      if( input%preconditioning_param /= 0 ) then
190 191 192 193 194 195 196 197 198
        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()
199 200 201 202 203
        call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
                     intfac, vacfac, resDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
      end if
    ! end of preconditioner

204 205

    !mixing of the densities
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 319 320 321 322 323 324 325 326 327 328
      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) 

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

329
    end if MPI0_c
330 331 332 333

  end subroutine mix

end module m_mix