Commit 34242da1 authored by Matteo Gerosa's avatar Matteo Gerosa
Browse files

Changes to fetch input and coulomb kernel (fft driver and igq map) consolidated and tested.

parent d55d9e8f
......@@ -24,6 +24,8 @@ MODULE class_coulomb
!
REAL(DP) :: div ! divergence
CHARACTER(LEN=7) :: singularity_removal_mode ! singularity_removal_mode
CHARACTER(LEN=5) :: cdriver ! FFT driver = "Wave", "Dense"
LOGICAL :: l_use_igq ! use igq map
INTEGER :: iq ! q-point
REAL(DP),ALLOCATABLE :: sqvc(:) ! square root of Coulomb potential in PW
!
......@@ -38,7 +40,7 @@ MODULE class_coulomb
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE sqvc_init(this,cdriver,singularity_removal_mode,iq)
SUBROUTINE sqvc_init(this,cdriver,l_use_igq,singularity_removal_mode,iq)
!-----------------------------------------------------------------------
!
! This routine computes results of applying sqVc to a potential
......@@ -47,8 +49,7 @@ MODULE class_coulomb
USE kinds, ONLY : DP
USE constants, ONLY : fpi, e2, eps8
USE cell_base, ONLY : at, tpiba2
USE gvect, ONLY : g
USE gvecs, ONLY : ngms
USE gvect, ONLY : g, ngm
USE westcom, ONLY : igq_q,npwqx,npwq
USE types_bz_grid, ONLY : q_grid
USE control_flags, ONLY : gamma_only
......@@ -59,6 +60,7 @@ MODULE class_coulomb
!
CLASS(coulomb) :: this
CHARACTER(LEN=*), INTENT(IN) :: cdriver
LOGICAL, INTENT(IN) :: l_use_igq
CHARACTER(LEN=*), INTENT(IN) :: singularity_removal_mode
INTEGER, INTENT(IN), OPTIONAL :: iq
!
......@@ -72,37 +74,54 @@ MODULE class_coulomb
!
CALL start_clock('sqvc_init')
!
SELECT CASE ( cdriver )
this%cdriver = cdriver
this%l_use_igq = l_use_igq
this%singularity_removal_mode = singularity_removal_mode
IF ( PRESENT(iq) ) THEN
this%iq = iq
ELSE
this%iq = 1 ! gamma-only
ENDIF
!
! ... Check compatibility between singularity removal mode and FFT driver
!
SELECT CASE ( this%singularity_removal_mode )
CASE ( 'default' )
CASE ( 'gb' )
IF ( this%cdriver == 'Wave' ) CALL errore("sqvc_init", "gb singularity removal mode requires Dense grid",1)
CASE DEFAULT
CALL errore( 'sqvc_init', 'singularity removal mode not supported, supported only default and gb', 1 )
END SELECT
!
! ... Check compatibility between FFT driver and use of igq map
!
SELECT CASE ( this%cdriver )
CASE ( 'Wave' )
IF (.NOT.gamma_only .AND. .NOT.this%l_use_igq) THEN
CALL errore("sqvc_init", "q-points case requires igq map when using Wave grid",1)
ELSEIF (gamma_only .AND. this%l_use_igq) THEN
CALL errore("sqvc_init", "igq map not needed in gamma-only case",1)
ENDIF
numg = npwq
numgx = npwqx
CASE ( 'Smooth' )
numg = ngms
numgx = ngms
CASE ( 'Dense' )
IF (this%l_use_igq) CALL errore("sqvc_init", "igq map not used with Dense grid",1)
numg = ngm
numgx = ngm
CASE DEFAULT
CALL errore("sqvc_init", "cdriver value not supported, supported only Wave and Smooth",1)
CALL errore("sqvc_init", "cdriver value not supported, supported only Wave and Dense",1)
END SELECT
!
IF( ALLOCATED(this%sqvc) ) DEALLOCATE( this%sqvc )
ALLOCATE( this%sqvc( numgx ) )
!
IF ( PRESENT(iq) ) THEN
this%iq = iq
ELSE
this%iq = 1 ! gamma-only
ENDIF
!
this%singularity_removal_mode = 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,numg
!
IF ( gamma_only .OR. cdriver == 'Smooth') THEN ! if smooth grid is used (X), no mapping is needed
qg(:) = g(:,ig) + q_grid%p_cart(:,this%iq)
ELSE
IF ( this%l_use_igq ) THEN
qg(:) = g(:,igq_q(ig,this%iq)) + q_grid%p_cart(:,this%iq)
ELSE
qg(:) = g(:,ig) + q_grid%p_cart(:,this%iq)
ENDIF
!
qgnorm2 = SUM( qg(:)**2 ) * tpiba2
......@@ -147,8 +166,7 @@ MODULE class_coulomb
USE control_flags, ONLY : gamma_only
USE gvecw, ONLY : ecutwfc
USE random_numbers, ONLY : randy
USE gvect, ONLY : g
USE gvecs, ONLY : ngms
USE gvect, ONLY : g, ngm
USE types_bz_grid, ONLY : q_grid
!
! I/O
......@@ -251,7 +269,7 @@ MODULE class_coulomb
!
DO iq = 1, q_grid%np
!
DO ig = 1,ngms
DO ig = 1,ngm
qg(:) = q_grid%p_cart(:,iq) + g(:,ig)
qgnorm2 = SUM( qg(:)**2 ) * tpiba2
on_double_grid = .TRUE.
......
......@@ -20,7 +20,6 @@ SUBROUTINE dfpt (m,dvg,dng,tr2)
USE fft_base, ONLY : dfftp,dffts
USE gvect, ONLY : nl,gstart,g,ngm
USE wavefunctions_module, ONLY : evc,psic
USE gvecs, ONLY : ngms
USE gvecw, ONLY : gcutw
USE mp, ONLY : mp_sum,mp_barrier,mp_bcast
USE mp_global, ONLY : inter_image_comm,inter_pool_comm,my_image_id
......@@ -340,7 +339,6 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
USE fft_base, ONLY : dfftp,dffts
USE gvect, ONLY : nl,nl,gstart,g,ngm
USE wavefunctions_module, ONLY : evc,psic
USE gvecs, ONLY : ngms,nls
USE gvecw, ONLY : gcutw
USE mp, ONLY : mp_sum,mp_barrier,mp_bcast
USE mp_global, ONLY : inter_image_comm,inter_pool_comm,my_image_id
......
......@@ -118,7 +118,7 @@ MODULE class_bz_grid
DO iq2 = 1, this%ngrid(2)
DO iq3 = 1, this%ngrid(3)
ip = ip + 1
IF ( ANY(qlist(:)) == ip ) THEN
IF ( ANY(qlist(:) == ip) ) THEN
iqlist = iqlist + 1
this%p_cryst(1,iqlist) = DBLE( iq1 - 1 ) / DBLE( this%ngrid(1) )
this%p_cryst(2,iqlist) = DBLE( iq2 - 1 ) / DBLE( this%ngrid(2) )
......
......@@ -11,7 +11,7 @@
! Marco Govoni
!
!-----------------------------------------------------------------------
SUBROUTINE fetch_input( driver, verbose )
SUBROUTINE fetch_input( num_drivers, driver, verbose )
!-----------------------------------------------------------------------
!
USE json_module, ONLY : json_file
......@@ -30,8 +30,8 @@ SUBROUTINE fetch_input( driver, verbose )
!
! I/O
!
!INTEGER, INTENT(IN) :: num_drivers
INTEGER, INTENT(IN) :: driver(:)
INTEGER, INTENT(IN) :: num_drivers
INTEGER, INTENT(IN) :: driver(num_drivers)
LOGICAL, INTENT(IN) :: verbose
!
! Workspace
......@@ -84,9 +84,10 @@ SUBROUTINE fetch_input( driver, verbose )
l_minimize_exx_if_active = .FALSE.
l_use_ecutrho = .FALSE.
IF( ALLOCATED(qlist) ) DEALLOCATE(qlist)
ALLOCATE( qlist(1) )
qlist(1) = 1
IF( .NOT. gamma_only ) THEN
IF ( gamma_only ) THEN
ALLOCATE( qlist(1) )
qlist(1) = 1
ELSE
ALLOCATE(qlist(nk1*nk2*nk3))
qlist = (/ (i, i=1,nk1*nk2*nk3,1) /)
ENDIF
......@@ -180,7 +181,7 @@ SUBROUTINE fetch_input( driver, verbose )
IF( found ) THEN
IF( ALLOCATED(qlist) ) DEALLOCATE(qlist)
ALLOCATE(qlist(SIZE(ivec)))
qlist = ivec
qlist(1:SIZE(ivec)) = ivec(1:SIZE(ivec))
ENDIF
ENDIF
!
......@@ -303,6 +304,11 @@ SUBROUTINE fetch_input( driver, verbose )
IF(tr2_dfpt<=0._DP) CALL errore('fetch_input','Err: tr2_dfpt<0.',1)
IF(trev_pdep<=0._DP) CALL errore('fetch_input','Err: trev_pdep<0.',1)
IF(trev_pdep_rel<=0._DP) CALL errore('fetch_input','Err: trev_pdep_rel<0.',1)
IF(gamma_only) THEN
IF (SIZE(qlist)/=1) CALL errore('fetch_input','Err: SIZE(qlist)/=1.',1)
ELSE
IF (SIZE(qlist)>nk1*nk2*nk3) CALL errore('fetch_input','Err: SIZE(qlist)>nk1*nk2*nk3.',1)
ENDIF
!
ENDIF
!
......@@ -484,7 +490,7 @@ SUBROUTINE fetch_input( driver, verbose )
CALL json%initialize()
CALL json%load_file(filename=TRIM(logfile))
!
CALL add_intput_parameters_to_json_file( driver, json )
CALL add_intput_parameters_to_json_file( num_drivers, driver, json )
!
OPEN( NEWUNIT=iunit, FILE=TRIM(logfile) )
CALL json%print_file( iunit )
......@@ -500,7 +506,7 @@ SUBROUTINE fetch_input( driver, verbose )
END SUBROUTINE
!
!-----------------------------------------------------------------------
SUBROUTINE add_intput_parameters_to_json_file( driver, json )
SUBROUTINE add_intput_parameters_to_json_file( num_drivers, driver, json )
!-----------------------------------------------------------------------
!
USE json_module, ONLY : json_file
......@@ -517,8 +523,8 @@ SUBROUTINE add_intput_parameters_to_json_file( driver, json )
!
! I/O
!
!INTEGER, INTENT(IN) :: num_drivers
INTEGER, INTENT(IN) :: driver(:)
INTEGER, INTENT(IN) :: num_drivers
INTEGER, INTENT(IN) :: driver(num_drivers)
TYPE(json_file), INTENT(INOUT) :: json
!
IF ( mpime == root ) THEN
......
......@@ -54,8 +54,8 @@ SUBROUTINE set_npwq()
!
ELSE
!
IF( l_use_ecutrho ) CALL errore("set_npwq", "Dense grid not implemented with q-points",1)
fftdriver = 'Wave'
! 'Dense' grid not yet implemented
!
npwqx = n_plane_waves( gcutw, q_grid%np, q_grid%p_cart, g, ngm )
!
......
......@@ -33,8 +33,7 @@ MODULE types_bz_grid
! ... if on exit ig0 == 0 --> G is not found
!
USE cell_base, ONLY : bg
USE gvecs, ONLY : ngms
USE gvect, ONLY : g
USE gvect, ONLY : g, ngm
USE constants, ONLY : eps8
!
IMPLICIT NONE
......@@ -61,7 +60,7 @@ MODULE types_bz_grid
! gtemp is in cart
!
ig0 = 0
DO ig = 1, ngms
DO ig = 1, ngm
IF ( ALL ( ABS( g(:,ig) - gtemp(:) ) < eps8 ) ) THEN
ig0 = ig
EXIT
......
......@@ -32,7 +32,6 @@ SUBROUTINE do_sxx ( )
USE fft_at_k, ONLY : single_invfft_k,single_fwfft_k
USE distribution_center, ONLY : pert
USE control_flags, ONLY : gamma_only
USE gvecs, ONLY : ngms
USE gvect, ONLY : g,nl,gstart,ngm_g,ngm
USE cell_base, ONLY : omega,at,alat
USE noncollin_module, ONLY : noncolin,npol
......
......@@ -69,9 +69,9 @@ SUBROUTINE dump_r ( auxr, fname )
ALLOCATE(auxr_(dffts%nnr))
auxr_ = CMPLX( auxr, 0.d0, KIND = DP)
IF( gamma_only ) THEN
CALL single_fwfft_gamma(dffts,ngm,ngm,auxr_,auxg,'Smooth')
CALL single_fwfft_gamma(dffts,ngm,ngm,auxr_,auxg,'Dense')
ELSE
CALL single_fwfft_k(dffts,ngm,ngm,auxr_,auxg,'Smooth')
CALL single_fwfft_k(dffts,ngm,ngm,auxr_,auxg,'Dense')
ENDIF
CALL write_wfc_spav ( 2005, TRIM(fname)//".spavr", auxg, westpp_r0, westpp_nr, westpp_rmax )
DEALLOCATE(auxg)
......
......@@ -42,7 +42,7 @@ SUBROUTINE westpp_setup
!
CALL set_npwq()
!
CALL pot3D%init('Wave','default')
CALL pot3D%init('Wave',.FALSE.,'default')
CALL pot3D%print_divergence()
!
CALL set_nbndocc()
......
......@@ -23,7 +23,6 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
USE io_files, ONLY : nwordwfc, iunwfc
USE scf, ONLY : rho, rho_core, rhog_core
USE gvect, ONLY : g,nl,gstart,ngm_g,ngm
USE gvecs, ONLY : ngms
USE gvecw, ONLY : gcutw
USE cell_base, ONLY : tpiba2,omega,tpiba,at,alat
USE fft_base, ONLY : dfftp,dffts
......@@ -71,14 +70,14 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
WRITE(stdout,'(5x,a)') '(X)-Sigma'
CALL io_push_bar()
!
ALLOCATE( pertg( ngms ) )
ALLOCATE( pertg( ngm ) )
IF(noncolin) THEN
ALLOCATE( pertr_nc( dffts%nnr, npol ) )
ELSE
ALLOCATE( pertr( dffts%nnr ) )
ENDIF
!
CALL pot3D%init('Smooth','gb')
CALL pot3D%init('Dense',.FALSE.,'gb')
!
! Set to zero
!
......@@ -156,14 +155,14 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
DO ir=1,dffts%nnr
pertr(ir)=psic(ir)*pertr(ir)
ENDDO
CALL single_fwfft_gamma(dffts,ngms,ngms,pertr,pertg,'Smooth')
CALL single_fwfft_gamma(dffts,ngm,ngm,pertr,pertg,'Dense')
ELSEIF(noncolin) THEN
CALL single_invfft_k(dffts,npw,npwx,evc(1 ,iv_glob),pertr_nc(1,1),'Wave',igk_k(1,current_k))
CALL single_invfft_k(dffts,npw,npwx,evc(1+npwx,iv_glob),pertr_nc(1,2),'Wave',igk_k(1,current_k))
DO ir=1,dffts%nnr
pertr_nc(ir,1)=DCONJG(psic_nc(ir,1))*pertr_nc(ir,1)+DCONJG(psic_nc(ir,2))*pertr_nc(ir,2)
ENDDO
CALL single_fwfft_k(dffts,ngms,ngms,pertr_nc(1,1),pertg,'Smooth') ! no igk
CALL single_fwfft_k(dffts,ngm,ngm,pertr_nc(1,1),pertg,'Dense') ! no igk
! ELSE
! CALL single_invfft_k(dffts,npw,npwx,evc(1,iv_glob),pertr,'Wave',igk_k(1,current_k))
! DO ir=1,dffts%nnr
......@@ -172,10 +171,10 @@ SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
! CALL single_fwfft_k(dffts,ngms,ngms,pertr,pertg,'Smooth') ! no igk
ENDIF
!
DO ig = 1,ngms
DO ig = 1,ngm
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
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - 2._DP * DDOT( 2*ngm, 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 ) - pot3D%div
!
......@@ -215,7 +214,6 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
USE io_files, ONLY : nwordwfc, iunwfc
USE scf, ONLY : rho, rho_core, rhog_core
USE gvect, ONLY : g,nl,gstart,ngm_g,ngm
USE gvecs, ONLY : ngms
USE gvecw, ONLY : gcutw
USE cell_base, ONLY : tpiba2,omega,tpiba,at,alat
USE fft_base, ONLY : dfftp,dffts
......@@ -270,7 +268,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
WRITE(stdout,'(5x,a)') '(X)-Sigma'
CALL io_push_bar()
!
ALLOCATE( pertg( ngms ) )
ALLOCATE( pertg( ngm ) )
ALLOCATE( phase(dffts%nnr) )
ALLOCATE( evckmq(npwx*npol,nbnd) )
IF(noncolin) THEN
......@@ -355,7 +353,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
!
CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
!
CALL pot3D%init('Smooth', 'gb', iq)
CALL pot3D%init('Dense',.FALSE.,'gb',iq)
!
CALL compute_phase( g0, 'cart', phase )
!
......@@ -375,19 +373,19 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
DO ir=1,dffts%nnr
pertr_nc(ir,1)=DCONJG(pertr_nc(ir,1)*phase(ir))*psic_nc(ir,1)+DCONJG(pertr_nc(ir,2)*phase(ir))*psic_nc(ir,2)
ENDDO
CALL single_fwfft_k(dffts,ngms,ngms,pertr_nc(1,1),pertg,'Smooth') ! no igk
CALL single_fwfft_k(dffts,ngm,ngm,pertr_nc(1,1),pertg,'Dense') ! no igk
ELSE
CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1,iv_glob),pertr,'Wave',igk_k(1,ikqs))
DO ir=1,dffts%nnr
pertr(ir)=DCONJG(pertr(ir)*phase(ir)) * psic(ir)
ENDDO
CALL single_fwfft_k(dffts,ngms,ngms,pertr,pertg,'Smooth') ! no igk
CALL single_fwfft_k(dffts,ngm,ngm,pertr,pertg,'Dense') ! no igk
ENDIF
!
DO ig = 1,ngms
DO ig = 1,ngm
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)
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - DDOT( 2*ngm, 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 ) - pot3D%div
!
......
......@@ -98,7 +98,7 @@ SUBROUTINE solve_gfreq_gamma(l_read_restart)
CALL deallocate_bec_type( becp )
CALL allocate_bec_type ( nkb, pert%nloc, becp ) ! I just need 2 becp at a time
!
CALL pot3D%init('Wave','default')
CALL pot3D%init('Wave',.FALSE.,'default')
!
IF(l_read_restart) THEN
CALL solvegfreq_restart_read( bks )
......@@ -501,7 +501,7 @@ SUBROUTINE solve_gfreq_k(l_read_restart)
!
! compute Coulomb potential
!
CALL pot3D%init('Wave', 'default', iq)
CALL pot3D%init('Wave',.TRUE.,'default',iq)
!
! The Hamiltonian is evaluated at k'
!
......
......@@ -163,7 +163,7 @@ SUBROUTINE solve_wfreq_gamma(l_read_restart,l_generate_plot)
CALL start_bar_type ( barra, 'wlanczos', barra_load )
ENDIF
!
CALL pot3D%init('Wave','default')
CALL pot3D%init('Wave',.FALSE.,'default')
!
! LOOP
!
......@@ -768,7 +768,7 @@ SUBROUTINE solve_wfreq_k(l_read_restart,l_generate_plot)
npwq = ngq(iq)
l_gammaq = q_grid%l_pIsGamma(iq)
!
CALL pot3D%init('Wave','default',iq)
CALL pot3D%init('Wave',.TRUE.,'default',iq)
!
DO iks = 1, k_grid%nps ! KPOINT-SPIN
!
......
......@@ -35,7 +35,7 @@ SUBROUTINE wfreq_readin()
!
! READ INPUT_WEST
!
CALL fetch_input((/1/),.TRUE.)
CALL fetch_input(1,(/1/),.TRUE.)
!
! read the input file produced by the pwscf program
! allocate memory and recalculate what is needed
......@@ -51,7 +51,7 @@ SUBROUTINE wfreq_readin()
!
! READ other sections of the input file
!
CALL fetch_input((/2,3/),.TRUE.)
CALL fetch_input(2,(/2,3/),.TRUE.)
!
CALL stop_clock('wfreq_readin')
!
......
......@@ -147,7 +147,7 @@ SUBROUTINE davidson_diago_gamma ( )
notcnv = nvec
dav_iter = -2
!
CALL pot3D%init('Wave','default')
CALL pot3D%init('Wave',.FALSE.,'default')
CALL pot3d%print_divergence()
!
! KIND OF CALCULATION
......@@ -583,7 +583,7 @@ SUBROUTINE davidson_diago_k ( )
!
! compute Coulomb potential
!
CALL pot3D%init('Wave','default',iq)
CALL pot3D%init('Wave',.TRUE.,'default',iq)
CALL pot3d%print_divergence()
!
IF ( q_grid%np > 1 ) THEN
......
......@@ -36,7 +36,7 @@ SUBROUTINE wstat_readin()
!
! READ INPUT_WEST
!
CALL fetch_input((/1/),.TRUE.)
CALL fetch_input(1,(/1/),.TRUE.)
!
! read the input file produced by the pwscf program
! allocate memory and recalculate what is needed
......@@ -52,7 +52,7 @@ SUBROUTINE wstat_readin()
!
! READ other sections of the input file
!
CALL fetch_input((/2/),.TRUE.)
CALL fetch_input(1,(/2/),.TRUE.)
!
CALL stop_clock('wstat_readin')
!
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment