mix.F90 12.2 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 19 20

SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
               hybrid,inDen,outDen,results)

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 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 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 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160

   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
   TYPE(t_potden),INTENT(IN)     :: inDen, outDen
   TYPE(t_results),INTENT(INOUT) :: results

   !Local type instances
   TYPE(t_potden) :: mixDen

   !Local Scalars
   REAL fix,intfac,vacfac
   INTEGER i,iter,imap,js,mit,irecl
   INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
   INTEGER iofl,archiveType,n_u_keep
   LOGICAL lexist,l_ldaU

   !Local Arrays
   REAL dist(6)
   REAL, ALLOCATABLE :: sm(:), fsm(:)
   CHARACTER(LEN=20) :: attributes(2)

   !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.
   IF (sym%invs) THEN
      intfac = 1.0
   ELSE
      intfac = 2.0
   END IF

   !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
   END IF

   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

   n_u_keep=atoms%n_u
   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
      END DO
      CLOSE (69)

      IF ( MOD(i,14*input%jspins) == 1 ) THEN     ! was already mixed in u_mix
         atoms%n_u = 0
      ELSE IF ( MOD(i,28*input%jspins)== 0 ) THEN ! mix here 
         atoms%n_u = i / (28 * input%jspins ) ! calculate number of U parameters
         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")
      END IF
   ELSE
      atoms%n_u=0
   END IF ! l_ldaU

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

   INQUIRE (file='broyd.'//CHAR(input%imix+48),exist=lexist)
   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
      IF ( .NOT.lexist) mit = 1
      WRITE (16,FMT='(a,f10.5)') 'BROYDEN FIRST MIXING',input%alpha
   ELSE IF (input%imix.EQ.5) THEN
      IF (.NOT.lexist) mit = 1
      WRITE (16,FMT='(a,f10.5)') 'BROYDEN SECOND MIXING',input%alpha
   ELSE IF (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

161 162
   !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
163 164 165 166 167
   CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
               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) 

168
   !put output charge density into array fsm
Gregor Michalicek's avatar
Gregor Michalicek committed
169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
   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)

   !store fsm - sm the difference on fsm
   DO imap = 1,nmap
      fsm(imap) = fsm(imap) - sm(imap)
   END DO

   ! 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',&
185
               recl=irecl,form='unformatted',status='unknown')
Gregor Michalicek's avatar
Gregor Michalicek committed
186 187 188
      ELSE
         OPEN (57,file='broyd',form='unformatted',status='unknown')
         OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
189
               recl=irecl,form='unformatted',status='unknown')
Gregor Michalicek's avatar
Gregor Michalicek committed
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
      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
   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))
      archiveType = CDN_ARCHIVE_TYPE_NOCO_const
   ELSE
      ALLOCATE (mixDen%cdom(1),mixDen%cdomvz(1,1),mixDen%cdomvxy(1,1,1))
      archiveType = CDN_ARCHIVE_TYPE_CDN1_const
   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) 

   !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
   iter = inDen%iter
   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
         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
   END DO
   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)
      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
         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

      !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,&
             mixDen%pw,mixDen%vacxy,mixDen%mt,mixDen%vacz,.FALSE.,.false., fix)

   !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.,iter,mixDen%mt,mixDen%pw,mixDen%vacz,&
                     mixDen%vacxy,mixDen%cdom,mixDen%cdomvz,mixDen%cdomvxy)

   IF (atoms%n_u > 0) THEN
      OPEN (69,file='n_mmp_mat',status='replace',form='formatted')
      WRITE (69,'(7f20.13)') mixDen%mmpMat(:,:,:,:)
      CLOSE (69)
   ENDIF

   IF (input%imix.GT.0) THEN
      CLOSE (57)
      CLOSE (59)
   END IF

   atoms%n_u=n_u_keep

   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