cleanup

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
 

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