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

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

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

   !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

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

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

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

   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

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

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

   !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

141
   !put input charge density into array sm
142

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

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

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

154
155
156
157
158
   l_pot = .FALSE.
   ! Apply metric w to fsm and store in fmMet:  w |fsm>
   CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
               mmap,nmaph,mapmt,mapvac2,fsm,fmMet,l_pot)

Gregor Michalicek's avatar
Gregor Michalicek committed
159
160
161
162
   !mixing of the densities
   IF (input%imix.EQ.0) THEN
      CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
   ELSE
163
164
      CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
                   hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
165
166
167

!     Replace the broyden call above by the commented metric and broyden2 calls
!     below to switch on the continuous restart of the Broyden method.
168
      ! Apply metric w to sm and store in smMet:  w |sm>
169
170
171
172
173
!      CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                  mmap,nmaph,mapmt,mapvac2,sm,smMet,l_pot)
!
!      CALL broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
!                    hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm,fmMet,smMet)
Gregor Michalicek's avatar
Gregor Michalicek committed
174
175
176
   END IF

   !initiatlize mixed density and extract it with brysh2 call
177
178
179
180
   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
181

182
   CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen) 
Gregor Michalicek's avatar
Gregor Michalicek committed
183
184
185
186

   !calculate the distance of charge densities...

   !induce metric in fsm use sm as an output array: |sm> = w |fsm>
187
188
!   CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!               mmap,nmaph,mapmt,mapvac2,fsm, sm)
Gregor Michalicek's avatar
Gregor Michalicek committed
189
190
191
192
193
194
195

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

Gregor Michalicek's avatar
Gregor Michalicek committed
197
   DO js = 1,input%jspins
198
      dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, fmMet(nmaph*(js-1)+1),1)
Gregor Michalicek's avatar
Gregor Michalicek committed
199
200
201
202
203
204

      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
205
206
         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
207
      ELSE
208
209
         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
210
211
      END IF
   END DO
212
   IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,fmMet(nmaph*2+1),1)
213
   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
214
215
216
217
218

   !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
219
      dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,fmMet(nmaph+1),1)
Gregor Michalicek's avatar
Gregor Michalicek committed
220
221
222
223
224
225
226
      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
227
228
229
230
         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
231
      ELSE
232
233
234
235
         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
236
237
238
239
240
241
242
      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)))
243
   DEALLOCATE (sm,fsm,smMet,fmMet)
Gregor Michalicek's avatar
Gregor Michalicek committed
244
245
246
247
   CALL closeXMLElement('densityConvergence')

   !fix charge of the new density
   CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
248
             inDen%pw,inDen%vacxy,inDen%mt,inDen%vacz,.FALSE.,.false., fix)
Gregor Michalicek's avatar
Gregor Michalicek committed
249

250
251
252
253
254
255
   IF(atoms%n_u.NE.n_u_keep) THEN
      inDen%mmpMat = n_mmpTemp
   END IF

   atoms%n_u=n_u_keep

256
257
258
259
260
261
262
263
264
   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
265
   IF (atoms%n_u > 0) THEN
266
267
      IF (.NOT.l_densityMatrixPresent) THEN
         inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
268
         CALL resetBroydenHistory()
269
      END IF
Gregor Michalicek's avatar
Gregor Michalicek committed
270
271
   ENDIF

272
273
274
275
   !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)

276
277
   inDen%iter = inDen%iter + 1

Gregor Michalicek's avatar
Gregor Michalicek committed
278
279
280
281
282
283
284
285
286
287
288
   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

289
END MODULE m_mix