2005-03-05 01:46:02 +08:00
|
|
|
!
|
|
|
|
|
! Copyright (C) 2001 PWSCF group
|
|
|
|
|
! This file is distributed under the terms of the
|
|
|
|
|
! GNU General Public License. See the file `License'
|
|
|
|
|
! in the root directory of the present distribution,
|
|
|
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
|
|
|
!
|
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
|
subroutine chi_test (dvscfs, chif, ik, depsi, auxr, auxg)
|
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
|
!
|
|
|
|
|
! This routine is just a debugging tool. Calculates the scalar product
|
|
|
|
|
! between the chi-wavefunction and Pc DH |psi> in two different ways.
|
|
|
|
|
! The results of the two procedures should be the same.
|
|
|
|
|
!
|
|
|
|
|
#include "f_defs.h"
|
|
|
|
|
|
2009-02-02 18:52:58 +08:00
|
|
|
USE kinds, ONLY : DP
|
|
|
|
|
USE wvfct, ONLY : npwx, nbnd
|
|
|
|
|
USE gsmooth, ONLY : nrxxs
|
|
|
|
|
use ramanm, ONLY : lrd2w, iud2w, jab
|
|
|
|
|
USE units_ph, ONLY : iuwfc, lrdwf, iudwf
|
|
|
|
|
USE qpoint, ONLY : npwq, nksq
|
|
|
|
|
USE eqv, ONLY : dpsi, dvpsi
|
|
|
|
|
USE control_ph, ONLY : nbnd_occ
|
2005-03-05 01:46:02 +08:00
|
|
|
USE wavefunctions_module, ONLY: evc
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
integer :: ik
|
2005-08-28 22:09:42 +08:00
|
|
|
complex(DP) :: dvscfs(nrxxs,3), chif(npwx,nbnd,6), depsi(npwx,nbnd,3), &
|
2005-03-05 01:46:02 +08:00
|
|
|
auxr(nrxxs), auxg(npwx)
|
|
|
|
|
|
2009-02-16 17:12:59 +08:00
|
|
|
complex(DP) :: tmp
|
|
|
|
|
complex(DP), EXTERNAL :: ZDOTC
|
2005-08-28 22:09:42 +08:00
|
|
|
complex(DP) , allocatable :: ps1(:,:,:), ps2(:,:,:), &
|
2005-03-05 01:46:02 +08:00
|
|
|
ps3(:,:,:,:), ps4(:,:,:), au2r(:)
|
|
|
|
|
integer :: ip, jp, ib, jb, ipa, ipb, nrec, ir, ig
|
|
|
|
|
|
|
|
|
|
allocate (ps1 (nbnd,3,6) )
|
|
|
|
|
allocate (ps2 (nbnd,3,6) )
|
|
|
|
|
allocate (ps3 (nbnd,3,nbnd,3) )
|
|
|
|
|
allocate (ps4 (nbnd,3,nbnd) )
|
|
|
|
|
allocate (au2r (nrxxs) )
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
!----------------------------------------------------------
|
|
|
|
|
!
|
|
|
|
|
do ip = 1, 3
|
|
|
|
|
nrec = (ip - 1) * nksq + ik
|
|
|
|
|
call davcio (depsi (1, 1, ip), lrdwf, iudwf, nrec, -1)
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
do jp = 1, 6
|
|
|
|
|
nrec = (jp - 1) * nksq + ik
|
|
|
|
|
call davcio (dpsi, lrd2w, iud2w, nrec, -1)
|
|
|
|
|
do ip = 1, 3
|
|
|
|
|
do ib = 1, nbnd_occ (ik)
|
|
|
|
|
ps1 (ib, ip, jp) = &
|
|
|
|
|
-ZDOTC (npwq, depsi (1, ib, ip), 1, &
|
|
|
|
|
dpsi (1, ib), 1)
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
do ip = 1, 3
|
|
|
|
|
do ib = 1, nbnd_occ (ik)
|
|
|
|
|
do jp = 1, 3
|
|
|
|
|
do jb = 1, nbnd_occ (ik)
|
|
|
|
|
ps3(ib, ip, jb, jp) = &
|
|
|
|
|
ZDOTC (npwq, depsi (1, ib, ip), 1, depsi (1, jb, jp), 1)
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
do ib = 1, nbnd_occ (ik)
|
|
|
|
|
call cft_wave ( evc (1, ib), au2r, +1 )
|
|
|
|
|
do ip = 1, 3
|
|
|
|
|
do ir = 1, nrxxs
|
|
|
|
|
auxr (ir) = au2r (ir) * dvscfs (ir, ip)
|
|
|
|
|
end do
|
|
|
|
|
auxg (:) = (0.d0, 0.d0)
|
|
|
|
|
call cft_wave (auxg, auxr, -1 )
|
|
|
|
|
do jb = 1, nbnd_occ (ik)
|
|
|
|
|
ps4 (ib, ip, jb) = &
|
|
|
|
|
ZDOTC (npwq, auxg, 1, evc (1, jb), 1)
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
do ip = 1, 3
|
|
|
|
|
do ib = 1, nbnd_occ (ik)
|
|
|
|
|
do ipa = 1, 3
|
|
|
|
|
do ipb = 1, 3
|
General cleanup of intrinsic functions:
conversion to real => DBLE
(including real part of a complex number)
conversion to complex => CMPLX
complex conjugate => CONJG
imaginary part => AIMAG
All functions are uppercase.
CMPLX is preprocessed by f_defs.h and performs an explicit cast:
#define CMPLX(a,b) cmplx(a,b,kind=DP)
This implies that 1) f_defs.h must be included whenever a CMPLX is present,
2) CMPLX should stay in a single line, 3) DP must be defined.
All occurrences of real, float, dreal, dfloat, dconjg, dimag, dcmplx
removed - please do not reintroduce any of them.
Tested only with ifc7 and g95 - beware unintended side effects
Maybe not the best solution (explicit casts everywhere would be better)
but it can be easily changed with a script if the need arises.
The following code might be used to test for possible trouble:
program test_intrinsic
implicit none
integer, parameter :: dp = selected_real_kind(14,200)
real (kind=dp) :: a = 0.123456789012345_dp
real (kind=dp) :: b = 0.987654321098765_dp
complex (kind=dp) :: c = ( 0.123456789012345_dp, 0.987654321098765_dp)
print *, ' A = ', a
print *, ' DBLE(A)= ', DBLE(a)
print *, ' C = ', c
print *, 'CONJG(C)= ', CONJG(c)
print *, 'DBLE(c),AIMAG(C) = ', DBLE(c), AIMAG(c)
print *, 'CMPLX(A,B,kind=dp)= ', CMPLX( a, b, kind=dp)
end program test_intrinsic
Note that CMPLX and REAL without a cast yield single precision numbers on
ifc7 and g95 !!!
git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2133 c92efa57-630b-4861-b058-cf58834340f0
2005-08-27 01:44:42 +08:00
|
|
|
tmp = CMPLX (0.d0, 0.d0)
|
2005-03-05 01:46:02 +08:00
|
|
|
do jb = 1, nbnd_occ (ik)
|
|
|
|
|
tmp = tmp + &
|
|
|
|
|
ps3 (ib, ip, jb, ipa) * ps4 (jb, ipb, ib)
|
|
|
|
|
enddo
|
|
|
|
|
if (ipa.eq.ipb) tmp = tmp * 2.d0
|
|
|
|
|
ps1 (ib, ip, jab (ipa, ipb)) = &
|
|
|
|
|
ps1 (ib, ip, jab (ipa, ipb)) - tmp
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
do ip = 1, 3
|
|
|
|
|
do ib = 1, nbnd_occ (ik)
|
|
|
|
|
call cft_wave (depsi (1, ib, ip), au2r, +1 )
|
|
|
|
|
do ipa = 1, 3
|
|
|
|
|
do ir = 1, nrxxs
|
|
|
|
|
auxr (ir) = au2r (ir) * dvscfs (ir, ipa)
|
|
|
|
|
enddo
|
|
|
|
|
auxg (:) = (0.d0, 0.d0)
|
|
|
|
|
call cft_wave (auxg, auxr, -1 )
|
|
|
|
|
do ipb = 1, 3
|
|
|
|
|
tmp = ZDOTC (npwq, auxg, 1, depsi (1, ib, ipb), 1)
|
|
|
|
|
if (ipa.eq.ipb) tmp = tmp * 2.d0
|
|
|
|
|
ps1 (ib, ip, jab (ipa, ipb)) = &
|
|
|
|
|
ps1 (ib, ip, jab (ipa, ipb)) + tmp
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
!
|
|
|
|
|
!----------------------------------------------------------
|
|
|
|
|
!
|
|
|
|
|
do ip = 1, 3
|
|
|
|
|
dpsi (:,:) = (0.d0, 0.d0)
|
|
|
|
|
call dvpsi_e (ik, ip)
|
|
|
|
|
|
|
|
|
|
do ib = 1, nbnd_occ (ik)
|
|
|
|
|
auxg (:) = (0.d0, 0.d0)
|
|
|
|
|
call DAXPY (2 * npwq, -1.d0, dvpsi (1,ib), 1, auxg, 1)
|
|
|
|
|
call cft_wave (evc (1, ib), auxr, +1 )
|
|
|
|
|
|
|
|
|
|
do ir = 1, nrxxs
|
|
|
|
|
auxr (ir) = auxr (ir) * dvscfs (ir, ip)
|
|
|
|
|
enddo
|
|
|
|
|
call cft_wave (auxg, auxr, -1 )
|
|
|
|
|
do jp = 1, 6
|
|
|
|
|
ps2 (ib, ip, jp) = &
|
|
|
|
|
ZDOTC (npwq, auxg, 1, chif (1, ib, jp), 1)
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
!
|
|
|
|
|
!----------------------------------------------------------
|
|
|
|
|
! If everything is ok, ps1 should be equal to ps2; if not
|
|
|
|
|
! there is a problem
|
|
|
|
|
!
|
|
|
|
|
do ib = 1, nbnd_occ (ik)
|
|
|
|
|
do jp = 1, 6
|
|
|
|
|
! write(6,9030) ib,jp, (ps1(ib,ip,jp),ip=1,3)
|
|
|
|
|
! write(6,9030) ib,jp, (ps2(ib,ip,jp),ip=1,3)
|
|
|
|
|
! write(6,'(/)')
|
|
|
|
|
write(6,9031) ib,jp, ( &
|
|
|
|
|
ps1 (ib, ip, jp) / ps2 (ib, ip, jp), ip = 1, 3)
|
|
|
|
|
enddo
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
9030 format(' bnd:',i5,' ip:',i5,6e12.6)
|
|
|
|
|
9031 format(' bnd:',i5,' ip:',i5,6f12.6)
|
|
|
|
|
|
|
|
|
|
deallocate (ps1 )
|
|
|
|
|
deallocate (ps2 )
|
|
|
|
|
deallocate (ps3 )
|
|
|
|
|
deallocate (ps4 )
|
|
|
|
|
deallocate (au2r )
|
|
|
|
|
|
|
|
|
|
return
|
|
|
|
|
end subroutine chi_test
|