cdngen.F90 8.49 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
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,&
15
                  inDen,outDen)
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

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

   USE m_constants
   USE m_prpqfftmap
   USE m_cdnval
   USE m_cdn_io
   USE m_wrtdop
   USE m_cdntot
   USE m_cdnovlp
   USE m_qfix
31
   USE m_genNewNocoInp
32 33
   USE m_types
   USE m_xmlOutput
34
   USE m_magMoms
35
   USE m_orbMagMoms
36
   USE m_cdncore
Daniel Wortmann's avatar
Daniel Wortmann committed
37
#ifdef CPP_MPI
38
   USE m_mpi_bc_potden
Daniel Wortmann's avatar
Daniel Wortmann committed
39
#endif
40

41
   IMPLICIT NONE
42

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

   ! Local type instances
   TYPE(t_noco) :: noco_new

   !Local Scalars
73
   REAL fix,qtot,scor,stot,sval,dummy
74
   INTEGER ivac,j,jspin,jspmax,k,iType
75
   LOGICAL l_enpara
76 77 78 79 80 81

   !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(:,:,:,:,:)
82

83 84
   !pk non-collinear (start)
   INTEGER igq2_fft(0:stars%kq1_fft*stars%kq2_fft-1)
Gregor Michalicek's avatar
Gregor Michalicek committed
85
   COMPLEX,ALLOCATABLE :: qa21(:), cdomvz(:,:)
86 87
   !pk non-collinear (end)

Gregor Michalicek's avatar
Gregor Michalicek committed
88
   CALL outDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

   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

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

   IF (mpi%irank.EQ.0) THEN
      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')
   END IF ! mpi%irank = 0

136 137
   CALL cdncore(results,mpi,dimension,oneD,sliceplot,input,vacuum,noco,sym,&
                stars,cell,sphhar,atoms,vTot,outDen,stdn,svdn)
138 139 140

   IF (mpi%irank.EQ.0) THEN
      !block 2 unnecessary for slicing: begin
141
      IF (.NOT.sliceplot%slice) THEN
142
         CALL openXMLElementNoAttributes('allElectronCharges')
Gregor Michalicek's avatar
Gregor Michalicek committed
143
         CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,outDen,.TRUE.,.true.,fix)
144
         CALL closeXMLElement('allElectronCharges')
145
         !pk non-collinear (start)
146
         IF (noco%l_noco) THEN
147
            !fix also the off-diagonal part of the density matrix
Gregor Michalicek's avatar
Gregor Michalicek committed
148
            outDen%pw(:stars%ng3,3) = fix*outDen%pw(:stars%ng3,3)
149
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
150
               outDen%vacz(:,:,3:4) = fix*outDen%vacz(:,:,3:4)
Gregor Michalicek's avatar
Gregor Michalicek committed
151
               outDen%vacxy(:,:,:,3) = fix*outDen%vacxy(:,:,:,3)
152
            END IF
153
         END IF
154
         !pk non-collinear (end)
155

156 157
         !spin densities at the nucleus
         !and magnetic moment in the spheres
158 159 160
         IF (input%jspins.EQ.2) THEN
            WRITE (6,FMT=8000)
            WRITE (16,FMT=8000)
161 162 163
            DO iType = 1,atoms%ntype
               sval = svdn(iType,1) - svdn(iType,input%jspins)
               stot = stdn(iType,1) - stdn(iType,input%jspins)
164
               scor = stot - sval
165 166
               WRITE (6,FMT=8010) iType,stot,sval,scor,svdn(iType,1),stdn(iType,1)
               WRITE (16,FMT=8010) iType,stot,sval,scor,svdn(iType,1),stdn(iType,1)
167
            END DO
168

169 170
            noco_new = noco

171
            CALL magMoms(dimension,input,atoms,noco_new,vTot,chmom,qa21)
172

173
            !Generate and save the new nocoinp file if the directions of the local
174
            !moments are relaxed or a constraint B-field is calculated.
175 176
            IF (ANY(noco%l_relax(:atoms%ntype)).OR.noco%l_constr) THEN
               CALL genNewNocoInp(input,atoms,jij,noco,noco_new)
177
            END IF
178

179
            IF (noco%l_soc) CALL orbMagMoms(dimension,atoms,noco,clmom)
180
         END IF
181 182 183 184 185 186 187 188
      !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)

189 190
      IF (sliceplot%slice) THEN
         OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
191
         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, outDen%iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
192
         IF (noco%l_noco) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
193
            WRITE (20) (outDen%pw(k,3),k=1,stars%ng3)
194
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
195 196 197 198 199 200 201
               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
202
               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
203
               DEALLOCATE(cdomvz)
204 205
            END IF
         END IF
206
         CLOSE(20) 
207 208 209 210 211 212 213 214 215 216
         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
217 218
      END IF

219 220 221 222 223 224 225 226 227 228
   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)
229

230
   RETURN
231

232
END SUBROUTINE cdngen
233

234
END MODULE m_cdngen