Merge branch 'myqe' into 'develop'

More I/O cleaning and workarounds

See merge request QEF/q-e!523
This commit is contained in:
giannozz 2019-07-12 09:44:29 +00:00
commit ac870c0433
8 changed files with 133 additions and 171 deletions

View File

@ -806,8 +806,8 @@ MODULE cp_restart_new
nbnd_ = nupdwn(1)
ALLOCATE( occ_(nbnd_, nspin), et_(nbnd_, nspin) )
CALL qexsd_copy_band_structure( output_obj%band_structure, lsda_, &
nk_, isk_, natomwfc, nbnd, nbnd_up, nbnd_dw, nelec_, wk_, &
occ_, ef, ef_up, ef_dw, et_ )
nk_, isk_, natomwfc, nbnd, nbnd_up, nbnd_dw, nelec_, xk, &
wk_, occ_, ef, ef_up, ef_dw, et_ )
! FIXME: in the call, the same array is passed as both occ0 and occm!
DO iss = 1, nspin
ib = iupdwn(iss)

View File

@ -23,7 +23,8 @@ MODULE qexsd_copy
qexsd_copy_atomic_species, qexsd_copy_atomic_structure, &
qexsd_copy_symmetry, qexsd_copy_algorithmic_info, &
qexsd_copy_basis_set, qexsd_copy_dft, qexsd_copy_band_structure, &
qexsd_copy_efield, qexsd_copy_magnetization
qexsd_copy_efield, qexsd_copy_magnetization, qexsd_copy_kpoints, &
qexsd_copy_efermi
!
CONTAINS
!-------------------------------------------------------------------------------
@ -445,7 +446,7 @@ CONTAINS
!
!------------------------------------------------------------------------
SUBROUTINE qexsd_copy_band_structure( band_struct_obj, lsda, nkstot, &
isk, natomwfc, nbnd, nbnd_up, nbnd_dw, nelec, wk, wg, &
isk, natomwfc, nbnd, nbnd_up, nbnd_dw, nelec, xk, wk, wg, &
ef, ef_up, ef_dw, et )
!------------------------------------------------------------------------
!
@ -458,14 +459,14 @@ CONTAINS
LOGICAL, INTENT(out) :: lsda
INTEGER, INTENT(out) :: nkstot, natomwfc, nbnd, nbnd_up, nbnd_dw, &
isk(:)
REAL(dp), INTENT(out):: nelec, ef, ef_up, ef_dw, wk(:)
REAL(dp), INTENT(out):: nelec, ef, ef_up, ef_dw, xk(:,:), wk(:)
REAL(dp), INTENT(inout), ALLOCATABLE :: wg(:,:), et(:,:)
!
LOGICAL :: two_fermi_energies
INTEGER :: ik
!
lsda = band_struct_obj%lsda
nkstot = band_struct_obj%nks
nelec = band_struct_obj%nelec
natomwfc = band_struct_obj%num_of_atomic_wfc
!
IF ( lsda) THEN
@ -506,26 +507,17 @@ CONTAINS
nbnd_dw = nbnd
isk(1:nkstot) = 1
END IF
!
IF ( band_struct_obj%fermi_energy_ispresent) THEN
ef = band_struct_obj%fermi_energy
ef_up = 0.d0
ef_dw = 0.d0
ELSE IF ( band_struct_obj%two_fermi_energies_ispresent ) THEN
ef = 0.d0
ef_up = band_struct_obj%two_fermi_energies(1)
ef_dw = band_struct_obj%two_fermi_energies(2)
ELSE
ef = 0.d0
ef_up = 0.d0
ef_dw = 0.d0
END IF
!
CALL qexsd_copy_efermi ( band_struct_obj, &
nelec, ef, two_fermi_energies, ef_up, ef_dw )
!
IF ( .NOT. ALLOCATED(et) ) ALLOCATE( et(nbnd,nkstot) )
IF ( .NOT. ALLOCATED(wg) ) ALLOCATE( wg(nbnd,nkstot) )
!
DO ik =1, band_struct_obj%ndim_ks_energies
IF ( band_struct_obj%lsda) THEN
xk(:,ik) = band_struct_obj%ks_energies(ik)%k_point%k_point(:)
xk(:,ik + band_struct_obj%ndim_ks_energies) = xk(:,ik)
wk(ik) = band_struct_obj%ks_energies(ik)%k_point%weight
wk(ik + band_struct_obj%ndim_ks_energies ) = wk(ik)
et(1:nbnd_up,ik) = band_struct_obj%ks_energies(ik)%eigenvalues%vector(1:nbnd_up)
@ -536,6 +528,7 @@ CONTAINS
wg(1:nbnd_dw,ik+band_struct_obj%ndim_ks_energies) = &
band_struct_obj%ks_energies(ik)%occupations%vector(nbnd_up+1:nbnd_up+nbnd_dw)*wk(ik)
ELSE
xk(:,ik) = band_struct_obj%ks_energies(ik)%k_point%k_point(:)
wk(ik) = band_struct_obj%ks_energies(ik)%k_point%weight
et (1:nbnd,ik) = band_struct_obj%ks_energies(ik)%eigenvalues%vector(1:nbnd)
wg (1:nbnd,ik) = band_struct_obj%ks_energies(ik)%occupations%vector(1:nbnd)*wk(ik)
@ -545,6 +538,34 @@ CONTAINS
!
END SUBROUTINE qexsd_copy_band_structure
!
SUBROUTINE qexsd_copy_efermi ( band_struct_obj, &
nelec, ef, two_fermi_energies, ef_up, ef_dw )
!------------------------------------------------------------------------
!
USE qes_types_module, ONLY : band_structure_type
!
IMPLICIT NONE
TYPE ( band_structure_type) :: band_struct_obj
LOGICAL, INTENT(out) :: two_fermi_energies
REAL(dp), INTENT(out):: nelec, ef, ef_up, ef_dw
!
nelec = band_struct_obj%nelec
two_fermi_energies = band_struct_obj%two_fermi_energies_ispresent
IF ( band_struct_obj%fermi_energy_ispresent) THEN
ef = band_struct_obj%fermi_energy
ef_up = 0.d0
ef_dw = 0.d0
ELSE IF ( two_fermi_energies ) THEN
ef = 0.d0
ef_up = band_struct_obj%two_fermi_energies(1)
ef_dw = band_struct_obj%two_fermi_energies(2)
ELSE
ef = 0.d0
ef_up = 0.d0
ef_dw = 0.d0
END IF
!
END SUBROUTINE qexsd_copy_efermi
!-----------------------------------------------------------------------
SUBROUTINE qexsd_copy_algorithmic_info ( algo_obj, &
real_space, tqr, okvan, okpaw )
@ -647,4 +668,56 @@ CONTAINS
!
END SUBROUTINE qexsd_copy_magnetization
!-----------------------------------------------------------------------
END MODULE qexsd_copy
!
!---------------------------------------------------------------------------
SUBROUTINE qexsd_copy_kpoints ( band_struct_obj, nks_start, xk_start,&
wk_start, nk1, nk2, nk3, k1, k2, k3 )
!---------------------------------------------------------------------------
!
USE qes_types_module, ONLY : band_structure_type
!
IMPLICIT NONE
!
TYPE ( band_structure_type ),INTENT(IN) :: band_struct_obj
INTEGER, INTENT(out) :: nks_start, nk1, nk2, nk3, k1, k2, k3
REAL(dp), ALLOCATABLE, INTENT(inout) :: xk_start(:,:), wk_start(:)
!
INTEGER :: ik
!
!
IF ( band_struct_obj%starting_k_points%monkhorst_pack_ispresent ) THEN
nks_start = 0
nk1 = band_struct_obj%starting_k_points%monkhorst_pack%nk1
nk2 = band_struct_obj%starting_k_points%monkhorst_pack%nk2
nk3 = band_struct_obj%starting_k_points%monkhorst_pack%nk3
k1 = band_struct_obj%starting_k_points%monkhorst_pack%k1
k2 = band_struct_obj%starting_k_points%monkhorst_pack%k2
k3 = band_struct_obj%starting_k_points%monkhorst_pack%k3
ELSE IF (band_struct_obj%starting_k_points%nk_ispresent ) THEN
nks_start = band_struct_obj%starting_k_points%nk
IF ( nks_start > 0 ) THEN
IF ( .NOT. ALLOCATED(xk_start) ) ALLOCATE (xk_start(3,nks_start))
IF ( .NOT. ALLOCATED(wk_start) ) ALLOCATE (wk_start(nks_start))
IF ( nks_start == size( band_struct_obj%starting_k_points%k_point ) ) THEN
DO ik =1, nks_start
xk_start(:,ik) = band_struct_obj%starting_k_points%k_point(ik)%k_point(:)
IF ( band_struct_obj%starting_k_points%k_point(ik)%weight_ispresent) THEN
wk_start(ik) = band_struct_obj%starting_k_points%k_point(ik)%weight
ELSE
wk_start(ik) = 0.d0
END IF
END DO
ELSE
CALL infomsg ( "qexsd_copy_kp: ", &
"actual number of start kpoint not equal to nks_start, set nks_start=0")
nks_start = 0
END IF
END IF
ELSE
CALL errore ("qexsd_copy_kp: ", &
" no information found for initializing brillouin zone information", 1)
END IF
!
END SUBROUTINE qexsd_copy_kpoints
!
END MODULE qexsd_copy

View File

@ -63,23 +63,23 @@ SUBROUTINE xc_gcx( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
!! length of the I/O arrays
INTEGER, INTENT(IN) :: ns
!! spin dimension for input
REAL(DP), INTENT(IN) :: rho(length,ns)
REAL(DP), INTENT(IN) :: rho(:,:)
!! Charge density
REAL(DP), INTENT(IN) :: grho(3,length,ns)
REAL(DP), INTENT(IN) :: grho(:,:,:)
!! gradient
REAL(DP), INTENT(OUT) :: ex(length)
REAL(DP), INTENT(OUT) :: ex(:)
!! exchange energy
REAL(DP), INTENT(OUT) :: ec(length)
REAL(DP), INTENT(OUT) :: ec(:)
!! correlation energy
REAL(DP), INTENT(OUT) :: v1x(length,ns)
REAL(DP), INTENT(OUT) :: v1x(:,:)
!! exchange potential (density part)
REAL(DP), INTENT(OUT) :: v2x(length,ns)
REAL(DP), INTENT(OUT) :: v2x(:,:)
!! exchange potential (gradient part)
REAL(DP), INTENT(OUT) :: v1c(length,ns)
REAL(DP), INTENT(OUT) :: v1c(:,:)
!! correlation potential (density part)
REAL(DP), INTENT(OUT) :: v2c(length,ns)
REAL(DP), INTENT(OUT) :: v2c(:,:)
!! correlation (gradient part)
REAL(DP), INTENT(OUT), OPTIONAL :: v2c_ud(length)
REAL(DP), INTENT(OUT), OPTIONAL :: v2c_ud(:)
!! correlation
!
! ... local variables

View File

@ -630,7 +630,7 @@ SUBROUTINE PAW_gcxc_potential(i, rho_lm,rho_core, v_lm, energy)
USE atom, ONLY : g => rgrid
USE constants, ONLY : sqrtpi, fpi,pi,e2
USE funct, ONLY : igcc_is_lyp
USE xc_gga, ONLY : xc_gcx ! gcxc, gcx_spin, gcc_spin, gcc_spin_more
USE xc_gga, ONLY : xc_gcx
USE mp, ONLY : mp_sum
!
TYPE(paw_info), INTENT(IN) :: i ! atom's minimal info
@ -659,7 +659,7 @@ SUBROUTINE PAW_gcxc_potential(i, rho_lm,rho_core, v_lm, energy)
!
!
!^^^
REAL(DP), ALLOCATABLE :: arho(:), grad2_v(:)
REAL(DP), ALLOCATABLE :: arho(:,:), grad2_v(:)
REAL(DP), ALLOCATABLE :: r_vec(:,:) !, rh(:), zeta(:) !, grhor(:,:), grhoud(:), grh2(:)
!
REAL(DP), DIMENSION(i%m,nspin_gga) :: v1x, v2x, v1c, v2c !workspace
@ -740,7 +740,7 @@ SUBROUTINE PAW_gcxc_potential(i, rho_lm,rho_core, v_lm, energy)
!
! GGA case
!
ALLOCATE( arho(i%m), grad2_v(i%m) )
ALLOCATE( arho(i%m,1), grad2_v(i%m) )
ALLOCATE( gradx(3,i%m,1))
!
!$omp do
@ -751,8 +751,8 @@ SUBROUTINE PAW_gcxc_potential(i, rho_lm,rho_core, v_lm, energy)
CALL PAW_gradient(i, ix, rho_lm, rho_rad, rho_core, grad2, grad)
!
DO k = 1, i%m
arho(k) = rho_rad(k,1)*g(i%t)%rm2(k) + rho_core(k)
arho(k) = ABS(arho(k))
arho(k,1) = rho_rad(k,1)*g(i%t)%rm2(k) + rho_core(k)
arho(k,1) = ABS(arho(k,1))
gradx(:,k,1) = grad(k,:,1)
ENDDO
!

View File

@ -40,7 +40,7 @@ MODULE pw_restart_new
PRIVATE
PUBLIC :: pw_write_schema, pw_write_binaries, pw_read_schema, &
read_collected_to_evc
PUBLIC :: readschema_ef, readschema_occupations, readschema_brillouin_zone
PUBLIC :: readschema_occupations
!
CONTAINS
!------------------------------------------------------------------------
@ -1003,91 +1003,19 @@ MODULE pw_restart_new
!
END SUBROUTINE pw_read_schema
!
!
!---------------------------------------------------------------------------
SUBROUTINE readschema_brillouin_zone( band_structure )
!---------------------------------------------------------------------------
!
USE lsda_mod, ONLY : lsda, isk
USE klist, ONLY : nkstot, xk, wk
USE start_k, ONLY : nks_start, xk_start, wk_start, &
nk1, nk2, nk3, k1, k2, k3
USE qes_types_module, ONLY : band_structure_type
!
IMPLICIT NONE
!
TYPE ( band_structure_type ),INTENT(IN) :: band_structure
INTEGER :: ik, isym, nks_
!
nks_ = band_structure%nks
nkstot = nks_
IF ( band_structure%lsda ) nkstot = nkstot * 2
!
!
DO ik = 1, nks_
xk(:,ik) = band_structure%ks_energies(ik)%k_point%k_point(:)
END DO
!! during lsda computations pw uses, for each k-point in the mesh, a distinct
!! k_point variable for the two spin channels, while in
!! the xml file only one k_point is present
IF ( band_structure%lsda ) THEN
DO ik = 1, nks_
xk(:,nks_+ik) = band_structure%ks_energies(ik)%k_point%k_point(:)
isk(ik) = 1
isk(ik+nks_) = 2
END DO
END IF
!
IF ( band_structure%starting_k_points%monkhorst_pack_ispresent ) THEN
nks_start = 0
nk1 = band_structure%starting_k_points%monkhorst_pack%nk1
nk2 = band_structure%starting_k_points%monkhorst_pack%nk2
nk3 = band_structure%starting_k_points%monkhorst_pack%nk3
k1 = band_structure%starting_k_points%monkhorst_pack%k1
k2 = band_structure%starting_k_points%monkhorst_pack%k2
k3 = band_structure%starting_k_points%monkhorst_pack%k3
ELSE IF (band_structure%starting_k_points%nk_ispresent ) THEN
nks_start = band_structure%starting_k_points%nk
IF ( nks_start > 0 ) THEN
IF ( .NOT. ALLOCATED(xk_start) ) ALLOCATE (xk_start(3,nks_start))
IF ( .NOT. ALLOCATED(wk_start) ) ALLOCATE (wk_start(nks_start))
IF ( nks_start == size( band_structure%starting_k_points%k_point ) ) THEN
DO ik =1, nks_start
xk_start(:,ik) = band_structure%starting_k_points%k_point(ik)%k_point(:)
IF ( band_structure%starting_k_points%k_point(ik)%weight_ispresent) THEN
wk_start(ik) = band_structure%starting_k_points%k_point(ik)%weight
ELSE
wk_start(ik) = 0.d0
END IF
END DO
ELSE
CALL infomsg ( "readschema_bz: ", &
"actual number of start kpoint not equal to nks_start, set nks_start=0")
nks_start = 0
END IF
END IF
ELSE
CALL errore ("readschema_bz: ", &
" no information found for initializing brillouin zone information", 1)
END IF
!
END SUBROUTINE readschema_brillouin_zone
!--------------------------------------------------------------------------------------------------
SUBROUTINE readschema_occupations( band_struct_obj )
!------------------------------------------------------------------------------------------------
!
USE lsda_mod, ONLY : lsda, nspin
USE fixed_occ, ONLY : tfixed_occ, f_inp
USE ktetra, ONLY : ntetra, tetra_type
USE klist, ONLY : ltetra, lgauss, ngauss, degauss, smearing
USE wvfct, ONLY : nbnd
USE input_parameters, ONLY : input_parameters_occupations => occupations
USE qes_types_module, ONLY : input_type, band_structure_type
USE qes_types_module, ONLY : band_structure_type
!
IMPLICIT NONE
!
TYPE ( band_structure_type ),INTENT(IN) :: band_struct_obj
INTEGER :: ispin, nk1, nk2, nk3, aux_dim1, aux_dim2
INTEGER :: nk1, nk2, nk3
!
lgauss = .FALSE.
ltetra = .FALSE.
@ -1270,27 +1198,5 @@ MODULE pw_restart_new
!
END SUBROUTINE read_collected_to_evc
!
!----------------------------------------------------------------------------------------
SUBROUTINE readschema_ef ( band_struct_obj )
!----------------------------------------------------------------------------------------
!
USE constants, ONLY : e2
USE ener, ONLY : ef, ef_up, ef_dw
USE klist, ONLY : two_fermi_energies, nelec
USE qes_types_module, ONLY : band_structure_type
!
IMPLICIT NONE
!
TYPE ( band_structure_type ),INTENT(IN) :: band_struct_obj
!
two_fermi_energies = band_struct_obj%two_fermi_energies_ispresent
nelec = band_struct_obj%nelec
IF ( two_fermi_energies) THEN
ef_up = band_struct_obj%two_fermi_energies(1)*e2
ef_dw = band_struct_obj%two_fermi_energies(2)*e2
ELSE IF ( band_struct_obj%fermi_energy_ispresent ) THEN
ef = band_struct_obj%fermi_energy*e2
END IF
END SUBROUTINE readschema_ef
!------------------------------------------------------------------------
END MODULE pw_restart_new

View File

@ -103,8 +103,10 @@ SUBROUTINE read_xml_file ( wfc_is_collected )
USE ions_base, ONLY : nat, nsp, ityp, amass, atm, tau, extfor
USE cell_base, ONLY : alat, at, bg, ibrav, celldm, omega
USE force_mod, ONLY : force
USE klist, ONLY : nks, nkstot, nelec, wk, tot_magnetization, &
nelup, neldw
USE klist, ONLY : nks, nkstot, xk, wk, tot_magnetization, &
nelec, nelup, neldw
USE start_k, ONLY : nks_start, xk_start, wk_start, &
nk1, nk2, nk3, k1, k2, k3
USE ener, ONLY : ef, ef_up, ef_dw
USE electrons_base, ONLY : nupdwn, set_nelup_neldw
USE wvfct, ONLY : npwx, nbnd, et, wg
@ -139,7 +141,7 @@ SUBROUTINE read_xml_file ( wfc_is_collected )
USE paw_variables, ONLY : okpaw
!
USE pw_restart_new, ONLY : pw_read_schema, &
readschema_occupations, readschema_brillouin_zone
readschema_occupations
USE qes_types_module,ONLY : output_type, parallel_info_type, &
general_info_type, input_type
USE qes_libs_module, ONLY : qes_reset
@ -147,8 +149,8 @@ SUBROUTINE read_xml_file ( wfc_is_collected )
qexsd_copy_algorithmic_info, qexsd_copy_atomic_species, &
qexsd_copy_atomic_structure, qexsd_copy_symmetry, &
qexsd_copy_basis_set, qexsd_copy_dft, qexsd_copy_efield, &
qexsd_copy_band_structure, qexsd_copy_magnetization
qexsd_copy_band_structure, qexsd_copy_magnetization, &
qexsd_copy_kpoints
#if defined(__BEOWULF)
USE qes_bcast_module,ONLY : qes_bcast
USE mp_images, ONLY : intra_image_comm
@ -219,6 +221,10 @@ SUBROUTINE read_xml_file ( wfc_is_collected )
ecutwfc = ecutwfc*e2
ecutrho = ecutrho*e2
dual = ecutrho/ecutwfc
! FIXME: next line ensures exact consistency between reciprocal and
! direct lattice vectors, preventing weird phonon symmetry errors
! (due to lousy algorithms, extraordinarily sensitive to tiny errors)
CALL recips ( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) )
!!
!! DFT section
CALL qexsd_copy_dft ( output_obj%dft, nsp, atm, &
@ -240,8 +246,8 @@ SUBROUTINE read_xml_file ( wfc_is_collected )
!! Band structure section
!! et and wg are allocated inside qexsd_copy_band_structure
CALL qexsd_copy_band_structure( output_obj%band_structure, lsda, &
nkstot, isk, natomwfc, nbnd, nupdwn(1), nupdwn(2), nelec, wk, wg, &
ef, ef_up, ef_dw, et )
nkstot, isk, natomwfc, nbnd, nupdwn(1), nupdwn(2), nelec, xk, &
wk, wg, ef, ef_up, ef_dw, et )
! convert to Ry
ef = ef*e2
ef_up = ef_up*e2
@ -273,7 +279,9 @@ SUBROUTINE read_xml_file ( wfc_is_collected )
END IF
!
CALL readschema_occupations( output_obj%band_structure )
CALL readschema_brillouin_zone( output_obj%band_structure )
!! Starting k-òoint information
CALL qexsd_copy_kpoints( output_obj%band_structure, nks_start, &
xk_start, wk_start, nk1, nk2, nk3, k1, k2, k3 )
!! Symmetry section
ALLOCATE ( irt(48,nat) )
IF ( lvalid_input ) THEN

View File

@ -49,7 +49,7 @@ SUBROUTINE setup()
tot_charge, tot_magnetization
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk, &
starting_magnetization
USE ener, ONLY : ef
USE ener, ONLY : ef, ef_up, ef_dw
USE electrons_base, ONLY : set_nelup_neldw
USE start_k, ONLY : nks_start, xk_start, wk_start, &
nk1, nk2, nk3, k1, k2, k3
@ -76,7 +76,8 @@ SUBROUTINE setup()
USE noncollin_module, ONLY : noncolin, npol, m_loc, i_cons, &
angle1, angle2, bfield, ux, nspin_lsda, &
nspin_gga, nspin_mag
USE pw_restart_new, ONLY : pw_read_schema, readschema_ef
USE pw_restart_new, ONLY : pw_read_schema
USE qexsd_copy, ONLY : qexsd_copy_efermi
USE qes_libs_module, ONLY : qes_reset
USE qes_types_module, ONLY : output_type, parallel_info_type, general_info_type
USE exx, ONLY : ecutfock, nbndproj
@ -95,8 +96,6 @@ SUBROUTINE setup()
LOGICAL, EXTERNAL :: check_para_diag
!
TYPE(output_type) :: output_obj
TYPE(parallel_info_type) :: parinfo_obj
TYPE(general_info_type) :: geninfo_obj
!
#if defined(__MPI)
LOGICAL :: lpara = .true.
@ -164,13 +163,12 @@ SUBROUTINE setup()
!
! ... in these cases, we need to read the Fermi energy
!
CALL pw_read_schema( ierr , output_obj, parinfo_obj, geninfo_obj )
CALL pw_read_schema( ierr , output_obj )
CALL errore( 'setup ', 'problem reading ef from file ' // &
& TRIM( tmp_dir ) // TRIM( prefix ) // '.save', ierr )
CALL readschema_ef ( output_obj%band_structure)
CALL qexsd_copy_efermi ( output_obj%band_structure, &
nelec, ef, two_fermi_energies, ef_up, ef_dw )
CALL qes_reset ( output_obj )
CALL qes_reset ( parinfo_obj )
CALL qes_reset ( geninfo_obj )
!
END IF
IF ( (lfcpopt .OR. lfcpdyn) .AND. restart ) THEN

View File

@ -135,29 +135,6 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
ENDDO
ENDDO
!
!
! DO k = 1, nrxx
! IF ( dft_is_meta() .AND. get_meta() /= 4 .AND. null_v(k) /= 0 ) THEN
! IF ( ABS(rhoaux(k,1))>epsr .AND. grho2(k,1)>epsg ) THEN
! !
! kedtau(k,1) = kedtau(k,1) / e2
! CALL tau_xc( rhoaux(k,1), grho2(k,1), kedtau(k,1), sx(k), sc(k), &
! v1x(k,1), v2x(k,1), v3x, v1c(k,1), v2c(k,1), v3c )
! kedtau(k,1) = kedtau(k,1) * e2
! !
! ENDIF
! ENDIF
! !
! DO l = 1, 3
! DO m = 1, l
! sigma_gradcorr(l,m) = sigma_gradcorr(l,m) + grho(l,k,1)*grho(m,k,1)* &
! e2 * (v2x(k,1) + v2c(k,1))
! ENDDO
! ENDDO
! !
! ENDDO
!
!
ELSEIF (nspin == 2) THEN
!
! This is the LSDA case