mix.F90 16.6 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 16 17 18 19 20 21 22 23 24 25 26
  !************************************************************************
CONTAINS
  SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym, cell, it, noco, oneD,hybrid)
    !
#include"cpp_double.h"
    USE m_loddop
    USE m_wrtdop
    USE m_brysh1
    USE m_stmix
    USE m_broyden
    USE m_brysh2
    USE m_metric
    USE m_qfix
27
    USE m_icorrkeys
28
    USE m_types
29
    USE m_xmlOutput
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
    IMPLICIT NONE
    TYPE(t_oneD),INTENT(IN)     :: oneD
    TYPE(t_hybrid),INTENT(INOUT):: hybrid !ddist is modified here
    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
    !     ..
    !     .. Scalar Arguments ..
    INTEGER :: nrhomfile=26
    INTEGER, INTENT (IN) :: it    


    !     ..
    !     .. Local Scalars ..
    REAL fix,intfac,vacfac
    INTEGER i,iter,imap,js,mit,nt,nt1,irecl
    INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
    INTEGER iq2,iq3,ivac,imz ,iofl
    INTEGER n_u_keep
    LOGICAL lexist,l_ldaU
    INTEGER d1,d10,asciioffset
    CHARACTER (len=5) cdnfile

    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: qpw(:,:),rhtxy(:,:,:,:)
    REAL, ALLOCATABLE :: rho(:,:,:,:),rht(:,:,:)
    REAL dist(6)
    REAL, ALLOCATABLE :: sm(:),fsm(:) 
    !---> off-diagonal part of the density matrix
    COMPLEX, ALLOCATABLE :: cdom(:),cdomvz(:,:),cdomvxy(:,:,:)
    COMPLEX, ALLOCATABLE :: n_mmp(:,:,:,:,:)
67
    CHARACTER(LEN=20)    :: attributes(2)
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 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 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
    !     ..
    !     .. 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
    mmaph = intfac*stars%n3d + atoms%ntypd*(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) mmap = mmap + 2*stars%n3d&
                           + 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
       IF ( MOD(i,14*input%jspins) == 1 ) THEN    ! was already mixed in u_mix
          atoms%n_u = 0
          ALLOCATE ( n_mmp(-3:-3,-3:-3,1,1,2) )
       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
          REWIND (69)
          ALLOCATE ( n_mmp(-3:3,-3:3,atoms%n_u,input%jspins,2) )
          READ (69,'(7f20.13)') n_mmp(:,:,:,:,:)
       ELSE
          CALL juDFT_error("strange n_mmp_mat-file...",calledby ="mix")
       ENDIF
       CLOSE (69)
    ELSE
       atoms%n_u=0
       ALLOCATE ( n_mmp(-3:-3,-3:-3,1,1,2) )
    ENDIF
    !
    ALLOCATE (sm(mmap),fsm(mmap))

    ALLOCATE (qpw(stars%n3d,input%jspins),rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,input%jspins),&
                   rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,input%jspins),rht(vacuum%nmzd,2,input%jspins) )

    IF (noco%l_noco) THEN
       ALLOCATE (cdom(stars%n3d),cdomvz(vacuum%nmzd,2), cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2))
    ELSE
       ALLOCATE ( cdom(1),cdomvz(1,1),cdomvxy(1,1,1) )
    ENDIF
    !
    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

    !--->    generate name of file to hold the results of this iteration
    !gs+  changed cdn filename gen. to cdn//tens//digits
    !gs   fixed cdn filehandle to 71, otherwise conflict after it>=18
    !gs   plus: you don't want to have open so many fh's
    !gs    construct 2char rep of iter --> '01'-->'Z9'
    d1 = MOD(it,10)
    d10 = INT( (it + 0.5)/10 )
    asciioffset = IACHAR('1')-1
    IF ( d10.GE.10 ) asciioffset = IACHAR('7')
    cdnfile = 'cdn'//ACHAR(d10+asciioffset)//ACHAR(d1+IACHAR('1')-1)
    !      WRITE (*,*) d1,d10,'>',cdnfile,'<' !gsdbg
    IF (.NOT.noco%l_noco) THEN
       nt1 = 72
       OPEN (nt1,file=cdnfile,form='unformatted',status='unknown')
       REWIND nt1
    ENDIF
    !gs-

    !---> reload densities of current iteration
    IF (noco%l_noco) THEN
       !---> initialize arrays for the off-diagonal part of the density matrix
       cdom(:) = CMPLX(0.0,0.0)
       IF (input%film) THEN
          cdomvz(:,:) = CMPLX(0.0,0.0)
          cdomvxy(:,:,:) = CMPLX(0.0,0.0)
       END IF

       !--->    reload input density matrix from file rhomat_inp
       OPEN (nrhomfile,FILE='rhomat_inp',FORM='unformatted', STATUS='unknown')
       !--->    first the diagonal elements of the density matrix
       CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
                        nrhomfile, iter,rho,qpw,rht,rhtxy)
       !--->    and then the off-diagonal part
       READ (nrhomfile,END=150,ERR=50) (cdom(iq3),iq3=1,stars%ng3)
       IF (input%film) THEN
          READ (nrhomfile,END=75,ERR=50) ((cdomvz(imz,ivac), imz=1,vacuum%nmz),ivac=1,vacuum%nvac)
          READ (nrhomfile,END=75,ERR=50) (((cdomvxy(imz,iq2-1,ivac), imz=1,vacuum%nmzxy),iq2=2,oneD%odi%nq2),ivac=1,vacuum%nvac)
       ENDIF
       GOTO 150
50     WRITE(6,*)'rhodirgen: ERROR: Problems while reading density'
       WRITE(6,*)'matrix from file rhomat_inp.'
       CALL juDFT_error("ERROR while reading file rhomat_inp",calledby ="mix")
75     WRITE(6,*)'rhomatdir: ERROR: reached end of file rhomat_inp'
       WRITE(6,*)'while reading the vacuum part of the off-diagonal'
       WRITE(6,*)'element of the desity matrix.'
       CALL juDFT_error("ERROR while reading file rhomat_inp",calledby ="mix")
150    CLOSE (nrhomfile)
    ELSE
       nt = 71                 !gs see above
       OPEN (nt,file='cdn1',form='unformatted',status='old')
       REWIND nt
       CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
            nt, iter,rho,qpw,rht,rhtxy)
       !--->      write density to file for storage
       CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
            nt1, iter,rho,qpw,rht,rhtxy)
    ENDIF
    !
    !--->  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,&
         intfac,vacfac,qpw,rho,rht,rhtxy,cdom,cdomvz,cdomvxy,n_mmp, nmap,nmaph,mapmt,mapvac,mapvac2,sm) 
    !     load output charge density
    !
    IF (noco%l_noco) THEN
       !--->    reload output density matrix from file rhomat_out
       OPEN (nrhomfile,FILE='rhomat_out',FORM='unformatted', STATUS='unknown')
       !--->    first the diagonal elements of the density matrix
       CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
                        nrhomfile, iter,rho,qpw,rht,rhtxy)
       !--->    and then the off-diagonal part
       READ (nrhomfile) (cdom(iq3),iq3=1,stars%ng3)
       IF (input%film) THEN
          READ (nrhomfile) ((cdomvz(imz,ivac),imz=1,vacuum%nmz),ivac=1,vacuum%nvac)
          READ (nrhomfile) (((cdomvxy(imz,iq2-1,ivac),imz=1,vacuum%nmzxy), iq2=2,oneD%odi%nq2),ivac=1,vacuum%nvac)
       ENDIF
       CLOSE (nrhomfile)
    ELSE
       CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
                        nt, iter,rho,qpw,rht,rhtxy)
       !--->      write density to file for storage
       CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
                        nt1, iter,rho,qpw,rht,rhtxy)
       CLOSE(nt1)
    ENDIF
    !
    !--->  put output charge density into arrays fsm 
    !
    CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD, intfac,vacfac,qpw,rho,rht,rhtxy,cdom,&
         cdomvz,cdomvxy,n_mmp, 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=(2*nmap+1)*8
    IF (input%imix.GE.3) THEN
       OPEN (57,file='broyd',form='unformatted',status='unknown')
       OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
            recl=irecl,form='unformatted',status='unknown')
    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
    !
    CALL brysh2(input,stars,atoms,sphhar, noco,vacuum, sym,sm, n_mmp,oneD,&
         qpw,rho,rht,rhtxy,cdom,cdomvz,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
    !
310
    IF(hybrid%l_calhf) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
311
       CALL openXMLElement('densityConvergence',(/'units  ','comment'/),(/'me/bohr^3','HF       '/))
312 313 314 315
    ELSE
       CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
    END IF
    DO js = 1,input%jspins
316 317 318 319
       dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, sm(nmaph*(js-1)+1),1)

       hybrid%ddist(js) = 1000*SQRT(ABS(dist(js)/cell%vol))

320
       attributes = ''
321 322 323
       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/)))
324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342
       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)
343 344
       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/)))
345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
       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
    DEALLOCATE (sm,fsm)
361
    CALL closeXMLElement('densityConvergence')
362 363 364 365 366
    !
    !----> output of mixed densities
    !
    !     ---> fix the new density
    CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
367
                   qpw,rhtxy,rho,rht,.FALSE., fix)
368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419

    iter = iter + 1
    IF (noco%l_noco) THEN
       !
       !--->    write mixed density to file rhomat_out
       !
       OPEN (nrhomfile,FILE='rhomat_inp',FORM='unformatted', STATUS='unknown')
       !--->    first the diagonal elements of the density matrix
       CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
                        nrhomfile, iter,rho,qpw,rht,rhtxy)
       !--->    and then the off-diagonal part
       WRITE (nrhomfile) (cdom(iq3),iq3=1,stars%ng3)
       IF (input%film) THEN
          WRITE (nrhomfile) ((cdomvz(imz,ivac),imz=1,vacuum%nmz),ivac=1,vacuum%nvac)
          WRITE (nrhomfile) (((cdomvxy(imz,iq2-1,ivac),imz=1,vacuum%nmzxy),&
               iq2=2,oneD%odi%nq2),ivac=1,vacuum%nvac)
       ENDIF
       CLOSE (nrhomfile)
    ELSE
       !
       !--->    write new density onto unit 71 (overwrite)
       !
       REWIND nt
       CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
                        nt, iter,rho,qpw,rht,rhtxy)
       CLOSE (nt)
    ENDIF
    DEALLOCATE ( cdom,cdomvz,cdomvxy )
    IF ( atoms%n_u > 0 ) THEN
       OPEN (69,file='n_mmp_mat',status='replace',form='formatted')
       WRITE (69,'(7f20.13)') n_mmp(:,:,:,:,1)
       CLOSE (69)
    ENDIF
    DEALLOCATE (n_mmp)

    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)

    DEALLOCATE (qpw,rhtxy,rho,rht)
    atoms%n_u=n_u_keep
  END SUBROUTINE mix
END MODULE m_mix