cdngen.F90 16.5 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 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
   END IF ! mpi%irank = 0
   IF (input%kcrel.EQ.0) THEN
      results%seigc = 0.

      ! 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
169
#ifdef CPP_MPI
170
         CALL mpi_bc_coreDen(mpi,atoms,input,dimension,rh,tec,qint)
171
#endif
172 173 174 175 176 177 178
      END IF

      DO jspin = 1,input%jspins
         IF ((input%jspins.EQ.2).AND.(mpi%irank.EQ.0)) THEN
            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
179 180
         END IF

181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
         !block 1 unnecessary for slicing: begin
         IF (.NOT.sliceplot%slice) THEN

            !add in core density
            IF (mpi%irank.EQ.0) THEN
               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
               IF (input%jspins.EQ.2) THEN
                  DO n = 1,atoms%ntype
                     stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
                  END DO
               END IF
            END IF  ! mpi%irank = 0

            !add core tail charge to outDen%pw
            IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
               !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
               IF (jspin .EQ. 2) THEN
203
                  DO ityp = 1,atoms%ntype
204 205 206 207 208 209 210
                     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
Gregor Michalicek's avatar
Gregor Michalicek committed
211
                     outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
212
                     !imaginary part rho_21
Gregor Michalicek's avatar
Gregor Michalicek committed
213
                     outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
214 215
                  END DO
               END IF
216 217 218 219 220 221 222 223 224
               !pk non-collinear (end)

            ELSE IF (input%ctail) THEN
               CALL cdnovlp(mpi,sphhar,stars,atoms,sym, dimension,vacuum,&
                            cell, input,oneD,l_st, jspin,rh(:,:,jspin),&
                            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
225
               END DO
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
            END IF
         !block 1 unnecessary for slicing: end
         END IF ! .NOT.sliceplot%slice

      END DO ! loop over spins
      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)
   ELSE ! input%kcrel.EQ.0
      !relativistic core implementation : kcrel.eq.1
      results%seigc = 0.
      IF ((input%jspins.EQ.2).AND.(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))
242
            END DO
243 244 245 246 247 248 249
         END DO
      END IF
      !block 1 unnecessary for slicing: begin
      IF (.NOT.sliceplot%slice) THEN
         !add in core density
         IF (mpi%irank.EQ.0) THEN
            CALL coredr(input,atoms,seig, outDen%mt,dimension,sphhar,vTot%mt(:,0,:,:),qint,rh)
250 251 252 253
            results%seigc = results%seigc + seig
            IF (input%jspins.EQ.2) THEN
               DO jspin = 1,input%jspins
                  DO n = 1,atoms%ntype
254
                     stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
255 256 257
                  END DO
               END DO
            END IF
258
         END IF
259

260 261 262 263 264 265 266 267 268 269 270 271 272
         IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
            !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
            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
Gregor Michalicek's avatar
Gregor Michalicek committed
273
               outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
274
               !imaginary part rho_21
Gregor Michalicek's avatar
Gregor Michalicek committed
275
               outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
276 277 278 279 280
            END DO
            !pk non-collinear (end)
         ELSE
            DO jspin = 1,input%jspins
               IF (input%ctail) THEN
281
!+gu hope this works as well
282 283 284 285 286 287 288 289 290
                  CALL cdnovlp(mpi,sphhar,stars,atoms,sym,dimension,vacuum,&
                               cell,input,oneD,l_st,jspin,rh(1,1,jspin),&
                               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
                  END DO
               END IF
            END DO
291
         END IF
292 293
      !block 1 unnecessary for slicing: end
      END IF ! .NOT.sliceplot%slice
294
! end relativistic core
295 296 297 298
   END IF ! input%kcrel.EQ.0

   IF (mpi%irank.EQ.0) THEN
      !block 2 unnecessary for slicing: begin
299
      IF (.NOT.sliceplot%slice) THEN
300
         CALL openXMLElementNoAttributes('allElectronCharges')
301 302
         CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
                   outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz,.TRUE.,.true.,fix)
303
         CALL closeXMLElement('allElectronCharges')
304
         !pk non-collinear (start)
305
         IF (noco%l_noco) THEN
306
            !fix also the off-diagonal part of the density matrix
Gregor Michalicek's avatar
Gregor Michalicek committed
307
            outDen%pw(:stars%ng3,3) = fix*outDen%pw(:stars%ng3,3)
308
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
309
               outDen%vacz(:,:,3:4) = fix*outDen%vacz(:,:,3:4)
Gregor Michalicek's avatar
Gregor Michalicek committed
310
               outDen%vacxy(:,:,:,3) = fix*outDen%vacxy(:,:,:,3)
311
            END IF
312
         END IF
313
         !pk non-collinear (end)
314

315 316
         !spin densities at the nucleus
         !and magnetic moment in the spheres
317 318 319
         IF (input%jspins.EQ.2) THEN
            WRITE (6,FMT=8000)
            WRITE (16,FMT=8000)
320
            DO n = 1,atoms%ntype
321 322 323 324 325
               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)
326
            END DO
327

328 329
            noco_new = noco

330
            CALL magMoms(dimension,input,atoms,noco_new,vTot,chmom,qa21)
331

332 333
            !save the new nocoinp file if the dierctions of the local
            !moments are relaxed or a constraint B-field is calculated.
334 335
            l_relax_any = .false.
            DO itype = 1,atoms%ntype
336 337
               l_relax_any = l_relax_any.OR.noco%l_relax(itype)
            END DO
338 339 340
            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")
341
               END IF
342
               iatom = 1
343 344
               DO itype = 1, atoms%ntype
                  IF (noco%l_ss) THEN
345 346 347 348
                     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
349 350 351 352 353 354 355 356 357
                     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
358
                  iatom= iatom + atoms%neq(n)
359
               END DO
360 361 362 363 364

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

367
            IF (noco%l_soc) CALL orbMagMoms(dimension,atoms,noco,clmom)
368
         END IF
369 370 371 372 373 374 375 376
      !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)

377 378
      IF (sliceplot%slice) THEN
         OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
379
         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, outDen%iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
380
         IF (noco%l_noco) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
381
            WRITE (20) (outDen%pw(k,3),k=1,stars%ng3)
382
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
383 384 385 386 387 388 389
               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
390
               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
391
               DEALLOCATE(cdomvz)
392 393
            END IF
         END IF
394
         CLOSE(20) 
395 396 397 398 399 400 401 402 403 404
         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
405 406
      END IF

407 408 409 410 411 412 413 414 415 416
   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)
417

418
   RETURN
419

420
END SUBROUTINE cdngen
421

422
END MODULE m_cdngen