cdngen.F90 20.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
#ifdef CPP_MPI
39 40
   USE m_mpi_bc_potden
   USE m_mpi_bc_coreden
41
#endif
42

43
   IMPLICIT NONE
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 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
   ! 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
   REAL slmom,slxmom,slymom,sum,thetai,phii,fermiEnergyTemp
   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)
91

92 93 94
   !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
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 327 328 329 330 331 332 333 334 335 336 337 338 339 340
            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
341 342
            WRITE (6,FMT=8020)
            WRITE (16,FMT=8020)
343
            noco_new = noco
344
            CALL openXMLElement('magneticMomentsInMTSpheres',(/'units'/),(/'muBohr'/))
345
            DO n = 1, atoms%ntype
346 347 348
               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)
349
               attributes = ''
350 351 352 353
               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)
354 355 356
               CALL writeXMLElementFormPoly('magneticMoment',(/'atomType      ','moment        ','spinUpCharge  ',&
                                                               'spinDownCharge'/),&
                                            attributes,reshape((/8,6,12,14,6,15,15,15/),(/4,2/)))
357
               IF (noco%l_mperp) THEN
358 359 360 361 362 363
                  !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
364
            CALL closeXMLElement('magneticMomentsInMTSpheres')
365

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

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

            IF (noco%l_soc) THEN
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432
               thetai = noco%theta
               phii   = noco%phi
               WRITE (6,FMT=9020)
               WRITE (16,FMT=9020)
               CALL openXMLElement('orbitalMagneticMomentsInMTSpheres',(/'units'/),(/'muBohr'/))
               DO n = 1, atoms%ntype
                  IF (noco%l_noco) THEN
                     thetai = noco%beta(n)
                     phii =   noco%alph(n)
                  END IF

                  ! magn. moment(-)
                  slxmom = clmom(1,n,1)+clmom(1,n,2)
                  slymom = clmom(2,n,1)+clmom(2,n,2)
                  slmom =  clmom(3,n,1)+clmom(3,n,2)

                  ! rotation: orbital moment || spin moment (extended to incude phi - hopefully)
                  slmom   = cos(thetai)*slmom + sin(thetai)* (cos(phii)*slxmom + sin(phii)*slymom)
                  clmom(3,n,1) = cos(thetai)*clmom(3,n,1) + &
                                 sin(thetai)*(cos(phii)*clmom(1,n,1) + sin(phii)*clmom(2,n,1))
                  clmom(3,n,2) = cos(thetai)*clmom(3,n,2) + &
                                 sin(thetai)*(cos(phii)*clmom(1,n,2) + sin(phii)*clmom(2,n,2))

                  WRITE (6,FMT=8030) n,slmom,(clmom(3,n,j),j=1,2)
                  WRITE (16,FMT=8030) n,slmom,(clmom(3,n,j),j=1,2)
                  attributes = ''
                  WRITE(attributes(1),'(i0)') n
                  WRITE(attributes(2),'(f15.10)') slmom
                  WRITE(attributes(3),'(f15.10)') clmom(3,n,1)
                  WRITE(attributes(4),'(f15.10)') clmom(3,n,2)
                  CALL writeXMLElementFormPoly('orbMagMoment',(/'atomType      ','moment        ','spinUpCharge  ',&
                                                                'spinDownCharge'/),&
                                               attributes,reshape((/8,6,12,14,6,15,15,15/),(/4,2/)))
               END DO
               CALL closeXMLElement('orbitalMagneticMomentsInMTSpheres')
433 434
            END IF
         END IF
435 436 437 438 439 440 441 442 443 444 445 446 447
      !block 2 unnecessary for slicing: end
      END IF ! .NOT.sliceplot%slice

      9020 FORMAT (/,/,10x,'orb. magnetic moments in the spheres:',/,10x,&
                   'type',t22,'moment',t33,'spin-up',t43,'spin-down')
      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)

448 449
      IF (sliceplot%slice) THEN
         OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
450
         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, outDen%iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
451
         IF (noco%l_noco) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
452
            WRITE (20) (outDen%pw(k,3),k=1,stars%ng3)
453
            IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
454 455 456 457 458 459 460
               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
461
               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
462
               DEALLOCATE(cdomvz)
463 464
            END IF
         END IF
465
         CLOSE(20) 
466 467 468 469 470 471 472 473 474 475
         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
476 477
      END IF

478 479 480 481 482 483 484 485 486 487
   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)
488

489
   RETURN
490

491
END SUBROUTINE cdngen
492

493
END MODULE m_cdngen