MODULE m_cdnval use m_juDFT CONTAINS SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,sliceplot,noco, input,banddos,cell,atoms,enpara,stars,& vacuum,dimension,sphhar,sym,obsolete,igq_fft,vr,vz,oneD,coreSpecInput,n_mmp,results, qpw,rhtxy,& rho,rht,cdom,cdomvz,cdomvxy,qvac,qvlay,qa21, chmom,clmom) ! ! *********************************************************** ! this subroutin is a modified version of cdnval.F. ! it calculates a layer charge distribution and an orbital ! composition for each state in a film. ! this information is written in the 'ek_orbcomp' file ! Yu.Koroteev 01.12.2003 ! *********************************************************** ! flapw7 valence density generator ! c.l.fu ! zelec used to calculate ne - 6.12.95 r.pentcheva ! ! changed subroutine to allow parallel writing of vacdos&dosinp ! used temporary direct access file 84,tmp_dos to store data used ! in cdninf ! call of cdninf changed, sympsi is called from cdnval now ! look for 'ifdef CPP_MPI' -blocks! d.wortmann 6.5.99 ! !******** ABBREVIATIONS ************************************************ ! nbands : number of bands in the energy window ! noccbd : number of occupied bands ! slice : when set to .true. the charge density of a enery range ! (slice) or a single state is calculated ! e1s,e2s : (used if slice) energy range for the slice. if both ! are set to 0.0 the charge density of the band nr. nnne ! is calculated ! pallst : (used if slice) if set to .true. bands above the ! Fermi-Energy are taken into account ! nnne : (used if slice) number of the band of which the charge ! density is calculated if e1s and e2s are 0.0 ! kk : (used if slice) if set to 0 all k-points are used to ! calculate the charge density of the slice, otherwise ! only k-points kk is taken into account ! nslibd : number of bands in slice ! ener : band energy averaged over all bands and k-points, ! wheighted with the l-like charge of each atom type ! sqal : l-like charge of each atom type. sum over all k-points ! and bands !*********************************************************************** ! USE m_eig66_io,ONLY: write_dos USE m_radfun USE m_radflo USE m_rhomt USE m_rhonmt USE m_rhomtlo USE m_rhonmtlo USE m_mcdinit USE m_sphpts USE m_points USE m_sympsi USE m_enpara, ONLY : w_enpara,mix_enpara USE m_eparas ! energy parameters and partial charges USE m_qal21 ! off-diagonal part of partial charges USE m_abcof USE m_topulay USE m_nmat ! calculate density matrix for LDA + U USE m_vacden USE m_nstm3 USE m_pwden USE m_forcea8 USE m_forcea12 USE m_forcea21 USE m_checkdop ! check continuity of density on MT radius R USE m_int21 ! integrate (spin) off-diagonal radial functions USE m_int21lo ! -"- for u_lo USE m_rhomt21 ! calculate (spin) off-diagonal MT-density coeff's USE m_rhonmt21 ! -"- non-MT-density coeff's USE m_cdnmt ! calculate the density and orbital moments etc. USE m_orbmom ! coeffd for orbital moments USE m_qmtsl ! These subroutines divide the input%film into vacuum%layers USE m_qintsl ! (slabs) and intergate the DOS in these vacuum%layers USE m_slabdim ! (mt + interstitial) USE m_slabgeom ! (written by Yu.Koroteev, 2003/2004) USE m_orbcomp ! calculate corbital composition (like p_x,p_y,p_z) USE m_Ekwritesl ! and write to file. USE m_abcrot2 USE m_doswrite USE m_cylpts USE m_cdnread, ONLY : cdn_read0, cdn_read USE m_corespec, only : l_cs ! calculation of core spectra (EELS) USE m_corespec_io, only : corespec_init USE m_corespec_eval, only : corespec_gaunt,corespec_rme,corespec_dos,corespec_ddscs #ifdef CPP_MPI USE m_mpi_col_den ! collect density data from parallel nodes #endif USE m_types USE m_xmlOutput IMPLICIT NONE 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_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 ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: eig_id,jspin ! .. Array Arguments .. COMPLEX, INTENT(INOUT) :: qpw(stars%ng3,dimension%jspd) COMPLEX, INTENT(INOUT) :: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,dimension%jspd) COMPLEX, INTENT(INOUT) :: cdom(stars%ng3) COMPLEX, INTENT(INOUT) :: cdomvz(vacuum%nmzd,2) COMPLEX, INTENT(INOUT) :: cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2) COMPLEX, INTENT(INOUT) :: qa21(atoms%ntype) INTEGER, INTENT (IN) :: igq_fft(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1) REAL, INTENT (IN) :: vz(vacuum%nmzd,2) REAL, INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) REAL, INTENT (OUT) :: chmom(atoms%ntype,dimension%jspd),clmom(3,atoms%ntype,dimension%jspd) REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) REAL, INTENT (INOUT) :: rht(vacuum%nmzd,2,dimension%jspd) REAL, INTENT (INOUT) :: qvac(dimension%neigd,2,kpts%nkpt,dimension%jspd) REAL, INTENT (INOUT) :: qvlay(dimension%neigd,vacuum%layerd,2,kpts%nkpt,dimension%jspd) COMPLEX, INTENT(INOUT) :: n_mmp(-3:3,-3:3,atoms%n_u) #ifdef CPP_MPI INCLUDE 'mpif.h' logical :: mpi_flag, mpi_status #endif ! .. Local Scalars .. TYPE(t_lapw):: lapw INTEGER :: llpd REAL wk,wronk,sign,emcd_lo,emcd_up INTEGER i,ie,iv,ivac,j,k,l,l1,lh ,n,ilo,isp,nat,& nbands,noded,nodeu,noccbd,nslibd,na,& ikpt,npd ,jsp_start,jsp_end,ispin INTEGER skip_t,skip_tt INTEGER n_size,i_rec,n_rank ,ncored,n_start,n_end,noccbd_l COMPLEX,parameter:: czero=(0.0,0.0) LOGICAL l_fmpl,l_mcd,l_evp,l_orbcomprot ! ...Local Arrays .. INTEGER n_bands(0:dimension%neigd),ncore(atoms%ntype) REAL cartk(3),xp(3,dimension%nspd),e_mcd(atoms%ntype,input%jspins,dimension%nstd) REAL ello(atoms%nlod,atoms%ntype,dimension%jspd),evac(2,dimension%jspd) REAL epar(0:atoms%lmaxd,atoms%ntype,dimension%jspd),evdu(2,dimension%jspd) REAL eig(dimension%neigd) REAL vz0(2) REAL uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype) REAL ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype) INTEGER, PARAMETER :: n2max_nstm3=13 INTEGER nsld,nsl ! INTEGER, ALLOCATABLE :: nmtsl(:,:),nslat(:,:) REAL, ALLOCATABLE :: zsl(:,:),volsl(:) REAL, ALLOCATABLE :: volintsl(:) REAL, ALLOCATABLE :: qintsl(:,:),qmtsl(:,:) REAL, ALLOCATABLE :: orbcomp(:,:,:),qmtp(:,:) REAL, ALLOCATABLE :: qis(:,:,:) !-new_sl !-dw INTEGER, ALLOCATABLE :: gvac1d(:),gvac2d(:) INTEGER, ALLOCATABLE :: jsym(:),ksym(:) REAL, ALLOCATABLE :: aclo(:,:,:),acnmt(:,:,:,:,:) REAL, ALLOCATABLE :: bclo(:,:,:),bcnmt(:,:,:,:,:) REAL, ALLOCATABLE :: cclo(:,:,:,:),ccnmt(:,:,:,:,:),we(:) REAL, ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) REAL, ALLOCATABLE :: uloulopn(:,:,:,:),uloulopn21(:,:,:) REAL, ALLOCATABLE :: uu(:,:,:),dd(:,:,:),du(:,:,:) REAL, ALLOCATABLE :: uunmt(:,:,:,:),ddnmt(:,:,:,:) REAL, ALLOCATABLE :: dunmt(:,:,:,:),udnmt(:,:,:,:),sqlo(:,:,:) REAL, ALLOCATABLE :: qal(:,:,:,:),sqal(:,:,:),ener(:,:,:) REAL, ALLOCATABLE :: svac(:,:),pvac(:,:),mcd(:,:,:) REAL, ALLOCATABLE :: enerlo(:,:,:),qmat(:,:,:,:) COMPLEX, ALLOCATABLE :: acof(:,:,:,:),bcof(:,:,:,:),ccof(:,:,:,:,:) COMPLEX, ALLOCATABLE :: acoflo(:,:,:,:),bcoflo(:,:,:,:) COMPLEX, ALLOCATABLE :: cveccof(:,:,:,:,:),f_a12(:,:) COMPLEX, ALLOCATABLE :: e1cof(:,:,:),e2cof(:,:,:),f_a21(:,:) COMPLEX, ALLOCATABLE :: f_b4(:,:),f_b8(:,:) COMPLEX, ALLOCATABLE :: aveccof(:,:,:,:),bveccof(:,:,:,:) COMPLEX, ALLOCATABLE :: uloulop21(:,:,:) COMPLEX, ALLOCATABLE :: uunmt21(:,:,:),ddnmt21(:,:,:) COMPLEX, ALLOCATABLE :: dunmt21(:,:,:),udnmt21(:,:,:) COMPLEX, ALLOCATABLE :: qstars(:,:,:,:),m_mcd(:,:,:,:) TYPE (t_orb), ALLOCATABLE :: orb(:,:,:,:) TYPE (t_orbl), ALLOCATABLE :: orbl(:,:,:,:) TYPE (t_orblo),ALLOCATABLE :: orblo(:,:,:,:,:) TYPE (t_mt21), ALLOCATABLE :: mt21(:,:) TYPE (t_lo21), ALLOCATABLE :: lo21(:,:) TYPE (t_usdus) :: usdus TYPE (t_zMat) :: zMat INTEGER :: nkpt_extended zmat%l_real=sym%invs.AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco) ! .. ! .. llpd=(atoms%lmaxd*(atoms%lmaxd+3))/2 !---> l_fmpl is meant as a switch to to a plot of the full magnet. !---> density without the atomic sphere approximation for the magnet. !---> density. It is not completely implemented (lo's missing). l_fmpl = .false. IF (noco%l_mperp) THEN !---> when the off-diag. part of the desinsity matrix, i.e. m_x and !---> m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval !---> is called only once. therefore, several spin loops have been !---> added. if l_mperp = F, these loops run only from jspin - jspin. jsp_start = 1 jsp_end = 2 ALLOCATE ( mt21(0:atoms%lmaxd,atoms%ntype),lo21(atoms%nlod,atoms%ntype) ) ! Deallocation at end of subroutine ALLOCATE ( uloulopn21(atoms%nlod,atoms%nlod,atoms%ntype) ) ALLOCATE ( uloulop21(atoms%nlod,atoms%nlod,atoms%ntype) ) ALLOCATE ( qmat(0:3,atoms%ntype,dimension%neigd,4) ) IF (l_fmpl) THEN ALLOCATE ( uunmt21((atoms%lmaxd+1)**2,sphhar%nlhd,atoms%ntype) ) ALLOCATE ( ddnmt21((atoms%lmaxd+1)**2,sphhar%nlhd,atoms%ntype) ) ALLOCATE ( dunmt21((atoms%lmaxd+1)**2,sphhar%nlhd,atoms%ntype) ) ALLOCATE ( udnmt21((atoms%lmaxd+1)**2,sphhar%nlhd,atoms%ntype) ) ELSE ALLOCATE ( uunmt21(1,1,1),ddnmt21(1,1,1) ) ALLOCATE ( dunmt21(1,1,1),udnmt21(1,1,1) ) ENDIF ELSE jsp_start = jspin jsp_end = jspin ALLOCATE ( mt21(1,1),lo21(1,1),uunmt21(1,1,1) ) ALLOCATE ( ddnmt21(1,1,1),dunmt21(1,1,1),udnmt21(1,1,1) ) ALLOCATE ( uloulopn21(1,1,1),uloulop21(1,1,1),qmat(1,1,1,1) ) ENDIF ! !---> if l_mperp = F, these variables are only needed for one spin !---> at a time, otherwise for both spins: ! ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) ) ! Deallocation before mpi_col_den ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) ) ALLOCATE ( usdus%us(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%uds(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%dus(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%duds(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%ddn(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( jsym(dimension%neigd),ksym(dimension%neigd) ) ALLOCATE ( gvac1d(dimension%nv2d),gvac2d(dimension%nv2d) ) ALLOCATE ( usdus%ulos(atoms%nlod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%dulos(atoms%nlod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%uulon(atoms%nlod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%dulon(atoms%nlod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( usdus%uloulopn(atoms%nlod,atoms%nlod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( uu(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( dd(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( du(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( uunmt(0:llpd,sphhar%nlhd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( ddnmt(0:llpd,sphhar%nlhd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( dunmt(0:llpd,sphhar%nlhd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( udnmt(0:llpd,sphhar%nlhd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( qal(0:3,atoms%ntype,dimension%neigd,jsp_start:jsp_end) ) ALLOCATE ( sqal(0:3,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( ener(0:3,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( sqlo(atoms%nlod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( enerlo(atoms%nlod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( svac(2,jsp_start:jsp_end) ) ALLOCATE ( pvac(2,jsp_start:jsp_end) ) ALLOCATE ( qstars(vacuum%nstars,dimension%neigd,vacuum%layerd,2) ) ! ! --> Initializations ! uu(:,:,:) = 0.0 ; dd(:,:,:) = 0.0 ; du(:,:,:) = 0.0 IF (noco%l_mperp) THEN mt21(:,:)%uu = czero ; mt21(:,:)%ud = czero mt21(:,:)%du = czero ; mt21(:,:)%dd = czero lo21(:,:)%uulo = czero ; lo21(:,:)%ulou = czero lo21(:,:)%dulo = czero ; lo21(:,:)%ulod = czero uloulop21(:,:,:) = czero ENDIF uunmt(:,:,:,:) = 0.0 ; ddnmt(:,:,:,:) = 0.0 udnmt(:,:,:,:) = 0.0 ; dunmt(:,:,:,:) = 0.0 IF (l_fmpl) THEN IF (.not.noco%l_mperp) CALL juDFT_error("for fmpl set noco%l_mperp = T!" ,calledby ="cdnval") uunmt21(:,:,:) = czero ; ddnmt21(:,:,:) = czero udnmt21(:,:,:) = czero ; dunmt21(:,:,:) = czero ENDIF svac(:,:) = 0.0 ; pvac(:,:) = 0.0 sqal(:,:,:) = 0.0 ; ener(:,:,:) = 0.0 !+soc IF (noco%l_soc) THEN ALLOCATE ( orb(0:atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( orbl(atoms%nlod,-atoms%llod:atoms%llod,atoms%ntype,jsp_start:jsp_end) ) ALLOCATE ( orblo(atoms%nlod,atoms%nlod,-atoms%llod:atoms%llod,atoms%ntype,jsp_start:jsp_end)) orb(:,:,:,:)%uu = 0.0 ; orb(:,:,:,:)%dd = 0.0 orb(:,:,:,:)%uum = czero ; orb(:,:,:,:)%uup = czero orb(:,:,:,:)%ddm = czero ; orb(:,:,:,:)%ddp = czero orbl(:,:,:,:)%uulo = 0.0 ; orbl(:,:,:,:)%dulo = 0.0 orbl(:,:,:,:)%uulom = czero ; orbl(:,:,:,:)%uulop = czero orbl(:,:,:,:)%dulom = czero ; orbl(:,:,:,:)%dulop = czero orblo(:,:,:,:,:)%z = 0.0 orblo(:,:,:,:,:)%p = czero ; orblo(:,:,:,:,:)%m = czero ELSE ALLOCATE ( orb(0:0,-atoms%lmaxd:-atoms%lmaxd,1,jsp_start:jsp_end) ) ALLOCATE ( orbl(1,-atoms%llod:-atoms%llod,1,jsp_start:jsp_end) ) ALLOCATE ( orblo(1,1,-atoms%llod:-atoms%llod,1,jsp_start:jsp_end) ) ENDIF !+for IF (input%l_f) THEN ALLOCATE ( f_a12(3,atoms%ntype),f_a21(3,atoms%ntype) ) ALLOCATE ( f_b4(3,atoms%ntype),f_b8(3,atoms%ntype) ) f_b4(:,:) = czero ; f_a12(:,:) = czero f_b8(:,:) = czero ; f_a21(:,:) = czero ELSE ALLOCATE ( f_b8(1,1) ) ENDIF ! INQUIRE (file='mcd_inp',exist=l_mcd) IF (l_mcd) THEN OPEN (23,file='mcd_inp',STATUS='old',FORM='formatted') READ (23,*) emcd_lo,emcd_up CLOSE (23) ALLOCATE ( m_mcd(dimension%nstd,(3+1)**2,3*atoms%ntype,2) ) ALLOCATE ( mcd(3*atoms%ntype,dimension%nstd,dimension%neigd) ) IF (.not.banddos%dos) WRITE (*,*) 'For mcd-spectra set banddos%dos=T!' ELSE ALLOCATE ( m_mcd(1,1,1,1),mcd(1,1,1) ) ENDIF ! calculation of core spectra (EELS) initializations -start- CALL corespec_init(input,atoms,coreSpecInput) IF(l_cs.AND.(mpi%isize.NE.1)) CALL juDFT_error('EELS + MPI not implemented', calledby = 'cdnval') IF(l_cs.AND.jspin.EQ.1) CALL corespec_gaunt() ! calculation of core spectra (EELS) initializations -end- IF (mpi%irank==0) THEN WRITE (6,FMT=8000) jspin WRITE (16,FMT=8000) jspin CALL openXMLElementPoly('mtCharges',(/'spin'/),(/jspin/)) END IF 8000 FORMAT (/,/,10x,'valence density: spin=',i2) CALL cdn_read0(& eig_id,& mpi%irank,mpi%isize,jspin,dimension%jspd,& noco%l_noco,& ello,evac,epar,wk,n_bands,n_size) #ifdef CPP_MPI ! Sinchronizes the RMA operations CALL MPI_BARRIER(mpi%mpi_comm,ie) #endif !+lo !---> if local orbitals are used, the eigenvector has a higher !---> dimension then nvd ALLOCATE ( aclo(atoms%nlod,atoms%ntype,jsp_start:jsp_end), & bclo(atoms%nlod,atoms%ntype,jsp_start:jsp_end),& cclo(atoms%nlod,atoms%nlod,atoms%ntype,jsp_start:jsp_end),& acnmt(0:atoms%lmaxd,atoms%nlod,sphhar%nlhd,atoms%ntype,jsp_start:jsp_end), & bcnmt(0:atoms%lmaxd,atoms%nlod,sphhar%nlhd,atoms%ntype,jsp_start:jsp_end), & ccnmt(atoms%nlod,atoms%nlod,sphhar%nlhd,atoms%ntype,jsp_start:jsp_end) ) aclo(:,:,:) = 0.0 ; bclo(:,:,:) = 0.0 ; ccnmt(:,:,:,:,:) = 0.0 acnmt(:,:,:,:,:)=0.0 ; bcnmt(:,:,:,:,:)=0.0 ; cclo(:,:,:,:)=0.0 ALLOCATE ( qis(dimension%neigd,kpts%nkpt,dimension%jspd)) skip_tt = dot_product(enpara%skiplo(:atoms%ntype,jspin),atoms%neq(:atoms%ntype)) IF (noco%l_soc.OR.noco%l_noco) skip_tt = 2 * skip_tt !-lo !---> set up l-dependent m.t. wavefunctions na = 1 ncored = 0 ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod,dimension%jspd) ) DO n = 1,atoms%ntype IF (input%cdinf.AND.mpi%irank==0) WRITE (6,FMT=8001) n DO l = 0,atoms%lmax(n) DO ispin =jsp_start,jsp_end CALL radfun(& l,n,ispin,epar(l,n,ispin),vr(1,0,n,ispin),atoms,& f(1,1,l,ispin),g(1,1,l,ispin),usdus,& nodeu,noded,wronk) IF (input%cdinf.AND.mpi%irank==0) WRITE (6,FMT=8002) l,& epar(l,n,ispin),usdus%us(l,n,ispin),usdus%dus(l,n,ispin),nodeu,& usdus%uds(l,n,ispin),usdus%duds(l,n,ispin),noded,usdus%ddn(l,n,ispin),& wronk END DO IF (noco%l_mperp) THEN CALL int_21(& f,g,atoms,n,l,& mt21(l,n)%uun,mt21(l,n)%udn,& mt21(l,n)%dun,mt21(l,n)%ddn) END IF END DO IF (l_mcd) THEN CALL mcd_init(& atoms,input,dimension,& vr(:,0,:,:),g,f,emcd_up,emcd_lo,n,jspin,& ncore,e_mcd,m_mcd) ncored = max(ncore(n),ncored) END IF IF(l_cs) CALL corespec_rme(atoms,input,n,dimension%nstd,& input%jspins,jspin,results%ef,& dimension%msh,vr(:,0,:,:),f,g) ! !---> generate the extra wavefunctions for the local orbitals, !---> if there are any. ! IF ( atoms%nlo(n) > 0 ) THEN DO ispin = jsp_start,jsp_end CALL radflo(atoms,n,ispin, ello(1,1,ispin),vr(:,0,n,ispin), f(1,1,0,ispin),& g(1,1,0,ispin),mpi, usdus, uuilon,duilon,ulouilopn, flo(:,:,:,ispin)) END DO END IF DO ilo = 1, atoms%nlo(n) IF (noco%l_mperp) THEN CALL int_21lo(f,g,atoms,n, flo,ilo,& lo21(ilo,n)%uulon,lo21(ilo,n)%dulon,& lo21(ilo,n)%uloun,lo21(ilo,n)%ulodn,& uloulopn21(1,1,n)) END IF END DO na = na + atoms%neq(n) END DO DEALLOCATE (flo) 8001 FORMAT (1x,/,/,' wavefunction parameters for atom type',i3,':',/,& t32,'radial function',t79,'energy derivative',/,t3,'l',t8,& 'energy',t26,'value',t39,'derivative',t53,'nodes',t68,& 'value',t81,'derivative',t95,'nodes',t107,'norm',t119,& 'wronskian') 8002 FORMAT (i3,f10.5,2 (5x,1p,2e16.7,i5),1p,2e16.7) IF (input%film) vz0(:) = vz(vacuum%nmz,:) nsld=1 !+q_sl IF ((banddos%ndir.EQ.-3).AND.banddos%dos) THEN IF (oneD%odi%d1) CALL juDFT_error("layer-resolved feature does not work with 1D",calledby ="cdnval") CALL slab_dim(atoms, nsld) ALLOCATE ( nmtsl(atoms%ntype,nsld),nslat(atoms%nat,nsld) ) ALLOCATE ( zsl(2,nsld),volsl(nsld) ) ALLOCATE ( volintsl(nsld) ) CALL slabgeom(& atoms,cell,nsld,& nsl,zsl,nmtsl,nslat,volsl,volintsl) ALLOCATE ( qintsl(nsld,dimension%neigd)) ALLOCATE ( qmtsl(nsld,dimension%neigd)) ALLOCATE ( orbcomp(dimension%neigd,23,atoms%nat) ) ALLOCATE ( qmtp(dimension%neigd,atoms%nat) ) IF (.NOT.input%film) qvac(:,:,:,jspin) = 0.0 ELSE ALLOCATE(nmtsl(1,1),nslat(1,1),zsl(1,1),volsl(1),volintsl(1)) ALLOCATE(qintsl(1,1),qmtsl(1,1),orbcomp(1,1,1),qmtp(1,1)) END IF !-q_sl ! !--> loop over k-points: each can be a separate task ! IF (kpts%nkpt < mpi%isize) THEN l_evp = .true. IF (l_mcd) THEN mcd(:,:,:) = 0.0 ENDIF ener(:,:,:) = 0.0 sqal(:,:,:) = 0.0 qal(:,:,:,:) = 0.0 enerlo(:,:,:) = 0.0 sqlo(:,:,:) = 0.0 ELSE l_evp = .false. END IF ALLOCATE ( we(dimension%neigd) ) i_rec = 0 ; n_rank = 0 ! For k-point paralelization: ! the number of iterations is adjusted to the number of MPI processes to synchronize RMA operations if (l_evp) then nkpt_extended = kpts%nkpt else nkpt_extended = (kpts%nkpt / mpi%isize + 1) * mpi%isize endif DO ikpt = 1,nkpt_extended i_rec = i_rec + 1 IF ((mod(i_rec-1,mpi%isize).EQ.mpi%irank).OR.l_evp) THEN IF ( ikpt < kpts%nkpt + 1) THEN !-t3e we=0.0 !---> determine number of occupied bands and set weights (we) noccbd = 0 DO i = 1,dimension%neigd ! nbands we(i) = results%w_iks(n_bands(n_rank)+i,ikpt,jspin) IF (noco%l_noco) we(i) = results%w_iks(i,ikpt,1) IF ((we(i).GE.1.e-8).OR.input%pallst) THEN noccbd = noccbd + 1 ELSE we(i)=0.0 END IF END DO ! uncomment this so that cdinf plots works for all states ! noccbd = neigd ! ! -> Gu test: distribute ev's among the processors... ! CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi) skip_t = skip_tt IF (l_evp.AND.(mpi%isize.GT.1)) THEN IF (banddos%dos) THEN noccbd_l = CEILING( real(n_bands(1)) / mpi%isize ) n_start = mpi%irank*noccbd_l + 1 n_end = min( (mpi%irank+1)*noccbd_l , n_bands(1) ) ELSE noccbd_l = CEILING( real(noccbd) / mpi%isize ) n_start = mpi%irank*noccbd_l + 1 n_end = min( (mpi%irank+1)*noccbd_l , noccbd ) END IF noccbd = n_end - n_start + 1 IF (noccbd<1) THEN noccbd=0 ELSE we(1:noccbd) = we(n_start:n_end) END IF IF (n_start > skip_tt) THEN skip_t = 0 END IF IF (n_end <= skip_tt) THEN skip_t = noccbd END IF IF ((n_start <= skip_tt).AND.(n_end > skip_tt)) THEN skip_t = mod(skip_tt,noccbd) END IF ELSE n_start = 1 IF (banddos%dos) THEN noccbd_l = n_bands(1) n_end = n_bands(1) noccbd = n_bands(1) ELSE noccbd_l = noccbd n_end = noccbd END IF END IF zMat%nbasfcn=lapw%nv(1)+atoms%nlotot IF (noco%l_noco) zMat%nbasfcn=zMat%nbasfcn+lapw%nv(2)+atoms%nlotot IF (zmat%l_real) THEN IF (.NOT.ALLOCATED(zMat%z_r)) THEN ALLOCATE (zMat%z_r(zmat%nbasfcn,dimension%neigd)) zMat%nbands = dimension%neigd END IF zMat%z_r = 0 ELSE IF (.NOT.ALLOCATED(zMat%z_c)) THEN ALLOCATE (zMat%z_c(zmat%nbasfcn,dimension%neigd)) zMat%nbands = dimension%neigd END IF zMat%z_c = 0 endif CALL cdn_read(& eig_id,dimension%nvd,dimension%jspd,mpi%irank,mpi%isize,& ikpt,jspin,zmat%nbasfcn,noco%l_ss,noco%l_noco,& noccbd,n_start,n_end,& ello,evdu,epar,& wk,nbands,eig,zMat) #ifdef CPP_MPI ! Sinchronizes the RMA operations if (l_evp) CALL MPI_BARRIER(mpi%mpi_comm,ie) #endif !IF (l_evp.AND.(isize.GT.1)) THEN ! eig(1:noccbd) = eig(n_start:n_end) !ENDIF ! IF (vacuum%nstm.EQ.3.AND.input%film) THEN CALL nstm3(& sym,atoms,vacuum,stars,ikpt,lapw%nv(jspin),& input,jspin,kpts,& cell,wk,lapw%k1(:,jspin),lapw%k2(:,jspin),& evac(1,jspin),vz,vz0,& gvac1d,gvac2d) END IF IF (noccbd.EQ.0) GO TO 199 ! !---> if slice, only a certain bands are taken into account !---> in order to do this the coresponding eigenvalues, eigenvectors !---> and weights have to be copied to the beginning of the arrays !---> eig, z and we and the number of occupied bands (noccbd) has to !---> changed IF (sliceplot%slice) THEN IF (mpi%irank==0) WRITE (16,FMT=*) 'NNNE',sliceplot%nnne IF (mpi%irank==0) WRITE (16,FMT=*) 'sliceplot%kk',sliceplot%kk nslibd = 0 IF (input%pallst) we(:nbands) = wk IF (sliceplot%kk.EQ.0) THEN IF (mpi%irank==0) THEN WRITE (16,FMT='(a)') 'ALL K-POINTS ARE TAKEN IN SLICE' WRITE (16,FMT='(a,i2)') ' sliceplot%slice: k-point nr.',ikpt END IF DO i = 1,nbands IF (eig(i).GE.sliceplot%e1s .AND. eig(i).LE.sliceplot%e2s) THEN nslibd = nslibd + 1 eig(nslibd) = eig(i) we(nslibd) = we(i) if (zmat%l_real) THEN zMat%z_r(:,nslibd) = zMat%z_r(:,i) else zMat%z_c(:,nslibd) = zMat%z_c(:,i) endif END IF END DO IF (mpi%irank==0) WRITE (16,'(a,i3)') ' eigenvalues in sliceplot%slice:',nslibd ELSE IF (sliceplot%kk.EQ.ikpt) THEN IF (mpi%irank==0) WRITE (16,FMT='(a,i2)') ' sliceplot%slice: k-point nr.',ikpt IF ((sliceplot%e1s.EQ.0.0) .AND. (sliceplot%e2s.EQ.0.0)) THEN IF (mpi%irank==0) WRITE (16,FMT='(a,i5,f10.5)') 'slice: eigenvalue nr.',& sliceplot%nnne,eig(sliceplot%nnne) nslibd = nslibd + 1 eig(nslibd) = eig(sliceplot%nnne) we(nslibd) = we(sliceplot%nnne) if (zmat%l_real) Then zMat%z_r(:,nslibd) = zMat%z_r(:,sliceplot%nnne) else zMat%z_c(:,nslibd) = zMat%z_c(:,sliceplot%nnne) endif ELSE DO i = 1,nbands IF (eig(i).GE.sliceplot%e1s .AND. eig(i).LE.sliceplot%e2s) THEN nslibd = nslibd + 1 eig(nslibd) = eig(i) we(nslibd) = we(i) if (zmat%l_real) THEN zMat%z_r(:,nslibd) = zMat%z_r(:,i) else zMat%z_c(:,nslibd) = zMat%z_c(:,i) endif END IF END DO IF (mpi%irank==0) WRITE (16,FMT='(a,i3)')' eigenvalues in sliceplot%slice:',nslibd END IF END IF noccbd = nslibd IF (nslibd.EQ.0) GO TO 199 !200 END IF ! sliceplot%slice !---> in normal iterations the charge density of the unoccupied !---> does not need to be calculated (in pwden, vacden and abcof) IF (banddos%dos.AND. .NOT.(l_evp.AND.(mpi%isize.GT.1)) ) THEN noccbd=nbands END IF ! ----> add in spin-doubling factor we(:noccbd) = 2.*we(:noccbd)/input%jspins !---> pk non-collinear !---> valence density in the interstitial and vacuum region !---> has to be called only once (if jspin=1) in the non-collinear !---> case ! ----> valence density in the interstitial region IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN CALL timestart("cdnval: pwden") CALL pwden(stars,kpts,banddos,oneD, input,mpi,noco,cell,atoms,sym,ikpt,& jspin,lapw,noccbd,igq_fft,we, eig,qpw,cdom,qis,results%force,f_b8,zMat) CALL timestop("cdnval: pwden") END IF !+new !---> charge of each valence state in this k-point of the SBZ !---> in the layer interstitial region of the film ! IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN CALL q_int_sl(jspin,stars,atoms,sym, volsl,volintsl,& cell,noccbd,lapw, nsl,zsl,nmtsl,oneD, qintsl(:,:),zMat) ! END IF END IF !-new c !---> valence density in the vacuum region IF (input%film) THEN IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN CALL timestart("cdnval: vacden") CALL vacden(vacuum,dimension,stars,oneD, kpts,input, cell,atoms,noco,banddos,& gvac1d,gvac2d, we,ikpt,jspin,vz,vz0, noccbd,lapw, evac,eig,& rhtxy,rht,qvac,qvlay, qstars,cdomvz,cdomvxy,zMat) CALL timestop("cdnval: vacden") END IF !---> perform Brillouin zone integration and summation over the !---> bands in order to determine the vacuum energy parameters. DO ispin = jsp_start,jsp_end DO ivac = 1,vacuum%nvac pvac(ivac,ispin)=pvac(ivac,ispin)+dot_product(eig(:noccbd)*qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd)) svac(ivac,ispin)=svac(ivac,ispin)+dot_product(qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd)) END DO END DO END IF !---> valence density in the atomic spheres !---> construct a(tilta) and b(tilta) IF (noco%l_mperp) THEN ALLOCATE ( acof(noccbd,0:dimension%lmd,atoms%nat,dimension%jspd),& ! Deallocated before call to sympsi bcof(noccbd,0:dimension%lmd,atoms%nat,dimension%jspd), & ccof(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat,dimension%jspd) ) ELSE ALLOCATE ( acof(noccbd,0:dimension%lmd,atoms%nat,jspin:jspin),& bcof(noccbd,0:dimension%lmd,atoms%nat,jspin:jspin),& ccof(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat,jspin:jspin) ) END IF DO ispin = jsp_start,jsp_end IF (input%l_f) THEN CALL timestart("cdnval: to_pulay") ALLOCATE (e1cof(noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),& ! Deallocated after call to force_a21 e2cof(noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),& acoflo(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat),& bcoflo(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat),& aveccof(3,noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),& bveccof(3,noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),& cveccof(3,-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat) ) CALL abcof(input,atoms,sym, cell,lapw,noccbd,usdus, noco,ispin,oneD,& acof(:,0:,:,ispin),bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin),zMat,& eig,acoflo,bcoflo,e1cof,e2cof,aveccof,bveccof,cveccof) ! CALL to_pulay(input,atoms,noccbd,sym, lapw, noco,cell,noccbd,eig,usdus,& ! ispin,oneD, acof(:,0:,:,ispin),bcof(:,0:,:,ispin),& ! e1cof,e2cof,aveccof,bveccof, ccof(-atoms%llod,1,1,1,ispin),acoflo,bcoflo,cveccof,zMat) CALL timestop("cdnval: to_pulay") ELSE CALL timestart("cdnval: abcof") CALL abcof(input,atoms,sym, cell,lapw,noccbd,usdus, noco,ispin,oneD,& acof(:,0:,:,ispin),bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin),zMat) CALL timestop("cdnval: abcof") END IF IF (atoms%n_u.GT.0) THEN CALL n_mat(atoms,sym,noccbd,usdus,ispin,we, acof(:,0:,:,ispin),bcof(:,0:,:,ispin),& ccof(-atoms%llod:,:,:,:,ispin), n_mmp) END IF ! !---> perform Brillouin zone integration and summation over the !---> bands in order to determine the energy parameters for each !---> atom and angular momentum ! IF (.not.sliceplot%slice) THEN CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,ccof,& skip_t,l_evp,acof(:,0:,:,ispin),bcof(:,0:,:,ispin),usdus,& ncore,l_mcd,m_mcd,& enerlo(1,1,ispin),sqlo(1,1,ispin),& ener(0,1,ispin),sqal(0,1,ispin),& qal(0:,:,:,ispin),mcd) IF (noco%l_mperp.AND.(ispin == jsp_end)) THEN CALL qal_21(atoms, input,noccbd,we,ccof,& noco,acof,bcof,mt21,lo21,uloulopn21,& qal,qmat) END IF END IF ! !+new !---> layer charge of each valence state in this k-point of the SBZ !---> from the mt-sphere region of the film ! IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN CALL q_mt_sl(ispin, atoms,noccbd,nsld, ikpt,noccbd,ccof(-atoms%llod,1,1,1,ispin),& skip_t,noccbd, acof(:,0:,:,ispin),bcof(:,0:,:,ispin),usdus,& nmtsl,nsl, qmtsl(:,:)) INQUIRE (file='orbcomprot',exist=l_orbcomprot) IF (l_orbcomprot) THEN ! rotate ab-coeffs CALL abcrot2(atoms, noccbd,& acof(:,0:,:,ispin),bcof(:,0:,:,ispin),& ccof(-atoms%llod:,:,:,:,ispin)) END IF CALL orb_comp(ispin,noccbd,atoms,noccbd,usdus,acof(1:,0:,1:,ispin),bcof(1:,0:,1:,ispin),& ccof(-atoms%llod:,1:,1:,1:,ispin), orbcomp, qmtp) END IF !-new !---> set up coefficients for the spherical and CALL timestart("cdnval: rhomt") CALL rhomt(atoms,we,noccbd, acof(:,0:,:,ispin),bcof(:,0:,:,ispin),& uu(0:,:,ispin),dd(0:,:,ispin),du(0:,:,ispin)) CALL timestop("cdnval: rhomt") !+soc IF (noco%l_soc) THEN CALL orbmom(atoms,noccbd, we,acof(:,0:,:,ispin),bcof(:,0:,:,ispin),& ccof(-atoms%llod:,:,:,:,ispin), orb(0:,-atoms%lmaxd:,:,ispin),orbl(:,-atoms%llod:,:,ispin),& orblo(:,:,-atoms%llod:,:,ispin) ) END IF ! -soc !---> non-spherical m.t. density CALL timestart("cdnval: rhonmt") CALL rhonmt(atoms,sphhar, we,noccbd,sym, acof(:,0:,:,ispin),bcof(:,0:,:,ispin),& uunmt(0:,:,:,ispin),ddnmt(0:,:,:,ispin), udnmt(0:,:,:,ispin),dunmt(0:,:,:,ispin)) CALL timestop("cdnval: rhonmt") !---> set up coefficients of the local orbitals and the !---> flapw - lo cross terms for the spherical and !---> non-spherical mt density CALL timestart("cdnval: rho(n)mtlo") CALL rhomtlo(atoms,& noccbd,we,acof(:,0:,:,ispin),bcof(:,0:,:,ispin),& ccof(-atoms%llod:,:,:,:,ispin),& aclo(:,:,ispin),bclo(:,:,ispin),cclo(:,:,:,ispin)) ! CALL rhonmtlo(& atoms,sphhar,& noccbd,we,acof(:,0:,:,ispin),& bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin),& acnmt(0:,:,:,:,ispin),bcnmt(0:,:,:,:,ispin),& ccnmt(:,:,:,:,ispin)) CALL timestop("cdnval: rho(n)mtlo") IF (input%l_f) THEN CALL timestart("cdnval: force_a12/21") IF (.not.input%l_useapw) THEN CALL force_a12(atoms,noccbd,sym, dimension,cell,oneD,& we,ispin,noccbd,usdus,acof(:,0:,:,ispin),& bcof(:,0:,:,ispin),e1cof,e2cof, acoflo,bcoflo, results,f_a12) ENDIF CALL force_a21(input,atoms,dimension,noccbd,sym,& oneD,cell,we,ispin,epar(0:,:,ispin),noccbd,eig,usdus,acof(:,0:,:,ispin),& bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin), aveccof,bveccof,cveccof,& results,f_a21,f_b4) DEALLOCATE (e1cof,e2cof,aveccof,bveccof) DEALLOCATE (acoflo,bcoflo,cveccof) CALL timestop("cdnval: force_a12/21") END IF IF(l_cs) THEN CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,& dimension%neigd,noccbd,results%ef,banddos%sig_dos,& eig,we,acof(1,0,1,ispin),bcof(1,0,1,ispin),& ccof(-atoms%llod,1,1,1,ispin)) END IF END DO !---> end loop over ispin IF (noco%l_mperp) THEN CALL rhomt21(atoms, we,noccbd,acof,bcof, ccof,& mt21,lo21,uloulop21) IF (l_fmpl) THEN CALL rhonmt21(atoms,llpd,sphhar, we,noccbd,sym, acof,bcof,& uunmt21,ddnmt21,udnmt21,dunmt21) END IF END IF DEALLOCATE (acof,bcof,ccof) ! 199 CONTINUE IF ((banddos%dos .OR. banddos%vacdos .OR. input%cdinf) ) THEN CALL timestart("cdnval: write_info") ! !---> calculate charge distribution of each state (l-character ...) !---> and write the information to the files dosinp and vacdos !---> for dos and bandstructure plots ! !--dw parallel writing of vacdos,dosinp.... ! write data to direct access file first, write to formated file later by PE 0 only! !--dw since z is no longer an argument of cdninf sympsi has to be called here! ! cartk=matmul(lapw%bkpt,cell%bmat) IF (banddos%ndir.GT.0) THEN CALL sympsi(lapw%bkpt,lapw%nv(jspin),lapw%k1(:,jspin),lapw%k2(:,jspin),& lapw%k3(:,jspin),sym,dimension,nbands,cell,eig,noco, ksym,jsym,zMat) END IF ! !--dw now write k-point data to tmp_dos ! CALL write_dos(eig_id,ikpt,jspin,qal(:,:,:,jspin),qvac(:,:,ikpt,jspin),qis(:,ikpt,jspin),& qvlay(:,:,:,ikpt,jspin),qstars,ksym,jsym,mcd,qintsl,& qmtsl(:,:),qmtp(:,:),orbcomp) CALL timestop("cdnval: write_info") !-new_sl END IF !---> end of loop over PE's IF (zmat%l_real) THEN DEALLOCATE (zMat%z_r) ELSE DEALLOCATE (zMat%z_c) END IF ELSE !(ikpt < nkpt + 1) #ifdef CPP_MPI ! Synchronizes the RMA operations CALL MPI_BARRIER(mpi%mpi_comm,ie) #endif END IF END IF ! --> end "IF ((mod(i_rec-1,mpi%isize).EQ.mpi%irank).OR.l_evp) THEN" END DO !---> end of k-point loop DEALLOCATE (we,f,g,usdus%us,usdus%dus,usdus%duds,usdus%uds,usdus%ddn) !+t3e #ifdef CPP_MPI CALL timestart("cdnval: mpi_col_den") DO ispin = jsp_start,jsp_end CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,& input,noco,l_fmpl,ispin,llpd, rhtxy(1,1,1,ispin),& rht(1,1,ispin),qpw(1,ispin), ener(0,1,ispin),sqal(0,1,ispin),& results,svac(1,ispin),pvac(1,ispin),uu(0,1,ispin),& dd(0,1,ispin),du(0,1,ispin),uunmt(0,1,1,ispin),ddnmt(0,1,1,ispin),& udnmt(0,1,1,ispin),dunmt(0,1,1,ispin),sqlo(1,1,ispin),& aclo(1,1,ispin),bclo(1,1,ispin),cclo(1,1,1,ispin),& acnmt(0,1,1,1,ispin),bcnmt(0,1,1,1,ispin),& ccnmt(1,1,1,1,ispin),enerlo(1,1,ispin),& orb(0,-atoms%lmaxd,1,ispin),orbl(1,-atoms%llod,1,ispin),& orblo(1,1,-atoms%llod,1,ispin),mt21,lo21,uloulop21,& uunmt21,ddnmt21,udnmt21,dunmt21,cdom,cdomvz,cdomvxy,n_mmp) END DO CALL timestop("cdnval: mpi_col_den") #endif IF(l_cs) CALL corespec_ddscs(jspin,input%jspins) IF (((jspin.eq.input%jspins).OR.noco%l_mperp) .AND. (banddos%dos.or.banddos%vacdos.or.input%cdinf) ) THEN CALL timestart("cdnval: dos") IF (mpi%irank==0) THEN CALL doswrite(& eig_id,dimension,kpts,atoms,vacuum,& input,banddos,& sliceplot,noco,sym,& cell,& l_mcd,ncored,ncore,e_mcd,& results%ef,results%bandgap,nsld,oneD) IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN CALL Ek_write_sl(& eig_id,dimension,kpts,atoms,vacuum,& nsld,input,jspin,& sym,cell,& nsl,nslat) END IF END IF #ifdef CPP_MPI CALL MPI_BARRIER(mpi%mpi_comm,ie) #endif CALL timestop("cdnval: dos") END IF IF (mpi%irank==0) THEN CALL cdnmt(& dimension%jspd,atoms,sphhar,llpd,& noco,l_fmpl,jsp_start,jsp_end,& epar,ello,vr(:,0,:,:),uu,du,dd,uunmt,udnmt,dunmt,ddnmt,& usdus,usdus%uloulopn,aclo,bclo,cclo,acnmt,bcnmt,ccnmt,& orb,orbl,orblo,mt21,lo21,uloulopn21,uloulop21,& uunmt21,ddnmt21,udnmt21,dunmt21,& chmom,clmom,& qa21,rho) DO ispin = jsp_start,jsp_end WRITE (6,*) 'Energy Parameters for spin:',ispin IF (.not.sliceplot%slice) THEN CALL mix_enpara(& ispin,atoms,vacuum,obsolete,input,& enpara,& vr(:,0,:,:),vz,pvac(1,ispin),& svac(1,ispin),& ener(0,1,ispin),sqal(0,1,ispin),& enerlo(1,1,ispin),& sqlo(1,1,ispin)) CALL w_enpara(& atoms,jspin,input%film,& enpara,16) END IF !---> check continuity of charge density IF (input%cdinf) THEN CALL timestart("cdnval: cdninf-stuff") WRITE (6,FMT=8210) ispin 8210 FORMAT (/,5x,'check continuity of cdn for spin=',i2) IF (input%film .AND. .NOT.oneD%odi%d1) THEN !---> vacuum boundaries npd = min(dimension%nspd,25) CALL points(xp,npd) DO ivac = 1,vacuum%nvac sign = 3. - 2.*ivac DO j = 1,npd xp(3,j) = sign*cell%z1/cell%amat(3,3) END DO CALL checkdop(& xp,npd,0,0,ivac,1,ispin,.true.,dimension,atoms,& sphhar,stars,sym,& vacuum,cell,oneD,& qpw,rho,rhtxy,rht) END DO ELSE IF (oneD%odi%d1) THEN !-odim npd = min(dimension%nspd,25) CALL cylpts(xp,npd,cell%z1) CALL checkdop(& xp,npd,0,0,ivac,1,ispin,.true.,dimension,atoms,& sphhar,stars,sym,& vacuum,cell,oneD,& qpw,rho,rhtxy,rht) !+odim END IF !---> m.t. boundaries nat = 1 DO n = 1, atoms%ntype CALL sphpts(xp,dimension%nspd,atoms%rmt(n),atoms%pos(1,atoms%nat)) CALL checkdop(& xp,dimension%nspd,n,nat,0,-1,ispin,.true.,& dimension,atoms,sphhar,stars,sym,& vacuum,cell,oneD,& qpw,rho,rhtxy,rht) nat = nat + atoms%neq(n) END DO CALL timestop("cdnval: cdninf-stuff") END IF !+for !---> forces of equ. A8 of Yu et al. IF ((input%l_f)) THEN CALL timestart("cdnval: force_a8") CALL force_a8(input,atoms,sphhar, ispin, vr(:,:,:,ispin),rho,& f_a12,f_a21,f_b4,f_b8,results%force) CALL timestop("cdnval: force_a8") END IF !-for END DO ! end of loop ispin = jsp_start,jsp_end CALL closeXMLElement('mtCharges') END IF ! end of (mpi%irank==0) !+t3e !Note: no deallocation anymore, we rely on Fortran08 :-) IF ((jsp_end.EQ.input%jspins)) THEN IF ((banddos%dos.OR.banddos%vacdos).AND.(banddos%ndir/=-2)) CALL juDFT_end("DOS OK",mpi%irank) IF (vacuum%nstm.EQ.3) CALL juDFT_end("VACWAVE OK",mpi%irank) END IF END SUBROUTINE cdnval END MODULE m_cdnval