mix.F90 16.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
MODULE m_mix
  USE m_juDFT
  !*************************************************************************
  !****  mixing of charge densities:                                    ****
  !****    IMIX= 0 : linear mixing                                     ****
  !****    IMIX = 3 : BROYDEN'S FIRST METHOD                            ****
  !****    IMIX = 5 : BROYDEN'S SECOND METHOD                           ****
  !****    IMIX = 7 : GENERALIZED ANDERSEN METHOD                       ****
  !****  implementation to flapw7 ..... R.Pentcheva, D.Vogtenhuber      ****
  !************************************************************************
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
    USE icorrkeys
    USE m_types
25
    USE m_xmlOutput
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
    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(:,:,:,:,:)
63
    CHARACTER(LEN=20)    :: attributes(2)
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
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
    !     ..
    !     .. 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
    !
306
    IF(hybrid%l_calhf) THEN
Daniel Wortmann's avatar
Bugfix    
Daniel Wortmann committed
307
       CALL openXMLElement('densityConvergence',(/'units  ','comment'/),(/'me/bohr^3','HF       '/))
308
309
310
311
    ELSE
       CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
    END IF
    DO js = 1,input%jspins
312
313
314
315
       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))

316
       attributes = ''
317
318
319
       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/)))
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
       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)
339
340
       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/)))
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
       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)
357
    CALL closeXMLElement('densityConvergence')
358
359
360
361
362
363
364
365
366
367
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
    !
    !----> output of mixed densities
    !
    !     ---> fix the new density
    CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
                   qpw,rhtxy,rho,rht, fix)

    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