cdngen.F90 7.66 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
                  archiveType,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
   ! 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
64
   TYPE(t_potden),INTENT(INOUT)     :: outDen
65 66

   !Scalar Arguments
67
   INTEGER, INTENT (IN)             :: eig_id, archiveType
68 69 70 71 72

   ! 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 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
   !pk non-collinear (end)

   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

   !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,&
121
                  sphhar,sym,obsolete,igq_fft,vTot,oneD,coreSpecInput,&
Gregor Michalicek's avatar
Gregor Michalicek committed
122
                  outDen,results,qvac,qvlay,qa21, chmom,clmom)
123 124 125 126 127 128 129 130 131
      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

132 133
   CALL cdncore(results,mpi,dimension,oneD,sliceplot,input,vacuum,noco,sym,&
                stars,cell,sphhar,atoms,vTot,outDen,stdn,svdn)
134 135 136

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

152 153
         !spin densities at the nucleus
         !and magnetic moment in the spheres
154 155 156
         IF (input%jspins.EQ.2) THEN
            WRITE (6,FMT=8000)
            WRITE (16,FMT=8000)
157 158 159
            DO iType = 1,atoms%ntype
               sval = svdn(iType,1) - svdn(iType,input%jspins)
               stot = stdn(iType,1) - stdn(iType,input%jspins)
160
               scor = stot - sval
161 162
               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)
163
            END DO
164

165 166
            noco_new = noco

167
            CALL magMoms(dimension,input,atoms,noco_new,vTot,chmom,qa21)
168

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

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

      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
192 193
      END IF

194 195 196 197 198
      IF (sliceplot%slice) THEN
         CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
                           1,-1.0,0.0,.FALSE.,outDen,'cdn_slice')
      END IF

199 200
   ENDIF ! mpi%irank.EQ.0

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

203 204 205 206 207 208 209 210
#ifdef CPP_MPI
   CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,outDen)
#endif

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

   RETURN
211

212
END SUBROUTINE cdngen
213

214
END MODULE m_cdngen