mix.F90 12.7 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
  !    IMIX = 3 : BROYDEN'S FIRST METHOD                            
  !    IMIX = 5 : BROYDEN'S SECOND METHOD                           
  !    IMIX = 7 : GENERALIZED ANDERSEN METHOD                       
14
  !************************************************************************
Gregor Michalicek's avatar
Gregor Michalicek committed
15

16
CONTAINS
Gregor Michalicek's avatar
Gregor Michalicek committed
17

18 19 20
SUBROUTINE mix( field, xcpot, dimension, obsolete, sliceplot, mpi, &
                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

   USE m_juDFT
   USE m_constants
   USE m_cdn_io
27
   USE m_broyd_io
Gregor Michalicek's avatar
Gregor Michalicek committed
28 29 30
   USE m_brysh1
   USE m_stmix
   USE m_broyden
31
   USE m_broyden2
Gregor Michalicek's avatar
Gregor Michalicek committed
32 33 34 35 36
   USE m_brysh2
   USE m_metric
   USE m_qfix
   USE m_types
   USE m_xmlOutput
37
   USE m_umix
38
   use m_vgen
Gregor Michalicek's avatar
Gregor Michalicek committed
39 40 41 42 43 44 45 46 47 48 49 50

   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
51 52 53 54 55 56
   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
Gregor Michalicek's avatar
Gregor Michalicek committed
57
   TYPE(t_atoms),INTENT(INOUT)   :: atoms !n_u is modified temporarily
58
   TYPE(t_potden),INTENT(INOUT)  :: outDen
Gregor Michalicek's avatar
Gregor Michalicek committed
59
   TYPE(t_results),INTENT(INOUT) :: results
60
   TYPE(t_potden),INTENT(INOUT)  :: inDen
61
   INTEGER, INTENT(IN)           :: archiveType
Gregor Michalicek's avatar
Gregor Michalicek committed
62 63 64

   !Local Scalars
   REAL fix,intfac,vacfac
65
   INTEGER i,imap,js, n, lh
Gregor Michalicek's avatar
Gregor Michalicek committed
66
   INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
67
   INTEGER iofl,n_u_keep
68
   LOGICAL l_exist,l_ldaU, l_densityMatrixPresent, l_pot
Gregor Michalicek's avatar
Gregor Michalicek committed
69 70 71

   !Local Arrays
   REAL dist(6)
72
   REAL, ALLOCATABLE :: sm(:), fsm(:), fmMet(:), smMet(:)
Gregor Michalicek's avatar
Gregor Michalicek committed
73
   CHARACTER(LEN=20) :: attributes(2)
74
   COMPLEX           :: n_mmpTemp(-3:3,-3:3,MAX(1,atoms%n_u),input%jspins)
Gregor Michalicek's avatar
Gregor Michalicek committed
75

76 77
   type(t_potden) :: resDen, resChM, vYukawa, vTot, vx

Gregor Michalicek's avatar
Gregor Michalicek committed
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
   l_densityMatrixPresent = ANY(inDen%mmpMat(:,:,:,:).NE.0.0)

Gregor Michalicek's avatar
Gregor Michalicek committed
91 92 93
   !In systems without inversions symmetry the interstitial star-
   !coefficients are complex. Thus twice as many numbers have to be
   !stored.
94 95
   intfac = 2.0
   IF (sym%invs) intfac = 1.0
Gregor Michalicek's avatar
Gregor Michalicek committed
96 97 98

   !The corresponding is true for the coeff. of the warping vacuum
   !density depending on the two dimensional inversion.
99 100
   vacfac = 2.0
   IF (sym%invs2) vacfac = 1.0
Gregor Michalicek's avatar
Gregor Michalicek committed
101 102 103 104 105 106 107 108 109 110 111 112

   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
   ! LDA+U (start)
114
   n_mmpTemp = inDen%mmpMat
Gregor Michalicek's avatar
Gregor Michalicek committed
115
   n_u_keep=atoms%n_u
116
   IF (atoms%n_u.GT.0) CALL u_mix(input,atoms,inDen%mmpMat,outDen%mmpMat)
117
   IF (l_densityMatrixPresent) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
118
      !In an LDA+U caclulation, also the density matrix is included in the
119
      !supervectors (sm,fsm) if no linear mixing is performed on it.
120
      IF (input%ldauLinMix) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
121 122
         atoms%n_u = 0
      ELSE
123
         mmap = mmap + 7 * 7 * 2 * atoms%n_u * input%jspins ! add 7*7 complex numbers per atoms%n_u and spin
Gregor Michalicek's avatar
Gregor Michalicek committed
124 125
      END IF
   ELSE
126 127 128
      atoms%n_u = 0
   END IF
   ! LDA+U (end)
Gregor Michalicek's avatar
Gregor Michalicek committed
129 130

   ALLOCATE (sm(mmap),fsm(mmap))
131
   ALLOCATE (smMet(mmap),fmMet(mmap))
132
   dist(:) = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

   !determine type of mixing:
   !imix=0:straight, imix=o broyden first, imix=5:broyden second
   !imix=:generalozed anderson mixing
   IF (input%imix.EQ.0) THEN
      WRITE (16,FMT='(a,2f10.5)') 'STRAIGHT MIXING',input%alpha
   ELSE IF (input%imix.EQ.3) THEN
      WRITE (16,FMT='(a,f10.5)') 'BROYDEN FIRST MIXING',input%alpha
   ELSE IF (input%imix.EQ.5) THEN
      WRITE (16,FMT='(a,f10.5)') 'BROYDEN SECOND MIXING',input%alpha
   ELSE IF (input%imix.EQ.7) THEN
      WRITE (16,FMT='(a,f10.5)') 'ANDERSON GENERALIZED',input%alpha
   ELSE
      CALL juDFT_error("mix: input%imix =/= 0,3,5,7 ",calledby ="mix")
   END IF

   IF (input%jspins.EQ.2.AND.input%imix.NE.0) THEN
      WRITE(6,'(''WARNING : for QUASI-NEWTON METHODS SPINF=1'')')
   END IF

153 154
   !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)
Gregor Michalicek's avatar
Gregor Michalicek committed
155
   CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
156
               intfac,vacfac,inDen,nmap,nmaph,mapmt,mapvac,mapvac2,sm) 
Gregor Michalicek's avatar
Gregor Michalicek committed
157

158
   !put output charge density into array fsm
Gregor Michalicek's avatar
Gregor Michalicek committed
159
   CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
160
               intfac,vacfac,outDen,nmap,nmaph,mapmt,mapvac,mapvac2,fsm)
Gregor Michalicek's avatar
Gregor Michalicek committed
161

162 163
   !store the difference fsm - sm in fsm
   fsm(:nmap) = fsm(:nmap) - sm(:nmap)
Gregor Michalicek's avatar
Gregor Michalicek committed
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
   ! preconditioning section
   if( input%preconditioning_param /= 0 ) then
     resDen = inDen 
     call brysh2( input, stars, atoms, sphhar, noco, vacuum, sym, fsm, oneD, resDen )   
     if( input%jspins == 2 ) call resDen%SpinsToChargeAndMagnetisation( resChM )
     call vYukawa%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 1 )
     call vTot%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 1 )
     call vx%init( stars, atoms, sphhar, vacuum, input%jspins, .false., 1 )
     call vgen( hybrid, field, input, xcpot, DIMENSION, atoms, sphhar, stars, vacuum, &
              sym, obsolete, cell, oneD, sliceplot, mpi, results, noco, .false., inDen, vTot, vx, &
              vYukawa )
     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
     end if
   ! end of preconditioning section

   if( input%jspins == 2 ) call resChM%ChargeAndMagnetisationToSpins( resDen )

189 190 191 192 193
   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)

Gregor Michalicek's avatar
Gregor Michalicek committed
194 195 196 197
   !mixing of the densities
   IF (input%imix.EQ.0) THEN
      CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
   ELSE
198 199
      CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
                   hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
200 201 202

!     Replace the broyden call above by the commented metric and broyden2 calls
!     below to switch on the continuous restart of the Broyden method.
203
      ! Apply metric w to sm and store in smMet:  w |sm>
204 205 206 207 208
!      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)
Gregor Michalicek's avatar
Gregor Michalicek committed
209 210 211
   END IF

   !initiatlize mixed density and extract it with brysh2 call
212
   inDen%mmpMat = CMPLX(0.0,0.0)
Gregor Michalicek's avatar
Gregor Michalicek committed
213

214
   CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen) 
Gregor Michalicek's avatar
Gregor Michalicek committed
215 216 217 218

   !calculate the distance of charge densities...

   !induce metric in fsm use sm as an output array: |sm> = w |fsm>
219 220
!   CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!               mmap,nmaph,mapmt,mapvac2,fsm, sm)
Gregor Michalicek's avatar
Gregor Michalicek committed
221 222 223 224 225 226 227

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

Gregor Michalicek's avatar
Gregor Michalicek committed
229
   DO js = 1,input%jspins
230
      dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, fmMet(nmaph*(js-1)+1),1)
Gregor Michalicek's avatar
Gregor Michalicek committed
231 232 233 234 235 236

      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
237 238
         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))
Gregor Michalicek's avatar
Gregor Michalicek committed
239
      ELSE
240 241
         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))
Gregor Michalicek's avatar
Gregor Michalicek committed
242 243
      END IF
   END DO
244
   IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,fmMet(nmaph*2+1),1)
245
   IF (noco%l_noco) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/cell%vol))
Gregor Michalicek's avatar
Gregor Michalicek committed
246 247 248 249 250

   !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
251
      dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,fmMet(nmaph+1),1)
Gregor Michalicek's avatar
Gregor Michalicek committed
252 253 254 255 256 257 258
      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
259 260 261 262
         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))
Gregor Michalicek's avatar
Gregor Michalicek committed
263
      ELSE
264 265 266 267
         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))
Gregor Michalicek's avatar
Gregor Michalicek committed
268 269 270 271 272 273 274
      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)))
275
   DEALLOCATE (sm,fsm,smMet,fmMet)
Gregor Michalicek's avatar
Gregor Michalicek committed
276 277 278
   CALL closeXMLElement('densityConvergence')

   !fix charge of the new density
279
   CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false., fix)
Gregor Michalicek's avatar
Gregor Michalicek committed
280

281 282 283 284 285 286
   IF(atoms%n_u.NE.n_u_keep) THEN
      inDen%mmpMat = n_mmpTemp
   END IF

   atoms%n_u=n_u_keep

287 288 289 290 291 292 293 294 295
   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

Gregor Michalicek's avatar
Gregor Michalicek committed
296
   IF (atoms%n_u > 0) THEN
297 298
      IF (.NOT.l_densityMatrixPresent) THEN
         inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
299
         CALL resetBroydenHistory()
300
      END IF
Gregor Michalicek's avatar
Gregor Michalicek committed
301 302
   ENDIF

303 304 305 306
   !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)

307 308
   inDen%iter = inDen%iter + 1

Gregor Michalicek's avatar
Gregor Michalicek committed
309 310 311 312 313 314 315 316 317 318 319
   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

320
END MODULE m_mix