mix.F90 11.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 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 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

   USE m_juDFT
   USE m_constants
   USE m_cdn_io
   USE m_brysh1
   USE m_stmix
   USE m_broyden
   USE m_brysh2
   USE m_metric
   USE m_qfix
   USE m_types
   USE m_xmlOutput

   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
47
   TYPE(t_potden),INTENT(IN)     :: outDen
Gregor Michalicek's avatar
Gregor Michalicek committed
48
   TYPE(t_results),INTENT(INOUT) :: results
49
   TYPE(t_potden),INTENT(INOUT)  :: inDen
50
   INTEGER, INTENT(IN)           :: archiveType
Gregor Michalicek's avatar
Gregor Michalicek committed
51 52 53

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

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

   !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

76 77
   l_densityMatrixPresent = ANY(inDen%mmpMat(:,:,:,:).NE.0.0)

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

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

   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

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

   ALLOCATE (sm(mmap),fsm(mmap))

118
   INQUIRE (file='broyd.'//CHAR(input%imix+48),exist=l_exist)
Gregor Michalicek's avatar
Gregor Michalicek committed
119 120 121 122 123 124 125 126 127 128 129
   DO i = 1,6
      dist(i) = 0.0
   END DO

   !determine type of mixing:
   !imix=0:straight, imix=o broyden first, imix=5:broyden second
   !imix=:generalozed anderson mixing
   mit = 0
   IF (input%imix.EQ.0) THEN
      WRITE (16,FMT='(a,2f10.5)') 'STRAIGHT MIXING',input%alpha
   ELSE IF (input%imix.EQ.3) THEN
130
      IF ( .NOT.l_exist) mit = 1
Gregor Michalicek's avatar
Gregor Michalicek committed
131 132
      WRITE (16,FMT='(a,f10.5)') 'BROYDEN FIRST MIXING',input%alpha
   ELSE IF (input%imix.EQ.5) THEN
133
      IF (.NOT.l_exist) mit = 1
Gregor Michalicek's avatar
Gregor Michalicek committed
134 135
      WRITE (16,FMT='(a,f10.5)') 'BROYDEN SECOND MIXING',input%alpha
   ELSE IF (input%imix.EQ.7) THEN
136
      IF (.NOT.l_exist) mit = 1
Gregor Michalicek's avatar
Gregor Michalicek committed
137 138 139 140 141 142 143 144 145
      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

146
   !put input charge density into array sm
147

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

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

156 157
   !store the difference fsm - sm in fsm
   fsm(:nmap) = fsm(:nmap) - sm(:nmap)
Gregor Michalicek's avatar
Gregor Michalicek committed
158 159 160 161 162 163 164

   ! open files for broyden
   irecl=(nmap+1)*8
   IF (input%imix.GE.3) THEN
      IF (hybrid%l_calhf) THEN
         OPEN (57,file='hf_broyd',form='unformatted',status='unknown')
         OPEN (59,file='hf_broyd.'//CHAR(input%imix+48),access='direct',&
165
               recl=irecl,form='unformatted',status='unknown')
Gregor Michalicek's avatar
Gregor Michalicek committed
166 167 168
      ELSE
         OPEN (57,file='broyd',form='unformatted',status='unknown')
         OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
169
               recl=irecl,form='unformatted',status='unknown')
Gregor Michalicek's avatar
Gregor Michalicek committed
170 171 172 173 174 175 176 177 178 179 180 181
      ENDIF
   END IF

   !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,&
                   mmap,nmaph,mapmt,mapvac2,nmap,fsm,mit,sm)
   END IF

   !initiatlize mixed density and extract it with brysh2 call
182 183 184 185
   inDen%cdom = CMPLX(0.0,0.0)
   inDen%cdomvz = CMPLX(0.0,0.0)
   inDen%cdomvxy = CMPLX(0.0,0.0)
   inDen%mmpMat = CMPLX(0.0,0.0)
Gregor Michalicek's avatar
Gregor Michalicek committed
186

187
   CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen) 
Gregor Michalicek's avatar
Gregor Michalicek committed
188 189 190 191 192 193 194 195 196 197 198 199 200

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

Gregor Michalicek's avatar
Gregor Michalicek committed
202 203 204 205 206 207 208 209
   DO js = 1,input%jspins
      dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, sm(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
210 211
         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
212
      ELSE
213 214
         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
215 216 217
      END IF
   END DO
   IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,sm(nmaph*2+1),1)
218
   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
219 220 221 222 223 224 225 226 227 228 229 230 231

   !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,sm(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
232 233 234 235
         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
236
      ELSE
237 238 239 240
         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
241 242 243 244 245 246 247 248 249 250 251 252
      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)
   CALL closeXMLElement('densityConvergence')

   !fix charge of the new density
   CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
253
             inDen%pw,inDen%vacxy,inDen%mt,inDen%vacz,.FALSE.,.false., fix)
Gregor Michalicek's avatar
Gregor Michalicek committed
254

255 256 257 258 259 260
   IF(atoms%n_u.NE.n_u_keep) THEN
      inDen%mmpMat = n_mmpTemp
   END IF

   atoms%n_u=n_u_keep

261 262 263 264 265 266 267 268 269
   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
270
   IF (atoms%n_u > 0) THEN
271 272
      IF (.NOT.l_densityMatrixPresent) THEN
         inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
273 274 275 276 277 278 279 280 281 282 283
         INQUIRE(file='broyd',exist=l_exist)
         IF (l_exist) THEN
            CALL system('rm broyd')
            PRINT *,"broyd file deleted because new density matrix is created."
         ENDIF
         INQUIRE(file='broyd.7',exist=l_exist)
         IF (l_exist) THEN
            CALL system('rm broyd.7')
            PRINT *,"broyd.7 file deleted because new density matrix is created."
         ENDIF
      END IF
Gregor Michalicek's avatar
Gregor Michalicek committed
284 285
   ENDIF

286 287 288 289
   !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)

290 291
   inDen%iter = inDen%iter + 1

Gregor Michalicek's avatar
Gregor Michalicek committed
292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307
   IF (input%imix.GT.0) THEN
      CLOSE (57)
      CLOSE (59)
   END IF

   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

308
END MODULE m_mix