mix.F90 11.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 10 11 12 13
  !  mixing of charge densities or potentials:
  !    IMIX= 0 : linear mixing                                     
  !    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

SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
19
               hybrid,archiveType,inDen,outDen,results)
Gregor Michalicek's avatar
Gregor Michalicek committed
20

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

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

   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_atoms),INTENT(INOUT)   :: atoms !n_u is modified temporarily
50
   TYPE(t_potden),INTENT(INOUT)  :: outDen
Gregor Michalicek's avatar
Gregor Michalicek committed
51
   TYPE(t_results),INTENT(INOUT) :: results
52
   TYPE(t_potden),INTENT(INOUT)  :: inDen
53
   INTEGER, INTENT(IN)           :: archiveType
Gregor Michalicek's avatar
Gregor Michalicek committed
54 55 56

   !Local Scalars
   REAL fix,intfac,vacfac
57
   INTEGER i,imap,js
Gregor Michalicek's avatar
Gregor Michalicek committed
58
   INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
59
   INTEGER iofl,n_u_keep
60
   LOGICAL l_exist,l_ldaU, l_densityMatrixPresent, l_pot
Gregor Michalicek's avatar
Gregor Michalicek committed
61 62 63

   !Local Arrays
   REAL dist(6)
64
   REAL, ALLOCATABLE :: sm(:), fsm(:), fmMet(:), smMet(:)
Gregor Michalicek's avatar
Gregor Michalicek committed
65
   CHARACTER(LEN=20) :: attributes(2)
66
   COMPLEX           :: n_mmpTemp(-3:3,-3:3,MAX(1,atoms%n_u),input%jspins)
Gregor Michalicek's avatar
Gregor Michalicek committed
67 68 69 70 71 72 73 74 75 76 77 78

   !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

79 80
   l_densityMatrixPresent = ANY(inDen%mmpMat(:,:,:,:).NE.0.0)

Gregor Michalicek's avatar
Gregor Michalicek committed
81 82 83
   !In systems without inversions symmetry the interstitial star-
   !coefficients are complex. Thus twice as many numbers have to be
   !stored.
84 85
   intfac = 2.0
   IF (sym%invs) intfac = 1.0
Gregor Michalicek's avatar
Gregor Michalicek committed
86 87 88

   !The corresponding is true for the coeff. of the warping vacuum
   !density depending on the two dimensional inversion.
89 90
   vacfac = 2.0
   IF (sym%invs2) vacfac = 1.0
Gregor Michalicek's avatar
Gregor Michalicek committed
91 92 93 94 95 96 97 98 99 100 101 102

   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

103
   ! LDA+U (start)
104
   n_mmpTemp = inDen%mmpMat
Gregor Michalicek's avatar
Gregor Michalicek committed
105
   n_u_keep=atoms%n_u
106
   IF (atoms%n_u.GT.0) CALL u_mix(input,atoms,inDen%mmpMat,outDen%mmpMat)
107
   IF (l_densityMatrixPresent) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
108
      !In an LDA+U caclulation, also the density matrix is included in the
109
      !supervectors (sm,fsm) if no linear mixing is performed on it.
110
      IF (input%ldauLinMix) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
111 112
         atoms%n_u = 0
      ELSE
113
         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
114 115
      END IF
   ELSE
116 117 118
      atoms%n_u = 0
   END IF
   ! LDA+U (end)
Gregor Michalicek's avatar
Gregor Michalicek committed
119 120

   ALLOCATE (sm(mmap),fsm(mmap))
121
   ALLOCATE (smMet(mmap),fmMet(mmap))
122
   dist(:) = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142

   !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

143
   !put input charge density into array sm
144

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

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

153 154
   !store the difference fsm - sm in fsm
   fsm(:nmap) = fsm(:nmap) - sm(:nmap)
Gregor Michalicek's avatar
Gregor Michalicek committed
155

156 157 158 159 160
   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
161 162 163 164
   !mixing of the densities
   IF (input%imix.EQ.0) THEN
      CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
   ELSE
165 166
      CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
                   hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
167 168 169

!     Replace the broyden call above by the commented metric and broyden2 calls
!     below to switch on the continuous restart of the Broyden method.
170
      ! Apply metric w to sm and store in smMet:  w |sm>
171 172 173 174 175
!      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
176 177 178
   END IF

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

181
   CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen) 
Gregor Michalicek's avatar
Gregor Michalicek committed
182 183 184 185

   !calculate the distance of charge densities...

   !induce metric in fsm use sm as an output array: |sm> = w |fsm>
186 187
!   CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!               mmap,nmaph,mapmt,mapvac2,fsm, sm)
Gregor Michalicek's avatar
Gregor Michalicek committed
188 189 190 191 192 193 194

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

Gregor Michalicek's avatar
Gregor Michalicek committed
196
   DO js = 1,input%jspins
197
      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
198 199 200 201 202 203

      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
204 205
         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
206
      ELSE
207 208
         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
209 210
      END IF
   END DO
211
   IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,fmMet(nmaph*2+1),1)
212
   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
213 214 215 216 217

   !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
218
      dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,fmMet(nmaph+1),1)
Gregor Michalicek's avatar
Gregor Michalicek committed
219 220 221 222 223 224 225
      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
226 227 228 229
         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
230
      ELSE
231 232 233 234
         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
235 236 237 238 239 240 241
      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)))
242
   DEALLOCATE (sm,fsm,smMet,fmMet)
Gregor Michalicek's avatar
Gregor Michalicek committed
243 244 245
   CALL closeXMLElement('densityConvergence')

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

248 249 250 251 252 253
   IF(atoms%n_u.NE.n_u_keep) THEN
      inDen%mmpMat = n_mmpTemp
   END IF

   atoms%n_u=n_u_keep

254 255 256 257 258 259 260 261 262
   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
263
   IF (atoms%n_u > 0) THEN
264 265
      IF (.NOT.l_densityMatrixPresent) THEN
         inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
266
         CALL resetBroydenHistory()
267
      END IF
Gregor Michalicek's avatar
Gregor Michalicek committed
268 269
   ENDIF

270 271 272 273
   !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)

274 275
   inDen%iter = inDen%iter + 1

Gregor Michalicek's avatar
Gregor Michalicek committed
276 277 278 279 280 281 282 283 284 285 286
   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

287
END MODULE m_mix