mix.F90 10.3 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 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
   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
48
   TYPE(t_potden),INTENT(IN)     :: outDen
Gregor Michalicek's avatar
Gregor Michalicek committed
49
   TYPE(t_results),INTENT(INOUT) :: results
50
   TYPE(t_potden),INTENT(INOUT)  :: inDen
51
   INTEGER, INTENT(IN)           :: archiveType
Gregor Michalicek's avatar
Gregor Michalicek committed
52 53 54

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

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

   !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

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

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

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

   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

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

   ALLOCATE (sm(mmap),fsm(mmap))
118
   dist(:) = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138

   !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

139
   !put input charge density into array sm
140

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

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

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

   !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,&
157
                   hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
Gregor Michalicek's avatar
Gregor Michalicek committed
158 159 160
   END IF

   !initiatlize mixed density and extract it with brysh2 call
161 162 163 164
   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
165

166
   CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen) 
Gregor Michalicek's avatar
Gregor Michalicek committed
167 168 169 170 171 172 173 174 175 176 177 178 179

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

Gregor Michalicek's avatar
Gregor Michalicek committed
181 182 183 184 185 186 187 188
   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
189 190
         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
191
      ELSE
192 193
         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
194 195 196
      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)
197
   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
198 199 200 201 202 203 204 205 206 207 208 209 210

   !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
211 212 213 214
         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
215
      ELSE
216 217 218 219
         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
220 221 222 223 224 225 226 227 228 229 230 231
      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,&
232
             inDen%pw,inDen%vacxy,inDen%mt,inDen%vacz,.FALSE.,.false., fix)
Gregor Michalicek's avatar
Gregor Michalicek committed
233

234 235 236 237 238 239
   IF(atoms%n_u.NE.n_u_keep) THEN
      inDen%mmpMat = n_mmpTemp
   END IF

   atoms%n_u=n_u_keep

240 241 242 243 244 245 246 247 248
   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
249
   IF (atoms%n_u > 0) THEN
250 251
      IF (.NOT.l_densityMatrixPresent) THEN
         inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
252
         CALL resetBroydenHistory()
253
      END IF
Gregor Michalicek's avatar
Gregor Michalicek committed
254 255
   ENDIF

256 257 258 259
   !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)

260 261
   inDen%iter = inDen%iter + 1

Gregor Michalicek's avatar
Gregor Michalicek committed
262 263 264 265 266 267 268 269 270 271 272
   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

273
END MODULE m_mix