Skip to content

GitLab

Commit cbe386bb authored by Matteo Gerosa's avatar Matteo Gerosa
Browse files

Completed coulomb class, all code made consistent with changes.

parent 415f03af
Hide whitespace changes
Inline Side-by-side
Coulomb_kernel/class_coulomb.f90
View file @ cbe386bb
......@@ -23,9 +23,10 @@ MODULE class_coulomb
TYPE, PUBLIC :: coulomb
!
REAL(DP) :: div ! divergece
CHARACTER(LEN=*) :: singularity_removal_mode ! singularity_removal_mode
CHARACTER(LEN=7) :: singularity_removal_mode ! singularity_removal_mode
INTEGER :: iq ! q-point
REAL(DP),ALLOCATABLE :: sqvc(:) ! square root of Coulomb potential in PW
INTEGER :: numg, numgx
!
CONTAINS
!
......@@ -38,16 +39,17 @@ MODULE class_coulomb
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE sqvc_init(this,iq,singularity_removal_mode)
SUBROUTINE sqvc_init(this,cdriver,singularity_removal_mode,iq)
!-----------------------------------------------------------------------
!
! This routine computes results of applying sqVc to a potential
! associated with vector q. Coulomb cutoff technique can be used
!
USE kinds, ONLY : DP
USE constants, ONLY : pi, tpi, fpi, e2, eps8
USE cell_base, ONLY : alat, omega, at, bg, tpiba, tpiba2
USE constants, ONLY : fpi, e2, eps8
USE cell_base, ONLY : at, tpiba2
USE gvect, ONLY : g
USE gvecs, ONLY : ngms
USE westcom, ONLY : igq_q,npwqx,npwq
USE types_bz_grid, ONLY : q_grid
USE control_flags, ONLY : gamma_only
......@@ -57,26 +59,46 @@ MODULE class_coulomb
! I/O
!
CLASS(coulomb) :: this
CHARACTER(LEN=*), INTENT(IN) :: cdriver
CHARACTER(LEN=*), INTENT(IN) :: singularity_removal_mode
INTEGER, INTENT(IN) :: iq
INTEGER, INTENT(IN), OPTIONAL :: iq
!
! Workspace
!
REAL(DP) :: qgnorm2,qg(3),x
INTEGER :: numg, numgx
INTEGER :: ig, ipol
LOGICAL :: on_double_grid, l_print
REAL(DP) :: grid_factor
!
CALL start_clock('storesqvc')
CALL start_clock('sqvc_init')
!
SELECT CASE ( cdriver )
CASE ( 'Wave' )
this%numg = npwq
this%numgx = npwqx
CASE ( 'Smooth' )
this%numg = ngms
this%numgx = ngms
CASE DEFAULT
CALL errore("sqvc_init", "cdriver value not supported, supported only Wave and Smooth",1)
END SELECT
!
IF( ALLOCATED(this%sqvc) ) DEALLOCATE( this%sqvc )
ALLOCATE( this%sqvc( npwqx ) )
ALLOCATE( this%sqvc( this%numgx ) )
!
IF ( PRESENT(iq) ) THEN
this%iq = iq
ELSE
this%iq = 1 ! gamma-only
ENDIF
!
this%iq = iq
this%singularity_removal_mode = singularity_removal_mode
this%singularity_removal_mode = TRIM(singularity_removal_mode)
IF (this%singularity_removal_mode /= "gb" .AND. this%singularity_removal_mode /= "default") &
& CALL errore( 'sqvc_init', 'singularity removal mode not supported, supported only default and gb', 1 )
!
this%sqvc = 0._DP
DO ig = 1,npwq
DO ig = 1,this%numg
!
IF ( gamma_only ) THEN
qg(:) = g(:,ig)
......@@ -109,32 +131,26 @@ MODULE class_coulomb
!
ENDDO
!
this%div = this%compute_divergence()
IF ( q_grid%l_pIsGamma(this%iq) ) this%div = this%compute_divergence()
!
CALL stop_clock('storesqvc')
CALL stop_clock('sqvc_init')
!
END SUBROUTINE
!
SUBROUTINE print_divergence( this )
!
! I/O
!
CLASS(coulomb) :: this
!
WRITE(stdout,"(5X,'Divergence = ',es14.6)") div
!
END SUBROUTINE
!
!
FUNCTION compute_divergence( this ) RESULT( div )
!
USE io_global, ONLY : stdout
USE constants, ONLY : pi, tpi, fpi, e2, eps8
USE cell_base, ONLY : omega, at, bg, tpiba2
USE mp, ONLY : mp_sum
USE mp_global, ONLY : intra_bgrp_comm
USE mp_world, ONLY : mpime, world_comm, nproc
USE control_flags, ONLY : gamma_only
USE gvecw, ONLY : ecutwfc
USE random_numbers, ONLY : randy
USE gvect, ONLY : g
USE westcom, ONLY : ngq,igq_q
USE types_bz_grid, ONLY : q_grid
!
! I/O
!
......@@ -149,13 +165,11 @@ MODULE class_coulomb
INTEGER :: i1, i2, i3, iq, ig, ipol
REAL(DP) :: prod(3,3), qhelp, edge(3), qbz(3), rand, qmo, vbz, vhelp, intcounter, x
!
div = 0
!
IF ( .NOT. q_grid%l_pIsGamma(this%iq) ) RETURN
div = 0._DP
!
SELECT CASE( this%singularity_removal_mode )
!
CASE("spherical")
CASE("default")
!
! In this case we use the spherical region
!
......@@ -238,8 +252,8 @@ MODULE class_coulomb
!
DO iq = 1, q_grid%np
!
DO ig = 1,ngq(iq) ! MATTEO CONTROLLA
qg(:) = q_grid%p_cart(:,iq) + g(:,igq_q(ig,iq)) ! MATTEO CONTROLLA
DO ig = 1,ngq(iq)
qg(:) = q_grid%p_cart(:,iq) + g(:,igq_q(ig,iq))
qgnorm2 = SUM( qg(:)**2 ) * tpiba2
on_double_grid = .TRUE.
DO ipol = 1,3
......@@ -265,7 +279,22 @@ MODULE class_coulomb
!
END SELECT
!
END SUBROUTINE
END FUNCTION
!
!
SUBROUTINE print_divergence( this )
!
USE io_global, ONLY : stdout
USE types_bz_grid, ONLY : q_grid
!
! I/O
!
CLASS(coulomb) :: this
!
IF ( .NOT. q_grid%l_pIsGamma(this%iq) ) RETURN
WRITE(stdout,"(5X,'Divergence = ',es14.6)") this%div
!
END SUBROUTINE
!
!
!
......
DFPT_kernel/Makefile
View file @ cbe386bb
......@@ -8,6 +8,7 @@ MODFLAGS= $(MOD_FLAG)../../iotk/src $(MOD_FLAG)../../Modules $(MOD_FLAG)../../LA
$(MOD_FLAG)../Modules \
$(MOD_FLAG)../Tools \
$(MOD_FLAG)../FFT_kernel \
$(MOD_FLAG)../Coulomb_kernel \
$(MOD_FLAG).
IFLAGS=
......
DFPT_kernel/dfpt.f90
View file @ cbe386bb
......@@ -36,9 +36,10 @@ SUBROUTINE dfpt (m,dvg,dng,tr2)
USE io_files, ONLY : tmp_dir, nwordwfc, iunwfc, diropn
USE uspp, ONLY : nkb, vkb, okvan
USE constants, ONLY : e2,fpi
USE westcom, ONLY : npwq,sqvc,nbnd_occ,iuwfc,lrwfc,npwqx,fftdriver
USE westcom, ONLY : npwq,nbnd_occ,iuwfc,lrwfc,npwqx,fftdriver
USE io_push, ONLY : io_push_title
USE mp_world, ONLY : mpime,world_comm
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -149,7 +150,7 @@ SUBROUTINE dfpt (m,dvg,dng,tr2)
aux_r = 0._DP
!
DO ig = 1, npwq ! perturbation acts only on body
aux_g(ig) = dvg(ig,ipert) * sqvc(ig)
aux_g(ig) = dvg(ig,ipert) * pot3D%sqvc(ig)
ENDDO
!
IF(gamma_only) THEN
......@@ -290,7 +291,7 @@ SUBROUTINE dfpt (m,dvg,dng,tr2)
ENDIF
!
DO ig=1,npwq ! pert acts only on body
dng(ig,ipert) = dng(ig,ipert) + 2._DP * wk(iks) * aux_g(ig) * sqvc(ig) / omega
dng(ig,ipert) = dng(ig,ipert) + 2._DP * wk(iks) * aux_g(ig) * pot3D%sqvc(ig) / omega
ENDDO
!
DEALLOCATE(aux_g)
......@@ -354,11 +355,12 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
USE control_flags, ONLY : gamma_only, io_level
USE io_files, ONLY : tmp_dir, nwordwfc, iunwfc, diropn
USE uspp, ONLY : nkb, vkb, okvan
USE westcom, ONLY : sqvc,nbnd_occ,iuwfc,lrwfc,npwqx,npwq,igq_q
USE westcom, ONLY : nbnd_occ,iuwfc,lrwfc,npwqx,npwq,igq_q
USE io_push, ONLY : io_push_title
USE mp_world, ONLY : mpime,world_comm
USE class_bz_grid, ONLY : bz_grid
USE types_bz_grid, ONLY : k_grid, q_grid, compute_phase
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -489,7 +491,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
aux_r = (0._DP,0._DP)
!
DO ig = 1, npwq ! perturbation acts only on body
aux_g(ig) = dvg(ig,ipert) * sqvc(ig)
aux_g(ig) = dvg(ig,ipert) * pot3D%sqvc(ig)
ENDDO
!
! inverse Fourier transform of the perturbation: (q+)G ---> R
......@@ -621,7 +623,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
CALL single_fwfft_k(dffts,npwq,npwqx,aux_r,aux_g,'Wave',igq_q(1,iq))
!
DO ig = 1, npwq ! pert acts only on body
dng(ig,ipert) = dng(ig,ipert) + 2._DP * wk(iks) * aux_g(ig) * sqvc(ig) / omega
dng(ig,ipert) = dng(ig,ipert) + 2._DP * k_grid%weight(iks) * aux_g(ig) * pot3D%sqvc(ig) / omega
ENDDO
!
DEALLOCATE( aux_g )
......
Modules/westcom.f90
View file @ cbe386bb
......@@ -19,7 +19,7 @@ MODULE scratch_area
SAVE
!
! COULOMB
REAL(DP),ALLOCATABLE :: sqvc(:)
! REAL(DP),ALLOCATABLE :: sqvc(:)
INTEGER :: npwq,npwqx,npwq_g
CHARACTER(LEN=6) :: fftdriver
INTEGER,ALLOCATABLE :: iks_l2g(:)
......@@ -70,7 +70,7 @@ MODULE scratch_area
! I/O
!INTEGER :: io_comm ! communicator for head of images (me_bgrp==0)
!
REAL(DP) :: isz
! REAL(DP) :: isz
!
END MODULE
!
......
Tools/Makefile
View file @ cbe386bb
......@@ -32,7 +32,6 @@ clean_scratchfiles.o \
set_nbndocc.o \
get_alpha_pv.o \
set_iks_l2g.o \
set_isz.o \
v_x.o \
gradcorr_x.o \
west_environment.o \
......
Westpp/do_sxx.f90
View file @ cbe386bb
......@@ -39,7 +39,8 @@ SUBROUTINE do_sxx ( )
USE mp_world, ONLY : mpime,root
USE constants, ONLY : rytoev
USE json_module, ONLY : json_file
USE coulomb, ONLY : store_sqvc
USE class_coulomb, ONLY : coulomb
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -47,9 +48,9 @@ SUBROUTINE do_sxx ( )
!
INTEGER :: ir, ip, ig, iks, ib, iv, ip_glob
COMPLEX(DP),ALLOCATABLE :: pertg(:),pertr(:),pertr_nc(:,:)
REAL(DP),ALLOCATABLE :: mysqvc(:)
! REAL(DP),ALLOCATABLE :: mysqvc(:)
TYPE(bar_type) :: barra
REAL(DP) :: mydiv
! REAL(DP) :: mydiv
REAL(DP),ALLOCATABLE :: sigma_exx( :, : )
REAL(DP),ALLOCATABLE :: sigma_sxx( :, : )
REAL(DP) :: peso
......@@ -61,8 +62,9 @@ SUBROUTINE do_sxx ( )
TYPE(json_file) :: json
INTEGER :: iunit
!
ALLOCATE( mysqvc(npwq) )
CALL store_sqvc(mysqvc,npwq,'spherical',1,.FALSE.,mydiv)
CALL pot3d%init('Wave','gb')
! ALLOCATE( mysqvc(npwq) )
! CALL store_sqvc(mysqvc,npwq,'spherical',1,.FALSE.,mydiv)
!CALL store_sqvc(mysqvc,npwq,1,mydiv)
!
CALL io_push_title("(S)creened eXact eXchange")
......@@ -155,11 +157,11 @@ SUBROUTINE do_sxx ( )
ENDIF
!
DO ig = 1,npwq
pertg(ig) = pertg(ig) * mysqvc(ig)
pertg(ig) = pertg(ig) * 3Dpot%sqvc(ig)
ENDDO
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - peso * DDOT( 2*npwq, pertg(1), 1, pertg(1), 1) / omega
!IF(gstart==2) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) + REAL( pertg(1), KIND = DP )**2 / omega
IF( ib == iv .AND. gstart == 2 ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - mydiv
IF( ib == iv .AND. gstart == 2 ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - pot3D%div
!
! -- < SXX >
IF( gamma_only ) THEN
......@@ -169,7 +171,7 @@ SUBROUTINE do_sxx ( )
ip_glob = pert%l2g(ip)
sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - dproj(1,ip)**2 * (ev(ip_glob)/(1._DP-ev(ip_glob))) / omega
ENDDO
IF( ib == iv ) sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - (1._DP/westpp_epsinfty-1._DP) * mydiv
IF( ib == iv ) sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - (1._DP/westpp_epsinfty-1._DP) * pot3D%div
ELSE
CALL glbrak_k( pertg, dvg, zproj, npwq, npwqx, 1, pert%nloc, 1, npol)
CALL mp_sum( zproj, intra_bgrp_comm )
......@@ -178,7 +180,7 @@ SUBROUTINE do_sxx ( )
sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - REAL(zproj(1,ip)*CONJG(zproj(1,ip)),KIND=DP) &
& * (ev(ip_glob)/(1._DP-ev(ip_glob))) / omega
ENDDO
IF( ib == iv ) sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - (1._DP/westpp_epsinfty-1._DP) * mydiv
IF( ib == iv ) sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - (1._DP/westpp_epsinfty-1._DP) * pot3D%div
ENDIF
! -- </ SXX >
!
......@@ -198,7 +200,7 @@ SUBROUTINE do_sxx ( )
sigma_sxx = sigma_exx + sigma_sxx
!
DEALLOCATE( pertg )
DEALLOCATE( mysqvc )
! DEALLOCATE( mysqvc )
IF( noncolin ) THEN
DEALLOCATE( pertr_nc )
ELSE
......
Westpp/westpp_setup.f90
View file @ cbe386bb
......@@ -14,9 +14,9 @@
SUBROUTINE westpp_setup
!-----------------------------------------------------------------------
!
USE westcom, ONLY : alphapv_dfpt,npwq,sqvc,west_prefix,westpp_save_dir,&
USE westcom, ONLY : alphapv_dfpt,npwq,west_prefix,westpp_save_dir,&
& n_imfreq,nbnd_occ,l_macropol,macropol_calculation,&
& n_refreq,isz,qp_bandrange,westpp_calculation,westpp_n_pdep_eigen_to_use
& n_refreq,qp_bandrange,westpp_calculation,westpp_n_pdep_eigen_to_use
USE mp, ONLY : mp_max
USE mp_global, ONLY : intra_bgrp_comm
USE pwcom, ONLY : nbnd
......@@ -28,7 +28,8 @@ SUBROUTINE westpp_setup
USE wavefunctions_module, ONLY : evc
USE mod_mpiio, ONLY : set_io_comm
USE pdep_db, ONLY : pdep_db_read
USE coulomb, ONLY : store_sqvc
USE class_coulomb, ONLY : coulomb
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -36,15 +37,16 @@ SUBROUTINE westpp_setup
REAL(DP) :: qq
COMPLEX(DP),EXTERNAL :: get_alpha_pv
INTEGER :: ig, i
LOGICAL :: l_printout_div = .TRUE.
!
CALL do_setup ( )
!
CALL set_npwq()
!
ALLOCATE(sqvc(npwq))
! ALLOCATE(sqvc(npwq))
!
CALL store_sqvc(sqvc,npwq,'spherical',1,.FALSE.,isz,l_printout_div)
CALL pot3D%init('Wave','default')
CALL pot3D%print_divergence()
! CALL store_sqvc(sqvc,npwq,'spherical',1,.FALSE.,isz,l_printout_div)
!CALL store_sqvc(sqvc,npwq,1,isz)
!
CALL set_nbndocc()
......
Wfreq/calc_corr.f90
View file @ cbe386bb
......@@ -24,13 +24,12 @@ SUBROUTINE calc_corr_gamma( sigma_corr, energy, l_verbose)
USE cell_base, ONLY : omega
USE constants, ONLY : tpi,fpi,rytoev,e2,pi
USE pwcom, ONLY : et,nks,current_spin,isk,xk,nbnd,lsda,g2kin,nspin,current_k,wk
USE westcom, ONLY : qp_bandrange,isz,&
& nbnd_occ,l_enable_lanczos,&
& n_lanczos,iks_l2g,l_macropol,&
USE westcom, ONLY : qp_bandrange,nbnd_occ,l_enable_lanczos,n_lanczos,iks_l2g,l_macropol,&
& d_head_ifr,z_head_rfr,d_body1_ifr,d_body2_ifr,d_diago,z_body_rfr
USE bar, ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
USE io_push, ONLY : io_push_bar,io_push_value,io_push_title
USE distribution_center, ONLY : pert,ifr,rfr,aband
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -129,7 +128,7 @@ SUBROUTINE calc_corr_gamma( sigma_corr, energy, l_verbose)
CALL mp_sum( partial_b, intra_bgrp_comm)
CALL mp_sum( partial_b, inter_image_comm)
!
sigma_corr(ib,iks) = sigma_corr(ib,iks) + CMPLX( partial_b/omega/pi + partial_h*isz/pi, 0._DP, KIND=DP )
sigma_corr(ib,iks) = sigma_corr(ib,iks) + CMPLX( partial_b/omega/pi + partial_h*pot3D%div/pi, 0._DP, KIND=DP )
!
IF(l_verbose) CALL update_bar_type( barra, 'sigmac_i', 1 )
!
......@@ -198,7 +197,7 @@ SUBROUTINE calc_corr_gamma( sigma_corr, energy, l_verbose)
CALL mp_sum( residues_b, intra_bgrp_comm )
CALL mp_sum( residues_b, inter_image_comm )
!
sigma_corr(ib,iks) = sigma_corr(ib,iks) + residues_b/omega + residues_h*isz
sigma_corr(ib,iks) = sigma_corr(ib,iks) + residues_b/omega + residues_h*pot3D%div
!
IF(l_verbose) CALL update_bar_type( barra, 'sigmac_r', 1 )
!
......@@ -225,15 +224,14 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
USE cell_base, ONLY : omega
USE constants, ONLY : tpi,fpi,rytoev,e2,pi
USE pwcom, ONLY : et,nks,current_spin,isk,xk,nbnd,lsda,g2kin,nspin,current_k,wk
USE westcom, ONLY : qp_bandrange,isz,&
& nbnd_occ,l_enable_lanczos,&
& n_lanczos,iks_l2g,l_macropol,&
USE westcom, ONLY : qp_bandrange,nbnd_occ,l_enable_lanczos,n_lanczos,iks_l2g,l_macropol,&
& z_head_ifr,z_head_rfr,z_body1_ifr_q,z_body2_ifr_q,d_diago_q,z_body_rfr_q
USE bar, ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
USE io_push, ONLY : io_push_bar,io_push_value,io_push_title
USE distribution_center, ONLY : pert,ifr,rfr,aband
USE class_bz_grid, ONLY : bz_grid
USE types_bz_grid, ONLY : k_grid, q_grid
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -350,7 +348,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
CALL mp_sum( partial_b, intra_bgrp_comm)
CALL mp_sum( partial_b, inter_image_comm)
!
sigma_corr(ib,iks) = sigma_corr(ib,iks) + partial_b/omega/pi + partial_h*isz/pi
sigma_corr(ib,iks) = sigma_corr(ib,iks) + partial_b/omega/pi + partial_h*pot3D%div/pi
!
IF(l_verbose) CALL update_bar_type( barra, 'sigmac_i', 1 )
!
......@@ -434,7 +432,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
CALL mp_sum( residues_b, intra_bgrp_comm )
CALL mp_sum( residues_b, inter_image_comm )
!
sigma_corr(ib,iks) = sigma_corr(ib,iks) + residues_b/omega + residues_h*isz
sigma_corr(ib,iks) = sigma_corr(ib,iks) + residues_b/omega + residues_h*pot3D%div
!
IF(l_verbose) CALL update_bar_type( barra, 'sigmac_r', 1 )
!
......
Wfreq/calc_exx2.f90
View file @ cbe386bb
......@@ -46,7 +46,8 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
USE class_idistribute, ONLY : idistribute
USE coulomb_vcut_module, ONLY : vcut_init, vcut_type, vcut_info, &
vcut_get, vcut_spheric_get, vcut_destroy
USE coulomb, ONLY : store_sqvc
USE class_coulomb, ONLY : coulomb
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -64,10 +65,8 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
TYPE(idistribute) :: vband
TYPE(bar_type) :: barra
INTEGER :: barra_load
REAL(DP),ALLOCATABLE :: mysqvc(:)
REAL(DP) :: ecutvcut
TYPE(vcut_type) :: vcut
REAL(DP) :: mydiv
! REAL(DP),ALLOCATABLE :: mysqvc(:)
! REAL(DP) :: mydiv
!
WRITE(stdout,'(5x,a)') ' '
CALL io_push_bar()
......@@ -75,14 +74,15 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
CALL io_push_bar()
!
ALLOCATE( pertg( ngms ) )
ALLOCATE( mysqvc(ngms) )
! ALLOCATE( mysqvc(ngms) )
IF(noncolin) THEN
ALLOCATE( pertr_nc( dffts%nnr, npol ) )
ELSE
ALLOCATE( pertr( dffts%nnr ) )
ENDIF
!
CALL store_sqvc( mysqvc, ngms, 'gb', 1, .FALSE., mydiv )
CALL pot3D%init('Smooth','gb')
! CALL store_sqvc( mysqvc, ngms, 'gb', 1, .FALSE., mydiv )
!CALL store_sqvc(mysqvc,ngms,div_kind_hf,mydiv)
!
! Set to zero
......@@ -178,11 +178,11 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
ENDIF
!
DO ig = 1,ngms
pertg(ig) = pertg(ig) * mysqvc(ig)
pertg(ig) = pertg(ig) * pot3D%sqvc(ig)
ENDDO
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - 2._DP * DDOT( 2*ngms, pertg(1), 1, pertg(1), 1) / omega
!IF(gstart==2) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) + REAL( pertg(1), KIND = DP )**2 / omega
IF( ib == iv_glob .AND. gstart == 2 ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - mydiv
IF( ib == iv_glob .AND. gstart == 2 ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - pot3D%div
!
ENDDO
!
......@@ -198,7 +198,7 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
CALL mp_sum( sigma_exx, inter_image_comm )
!
DEALLOCATE( pertg )
DEALLOCATE( mysqvc )
! DEALLOCATE( mysqvc )
IF( noncolin ) THEN
DEALLOCATE( pertr_nc )
ELSE
......@@ -230,7 +230,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
USE fft_at_gamma, ONLY : single_invfft_gamma,single_fwfft_gamma
USE fft_at_k, ONLY : single_invfft_k,single_fwfft_k
USE wavefunctions_module, ONLY : evc,psic,psic_nc
USE westcom, ONLY : iuwfc,lrwfc,npwq,nbnd_occ,div_kind_hf
USE westcom, ONLY : iuwfc,lrwfc,npwq,nbnd_occ
USE control_flags, ONLY : gamma_only
USE noncollin_module, ONLY : noncolin,npol
USE buffers, ONLY : get_buffer
......@@ -245,7 +245,8 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
vcut_get, vcut_spheric_get, vcut_destroy
USE class_bz_grid, ONLY : bz_grid
USE types_bz_grid, ONLY : k_grid, q_grid, compute_phase
USE coulomb, ONLY : store_sqvc
USE class_coulomb, ONLY : coulomb
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
......@@ -266,11 +267,9 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
TYPE(idistribute) :: vband
TYPE(bar_type) :: barra
INTEGER :: barra_load
REAL(DP),ALLOCATABLE :: mysqvc(:)
! REAL(DP),ALLOCATABLE :: mysqvc(:)
LOGICAL :: l_gammaq
REAL(DP) :: ecutvcut
TYPE(vcut_type) :: vcut
REAL(DP) :: mydiv
! REAL(DP) :: mydiv
REAL(DP) :: g0(3)
!
WRITE(stdout,'(5x,a)') ' '
......@@ -279,7 +278,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
CALL io_push_bar()
!
ALLOCATE( pertg( ngms ) )
ALLOCATE( mysqvc( ngms ) )
! ALLOCATE( mysqvc( ngms ) )
ALLOCATE( phase(dffts%nnr) )
ALLOCATE( evckmq(npwx*npol,nbnd) )
IF(noncolin) THEN
......@@ -369,7 +368,8 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
!ikqs = k_grid%find( kmq, 'cart' )
!!ikqs = kmq_grid%index_kq(iks,iq)
!
CALL store_sqvc( mysqvc, ngms, 'gb', iq, .FALSE., mydiv )
CALL pot3D%init('Smooth', 'gb', iq)
! CALL store_sqvc( mysqvc, ngms, 'gb', iq, .FALSE., mydiv )
!IF ( l_gammaq ) THEN
! CALL store_sqvc(mysqvc,ngms,div_kind_hf,mydiv)
!ELSE
......@@ -406,11 +406,11 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
ENDIF
!
DO ig = 1,ngms
pertg(ig) = pertg(ig) * mysqvc(ig)
pertg(ig) = pertg(ig) * pot3D%sqvc(ig)
ENDDO
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - DDOT( 2*ngms, pertg(1), 1, pertg(1), 1)/omega*q_grid%weight(iq)
!IF(gstart==2) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) + REAL( pertg(1), KIND = DP )**2 / omega
IF( ib == iv_glob .AND. gstart == 2 .AND. l_gammaq ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - mydiv
IF( ib == iv_glob .AND. gstart == 2 .AND. l_gammaq ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - pot3D%div
!
ENDDO ! iv
!
......@@ -428,7 +428,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )