mix.F90 13.8 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 7 8
MODULE m_mix
  USE m_juDFT
  !*************************************************************************
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 15
  !************************************************************************
CONTAINS
16 17
  SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,it,noco,oneD,&
                 hybrid,inDen,outDen,results)
18 19
    !
#include"cpp_double.h"
20
    USE m_constants
21
    USE m_cdn_io
22 23 24 25 26 27 28
    USE m_brysh1
    USE m_stmix
    USE m_broyden
    USE m_brysh2
    USE m_metric
    USE m_qfix
    USE m_types
29
    USE m_xmlOutput
30 31
    IMPLICIT NONE
    TYPE(t_oneD),INTENT(IN)     :: oneD
Daniel Wortmann's avatar
Daniel Wortmann committed
32
    TYPE(t_hybrid),INTENT(IN)   :: hybrid 
33 34 35 36 37 38 39 40
    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
41
    TYPE(t_potden),INTENT(IN)   :: inDen, outDen
42
    TYPE(t_results),INTENT(INOUT)::results
43 44 45 46 47
    !     ..
    !     .. Scalar Arguments ..
    INTEGER :: nrhomfile=26
    INTEGER, INTENT (IN) :: it    

48 49
    ! Local type instances
    TYPE(t_potden) :: mixDen
50 51 52

    !     ..
    !     .. Local Scalars ..
53
    REAL fix,intfac,vacfac, fermiEnergyTemp
54
    INTEGER i,iter,imap,js,mit,nt,irecl
55
    INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
56
    INTEGER iq2,iq3,ivac,imz ,iofl, archiveType
57
    INTEGER n_u_keep
58
    LOGICAL lexist,l_ldaU, l_qfix
59 60 61 62 63 64

    !     ..
    !     .. Local Arrays ..
    REAL dist(6)
    REAL, ALLOCATABLE :: sm(:),fsm(:) 
    !---> off-diagonal part of the density matrix
65
    CHARACTER(LEN=20)    :: attributes(2)
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
    !     ..
    !     .. Intrinsic Functions ..
    INTRINSIC char,sqrt
    !     .. 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
    !---> In systems without inversions symmetry the interstitial star-
    !---> coefficients are complex. Thus twice as many numbers have to be
    !---> stored.
    n_u_keep=atoms%n_u
    IF (sym%invs) THEN
       intfac = 1.0
    ELSE
       intfac = 2.0
    ENDIF
    !---> The corresponding is true for the coeff. of the warping vacuum
    !---> density depending on the two dimensional inversion.
    IF (sym%invs2) THEN
       vacfac = 1.0
    ELSE
       vacfac = 2.0
    ENDIF
97
    mmaph = intfac*stars%ng3 + atoms%ntype*(sphhar%nlhd+1)*atoms%jmtd&
98 99 100 101 102
               + 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.
103
    IF (noco%l_noco) mmap = mmap + 2*stars%ng3&
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
                           + 2*vacuum%nmzxyd*(oneD%odi%n2d-1)*vacuum%nvac + 2*vacuum%nmzd*vacuum%nvac

    INQUIRE (file='n_mmp_mat',exist=l_ldaU) 
    IF (l_ldaU) THEN
       !
       ! In an LDA+U caclulation, also the density matrix is included in the
       ! supervectors (sm,fsm) if no mixing factors are in the n_mmp_mat-file
       !
       OPEN (69,file='n_mmp_mat',status='old',form='formatted')
       i = 0 
       DO
          READ (69,*,iostat=iofl)
          IF (iofl < 0) EXIT
          i = i + 1
       ENDDO
119
       CLOSE (69)
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
       IF ( MOD(i,14*input%jspins) == 1 ) THEN    ! was already mixed in u_mix
          atoms%n_u = 0
       ELSEIF ( MOD(i,28*input%jspins)== 0 ) THEN ! mix here 
          atoms%n_u = i / (28 * input%jspins )    ! atoms%n_u atoms have lda+u applied
          mmap = mmap + 7 * i / 2     ! add 7*7 complex numbers per atoms%n_u
       ELSE
          CALL juDFT_error("strange n_mmp_mat-file...",calledby ="mix")
       ENDIF
    ELSE
       atoms%n_u=0
    ENDIF
    !
    ALLOCATE (sm(mmap),fsm(mmap))

    IF (noco%l_noco) THEN
135
       archiveType = CDN_ARCHIVE_TYPE_NOCO_const
136
    ELSE
137
       archiveType = CDN_ARCHIVE_TYPE_CDN1_const
138
    ENDIF
139

140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
    !
    INQUIRE (file='broyd.'//CHAR(input%imix+48),exist=lexist)
    DO i = 1,6
       dist(i) = 0.0
    END DO
    mit = 0
    !---> detremine 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
    ELSEIF (input%imix.EQ.3) THEN
       IF ( .NOT.lexist) mit = 1
       WRITE (16,FMT='(a,f10.5)') 'BROYDEN FIRST MIXING',input%alpha
    ELSEIF (input%imix.EQ.5) THEN
       IF (.NOT.lexist) mit = 1
       WRITE (16,FMT='(a,f10.5)') 'BROYDEN SECOND MIXING',input%alpha
    ELSEIF (input%imix.EQ.7) THEN
       IF (.NOT.lexist) mit = 1
       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

    !---> reload densities of current iteration
169 170
!    CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
!                     CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
171 172 173 174 175 176 177

    !
    !--->  put input charge density into arrays sm 
    !      in the spin polarized case the arrays consist of 
    !      spin up and spin down densities

    CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
178 179 180 181 182 183
                intfac,vacfac,inDen%pw,inDen%mt,inDen%vacz,inDen%vacxy,inDen%cdom,&
                inDen%cdomvz,inDen%cdomvxy,inDen%mmpMat(-lmaxU_const,-lmaxU_const,1,1),&
                nmap,nmaph,mapmt,mapvac,mapvac2,sm) 

!    CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
!         intfac,vacfac,qpw,rho,rht,rhtxy,cdom,cdomvz,cdomvxy,n_mmp(-3,-3,1,1,1), nmap,nmaph,mapmt,mapvac,mapvac2,sm) 
184

185
    !     load output charge density
186 187
!    CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
!                     CDN_OUTPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
188

189 190 191
    !
    !--->  put output charge density into arrays fsm 
    !
192 193 194 195 196 197 198 199

    CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
                intfac,vacfac,outDen%pw,outDen%mt,outDen%vacz,outDen%vacxy,outDen%cdom,&
                outDen%cdomvz,outDen%cdomvxy,outDen%mmpMat(-lmaxU_const,-lmaxU_const,1,1),&
                nmap,nmaph,mapmt,mapvac,mapvac2,fsm)

!    CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD, intfac,vacfac,qpw,rho,rht,rhtxy,cdom,&
!         cdomvz,cdomvxy,n_mmp(-3,-3,1,1,2), nmap,nmaph,mapmt,mapvac,mapvac2,fsm)
200 201 202 203 204 205 206 207 208
    !
    ! --> store fsm - sm the difference on fsm
    !
    DO imap = 1,nmap
       fsm(imap) = fsm(imap) - sm(imap)
    END DO
    !
    ! open files for broyden
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
209
    irecl=(nmap+1)*8
210
    IF (input%imix.GE.3) THEN
211 212 213 214 215 216 217 218 219
       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',&
               recl=irecl,form='unformatted',status='unknown')
       ELSE
          OPEN (57,file='broyd',form='unformatted',status='unknown')
          OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
               recl=irecl,form='unformatted',status='unknown')
       ENDIF
220 221 222 223 224 225 226 227 228 229 230 231 232
    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
    !     call timestamp(-3)
    !----->  load output densities
    !
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254

    CALL mixDen%init(stars,atoms,sphhar,vacuum,oneD,input%jspins,.FALSE.)
    IF (noco%l_noco) THEN
       ALLOCATE (mixDen%cdom(stars%ng3),mixDen%cdomvz(vacuum%nmzd,2))
       ALLOCATE (mixDen%cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2))
    ELSE
       ALLOCATE (mixDen%cdom(1),mixDen%cdomvz(1,1),mixDen%cdomvxy(1,1,1))
    ENDIF
    IF (atoms%n_u.GT.0) THEN
       ALLOCATE (mixDen%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u,input%jspins))
    ELSE
       ALLOCATE (mixDen%mmpMat(-lmaxU_const:-lmaxU_const,-lmaxU_const:-lmaxU_const,1,input%jspins))
    END IF
    mixDen%cdom = CMPLX(0.0,0.0)
    mixDen%cdomvz = CMPLX(0.0,0.0)
    mixDen%cdomvxy = CMPLX(0.0,0.0)
    mixDen%mmpMat = CMPLX(0.0,0.0)


    CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,mixDen%mmpMat,oneD,&
                mixDen%pw,mixDen%mt,mixDen%vacz,mixDen%vacxy,mixDen%cdom,&
                mixDen%cdomvz,mixDen%cdomvxy) 
255 256 257 258 259 260 261 262 263 264 265 266

    !
    !----->  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
    !
267
    IF(hybrid%l_calhf) THEN
Daniel Wortmann's avatar
Bugfix  
Daniel Wortmann committed
268
       CALL openXMLElement('densityConvergence',(/'units  ','comment'/),(/'me/bohr^3','HF       '/))
269 270 271
    ELSE
       CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
    END IF
272
    iter = inDen%iter
273
    DO js = 1,input%jspins
274 275
       dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, sm(nmaph*(js-1)+1),1)

276
       attributes = ''
277 278 279
       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/)))
280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
       IF( hybrid%l_calhf ) THEN
          WRITE (16,FMT=7901) js,iter,1000*SQRT(ABS(dist(js)/cell%vol))
          WRITE ( 6,FMT=7901) js,iter,1000*SQRT(ABS(dist(js)/cell%vol))
       ELSE
          WRITE (16,FMT=7900) js,iter,1000*SQRT(ABS(dist(js)/cell%vol))
          WRITE ( 6,FMT=7900) js,iter,1000*SQRT(ABS(dist(js)/cell%vol))
       END IF
    ENDDO
    IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,sm(nmaph*2+1),1)
    IF (noco%l_noco) WRITE (6,FMT=7900) 3,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,sm(nmaph+1),1)
       dist(4) = dist(1) + dist(2) + 2.0e0*dist(3)
       dist(5) = dist(1) + dist(2) - 2.0e0*dist(3)
299 300
       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/)))
301 302 303 304 305 306 307 308 309 310 311
       IF( hybrid%l_calhf ) THEN
          WRITE (16,FMT=8001) iter,1000*SQRT(ABS(dist(4)/cell%vol))
          WRITE (16,FMT=8011) iter,1000*SQRT(ABS(dist(5)/cell%vol))
          WRITE ( 6,FMT=8001) iter,1000*SQRT(ABS(dist(4)/cell%vol))
          WRITE ( 6,FMT=8011) iter,1000*SQRT(ABS(dist(5)/cell%vol))
       ELSE
          WRITE (16,FMT=8000) iter,1000*SQRT(ABS(dist(4)/cell%vol))
          WRITE (16,FMT=8010) iter,1000*SQRT(ABS(dist(5)/cell%vol))
          WRITE ( 6,FMT=8000) iter,1000*SQRT(ABS(dist(4)/cell%vol))
          WRITE ( 6,FMT=8010) iter,1000*SQRT(ABS(dist(5)/cell%vol))
       END IF
312
      
313 314 315 316
       ! 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
317
    results%last_distance=maxval(1000*SQRT(ABS(dist/cell%vol)))
318
    DEALLOCATE (sm,fsm)
319
    CALL closeXMLElement('densityConvergence')
320 321 322 323 324
    !
    !----> output of mixed densities
    !
    !     ---> fix the new density
    CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
325
              mixDen%pw,mixDen%vacxy,mixDen%mt,mixDen%vacz,.FALSE.,.false., fix)
326

327
    CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
328 329
                      1,results%last_distance,results%ef,.TRUE.,iter,mixDen%mt,mixDen%pw,mixDen%vacz,&
                      mixDen%vacxy,mixDen%cdom,mixDen%cdomvz,mixDen%cdomvxy)
330

331 332
    IF ( atoms%n_u > 0 ) THEN
       OPEN (69,file='n_mmp_mat',status='replace',form='formatted')
333
       WRITE (69,'(7f20.13)') mixDen%mmpMat(:,:,:,:)
334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352
       CLOSE (69)
    ENDIF

    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)

    atoms%n_u=n_u_keep
  END SUBROUTINE mix
END MODULE m_mix