diff --git a/main/cdngen.F90 b/main/cdngen.F90 index da478a7d2852087285e9e5e886a063edc95c14a9..e859d55eda7788fee4ec0f1827405d57d3cd2c2d 100644 --- a/main/cdngen.F90 +++ b/main/cdngen.F90 @@ -231,7 +231,6 @@ subroutine save_kinED(xcpot, input, noco, stars, cell, sym) use m_types use m_pw_tofrom_grid use m_judft_stop - use m_metagga, only: give_stats use m_npy implicit none diff --git a/xc-pot/metagga.F90 b/xc-pot/metagga.F90 index f0db234eb6172daf9b80fdc54d18b9c383a357c7..2ba6237efd9a299d581b4af11ff661ea537b638e 100644 --- a/xc-pot/metagga.F90 +++ b/xc-pot/metagga.F90 @@ -52,19 +52,6 @@ CONTAINS #endif END SUBROUTINE calc_kinEnergyDen - SUBROUTINE give_stats(array, array_name) - implicit none - real, intent(in) :: array(:,:) - character(len=*), optional :: array_name - - write (*,*) "#############################" - if(present(array_name)) write (*,*) "Array = ", array_name - write (*,'(A,ES17.10)') "min = ", minval(array) - write (*,'(A,ES17.10)') "max = ", maxval(array) - write (*,'(A,ES17.10)') "mean = ", sum(array) / size(array) - write (*,*) "#############################" - END SUBROUTINE give_stats - SUBROUTINE dump_array(array, filename) implicit none real, intent(in) :: array(:,:) @@ -160,7 +147,6 @@ CONTAINS INTEGER, INTENT(in) :: jspin TYPE(t_kpts), INTENT(in) :: kpts REAL, INTENT(inout) :: f_ik(:,:) ! f_ik(band_idx, kpt_idx) - REAL, allocatable :: eigenvalues(:,:) ! local vars REAL :: e_i(SIZE(f_ik,dim=1)) @@ -168,7 +154,6 @@ CONTAINS DO ikpt = 1,kpts%nkpt CALL read_eig(eig_id,ikpt,jspin, eig=e_i) - eigenvalues(:,ikpt) = e_i f_ik(:,ikpt) = f_ik(:,ikpt) * e_i ENDDO END SUBROUTINE calc_EnergyDen_auxillary_weights @@ -307,20 +292,20 @@ CONTAINS IF (banddos%l_mcd.AND..NOT.PRESENT(mcd)) CALL juDFT_error("mcd is missing",calledby ="cdnval") ! calculation of core spectra (EELS) initializations -start- - l_coreSpec = .FALSE. - IF (PRESENT(coreSpecInput)) THEN - 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() - l_coreSpec = l_cs - END IF + !l_coreSpec = .FALSE. + !IF (PRESENT(coreSpecInput)) THEN + !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() + !l_coreSpec = l_cs + !END IF ! 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 + !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) DO iType = 1, atoms%ntype @@ -380,59 +365,58 @@ CONTAINS IF (.NOT.((jspin.EQ.2).AND.noco%l_noco)) THEN ! valence density in the interstitial region CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,& - jspin,lapw,noccbd,we,eig,kinED,results,force%f_b8,zMat,dos) + jspin,lapw,noccbd,we,eig,kinED,results,force%f_b8,zPrime,dos) ! charge of each valence state in this k-point of the SBZ in the layer interstitial region of the film - IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat) - ! valence density in the vacuum region - IF (input%film) THEN - CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,& - gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac,eig,kinED,zMat,dos) - END IF + !IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat) + !! valence density in the vacuum region + !IF (input%film) THEN + !CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,& + !gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac,eig,kinED,zMat,dos) + !END IF END IF - IF (input%film) CALL regCharges%sumBandsVac(vacuum,dos,noccbd,ikpt,jsp_start,jsp_end,eig,we) ! valence density in the atomic spheres - CALL eigVecCoeffs%init(input,DIMENSION,atoms,noco,jspin,noccbd) - DO ispin = jsp_start, jsp_end - IF (input%l_f) CALL force%init2(noccbd,input,atoms) - CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,& - eigVecCoeffs%acof(:,0:,:,ispin),eigVecCoeffs%bcof(:,0:,:,ispin),& - eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force) - IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,kinED%mmpMat(:,:,:,jspin)) - - ! perform Brillouin zone integration and summation over the - ! bands in order to determine the energy parameters for each atom and angular momentum - CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,& - skip_t,cdnvalJob%l_evp,eigVecCoeffs,usdus,regCharges,dos,banddos%l_mcd,mcd) - - IF (noco%l_mperp.AND.(ispin==jsp_end)) THEN - CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos) - ENDIF - - ! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film - IF (l_dosNdir) THEN - IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab) - - INQUIRE (file='orbcomprot',exist=l_orbcomprot) - IF (l_orbcomprot) CALL abcrot2(atoms,noccbd,eigVecCoeffs,ispin) ! rotate ab-coeffs - - IF (PRESENT(orbcomp)) CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp) - END IF - - CALL calcDenCoeffs(atoms,sphhar,sym,we,noccbd,eigVecCoeffs,ispin,denCoeffs) - - IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb) - IF (input%l_f) CALL force%addContribsA21A12(input,atoms,dimension,sym,cell,oneD,enpara,& - usdus,eigVecCoeffs,noccbd,ispin,eig,we,results) - IF(l_coreSpec) CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,& - noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs) - END DO ! end loop over ispin - IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd) - - IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf).AND.(banddos%ndir.GT.0)) THEN - ! since z is no longer an argument of cdninf sympsi has to be called here! - CALL sympsi(lapw,jspin,sym,dimension,nbands,cell,eig,noco,dos%ksym(:,ikpt,jspin),dos%jsym(:,ikpt,jspin),zMat) - END IF + !CALL eigVecCoeffs%init(input,DIMENSION,atoms,noco,jspin,noccbd) + !DO ispin = jsp_start, jsp_end + !IF (input%l_f) CALL force%init2(noccbd,input,atoms) + !CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,& + !eigVecCoeffs%acof(:,0:,:,ispin),eigVecCoeffs%bcof(:,0:,:,ispin),& + !eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force) + !IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,kinED%mmpMat(:,:,:,jspin)) + + !! perform Brillouin zone integration and summation over the + !! bands in order to determine the energy parameters for each atom and angular momentum + !CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,& + !skip_t,cdnvalJob%l_evp,eigVecCoeffs,usdus,regCharges,dos,banddos%l_mcd,mcd) + + !IF (noco%l_mperp.AND.(ispin==jsp_end)) THEN + !CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos) + !ENDIF + + !! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film + !IF (l_dosNdir) THEN + !IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab) + + !INQUIRE (file='orbcomprot',exist=l_orbcomprot) + !IF (l_orbcomprot) CALL abcrot2(atoms,noccbd,eigVecCoeffs,ispin) ! rotate ab-coeffs + + !IF (PRESENT(orbcomp)) CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp) + !END IF + + !CALL calcDenCoeffs(atoms,sphhar,sym,we,noccbd,eigVecCoeffs,ispin,denCoeffs) + + !IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb) + !IF (input%l_f) CALL force%addContribsA21A12(input,atoms,dimension,sym,cell,oneD,enpara,& + !usdus,eigVecCoeffs,noccbd,ispin,eig,we,results) + !IF(l_coreSpec) CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,& + !noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs) + !END DO ! end loop over ispin + !IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd) + + !IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf).AND.(banddos%ndir.GT.0)) THEN + !! since z is no longer an argument of cdninf sympsi has to be called here! + !CALL sympsi(lapw,jspin,sym,dimension,nbands,cell,eig,noco,dos%ksym(:,ikpt,jspin),dos%jsym(:,ikpt,jspin),zMat) + !END IF END DO ! end of k-point loop #ifdef CPP_MPI @@ -442,26 +426,27 @@ CONTAINS END DO #endif - IF (mpi%irank==0) THEN - CALL cdnmt(mpi, input%jspins,atoms,sphhar,noco,jsp_start,jsp_end,& - enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,kinED%mt) - IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins) - DO ispin = jsp_start,jsp_end - IF (input%cdinf) THEN - WRITE (6,FMT=8210) ispin -8210 FORMAT (/,5x,'check continuity of cdn for spin=',i2) - CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,kinED,ispin) - END IF - IF (input%l_f) CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),kinED%mt,force,results) - END DO - CALL closeXMLElement('mtCharges') - END IF + !IF (mpi%irank==0) THEN + !CALL cdnmt(mpi, input%jspins,atoms,sphhar,noco,jsp_start,jsp_end,& + !enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,kinED%mt) + !IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins) + !DO ispin = jsp_start,jsp_end + !IF (input%cdinf) THEN + !WRITE (6,FMT=8210) ispin +!8210 FORMAT (/,5x,'check continuity of cdn for spin=',i2) + !CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,kinED,ispin) + !END IF + !IF (input%l_f) CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),kinED%mt,force,results) + !END DO + !CALL closeXMLElement('mtCharges') + !END IF CALL timestop("cdnval") END SUBROUTINE calc_kinED_pw subroutine set_zPrime(dim_idx, zMat, kpt, lapw, cell, zPrime) USE m_types + USE m_constants implicit none INTEGER, intent(in) :: dim_idx TYPE (t_mat), intent(in) :: zMat @@ -482,9 +467,9 @@ CONTAINS fac = k_plus_g(dim_idx) if(zPrime%l_real) then - zPrime%data_r(basis_idx,:) = fac * zPrime%data_r(basis_idx,:) + zPrime%data_r(basis_idx,:) = fac * zMat%data_r(basis_idx,:) else - zPrime%data_c(basis_idx,:) = fac * zPrime%data_c(basis_idx,:) + zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:) endif enddo end subroutine set_zPrime