Commit 4ea2d3b1 authored by Matteo Gerosa's avatar Matteo Gerosa
Browse files

Revised implementation of the findp method in class bz_grid. Search of G0 is...

Revised implementation of the findp method in class bz_grid. Search of G0 is now parallel. Propagated changes to all code.
parent 5dbc5432
......@@ -373,7 +373,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
!
! Workspace
!
INTEGER :: ipert, ig, ir, ibnd, iks, ikqs, ik
INTEGER :: ipert, ig, ir, ibnd, iks, ikqs, ik, is
! Counter on perturbations
! Counter on plane waves
! Counter on real space grids
......@@ -428,13 +428,15 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
!
DO iks = 1, k_grid%nps ! KPOINT-SPIN LOOP
!
ik = k_grid%ip(iks)
ik = k_grid%ip(iks)
is = k_grid%is(iks)
!
current_k = iks
!
CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, "cart" )
ikqs = k_grid%find( kmq, 'cart' )
!ikqs = kmq_grid%index_kq(iks,iq)
CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
!CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
!ikqs = k_grid%find( kmq, 'cart' )
!!ikqs = kmq_grid%index_kq(iks,iq)
!
! computes the phase needed to bring the wavefunction at k-q
! to the equivalent [k-q] point in the first BZ
......
......@@ -43,7 +43,7 @@ MODULE class_bz_grid
!
PROCEDURE :: init => k_or_q_grid_init
PROCEDURE :: find => findp
PROCEDURE :: add => addp
! PROCEDURE :: add => addp
! PROCEDURE :: init_kq => kq_grid_init
! PROCEDURE :: init_q => q_grid_init
! PROCEDURE :: get_phase => get_phase
......@@ -100,6 +100,7 @@ MODULE class_bz_grid
!
! set weights
!
ALLOCATE ( this%weight (this%nps) )
DO ips = 1, this%nps
this%weight(ips) = wk(ips)
ENDDO
......@@ -165,91 +166,123 @@ MODULE class_bz_grid
END SUBROUTINE
!
!
FUNCTION findp(this,p,unit_type) RESULT(ip)
!
!FUNCTION findp(this,p,unit_type) RESULT(ip)
SUBROUTINE findp( this, p, is, unit_type, ip, g0 )
!
! ... ip is the index of p (unit_type = [ "cryst", "cart"])
! ... if on exit ip == 0 --> p is not commensurate with this grid
! ... g0 relates p to an equivalent vector inside the 1BZ
!
USE constants, ONLY : eps8
USE cell_base, ONLY : at, bg
!
IMPLICIT NONE
!
! I/O
!
CLASS(bz_grid),INTENT(IN) :: this
REAL(DP),INTENT(IN) :: p(3)
CHARACTER(LEN=*),INTENT(IN) :: unit_type
INTEGER :: ip
CLASS(bz_grid), INTENT(IN) :: this
REAL(DP), INTENT(IN) :: p(3)
INTEGER, INTENT(IN) :: is
CHARACTER(LEN=*), INTENT(IN) :: unit_type
INTEGER, INTENT(OUT) :: ip
REAL(DP), INTENT(OUT) :: g0(3)
!
! Workspace
!
INTEGER :: i
!
ip = 0
SELECT CASE( unit_type )
CASE("cryst")
DO i = 1, this%np
IF( ( ALL( ABS ( p(:) - this%p_cryst(:,i) ) .LT. eps8 ) ) ) THEN
ip = i
EXIT
ENDIF
ENDDO
CASE("cart")
DO i = 1, this%np
IF( ( ALL( ABS ( p(:) - this%p_cart(:,i) ) .LT. eps8 ) ) ) THEN
ip = i
EXIT
ENDIF
ENDDO
CASE DEFAULT
CALL errore( 1, "unit_type not supported, supported only cryst or cart" )
END SELECT
!
END FUNCTION
!
!
SUBROUTINE addp( this, pin1, pin2, pout, g0, unit_type )
!
! ... out : pout and g0
! ... pout = pin1 + pin2 - g0 ( g0 makes sure that pout is in 1BZ )
! ... unit_type determines the units of pin1, pin2 and pout, g0
!
USE cell_base, ONLY : at, bg
!
IMPLICIT NONE
!
! I/O
!
CLASS(bz_grid), INTENT(IN) :: this
REAL(DP), INTENT(IN) :: pin1(3), pin2(3)
REAL(DP), INTENT(OUT) :: pout(3)
REAL(DP), INTENT(OUT) :: g0(3)
CHARACTER(LEN=*),INTENT(IN) :: unit_type
!
! Workspace
!
REAL(DP) :: ptemp(3)
REAL(DP) :: dp(3)
!
SELECT CASE(unit_type)
CASE("cryst","cart")
CASE DEFAULT
CALL errore( 1, "unit_type not supported, supported only cryst or cart" )
CALL errore( "types_bz_grid", "unit_type not supported, supported only cryst or cart", 1 )
END SELECT
!
ptemp = pin1 + pin2
IF( unit_type == "cart" ) CALL cryst_to_cart( 1, ptemp, at, -1 )
!
! ptemp is now in cryst
! The search must be performed in crystalline coordinates
!
g0 = NINT( ptemp ) ! in cryst
pout = ptemp - g0 ! in cryst
IF ( unit_type == "cart" ) CALL cryst_to_cart( 1, p, at, -1 )
!
IF( unit_type == "cart" ) THEN
CALL cryst_to_cart( 1, pout , bg, 1 )
CALL cryst_to_cart( 1, g0 , bg, 1 )
ENDIF
ip = 0
DO i = 1, this%np
dp(:) = p(:) - this%p_cryst(:,i) - NINT( p(:) - this%p_cryst(:,i) )
IF ( ALL ( ABS ( dp ) .LT. eps8 ) ) THEN
ip = i + (is-1) * this%np
g0(:) = p(:) - this%p_cryst(:,ip)
EXIT
ENDIF
ENDDO
!
! Tranform g0 back to cartesian coordinates if needed
!
IF ( unit_type == "cart" ) CALL cryst_to_cart( 1, g0, bg, 1 )
!
!ip = 0
!SELECT CASE( unit_type )
!CASE("cryst")
! DO i = 1, this%np
! IF( ( ALL( ABS ( p(:) - this%p_cryst(:,i) ) .LT. eps8 ) ) ) THEN
! ip = i
! EXIT
! ENDIF
! ENDDO
!CASE("cart")
! DO i = 1, this%np
! IF( ( ALL( ABS ( p(:) - this%p_cart(:,i) ) .LT. eps8 ) ) ) THEN
! ip = i
! EXIT
! ENDIF
! ENDDO
!CASE DEFAULT
! CALL errore( "class_bz_grid", "unit_type not supported, supported only cryst or cart", 1 )
!END SELECT
!
!END FUNCTION
END SUBROUTINE
!
!
!SUBROUTINE addp( this, pin1, pin2, pout, g0, unit_type )
! !
! ! ... out : pout and g0
! ! ... pout = pin1 + pin2 - g0 ( g0 makes sure that pout is in 1BZ )
! ! ... unit_type determines the units of pin1, pin2 and pout, g0
! !
! USE cell_base, ONLY : at, bg
! USE constants, ONLY : eps8
! !
! IMPLICIT NONE
! !
! ! I/O
! !
! CLASS(bz_grid), INTENT(IN) :: this
! REAL(DP), INTENT(IN) :: pin1(3), pin2(3)
! REAL(DP), INTENT(OUT) :: pout(3)
! REAL(DP), INTENT(OUT) :: g0(3)
! CHARACTER(LEN=*),INTENT(IN) :: unit_type
! !
! ! Workspace
! !
! REAL(DP) :: ptemp(3)
! !
! SELECT CASE(unit_type)
! CASE("cryst","cart")
! CASE DEFAULT
! CALL errore( "types_bz_grid", "unit_type not supported, supported only cryst or cart", 1 )
! END SELECT
! !
! ptemp = pin1 + pin2
! IF( unit_type == "cart" ) CALL cryst_to_cart( 1, ptemp, at, -1 )
! !
! ! ptemp is now in cryst
! !
! g0 = NINT( ptemp ) ! in cryst
! pout = ptemp - g0 ! in cryst
! !
! IF( unit_type == "cart" ) THEN
! CALL cryst_to_cart( 1, pout , bg, 1 )
! CALL cryst_to_cart( 1, g0 , bg, 1 )
! ENDIF
! !
!END SUBROUTINE
! !
! IF ( sig == +1 ) THEN
! csig = '+'
......
......@@ -29,7 +29,9 @@ MODULE types_bz_grid
CONTAINS
!
!
FUNCTION findG(g0, unit_type) RESULT(ig0)
!
!
FUNCTION findG(g0, unit_type) RESULT(ig0)
!
! ... ig0 is the index of G (unit_type = [ "cryst", "cart"])
! ... if on exit ig0 == 0 --> G is not found
......@@ -54,10 +56,10 @@ MODULE types_bz_grid
SELECT CASE(unit_type)
CASE("cryst","cart")
CASE DEFAULT
CALL errore( 1, "unit_type not supported, supported only cryst or cart" )
CALL errore( "types_bz_grid", "unit_type not supported, supported only cryst or cart", 1 )
END SELECT
!
gtemp = g0
gtemp = g0
IF( unit_type == "cryst" ) CALL cryst_to_cart( 1, gtemp, bg, 1)
!
! gtemp is in cart
......@@ -66,6 +68,7 @@ MODULE types_bz_grid
DO ig = 1, ngms
IF ( ALL ( ABS( g(:,ig) - gtemp(:) ) < eps8 ) ) THEN
ig0 = ig
EXIT
ENDIF
ENDDO
!
......@@ -79,6 +82,8 @@ MODULE types_bz_grid
USE gvecs, ONLY : nls
USE fft_base, ONLY : dffts
USE fft_interfaces, ONLY : invfft
USE mp, ONLY : mp_max
USE mp_bands, ONLY : intra_bgrp_comm
!
IMPLICIT NONE
!
......@@ -90,22 +95,27 @@ MODULE types_bz_grid
!
! Workspace
!
INTEGER :: ipol, ig, ig0
INTEGER :: ipol, ig, ig0, glob_ig0
!
SELECT CASE(unit_type)
CASE("cryst","cart")
CASE DEFAULT
CALL errore( 1, "unit_type not supported, supported only cryst or cart" )
CALL errore( "types_bz_grid", "unit_type not supported, supported only cryst or cart", 1 )
END SELECT
!
ig0 = findG(g0,unit_type)
IF( ig0 == 0 ) CALL errore( 1, "G0 not found" )
glob_ig0 = ig0
CALL mp_max( glob_ig0, intra_bgrp_comm )
!
phase = (0._DP, 0._DP)
phase( nls(ig0) ) = (1._DP, 0._DP)
IF( glob_ig0 == 0 ) CALL errore( "types_bz_grid", "G0 not found", 1 )
!
! phase(r) = exp(-iG_0*r)
!
phase = (0._DP, 0._DP)
!
IF ( ig0 /= 0 ) THEN
phase( nls(ig0) ) = (1._DP, 0._DP)
ENDIF
CALL invfft( 'Wave', phase, dffts )
phase(1:dffts%nnr) = DCONJG( phase(1:dffts%nnr) )
!
......
......@@ -246,7 +246,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
!
! Workspace
!
INTEGER :: ik,iks,iq,ikqs,ib,ifreq,glob_ifreq,il,im,glob_im,ip
INTEGER :: ik,is,iks,iq,ikqs,ib,ifreq,glob_ifreq,il,im,glob_im,ip
INTEGER :: nbndval
!
REAL(DP),EXTERNAL :: integrate_imfreq
......@@ -286,7 +286,9 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
! LOOP
!
DO iks = 1, k_grid%nps ! KPOINT-SPIN (MATRIX ELEMENT)
!
ik = k_grid%ip(iks)
is = k_grid%is(iks)
!
DO ib = qp_bandrange(1), qp_bandrange(2)
!
......@@ -295,8 +297,9 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
!
DO iq = 1, q_grid%np ! Q-POINT
!
CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
ikqs = k_grid%find( kmq, 'cart' )
CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
!CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
!ikqs = k_grid%find( kmq, 'cart' )
!ikqs = kmq_grid%index_kq(iks,iq)
l_gammaq = q_grid%l_pIsGamma(iq)
nbndval = nbnd_occ(ikqs)
......@@ -369,6 +372,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
DO iks = 1, k_grid%nps
!
ik = k_grid%ip(iks)
is = k_grid%is(iks)
!
DO ib = qp_bandrange(1), qp_bandrange(2)
!
......@@ -377,11 +381,12 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
residues_b = 0._DP
residues_h = 0._DP
!
DO iq = 1, q_grid%np
DO iq = 1, q_grid%np ! Q-POINT
!
CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
ikqs = k_grid%find( kmq, 'cart' )
!ikqs = kmq_grid%index_kq(iks,iq)
CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
!CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
!ikqs = k_grid%find( kmq, 'cart' )
!!ikqs = kmq_grid%index_kq(iks,iq)
l_gammaq = q_grid%l_pIsGamma(iq)
nbndval = nbnd_occ(ikqs)
!
......
......@@ -266,7 +266,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
COMPLEX(DP), ALLOCATABLE :: evckmq(:,:), phase(:)
REAL(DP), EXTERNAL :: DDOT
COMPLEX(DP), EXTERNAL :: ZDOTC
INTEGER :: ib,iv,i1,ir,iks,ik,ig,iv_glob,iq,ikqs
INTEGER :: ib,iv,i1,ir,iks,ik,is,ig,iv_glob,iq,ikqs
INTEGER :: nbndval
INTEGER :: npwkq
TYPE(idistribute) :: vband
......@@ -309,6 +309,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
DO iks = 1, k_grid%nps ! KPOINT-SPIN
!
ik = k_grid%ip(iks)
is = k_grid%is(iks)
!
! ... Set k-point, spin, kinetic energy, needed by Hpsi
!
......@@ -373,9 +374,10 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
!
l_gammaq = q_grid%l_pIsGamma(iq)
!
CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
ikqs = k_grid%find( kmq, 'cart' )
!ikqs = kmq_grid%index_kq(iks,iq)
CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
!CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
!ikqs = k_grid%find( kmq, 'cart' )
!!ikqs = kmq_grid%index_kq(iks,iq)
!
IF ( l_gammaq ) THEN
CALL store_sqvc(mysqvc,ngms,div_kind_hf,mydiv)
......
......@@ -463,6 +463,7 @@ SUBROUTINE solve_gfreq_k(l_read_restart)
!
DO ikks = 1, k_grid%nps ! KPOINT-SPIN (MATRIX ELEMENT)
IF(ikks<bks%lastdone_ks) CYCLE
!
ikk = k_grid%ip(ikks)
!
npwk = ngk(ikks)
......@@ -493,13 +494,15 @@ SUBROUTINE solve_gfreq_k(l_read_restart)
!
DO iks = 1, k_grid%nps ! KPOINT-SPIN (INTEGRAL OVER K')
IF (ikks==bks%lastdone_ks .AND. ib==bks%lastdone_band .AND. iks <= bks%lastdone_ki) CYCLE
!
ik = k_grid%ip(iks)
!
time_spent(1) = get_clock( 'glanczos' )
!
CALL k_grid%add( k_grid%p_cart(:,ikk), -k_grid%p_cart(:,ik), q, g0, 'cart' )
iq = q_grid%find( q, 'cart' )
!iq = q_grid_aux%index_q(ikks,iks)
CALL q_grid%find( k_grid%p_cart(:,ikk) - k_grid%p_cart(:,ik), 1, 'cart', iq, g0 )
!CALL k_grid%add( k_grid%p_cart(:,ikk), -k_grid%p_cart(:,ik), q, g0, 'cart' )
!iq = q_grid%find( q, 'cart' )
!!iq = q_grid_aux%index_q(ikks,iks)
!
CALL preallocate_solvegfreq_q( iks_l2g(ikks), iks_l2g(iks), qp_bandrange(1), qp_bandrange(2), pert)
!
......
......@@ -692,6 +692,7 @@ SUBROUTINE solve_qp_k(l_secant,l_generate_plot)
! WARNING: iks and ikks are switched w.r.t. solve_gfreq_k
!
DO iks = 1, k_grid%nps ! KPOINT-SPIN (MATRIX ELEMENT)
!
ik = k_grid%ip(iks)
!
nbndval = nbnd_occ(iks)
......@@ -699,11 +700,13 @@ SUBROUTINE solve_qp_k(l_secant,l_generate_plot)
DO ib = qp_bandrange(1), qp_bandrange(2)
!
DO ikks = 1, k_grid%nps ! KPOINT-SPIN (INTEGRAL OVER K')
!
ikk = k_grid%ip(ikks)
!
CALL k_grid%add( k_grid%p_cart(:,ik), -k_grid%p_cart(:,ikk), q, g0, 'cart' )
iq = q_grid%find( q, 'cart' )
!iq = q_grid_aux%index_q(iks,ikks)
CALL q_grid%find( k_grid%p_cart(:,ik) - k_grid%p_cart(:,ikk), 1, 'cart', iq, g0 )
!CALL k_grid%add( k_grid%p_cart(:,ik), -k_grid%p_cart(:,ikk), q, g0, 'cart' )
!iq = q_grid%find( q, 'cart' )
!!iq = q_grid_aux%index_q(iks,ikks)
!
IF(ALLOCATED(overlap)) DEALLOCATE(overlap)
ALLOCATE(overlap(pert%nglob, nbnd ) )
......
......@@ -648,7 +648,7 @@ SUBROUTINE solve_wfreq_k(l_read_restart,l_generate_plot)
!
! Workspace
!
INTEGER :: i1,i2,i3,im,ip,ig,glob_ip,ir,iv,iks,ik,iq,ikqs,ipol,m
INTEGER :: i1,i2,i3,im,ip,ig,glob_ip,ir,iv,iks,ik,is,iq,ikqs,ipol,m
CHARACTER(LEN=512) :: fname
CHARACTER(LEN=6) :: my_label_b
CHARACTER(LEN=5) :: my_label_q
......@@ -771,7 +771,10 @@ SUBROUTINE solve_wfreq_k(l_read_restart,l_generate_plot)
ENDIF
!
DO iks = 1, k_grid%nps ! KPOINT-SPIN
!
ik = k_grid%ip(iks)
is = k_grid%is(iks)
!
IF(iq==bks%lastdone_q .AND. iks<bks%lastdone_ks) CYCLE
!
! ... Set k-point, spin, kinetic energy, needed by Hpsi
......@@ -816,8 +819,9 @@ SUBROUTINE solve_wfreq_k(l_read_restart,l_generate_plot)
!ikqs = kpq_grid%index_kq(iks,iq)
!npwkq = ngk(ikqs)
!
CALL k_grid%add( k_grid%p_cart(:,ik), q_grid%p_cart(:,iq), kpq, g0, 'cart' )
ikqs = k_grid%find( kpq, 'cart' )
CALL k_grid%find( k_grid%p_cart(:,ik) + q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
!CALL k_grid%add( k_grid%p_cart(:,ik), q_grid%p_cart(:,iq), kpq, g0, 'cart' )
!ikqs = k_grid%find( kpq, 'cart' )
!
npwkq = ngk(ikqs)
!
......
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