cdngen.F90 21.6 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 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 38
MODULE m_cdngen

USE m_juDFT

CONTAINS

SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
                  dimension,kpts,atoms,sphhar,stars,sym,obsolete,&
                  enpara,cell,noco,jij,results,oneD,coreSpecInput,&
                  inIter,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_pot_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
39
#ifdef CPP_MPI
40 41 42
   USE m_mpi_bc_pot
   USE m_mpi_bc_potden
   USE m_mpi_bc_coreden
43
#endif
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 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
   IMPLICIT NONE

   ! 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(INOUT)     :: 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,archiveType
   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 :: vz(:,:,:),vr(:,:,:,:)
   REAL   ,ALLOCATABLE :: qvac(:,:,:,:),qvlay(:,:,:,:,:)
   COMPLEX,ALLOCATABLE :: vpw(:,:),vzxy(:,:,:,:)
   CHARACTER(LEN=20)   :: attributes(4)

   !pk non-collinear (start)
   REAL    rhoint,momint,alphdiff(atoms%ntype)
   INTEGER igq2_fft(0:stars%kq1_fft*stars%kq2_fft-1)
   COMPLEX,ALLOCATABLE :: qa21(:)
   !pk non-collinear (end)

   !Read Potential and keep only vr(:,0,:,:) and vz
   ALLOCATE(vpw(stars%ng3,dimension%jspd),vzxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,dimension%jspd))
   ALLOCATE(vz(vacuum%nmzd,2,dimension%jspd),vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd))
   IF (mpi%irank.EQ.0) THEN
      CALL readPotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_TOT_const,&
                         iter,vr,vpw,vz,vzxy)
   END IF
108
#ifdef CPP_MPI
109 110
   CALL mpi_bc_pot(mpi,stars,sphhar,atoms,input,vacuum,&
                   iter,vr,vpw,vz,vzxy)
111
#endif
112
   DEALLOCATE ( vpw,vzxy )
113

114 115 116 117
   iter = inIter
   CALL outDen%init(stars,atoms,sphhar,vacuum,oneD,DIMENSION%jspd,.FALSE.)
   archiveType = CDN_ARCHIVE_TYPE_CDN1_const
   IF(noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
118

119 120 121 122 123 124 125 126 127 128
   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 (outDen%cdom(stars%ng3),outDen%cdomvz(vacuum%nmzd,2))
   ALLOCATE (outDen%cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2))
   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))
129

130 131 132 133 134 135 136 137 138 139 140 141 142 143
   !initialize density arrays with zero
   qa21(:) = cmplx(0.0,0.0)
   qvac(:,:,:,:) = 0.0 
   qvlay(:,:,:,:,:) = 0.0
   outDen%mt(:,:,:,:) = 0.0
   outDen%pw(:,:) = cmplx(0.0,0.0)
   outDen%cdom(:) =  cmplx(0.0,0.0)
   IF (input%film) THEN
      outDen%vacz(:,:,:) = 0.0
      outDen%cdomvz(:,:) = cmplx(0.0,0.0)
      outDen%vacxy(:,:,:,:) = cmplx(0.0,0.0)
      outDen%cdomvxy(:,:,:) = cmplx(0.0,0.0)
   END IF
   outDen%iter = iter
144
        
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
   !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)

   !LDA+U: initialise density-matrix if needed
   IF (atoms%n_u.GT.0) THEN
      ALLOCATE (outDen%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u,input%jspins))
   ELSE
      ALLOCATE (outDen%mmpMat(-lmaxU_const:-lmaxU_const,-lmaxU_const:-lmaxU_const,1,input%jspins))
   END IF
   outDen%mmpMat(:,:,:,:) = CMPLX(0.0,0.0)


   !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,vr,vz(:,:,jspin),oneD,coreSpecInput,&
                  outDen%mmpMat(-lmaxU_const:,-lmaxU_const:,:,jspin),results, outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz,&
                  outDen%cdom,outDen%cdomvz,outDen%cdomvxy,qvac,qvlay,qa21, chmom,clmom)
      CALL timestop("cdngen: cdnval")
177
!-fo
178 179 180 181 182
   END DO

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

183
!+t3e
184
   IF (mpi%irank.EQ.0) THEN
185
!-t3e
186 187 188
      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')
189
!---> changes
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
   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
207
#ifdef CPP_MPI
208
         CALL mpi_bc_coreDen(mpi,atoms,input,dimension,rh,tec,qint)
209
#endif
210 211 212 213 214 215 216
      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
217 218
         END IF

219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
         !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, vr(:,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
241
                  DO ityp = 1,atoms%ntype
242 243 244 245 246 247 248 249 250 251
                     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
                     outDen%cdom(1) = outDen%cdom(1) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
                     !imaginary part rho_21
                     outDen%cdom(1) = outDen%cdom(1) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
252 253
                  END DO
               END IF
254 255 256 257 258 259 260 261 262
               !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
263
               END DO
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279
            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))
280
            END DO
281 282 283 284 285 286
         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
287
            CALL coredr(input,atoms,seig, outDen%mt,dimension,sphhar,vr(:,0,:,:),qint,rh)
288 289 290 291
            results%seigc = results%seigc + seig
            IF (input%jspins.EQ.2) THEN
               DO jspin = 1,input%jspins
                  DO n = 1,atoms%ntype
292
                     stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
293 294 295
                  END DO
               END DO
            END IF
296
         END IF
297

298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318
         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
               outDen%cdom(1) = outDen%cdom(1) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
               !imaginary part rho_21
               outDen%cdom(1) = outDen%cdom(1) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
            END DO
            !pk non-collinear (end)
         ELSE
            DO jspin = 1,input%jspins
               IF (input%ctail) THEN
319
!+gu hope this works as well
320 321 322 323 324 325 326 327 328
                  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
329
         END IF
330 331
      !block 1 unnecessary for slicing: end
      END IF ! .NOT.sliceplot%slice
332
! end relativistic core
333 334 335 336
   END IF ! input%kcrel.EQ.0

   IF (mpi%irank.EQ.0) THEN
      !block 2 unnecessary for slicing: begin
337
      IF (.NOT.sliceplot%slice) THEN
338
         CALL openXMLElementNoAttributes('allElectronCharges')
339 340
         CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
                   outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz,.TRUE.,.true.,fix)
341
         CALL closeXMLElement('allElectronCharges')
342
         !pk non-collinear (start)
343
         IF (noco%l_noco) THEN
344
            !fix also the off-diagonal part of the density matrix
345
            outDen%cdom(:stars%ng3) = fix*outDen%cdom(:stars%ng3)
346 347 348 349
            IF (input%film) THEN
               outDen%cdomvz(:,:) = fix*outDen%cdomvz(:,:)
               outDen%cdomvxy(:,:,:) = fix*outDen%cdomvxy(:,:,:)
            END IF
350
         END IF
351
         !pk non-collinear (end)
352

353 354
         !spin densities at the nucleus
         !and magnetic moment in the spheres
355 356 357
         IF (input%jspins.EQ.2) THEN
            WRITE (6,FMT=8000)
            WRITE (16,FMT=8000)
358
            DO n = 1,atoms%ntype
359 360 361 362 363
               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)
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
            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
379 380
            WRITE (6,FMT=8020)
            WRITE (16,FMT=8020)
381
            noco_new = noco
382
            CALL openXMLElement('magneticMomentsInMTSpheres',(/'units'/),(/'muBohr'/))
383
            DO n = 1, atoms%ntype
384 385 386
               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)
387
               attributes = ''
388 389 390 391
               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)
392 393 394
               CALL writeXMLElementFormPoly('magneticMoment',(/'atomType      ','moment        ','spinUpCharge  ',&
                                                               'spinDownCharge'/),&
                                            attributes,reshape((/8,6,12,14,6,15,15,15/),(/4,2/)))
395
               IF (noco%l_mperp) THEN
396 397 398
                  !calculate the perpendicular part of the local moment
                  !and relax the angle of the local moment or calculate
                  !the constraint B-field.
399
                  CALL m_perp(atoms,n,noco_new,vr(:,0,:,:),chmom,qa21,alphdiff)
400 401
               END IF
            END DO
402
            CALL closeXMLElement('magneticMomentsInMTSpheres')
403

404 405
            !save the new nocoinp file if the dierctions of the local
            !moments are relaxed or a constraint B-field is calculated.
406 407 408
            l_relax_any = .false.
            iatom = 1
            DO itype = 1,atoms%ntype
409 410
               l_relax_any = l_relax_any.OR.noco%l_relax(itype)
            END DO
411 412 413
            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")
414 415 416 417 418 419 420 421 422 423 424 425 426 427
               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
428 429 430 431 432

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

            IF (noco%l_soc) THEN
436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470
               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')
471 472
            END IF
         END IF
473 474 475 476 477 478 479 480 481 482 483 484 485
      !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)

486 487
      IF (sliceplot%slice) THEN
         OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
488
         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
489
         IF (noco%l_noco) THEN
490
            WRITE (20) (outDen%cdom(k),k=1,stars%ng3)
491
            IF (input%film) THEN
492 493
               WRITE (20) ((outDen%cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac)
               WRITE (20) (((outDen%cdomvxy(j,k-1,ivac),j=1,vacuum%nmzxy),k=2,oneD%odi%nq2) ,ivac=1,vacuum%nvac)
494 495
            END IF
         END IF
496
         CLOSE(20) 
497
         CALL juDFT_end("slice OK")
498
      END IF
499

500
      CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
501
                        CDN_OUTPUT_DEN_const,0,results%last_distance,results%ef,.FALSE.,iter,&
502
                        outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy,outDen%cdom,outDen%cdomvz,outDen%cdomvxy)
503 504 505 506 507 508 509 510 511 512
   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)
513

514
   RETURN
515

516
END SUBROUTINE cdngen
517

518
END MODULE m_cdngen