TABLE OF CONTENTS
ABINIT/meanvalue_g [ Functions ]
NAME
meanvalue_g
FUNCTION
Compute the mean value of one wavefunction, in reciprocal space, for an operator that is real, diagonal in reciprocal space: <wf|op|wf> For the time being, only spin-independent operators are treated.
COPYRIGHT
Copyright (C) 2003-2018 ABINIT group (XG,BA) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
diag(npw)=diagonal operator (real, spin-independent !)
filter= if 1, need to filter on the value of diag, that must be less than huge(0.0d0)*1.d-11
otherwise, should be 0
istwf_k=storage mode of the vectors
npw=number of planewaves of the vector
nspinor=number of spinor components
vect(2,npw*nspinor)=vector
vect1(2,npw*nspinor*use_ndo)=vector1 (=vector in most of the cases)
use_ndo = says if vect=/vect1
OUTPUT
ar=mean value
PARENTS
dfpt_vtowfk,energy,forstrnps,vtowfk
CHILDREN
xmpi_sum
SOURCE
40 #if defined HAVE_CONFIG_H 41 #include "config.h" 42 #endif 43 44 #include "abi_common.h" 45 46 47 subroutine meanvalue_g(ar,diag,filter,istwf_k,mpi_enreg,npw,nspinor,vect,vect1,use_ndo,ar_im) 48 49 use defs_basis 50 use defs_abitypes 51 use m_profiling_abi 52 use m_errors 53 use m_xmpi 54 55 !This section has been created automatically by the script Abilint (TD). 56 !Do not modify the following lines by hand. 57 #undef ABI_FUNC 58 #define ABI_FUNC 'meanvalue_g' 59 !End of the abilint section 60 61 implicit none 62 63 !Arguments ------------------------------------ 64 !scalars 65 integer,intent(in) :: filter,istwf_k,npw,nspinor,use_ndo 66 real(dp),intent(out) :: ar 67 real(dp),intent(out),optional :: ar_im 68 type(MPI_type),intent(in) :: mpi_enreg 69 !arrays 70 real(dp),intent(in) :: diag(npw),vect(2,npw*nspinor) 71 real(dp),intent(in) :: vect1(2,npw*nspinor) 72 73 !Local variables------------------------------- 74 !scalars 75 integer :: i1,ierr,ipw,jpw,me_g0 76 character(len=500) :: message 77 !arrays 78 79 ! ************************************************************************* 80 me_g0 = mpi_enreg%me_g0 81 82 DBG_CHECK(ANY(filter==(/0,1/)),"Wrong filter") 83 DBG_CHECK(ANY(nspinor==(/1,2/)),"Wrong nspinor") 84 DBG_CHECK(ANY(istwf_k==(/(ipw,ipw=1,9)/)),"Wrong istwf_k") 85 86 if(nspinor==2 .and. istwf_k/=1)then 87 write(message,'(a,a,a,i6,a,i6)')& 88 & ' When istwf_k/=1, nspinor must be 1,',ch10,& 89 & ' however, nspinor=',nspinor,', and istwf_k=',istwf_k 90 MSG_BUG(message) 91 end if 92 93 if(use_ndo==1 .and. (istwf_k==2 .and.me_g0==1)) then 94 write(message,'(a,a)') ch10,' use_ndo==1, not tested' 95 MSG_BUG(message) 96 end if 97 98 ar=zero 99 if(present(ar_im)) ar_im=zero 100 101 !Normal storage mode 102 if(istwf_k==1)then 103 104 ! No filter 105 if(filter==0)then 106 !$OMP PARALLEL DO REDUCTION(+:ar) 107 do ipw=1,npw 108 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 109 end do 110 if(nspinor==2)then 111 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw) 112 do ipw=1+npw,2*npw 113 jpw=ipw-npw 114 ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 115 end do 116 end if 117 if(use_ndo==1 .and. nspinor==2)then 118 !$OMP PARALLEL DO REDUCTION(+:ar_im) 119 do ipw=1,npw 120 ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 121 end do 122 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw) 123 do ipw=1+npw,2*npw 124 jpw=ipw-npw 125 ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 126 end do 127 end if 128 129 ! !$OMP PARALLEL DO REDUCTION(+:ar,ar_im) 130 ! do ipw=1,npw 131 ! ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 132 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 133 ! end do 134 ! if(nspinor==2)then 135 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 136 ! do ipw=1+npw,2*npw 137 ! ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 138 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 139 ! end do 140 ! end if 141 else ! will filter 142 143 !$OMP PARALLEL DO REDUCTION(+:ar) 144 do ipw=1,npw 145 if(diag(ipw)<huge(0.0d0)*1.d-11)then 146 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 147 end if 148 end do 149 if(nspinor==2)then 150 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw) 151 do ipw=1+npw,2*npw 152 jpw=ipw-npw 153 if(diag(jpw)<huge(0.0d0)*1.d-11)then 154 ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 155 end if 156 end do 157 end if 158 if(use_ndo==1 .and. nspinor==2)then 159 !$OMP PARALLEL DO REDUCTION(+:ar_im) 160 do ipw=1,npw 161 if(diag(ipw)<huge(0.0d0)*1.d-11)then 162 ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 163 end if 164 end do 165 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw) 166 do ipw=1+npw,2*npw 167 jpw=ipw-npw 168 if(diag(jpw)<huge(0.0d0)*1.d-11)then 169 ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 170 end if 171 end do 172 end if 173 174 175 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 176 ! do ipw=1,npw 177 ! if(diag(ipw)<huge(0.0d0)*1.d-11)then 178 ! ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 179 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 180 ! end if 181 ! end do 182 ! if(nspinor==2)then 183 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 184 ! do ipw=1+npw,2*npw 185 ! if(diag(ipw-npw)<huge(0.0d0)*1.d-11)then 186 ! ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 187 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 188 ! end if 189 ! end do 190 ! end if ! nspinor==2 191 192 end if ! filter==0 193 194 else if(istwf_k>=2)then 195 196 if(filter==0)then 197 i1=1 198 if(istwf_k==2 .and. me_g0==1)then ! MPIWF need to know which proc has G=0 199 ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2 200 end if 201 202 !$OMP PARALLEL DO REDUCTION(+:ar) 203 do ipw=i1,npw 204 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 205 end do 206 207 else ! filter/=0 208 i1=1 209 if(istwf_k==2 .and. me_g0==1)then 210 if(diag(1)<huge(0.0d0)*1.d-11)then 211 ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2 212 end if 213 end if 214 215 !$OMP PARALLEL DO REDUCTION(+:ar) 216 do ipw=i1,npw 217 if(diag(ipw)<huge(0.0d0)*1.d-11)then 218 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 219 end if 220 end do 221 end if ! filter==0 222 223 ar=two*ar 224 225 end if ! istwf_k 226 227 !MPIWF need to make reduction on ar and ai . 228 if(mpi_enreg%paral_kgb==1)then 229 call xmpi_sum(ar,mpi_enreg%comm_bandspinorfft ,ierr) 230 if(present(ar_im))then 231 call xmpi_sum(ar_im,mpi_enreg%comm_bandspinorfft,ierr) 232 end if 233 end if 234 235 end subroutine meanvalue_g