cdngen.F90 18.3 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 37
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_m_perp
   USE m_types
   USE m_xmlOutput
38
   USE m_orbMagMoms
Daniel Wortmann's avatar
Daniel Wortmann committed
39
#ifdef CPP_MPI
40 41
   USE m_mpi_bc_potden
   USE m_mpi_bc_coreden
Daniel Wortmann's avatar
Daniel Wortmann committed
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 77
   ! 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
   REAL fix,qtot,scor,seig,smom,stot,sval,dummy
78
   REAL sum,fermiEnergyTemp
79 80 81 82 83 84 85 86 87 88 89 90 91
   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(:,:,:,:,:)
   CHARACTER(LEN=20)   :: attributes(4)
92

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

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

   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
117
        
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
   !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
138
                  outDen,results,qvac,qvlay,qa21, chmom,clmom)
139
      CALL timestop("cdngen: cdnval")
140
!-fo
141 142 143 144 145
   END DO

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

146
!+t3e
147
   IF (mpi%irank.EQ.0) THEN
148
!-t3e
149 150 151
      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
152
!---> changes
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
   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
170
#ifdef CPP_MPI
171
         CALL mpi_bc_coreDen(mpi,atoms,input,dimension,rh,tec,qint)
172
#endif
173 174 175 176 177 178 179
      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
180 181
         END IF

182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
         !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
204
                  DO ityp = 1,atoms%ntype
205 206 207 208 209 210 211
                     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
212
                     outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
213
                     !imaginary part rho_21
Gregor Michalicek's avatar
Gregor Michalicek committed
214
                     outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
215 216
                  END DO
               END IF
217 218 219 220 221 222 223 224 225
               !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
Daniel Wortmann's avatar
Daniel Wortmann committed
226
               END DO
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242
            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))
Daniel Wortmann's avatar
Daniel Wortmann committed
243
            END DO
244 245 246 247 248 249 250
         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)
251 252 253 254
            results%seigc = results%seigc + seig
            IF (input%jspins.EQ.2) THEN
               DO jspin = 1,input%jspins
                  DO n = 1,atoms%ntype
255
                     stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
256 257 258
                  END DO
               END DO
            END IF
259
         END IF
260

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

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

316 317
         !spin densities at the nucleus
         !and magnetic moment in the spheres
318 319 320
         IF (input%jspins.EQ.2) THEN
            WRITE (6,FMT=8000)
            WRITE (16,FMT=8000)
321
            DO n = 1,atoms%ntype
322 323 324 325 326
               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)
327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
            END DO
            IF (noco%l_mperp) THEN
               ! angles in nocoinp file are (alph-alphdiff)
               iatom = 1
               DO n = 1,atoms%ntype
                  IF (noco%l_ss) THEN
                     alphdiff(n)= 2.*pi_const*(noco%qss(1)*atoms%taual(1,iatom) + &
                                               noco%qss(2)*atoms%taual(2,iatom) + &
                                               noco%qss(3)*atoms%taual(3,iatom) )
                  ELSE
                     alphdiff(n)= 0.
                  END IF
                  iatom= iatom + atoms%neq(n)
               END DO
            END IF
342 343
            WRITE (6,FMT=8020)
            WRITE (16,FMT=8020)
344
            noco_new = noco
345
            CALL openXMLElement('magneticMomentsInMTSpheres',(/'units'/),(/'muBohr'/))
346
            DO n = 1, atoms%ntype
347 348 349
               smom = chmom(n,1) - chmom(n,input%jspins)
               WRITE (6,FMT=8030) n,smom, (chmom(n,j),j=1,input%jspins)
               WRITE (16,FMT=8030) n,smom, (chmom(n,j),j=1,input%jspins)
350
               attributes = ''
351 352 353 354
               WRITE(attributes(1),'(i0)') n
               WRITE(attributes(2),'(f15.10)') smom
               WRITE(attributes(3),'(f15.10)') chmom(n,1)
               WRITE(attributes(4),'(f15.10)') chmom(n,2)
355 356 357
               CALL writeXMLElementFormPoly('magneticMoment',(/'atomType      ','moment        ','spinUpCharge  ',&
                                                               'spinDownCharge'/),&
                                            attributes,reshape((/8,6,12,14,6,15,15,15/),(/4,2/)))
358
               IF (noco%l_mperp) THEN
359 360 361 362 363 364
                  !calculate the perpendicular part of the local moment
                  !and relax the angle of the local moment or calculate
                  !the constraint B-field.
                  CALL m_perp(atoms,n,noco_new,vTot%mt(:,0,:,:),chmom,qa21,alphdiff)
               END IF
            END DO
365
            CALL closeXMLElement('magneticMomentsInMTSpheres')
366

367 368
            !save the new nocoinp file if the dierctions of the local
            !moments are relaxed or a constraint B-field is calculated.
369 370 371
            l_relax_any = .false.
            iatom = 1
            DO itype = 1,atoms%ntype
372 373
               l_relax_any = l_relax_any.OR.noco%l_relax(itype)
            END DO
374 375 376
            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")
377 378 379 380 381 382 383 384 385 386 387 388 389 390
               END IF
               DO itype = 1, atoms%ntype
                  IF (noco%l_ss) THEN
                     noco_new%alph(itype) = noco%alph(itype) - alphdiff(itype)
                     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
               END DO
391 392 393 394 395

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

398
            IF (noco%l_soc) CALL orbMagMoms(dimension,atoms,noco,clmom)
399
         END IF
400 401 402 403 404 405 406 407 408 409 410
      !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)
      8020 FORMAT (/,/,2x,'-->  magnetic moments in the spheres:',/,2x,&
                   'mm -->   type',t22,'moment',t33,'spin-up',t43,'spin-down')
      8030 FORMAT (2x,'--> mm',i8,2x,3f12.5)

411 412
      IF (sliceplot%slice) THEN
         OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
413
         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, outDen%iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
414
         IF (noco%l_noco) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
415
            WRITE (20) (outDen%pw(k,3),k=1,stars%ng3)
416
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
417 418 419 420 421 422 423
               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
424
               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
425
               DEALLOCATE(cdomvz)
426 427
            END IF
         END IF
428
         CLOSE(20) 
429 430 431 432 433 434 435 436 437 438
         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
439 440
      END IF

441 442 443 444 445 446 447 448 449 450
   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)
451

452
   RETURN
453

454
END SUBROUTINE cdngen
455

456
END MODULE m_cdngen