cdngen.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
MODULE m_cdngen
7

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
USE m_juDFT

CONTAINS

SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
                  dimension,kpts,atoms,sphhar,stars,sym,obsolete,&
                  enpara,cell,noco,jij,vTot,results,oneD,coreSpecInput,&
                  inIter,inDen,outDen)

   !*****************************************************
   !    Charge density generator
   !    calls cdnval to generate the valence charge and the
   !    core routines for the core contribution
   !*****************************************************

   USE m_constants
   USE m_umix
   USE m_prpqfftmap
   USE m_cdnval
   USE m_cdn_io
   USE m_wrtdop
   USE m_cdntot
   USE m_cdnovlp
   USE m_qfix
   USE m_rwnoco
   USE m_cored
   USE m_coredr
   USE m_types
   USE m_xmlOutput
37
   USE m_magMoms
38
   USE m_orbMagMoms
39
#ifdef CPP_MPI
40 41
   USE m_mpi_bc_potden
   USE m_mpi_bc_coreden
42
#endif
43

44
   IMPLICIT NONE
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
   ! Type instance arguments
   TYPE(t_results),INTENT(INOUT)    :: results
   TYPE(t_mpi),INTENT(IN)           :: mpi
   TYPE(t_dimension),INTENT(IN)     :: dimension
   TYPE(t_oneD),INTENT(IN)          :: oneD
   TYPE(t_enpara),INTENT(INOUT)     :: enpara
   TYPE(t_obsolete),INTENT(IN)      :: obsolete
   TYPE(t_banddos),INTENT(IN)       :: banddos
   TYPE(t_sliceplot),INTENT(IN)     :: sliceplot
   TYPE(t_input),INTENT(IN)         :: input
   TYPE(t_vacuum),INTENT(IN)        :: vacuum
   TYPE(t_noco),INTENT(IN)          :: noco
   TYPE(t_jij),INTENT(IN)           :: jij
   TYPE(t_sym),INTENT(IN)           :: sym
   TYPE(t_stars),INTENT(IN)         :: stars
   TYPE(t_cell),INTENT(IN)          :: cell
   TYPE(t_kpts),INTENT(IN)          :: kpts
   TYPE(t_sphhar),INTENT(IN)        :: sphhar
   TYPE(t_atoms),INTENT(IN)         :: atoms
   TYPE(t_coreSpecInput),INTENT(IN) :: coreSpecInput
   TYPE(t_potden),INTENT(IN)        :: vTot
   TYPE(t_potden),INTENT(INOUT)     :: inDen,outDen

   !Scalar Arguments
   INTEGER, INTENT (IN)             :: eig_id
   INTEGER, INTENT (IN)             :: inIter

   ! Local type instances
   TYPE(t_noco) :: noco_new

   !Local Scalars
77
   REAL fix,qtot,scor,seig,stot,sval,dummy
78
   REAL sum,fermiEnergyTemp
79 80 81 82 83 84 85 86 87 88 89 90
   INTEGER iter,ivac,j,jspin,jspmax,k,n,nt,ieig,ikpt
   INTEGER  ityp,ilayer,urec,itype,iatom
   LOGICAL l_relax_any,exst,n_exist,l_qfix, l_enpara
   LOGICAL, PARAMETER :: l_st=.FALSE.

   !Local Arrays
   REAL stdn(atoms%ntype,dimension%jspd),svdn(atoms%ntype,dimension%jspd),alpha_l(atoms%ntype)
   REAL rh(dimension%msh,atoms%ntype,dimension%jspd),qint(atoms%ntype,dimension%jspd)
   REAL tec(atoms%ntype,DIMENSION%jspd),rhTemp(dimension%msh,atoms%ntype,dimension%jspd)
   REAL chmom(atoms%ntype,dimension%jspd),clmom(3,atoms%ntype,dimension%jspd)
   INTEGER,ALLOCATABLE :: igq_fft(:)
   REAL   ,ALLOCATABLE :: qvac(:,:,:,:),qvlay(:,:,:,:,:)
91

92
   !pk non-collinear (start)
93
   REAL    rhoint,momint,alphdiff
94
   INTEGER igq2_fft(0:stars%kq1_fft*stars%kq2_fft-1)
Gregor Michalicek's avatar
Gregor Michalicek committed
95
   COMPLEX,ALLOCATABLE :: qa21(:), cdomvz(:,:)
96 97 98
   !pk non-collinear (end)

   iter = inIter
Gregor Michalicek's avatar
Gregor Michalicek committed
99
   CALL outDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115

   IF (mpi%irank.EQ.0) THEN
      INQUIRE(file='enpara',exist=l_enpara)
      IF (l_enpara) OPEN (40,file ='enpara',form = 'formatted',status ='unknown')
   ENDIF
   ALLOCATE (qa21(atoms%ntype))
   ALLOCATE (qvac(dimension%neigd,2,kpts%nkpt,dimension%jspd))
   ALLOCATE (qvlay(dimension%neigd,vacuum%layerd,2,kpts%nkpt,dimension%jspd))
   ALLOCATE (igq_fft(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1))

   !initialize density arrays with zero
   qa21(:) = cmplx(0.0,0.0)
   qvac(:,:,:,:) = 0.0 
   qvlay(:,:,:,:,:) = 0.0

   outDen%iter = iter
116
        
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
   !Set up pointer for backtransformation of from g-vector in
   !positive domain fof carge density fftibox into stars
   !In principle this can also be done in main program once.
   !It is done here to save memory.
   CALL prp_qfft_map(stars,sym, input, igq2_fft,igq_fft)

   !in a non-collinear calcuation where the off-diagonal part of
   !density matrix in the muffin-tins is calculated, the a- and
   !b-coef. for both spins are needed at once. Thus, cdnval is only
   !called once and both spin directions are calculated in a single
   !go.
   IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('valenceDensity')

   jspmax = input%jspins
   IF (noco%l_mperp) jspmax = 1
   DO jspin = 1,jspmax
      CALL timestart("cdngen: cdnval")
      CALL cdnval(eig_id,&
                  mpi,kpts,jspin,sliceplot,noco, input,banddos,cell,atoms,enpara,stars, vacuum,dimension,&
                  sphhar,sym,obsolete,igq_fft,vTot%mt,vTot%vacz(:,:,jspin),oneD,coreSpecInput,&
Gregor Michalicek's avatar
Gregor Michalicek committed
137
                  outDen,results,qvac,qvlay,qa21, chmom,clmom)
138
      CALL timestop("cdngen: cdnval")
139
!-fo
140 141 142 143 144
   END DO

   ! lda+u
   IF ((atoms%n_u.GT.0).and.(mpi%irank.EQ.0)) CALL u_mix(input,atoms,inDen%mmpMat,outDen%mmpMat)

145
!+t3e
146
   IF (mpi%irank.EQ.0) THEN
147
!-t3e
148 149 150
      IF (l_enpara) CLOSE (40)
      CALL cdntot(stars,atoms,sym, vacuum,input,cell,oneD, outDen%pw,outDen%mt,outDen%vacz,.TRUE., qtot,dummy)
      CALL closeXMLElement('valenceDensity')
151
!---> changes
152 153
   END IF ! mpi%irank = 0

154 155 156 157 158 159 160 161 162 163
   results%seigc = 0.0
   IF (mpi%irank.EQ.0) THEN
      DO jspin = 1,input%jspins
         DO n = 1,atoms%ntype
            svdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
         END DO
      END DO
   END IF

   IF (input%kcrel.EQ.0) THEN
164 165 166 167 168 169 170 171 172 173 174 175 176
      ! Generate input file ecore for subsequent GW calculation
      ! 11.2.2004 Arno Schindlmayr
      IF ((input%gw.eq.1 .or. input%gw.eq.3).AND.(mpi%irank.EQ.0)) THEN
         OPEN (15,file='ecore',status='unknown', action='write',form='unformatted')
      END IF

      rh = 0.0
      tec = 0.0
      qint = 0.0
      IF (input%frcor) THEN
         IF (mpi%irank.EQ.0) THEN
            CALL readCoreDensity(input,atoms,dimension,rh,tec,qint)
         END IF
177
#ifdef CPP_MPI
178
         CALL mpi_bc_coreDen(mpi,atoms,input,dimension,rh,tec,qint)
179
#endif
180
      END IF
181
   END IF
182

183 184 185 186 187 188 189 190 191 192 193 194 195 196
   IF (.NOT.sliceplot%slice) THEN
      !add in core density
      IF (mpi%irank.EQ.0) THEN
         IF (input%kcrel.EQ.0) THEN
            DO jspin = 1,input%jspins
               CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh,tec,seig)
               rhTemp(:,:,jspin) = rh(:,:,jspin)
               results%seigc = results%seigc + seig
            END DO
         ELSE
            CALL coredr(input,atoms,seig, outDen%mt,dimension,sphhar,vTot%mt(:,0,:,:),qint,rh)
            results%seigc = results%seigc + seig
         END IF
      END IF
197
      DO jspin = 1,input%jspins
198
         IF (mpi%irank.EQ.0) THEN
199
            DO n = 1,atoms%ntype
200
               stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
201
            END DO
202
         END IF
203 204
         IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
            IF (jspin.EQ.2) THEN
205 206 207 208
               !pk non-collinear (start)
               !add the coretail-charge to the constant interstitial
               !charge (star 0), taking into account the direction of
               !magnetisation of this atom
209 210 211 212 213 214 215 216 217 218 219 220
               DO ityp = 1,atoms%ntype
                  rhoint = (qint(ityp,1) + qint(ityp,2)) /cell%volint/input%jspins/2.0
                  momint = (qint(ityp,1) - qint(ityp,2)) /cell%volint/input%jspins/2.0
                  !rho_11
                  outDen%pw(1,1) = outDen%pw(1,1) + rhoint + momint*cos(noco%beta(ityp))
                  !rho_22
                  outDen%pw(1,2) = outDen%pw(1,2) + rhoint - momint*cos(noco%beta(ityp))
                  !real part rho_21
                  outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
                  !imaginary part rho_21
                  outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
               END DO
221
               !pk non-collinear (end)
222 223 224 225 226 227
            END IF
         ELSE
            IF (input%ctail) THEN
               !+gu hope this works as well
               CALL cdnovlp(mpi,sphhar,stars,atoms,sym,dimension,vacuum,&
                            cell,input,oneD,l_st,jspin,rh(:,:,jspin),&
228 229 230 231
                            outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz)
            ELSE IF (mpi%irank.EQ.0) THEN
               DO ityp = 1,atoms%ntype
                  outDen%pw(1,jspin) = outDen%pw(1,jspin) + qint(ityp,jspin)/input%jspins/cell%volint
232
               END DO
233
            END IF
234 235 236
         END IF
      END DO
   END IF
237

238
   IF (input%kcrel.EQ.0) THEN
239 240 241 242
      IF (mpi%irank.EQ.0) THEN
         CALL writeCoreDensity(input,atoms,dimension,rhTemp,tec,qint)
      END IF
      IF ((input%gw.eq.1 .or. input%gw.eq.3).AND.(mpi%irank.EQ.0)) CLOSE(15)
243
   END IF
244 245 246

   IF (mpi%irank.EQ.0) THEN
      !block 2 unnecessary for slicing: begin
247
      IF (.NOT.sliceplot%slice) THEN
248
         CALL openXMLElementNoAttributes('allElectronCharges')
249 250
         CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
                   outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz,.TRUE.,.true.,fix)
251
         CALL closeXMLElement('allElectronCharges')
252
         !pk non-collinear (start)
253
         IF (noco%l_noco) THEN
254
            !fix also the off-diagonal part of the density matrix
Gregor Michalicek's avatar
Gregor Michalicek committed
255
            outDen%pw(:stars%ng3,3) = fix*outDen%pw(:stars%ng3,3)
256
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
257
               outDen%vacz(:,:,3:4) = fix*outDen%vacz(:,:,3:4)
Gregor Michalicek's avatar
Gregor Michalicek committed
258
               outDen%vacxy(:,:,:,3) = fix*outDen%vacxy(:,:,:,3)
259
            END IF
260
         END IF
261
         !pk non-collinear (end)
262

263 264
         !spin densities at the nucleus
         !and magnetic moment in the spheres
265 266 267
         IF (input%jspins.EQ.2) THEN
            WRITE (6,FMT=8000)
            WRITE (16,FMT=8000)
268
            DO n = 1,atoms%ntype
269 270 271 272 273
               sval = svdn(n,1) - svdn(n,input%jspins)
               stot = stdn(n,1) - stdn(n,input%jspins)
               scor = stot - sval
               WRITE (6,FMT=8010) n,stot,sval,scor,svdn(n,1),stdn(n,1)
               WRITE (16,FMT=8010) n,stot,sval,scor,svdn(n,1),stdn(n,1)
274
            END DO
275

276 277
            noco_new = noco

278
            CALL magMoms(dimension,input,atoms,noco_new,vTot,chmom,qa21)
279

280 281
            !save the new nocoinp file if the dierctions of the local
            !moments are relaxed or a constraint B-field is calculated.
282 283
            l_relax_any = .false.
            DO itype = 1,atoms%ntype
284 285
               l_relax_any = l_relax_any.OR.noco%l_relax(itype)
            END DO
286 287 288
            IF (l_relax_any.OR.noco%l_constr) THEN
               IF (.not. noco%l_mperp) THEN
                  CALL juDFT_error ("(l_relax_any.OR.noco).AND.(.NOT. )" ,calledby ="cdngen")
289
               END IF
290
               iatom = 1
291 292
               DO itype = 1, atoms%ntype
                  IF (noco%l_ss) THEN
293 294 295 296
                     alphdiff = 2.0*pi_const*(noco%qss(1)*atoms%taual(1,iatom) + &
                                              noco%qss(2)*atoms%taual(2,iatom) + &
                                              noco%qss(3)*atoms%taual(3,iatom) )
                     noco_new%alph(itype) = noco%alph(itype) - alphdiff
297 298 299 300 301 302 303 304 305
                     DO WHILE (noco_new%alph(n) > +pi_const)
                        noco_new%alph(n)= noco_new%alph(n) - 2.*pi_const
                     END DO
                     DO WHILE (noco_new%alph(n) < -pi_const)
                        noco_new%alph(n)= noco_new%alph(n) + 2.*pi_const
                     END DO
                  ELSE
                     noco_new%alph(itype) = noco%alph(itype)
                  END IF
306
                  iatom= iatom + atoms%neq(n)
307
               END DO
308 309 310 311 312

               OPEN (24,file='nocoinp',form='formatted', status='old')
               REWIND (24)
               CALL rw_noco_write(atoms,jij,noco_new, input)
               CLOSE (24)
313
            END IF
314

315
            IF (noco%l_soc) CALL orbMagMoms(dimension,atoms,noco,clmom)
316
         END IF
317 318 319 320 321 322 323 324
      !block 2 unnecessary for slicing: end
      END IF ! .NOT.sliceplot%slice

      8000 FORMAT (/,/,10x,'spin density at the nucleus:',/,10x,'type',t25,&
                   'input%total',t42,'valence',t65,'core',t90,&
                   'majority valence and input%total density',/)
      8010 FORMAT (i13,2x,3e20.8,5x,2e20.8)

325 326
      IF (sliceplot%slice) THEN
         OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
327
         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, outDen%iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
328
         IF (noco%l_noco) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
329
            WRITE (20) (outDen%pw(k,3),k=1,stars%ng3)
330
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
331 332 333 334 335 336 337
               ALLOCATE(cdomvz(vacuum%nmz,vacuum%nvac))
               DO ivac = 1, vacuum%nvac
                  DO j = 1, vacuum%nmz
                     cdomvz(j,ivac) = CMPLX(outDen%vacz(j,ivac,3),outDen%vacz(j,ivac,4))
                  END DO
               END DO
               WRITE (20) ((cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac)
Gregor Michalicek's avatar
Gregor Michalicek committed
338
               WRITE (20) (((outDen%vacxy(j,k-1,ivac,3),j=1,vacuum%nmzxy),k=2,oneD%odi%nq2) ,ivac=1,vacuum%nvac)
Gregor Michalicek's avatar
Gregor Michalicek committed
339
               DEALLOCATE(cdomvz)
340 341
            END IF
         END IF
342
         CLOSE(20) 
343 344 345 346 347 348 349 350 351 352
         CALL juDFT_end("slice OK")
      END IF

      IF(vacuum%nvac.EQ.1) THEN
         outDen%vacz(:,2,:) = outDen%vacz(:,1,:)
         IF (sym%invs) THEN
            outDen%vacxy(:,:,2,:) = CONJG(outDen%vacxy(:,:,1,:))
         ELSE
            outDen%vacxy(:,:,2,:) = outDen%vacxy(:,:,1,:)
         END IF
353 354
      END IF

355 356 357 358 359 360 361 362 363 364
   ENDIF ! mpi%irank.EQ.0

#ifdef CPP_MPI
   CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,outDen)
#endif

   DEALLOCATE (qvac,qvlay,qa21)
   DEALLOCATE (igq_fft)

   IF (sliceplot%slice) CALL juDFT_end("sliceplot%slice OK",mpi%irank)
365

366
   RETURN
367

368
END SUBROUTINE cdngen
369

370
END MODULE m_cdngen