cdngen.F90 10.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_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
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_types
   USE m_xmlOutput
35
   USE m_magMoms
36
   USE m_orbMagMoms
37
   USE m_cdncore
Daniel Wortmann's avatar
Daniel Wortmann committed
38
#ifdef CPP_MPI
39
   USE m_mpi_bc_potden
Daniel Wortmann's avatar
Daniel Wortmann committed
40
#endif
41

42
   IMPLICIT NONE
43

44 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
   ! 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
75 76 77 78
   REAL fix,qtot,scor,stot,sval,dummy
   INTEGER iter,ivac,j,jspin,jspmax,k,n
   INTEGER itype,iatom
   LOGICAL l_relax_any,exst, l_enpara
79 80 81 82 83 84

   !Local Arrays
   REAL stdn(atoms%ntype,dimension%jspd),svdn(atoms%ntype,dimension%jspd),alpha_l(atoms%ntype)
   REAL chmom(atoms%ntype,dimension%jspd),clmom(3,atoms%ntype,dimension%jspd)
   INTEGER,ALLOCATABLE :: igq_fft(:)
   REAL   ,ALLOCATABLE :: qvac(:,:,:,:),qvlay(:,:,:,:,:)
85

86
   !pk non-collinear (start)
87
   REAL    alphdiff
88
   INTEGER igq2_fft(0:stars%kq1_fft*stars%kq2_fft-1)
Gregor Michalicek's avatar
Gregor Michalicek committed
89
   COMPLEX,ALLOCATABLE :: qa21(:), cdomvz(:,:)
90 91 92
   !pk non-collinear (end)

   iter = inIter
Gregor Michalicek's avatar
Gregor Michalicek committed
93
   CALL outDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109

   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
Daniel Wortmann's avatar
Daniel Wortmann committed
110
        
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
   !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,&
130
                  sphhar,sym,obsolete,igq_fft,vTot,oneD,coreSpecInput,&
Gregor Michalicek's avatar
Gregor Michalicek committed
131
                  outDen,results,qvac,qvlay,qa21, chmom,clmom)
132
      CALL timestop("cdngen: cdnval")
133
!-fo
134 135 136 137 138
   END DO

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

139
!+t3e
140
   IF (mpi%irank.EQ.0) THEN
141
!-t3e
142 143 144
      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')
Daniel Wortmann's avatar
Daniel Wortmann committed
145
!---> changes
146 147
   END IF ! mpi%irank = 0

148 149
   CALL cdncore(results,mpi,dimension,oneD,sliceplot,input,vacuum,noco,sym,&
                stars,cell,sphhar,atoms,vTot,outDen,stdn,svdn)
150 151 152

   IF (mpi%irank.EQ.0) THEN
      !block 2 unnecessary for slicing: begin
153
      IF (.NOT.sliceplot%slice) THEN
154
         CALL openXMLElementNoAttributes('allElectronCharges')
155 156
         CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
                   outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz,.TRUE.,.true.,fix)
157
         CALL closeXMLElement('allElectronCharges')
158
         !pk non-collinear (start)
159
         IF (noco%l_noco) THEN
160
            !fix also the off-diagonal part of the density matrix
Gregor Michalicek's avatar
Gregor Michalicek committed
161
            outDen%pw(:stars%ng3,3) = fix*outDen%pw(:stars%ng3,3)
162
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
163
               outDen%vacz(:,:,3:4) = fix*outDen%vacz(:,:,3:4)
Gregor Michalicek's avatar
Gregor Michalicek committed
164
               outDen%vacxy(:,:,:,3) = fix*outDen%vacxy(:,:,:,3)
165
            END IF
166
         END IF
167
         !pk non-collinear (end)
168

169 170
         !spin densities at the nucleus
         !and magnetic moment in the spheres
171 172 173
         IF (input%jspins.EQ.2) THEN
            WRITE (6,FMT=8000)
            WRITE (16,FMT=8000)
174
            DO n = 1,atoms%ntype
175 176 177 178 179
               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)
180
            END DO
181

182 183
            noco_new = noco

184
            CALL magMoms(dimension,input,atoms,noco_new,vTot,chmom,qa21)
185

186 187
            !save the new nocoinp file if the dierctions of the local
            !moments are relaxed or a constraint B-field is calculated.
188 189
            l_relax_any = .false.
            DO itype = 1,atoms%ntype
190 191
               l_relax_any = l_relax_any.OR.noco%l_relax(itype)
            END DO
192 193 194
            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")
195
               END IF
196
               iatom = 1
197 198
               DO itype = 1, atoms%ntype
                  IF (noco%l_ss) THEN
199 200 201 202
                     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
203 204 205 206 207 208 209 210 211
                     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
212
                  iatom= iatom + atoms%neq(n)
213
               END DO
214 215 216 217 218

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

221
            IF (noco%l_soc) CALL orbMagMoms(dimension,atoms,noco,clmom)
222
         END IF
223 224 225 226 227 228 229 230
      !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)

231 232
      IF (sliceplot%slice) THEN
         OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
233
         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, outDen%iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
234
         IF (noco%l_noco) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
235
            WRITE (20) (outDen%pw(k,3),k=1,stars%ng3)
236
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
237 238 239 240 241 242 243
               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
244
               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
245
               DEALLOCATE(cdomvz)
246 247
            END IF
         END IF
248
         CLOSE(20) 
249 250 251 252 253 254 255 256 257 258
         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
259 260
      END IF

261 262 263 264 265 266 267 268 269 270
   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)
271

272
   RETURN
273

274
END SUBROUTINE cdngen
275

276
END MODULE m_cdngen