mirror of https://gitlab.com/QEF/q-e.git
Initial release (previously developed against tag QE-3-1-1 in directory
NMR_new). git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3662 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
seitsonen
parent
d9a096994a
commit
f926db4ce1
|
|
@ -0,0 +1,107 @@
|
|||
# Makefile for NMR_new
|
||||
|
||||
include ../make.sys
|
||||
|
||||
NMROBJS = \
|
||||
nmr_module.o \
|
||||
magn_main.o \
|
||||
stop_code.o \
|
||||
apply_p.o \
|
||||
apply_vel.o \
|
||||
greenfunction.o \
|
||||
h_psiq.o \
|
||||
ch_psi_all.o \
|
||||
sym_cart_tensor.o \
|
||||
symmetrize_field.o \
|
||||
test_sum_rule.o \
|
||||
compute_u_kq.o \
|
||||
suscept_crystal.o \
|
||||
j_para.o \
|
||||
biot_savart.o \
|
||||
compute_sigma_bare.o \
|
||||
init_us_2_no_phase.o \
|
||||
init_paw_2_no_phase.o \
|
||||
g_tensor_crystal.o \
|
||||
write_tensor_field.o
|
||||
|
||||
MODULES = \
|
||||
../Modules/atom.o \
|
||||
../Modules/autopilot.o \
|
||||
../Modules/basic_algebra_routines.o \
|
||||
../Modules/bfgs_module.o \
|
||||
../Modules/berry_phase.o \
|
||||
../Modules/cell_base.o \
|
||||
../Modules/check_stop.o \
|
||||
../Modules/clocks.o \
|
||||
../Modules/constants.o \
|
||||
../Modules/constraints_module.o \
|
||||
../Modules/control_flags.o \
|
||||
../Modules/electrons_base.o \
|
||||
../Modules/fft_base.o \
|
||||
../Modules/fft_scalar.o \
|
||||
../Modules/fft_types.o \
|
||||
../Modules/functionals.o \
|
||||
../Modules/input_parameters.o \
|
||||
../Modules/io_files.o \
|
||||
../Modules/io_global.o \
|
||||
../Modules/ions_base.o \
|
||||
../Modules/ions_nose.o \
|
||||
../Modules/kind.o \
|
||||
../Modules/metadyn_base.o \
|
||||
../Modules/metadyn_io.o \
|
||||
../Modules/metadyn_vars.o \
|
||||
../Modules/mp_global.o \
|
||||
../Modules/mp_wave.o \
|
||||
../Modules/mp.o \
|
||||
../Modules/parallel_include.o \
|
||||
../Modules/parameters.o \
|
||||
../Modules/parser.o \
|
||||
../Modules/path_base.o \
|
||||
../Modules/path_formats.o \
|
||||
../Modules/path_variables.o \
|
||||
../Modules/path_opt_routines.o \
|
||||
../Modules/path_io_routines.o \
|
||||
../Modules/path_reparametrisation.o \
|
||||
../Modules/printout_base.o \
|
||||
../Modules/pseudo_types.o \
|
||||
../Modules/ptoolkit.o \
|
||||
../Modules/random_numbers.o \
|
||||
../Modules/read_cards.o \
|
||||
../Modules/read_namelists.o \
|
||||
../Modules/read_uspp.o \
|
||||
../Modules/read_upf.o \
|
||||
../Modules/recvec.o \
|
||||
../Modules/splinelib.o \
|
||||
../Modules/stick_base.o \
|
||||
../Modules/shmem_include.o \
|
||||
../Modules/task_groups.o \
|
||||
../Modules/timestep.o \
|
||||
../Modules/uspp.o \
|
||||
../Modules/version.o \
|
||||
../Modules/wavefunctions.o \
|
||||
../Modules/wave_base.o \
|
||||
../Modules/xml_io_base.o
|
||||
|
||||
PWOBJS = ../PW/libpw.a
|
||||
|
||||
PHOBJS = ../PH/libph.a
|
||||
|
||||
TLDEPS=bindir ph pw mods libs libiotk
|
||||
|
||||
all : tldeps magn.x
|
||||
|
||||
magn.x : $(NMROBJS) $(PWOBJS) $(PHOBJS) $(LIBOBJS)
|
||||
$(MPIF90) $(LDFLAGS) -o $@ \
|
||||
$(NMROBJS) $(MODULES) $(PHOBJS) $(PWOBJS) $(LIBOBJS) $(LIBS)
|
||||
- ( cd ../bin; ln -fs ../NMR_new/$@ . )
|
||||
|
||||
tldeps:
|
||||
test -n "$(TLDEPS)" && ( cd .. ; $(MAKE) $(MFLAGS) $(TLDEPS) || exit 1) || :
|
||||
|
||||
clean :
|
||||
- /bin/rm -f magn.x *.o *~ *.F90 *.d *.mod *.i work.pc
|
||||
|
||||
include make.depend
|
||||
|
||||
# DO NOT DELETE
|
||||
|
||||
|
|
@ -0,0 +1,45 @@
|
|||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 apply_p(psi, p_psi, ik, ipol, q)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Apply the kinetic part of the velocity operator
|
||||
! ... |p_psi> = (G+k+q/2)_{ipol} |psi>
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE klist, ONLY : xk
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- parameters ---------------------------------------------------------
|
||||
IMPLICIT none
|
||||
INTEGER, INTENT(IN) :: ik ! k-point
|
||||
INTEGER, INTENT(IN) :: ipol ! cartesian direction (1..3)
|
||||
REAL(DP), INTENT(IN) :: q(3)
|
||||
COMPLEX(DP), INTENT(IN) :: psi(npwx,nbnd)
|
||||
COMPLEX(DP), INTENT(OUT) :: p_psi(npwx,nbnd)
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
REAL(DP) :: gk
|
||||
INTEGER :: ig, ibnd
|
||||
|
||||
call start_clock('apply_p')
|
||||
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
do ig = 1, npw
|
||||
gk = xk(ipol,ik) + g(ipol,igk(ig)) + q(ipol)
|
||||
p_psi(ig,ibnd) = gk * tpiba * psi(ig,ibnd)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call stop_clock('apply_p')
|
||||
|
||||
END SUBROUTINE apply_p
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,293 @@
|
|||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE apply_vel(psi, vel_psi, ik, ipol, q)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Apply the velocity operator
|
||||
! ... v = p + dV^{NL}_{k+q,k}/dk
|
||||
! ...
|
||||
! ... Here we use Hartree atomic units, so that:
|
||||
! ... V^{NL} => V^{NL} * ryd_to_hartree
|
||||
!-----------------------------------------------------------------------
|
||||
USE kinds, ONLY : DP
|
||||
USE ions_base, ONLY : nat, ityp, ntyp => nsp
|
||||
USE klist, ONLY : xk
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk
|
||||
USE becmod, ONLY : becp
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- paramters ----------------------------------------------------------
|
||||
IMPLICIT NONE
|
||||
INTEGER, INTENT(IN) :: ipol ! cartesian direction (1..3)
|
||||
INTEGER, INTENT(IN) :: ik ! k-point
|
||||
COMPLEX(DP), INTENT(IN) :: psi(npwx,nbnd)
|
||||
COMPLEX(DP), INTENT(OUT) :: vel_psi(npwx,nbnd)
|
||||
REAL(DP), INTENT(IN) :: q(3)
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
real(dp), parameter :: ryd_to_hartree = 0.5d0
|
||||
complex(dp), allocatable :: aux(:,:), vkb_save(:,:)
|
||||
real(dp) :: dk, dxk(3)
|
||||
integer :: isign
|
||||
logical :: q_is_zero
|
||||
|
||||
|
||||
! first apply p
|
||||
call apply_p(psi, vel_psi, ik, ipol, q)
|
||||
|
||||
call start_clock('apply_vel')
|
||||
|
||||
! set dk
|
||||
dk = q_nmr/2.d0
|
||||
|
||||
! if no projectors, return
|
||||
if (nkb == 0) return
|
||||
|
||||
! check if |q| is zero
|
||||
q_is_zero = .false.
|
||||
if (sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3)) < 1d-8) q_is_zero = .true.
|
||||
|
||||
! allocate temporary arrays, save old NL-potential
|
||||
allocate(aux(npwx,nbnd), vkb_save(npwx,nkb))
|
||||
vkb_save = vkb
|
||||
|
||||
#if 1
|
||||
!====================================================================
|
||||
! compute (1/2|dk|) ( V^{NL}_{k+dk+q,k+dk} |psi> -
|
||||
! V^{NL}_{k-dk+q,k-dk} |psi> )
|
||||
!====================================================================
|
||||
do isign = -1,1,2
|
||||
dxk(:) = xk(:,ik)
|
||||
dxk(ipol) = dxk(ipol) + isign * dk ! k + dk
|
||||
|
||||
! compute <\beta(k \pm dk)| and project on |psi>
|
||||
call init_us_2_no_phase(npw, igk, dxk, vkb)
|
||||
call ccalbec (nkb, npwx, npw, nbnd_occ(ik), becp, vkb, psi)
|
||||
|
||||
! |q|!=0 => compute |\beta(k \pm dk + q)>
|
||||
if (.not. q_is_zero) then
|
||||
dxk(:) = dxk(:) + q(:)
|
||||
call init_us_2_no_phase(npw, igk, dxk, vkb)
|
||||
endif
|
||||
|
||||
! apply |\beta(k \pm dk+q)>D<\beta(k \pm dk)| to |psi>
|
||||
aux = (0.d0,0.d0)
|
||||
call add_vuspsi(npwx, npw, nbnd_occ(ik), psi, aux)
|
||||
vel_psi = vel_psi + dble(isign) * ryd_to_hartree * aux/(2.d0*dk*tpiba)
|
||||
enddo
|
||||
#else
|
||||
|
||||
do isign = -1,1,2
|
||||
! compute <\beta(k)| and project on |psi>
|
||||
call init_us_2_no_phase(npw, igk, xk(1,ik), vkb)
|
||||
call ccalbec (nkb, npwx, npw, nbnd_occ(ik), becp, vkb, psi)
|
||||
|
||||
dxk(ipol) = xk(ipol,ik) + isign * dk ! k + dk
|
||||
call init_us_2_no_phase(npw, igk, dxk, vkb)
|
||||
|
||||
! apply |\beta(k \pm dk+q)>D<\beta(k)| to |psi>
|
||||
aux = (0.d0,0.d0)
|
||||
call add_vuspsi(npwx, npw, nbnd_occ(ik), psi, aux)
|
||||
vel_psi = vel_psi + 0.5*dble(isign) * ryd_to_hartree * aux/(2.d0*dk*tpiba)
|
||||
|
||||
|
||||
! compute <\beta(k \pm dk)| and project on |psi>
|
||||
dxk(:) = xk(:,ik)
|
||||
dxk(ipol) = dxk(ipol) + isign * dk ! k + dk
|
||||
call init_us_2_no_phase(npw, igk, xk(1,ik), vkb)
|
||||
call ccalbec (nkb, npwx, npw, nbnd_occ(ik), becp, vkb, psi)
|
||||
|
||||
dxk(:) = xk(:,ik)
|
||||
call init_us_2_no_phase(npw, igk, dxk, vkb)
|
||||
|
||||
! apply |\beta(k+q)>D<\beta(k \pm dk)| to |psi>
|
||||
aux = (0.d0,0.d0)
|
||||
call add_vuspsi(npwx, npw, nbnd_occ(ik), psi, aux)
|
||||
vel_psi = vel_psi + 0.5*dble(isign) * ryd_to_hartree * aux/(2.d0*dk*tpiba)
|
||||
enddo
|
||||
#endif
|
||||
|
||||
! restore NL-potential at k
|
||||
vkb = vkb_save
|
||||
|
||||
! free memory
|
||||
deallocate(aux, vkb_save)
|
||||
|
||||
call stop_clock('apply_vel')
|
||||
|
||||
END SUBROUTINE apply_vel
|
||||
|
||||
|
||||
|
||||
|
||||
! Here are some previous and testing version of the subroutine
|
||||
! it looks like that the analytic derivative has some troubles.
|
||||
!
|
||||
! The numerical one (that computes just |dbeta/dk>) is just fine
|
||||
! but is is more involved that the routine above.
|
||||
|
||||
#if 0
|
||||
real(DP), allocatable :: gk (:,:)
|
||||
complex(DP), allocatable :: ps2(:,:,:), dvkb (:,:), dvkb1 (:,:), &
|
||||
work (:,:), becp2(:,:), becp1(:,:), aux(:,:), vkb_save(:,:)
|
||||
complex(DP), external :: ZDOTC
|
||||
integer :: ig, na, ibnd, ikb, jkb, nt, ih, jh, ijkb0
|
||||
integer :: isign
|
||||
real(DP) :: axis(3)
|
||||
|
||||
|
||||
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
print*, ibnd
|
||||
print '(2(F12.6,2X))', vel_psi(1:5,ibnd)
|
||||
print*
|
||||
enddo
|
||||
|
||||
!!!vl_psi = 0.d0
|
||||
|
||||
! this is the old routine, applying |beta><dbeta/dk| + |dbeta/dk><beta|
|
||||
! compute (k+G)/|k+G|
|
||||
allocate (gk(3,npwx))
|
||||
do ig = 1, npw
|
||||
gk (1, ig) = (xk (1, ik) + g (1, igk (ig) ) ) * tpiba
|
||||
gk (2, ig) = (xk (2, ik) + g (2, igk (ig) ) ) * tpiba
|
||||
gk (3, ig) = (xk (3, ik) + g (3, igk (ig) ) ) * tpiba
|
||||
g2kin (ig) = gk (1, ig) **2 + gk (2, ig) **2 + gk (3, ig) **2
|
||||
|
||||
if (g2kin (ig) < 1.0d-10) then
|
||||
gk (1, ig) = 0.d0
|
||||
gk (2, ig) = 0.d0
|
||||
gk (3, ig) = 0.d0
|
||||
else
|
||||
gk (1, ig) = gk (1, ig) / sqrt (g2kin (ig) )
|
||||
gk (2, ig) = gk (2, ig) / sqrt (g2kin (ig) )
|
||||
gk (3, ig) = gk (3, ig) / sqrt (g2kin (ig) )
|
||||
endif
|
||||
enddo
|
||||
|
||||
! allocate variable for the derivative of the beta's
|
||||
allocate(dvkb(npwx,nkb), dvkb1(npwx,nkb))
|
||||
dvkb = (0.d0,0.d0)
|
||||
dvkb1 = (0.d0,0.d0)
|
||||
work = (0.d0,0.d0)
|
||||
|
||||
! and this is the contribution from nonlocal pseudopotentials
|
||||
axis = 0.d0; axis(ipol) = 1.d0
|
||||
!!!call gen_us_dy (ik, at (1, ipol), dvkb1)
|
||||
call gen_us_dy (ik, axis(1), dvkb1)
|
||||
call gen_us_dj (ik, dvkb)
|
||||
print*, dvkb1(2,:)
|
||||
|
||||
jkb = 0
|
||||
do nt = 1, ntyp
|
||||
do na = 1, nat
|
||||
if (nt == ityp (na)) then
|
||||
do ikb = 1, nh (nt)
|
||||
jkb = jkb + 1
|
||||
do ig = 1, npw
|
||||
!!!work (ig,jkb) = dvkb1 (ig, jkb) + dvkb (ig, jkb) * &
|
||||
!!! (at (1, ipol) * gk (1, ig) + &
|
||||
!!! at (2, ipol) * gk (2, ig) + &
|
||||
!!! at (3, ipol) * gk (3, ig) )
|
||||
work (ig,jkb) = dvkb1 (ig, jkb) + dvkb (ig, jkb) * gk(ipol,ig)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
work(1,3) = (0.d0,-0.133684d0)
|
||||
jkb = 3
|
||||
write(98,'(2F12.6,2X))') work(:,jkb)
|
||||
print*
|
||||
|
||||
! this was the previous numerical derivative
|
||||
!allocate(dvkb(npwx,nkb), dvkb1(npwx,nkb))
|
||||
dvkb = (0.d0,0.d0)
|
||||
dvkb1 = (0.d0,0.d0)
|
||||
work = (0.d0,0.d0)
|
||||
vel_psi = 0.d0
|
||||
|
||||
! compute it by finite differences
|
||||
work = (0.d0,0.d0)
|
||||
dxk(:) = xk(:,ik); dxk(ipol) = dxk(ipol) + dk
|
||||
call init_us_2(npw,igk,dxk,work)
|
||||
|
||||
dxk(:) = xk(:,ik); dxk(ipol) = dxk(ipol) - dk
|
||||
call init_us_2(npw,igk,dxk,dvkb)
|
||||
|
||||
work(:,:) = work(:,:) - dvkb(:,:)
|
||||
work = 0.5d0 * work / (dk * tpiba)
|
||||
write(99,'(2F12.6,2X))') work(:,jkb)
|
||||
print*
|
||||
|
||||
|
||||
allocate( becp2(nkb,nbnd), becp1(nkb,nbnd))
|
||||
becp2(:,:) = (0.d0,0.d0)
|
||||
becp1(:,:) = (0.d0,0.d0)
|
||||
call ccalbec (nkb, npwx, npw, nbnd, becp2, work, psi)
|
||||
call init_us_2(npw,igk,xk(1,ik),vkb)
|
||||
call ccalbec (nkb, npwx, npw, nbnd, becp1, vkb, psi)
|
||||
|
||||
ijkb0 = 0
|
||||
allocate (ps2 ( nkb, nbnd, 2))
|
||||
ps2(:,:,:)=(0.d0,0.d0)
|
||||
do nt = 1, ntyp
|
||||
do na = 1, nat
|
||||
if (nt == ityp (na)) then
|
||||
do ih = 1, nh (nt)
|
||||
ikb = ijkb0 + ih
|
||||
do jh = 1, nh (nt)
|
||||
jkb = ijkb0 + jh
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
!! CHECK factor (-1.d0,0.d0)
|
||||
!!!ps2(ikb,ibnd,1) = ps2(ikb,ibnd,1)+ becp2(jkb,ibnd)* &
|
||||
!!! (-1.d0,0.d0) * (deeq(ih,jh,na,current_spin) &
|
||||
!!! -et(ibnd,ik)*qq(ih,jh,nt))
|
||||
!!!ps2(ikb,ibnd,2) = ps2(ikb,ibnd,2) +becp1(jkb,ibnd,ik) * &
|
||||
!!! (-1.d0,0.d0) * (deeq(ih,jh,na,current_spin)&
|
||||
!!! -et(ibnd,ik)*qq(ih,jh,nt))
|
||||
ps2(ikb,ibnd,1) = ps2(ikb,ibnd,1) + becp2(jkb,ibnd)* &
|
||||
(1.d0,0.d0) * deeq(ih,jh,na,current_spin)
|
||||
ps2(ikb,ibnd,2) = ps2(ikb,ibnd,2) + becp1(jkb,ibnd) * &
|
||||
(1.d0,0.d0) * deeq(ih,jh,na,current_spin)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
ijkb0=ijkb0+nh(nt)
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
if (ikb /= nkb .OR. jkb /= nkb) call errore ('apply_vel', 'unexpected error',1)
|
||||
!print*, 'ps2(:,:,1)='
|
||||
!print*, ps2(:,:,1)
|
||||
!print*, 'ps2(:,:,2)='
|
||||
!print*, ps2(:,:,2)
|
||||
|
||||
CALL ZGEMM( 'N', 'N', npw, nbnd_occ(ik), nkb, &
|
||||
(1.d0,0.d0), vkb(1,1), npwx, ps2(1,1,1), nkb, (1.d0,0.d0), &
|
||||
vel_psi(1,1), npwx )
|
||||
|
||||
CALL ZGEMM( 'N', 'N', npw, nbnd_occ(ik), nkb, &
|
||||
(1.d0,0.d0), work(1,1), npwx, ps2(1,1,2), nkb, (1.d0,0.d0), &
|
||||
vel_psi(1,1), npwx )
|
||||
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
print*, ibnd
|
||||
print '(2(F12.6,2X))', vel_psi(1:5,ibnd)
|
||||
print*
|
||||
enddo
|
||||
STOP
|
||||
|
||||
deallocate (ps2)
|
||||
deallocate (gk)
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,95 @@
|
|||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 biot_savart(jpol)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Compute the induced magentic field via the Biot-Savart law
|
||||
! ... B_ind(r) = (1/c) \int d^3r' j(r') (r-r')/|r-r'|
|
||||
! ... which in reciprocal space reads:
|
||||
! ... B_ind(G) = (4\pi/c) (i G \times j(G))/G^2
|
||||
! ... the G=0 is not computed here and is given by chi_bare
|
||||
USE kinds, ONLY : DP
|
||||
USE constants, ONLY : fpi
|
||||
USE klist, ONLY : xk
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk
|
||||
USE gvect, ONLY : ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, &
|
||||
nrx3, nrxx, nl, nlm, g, gg, ecutwfc, gcutm
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- parameters ---------------------------------------------------------
|
||||
IMPLICIT none
|
||||
INTEGER, INTENT(IN) :: jpol
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
COMPLEX(DP), allocatable :: aux(:), j_of_g(:,:)
|
||||
REAL(DP) :: gk
|
||||
complex(dp) :: fact
|
||||
INTEGER :: ig, ipol, ispin
|
||||
|
||||
call start_clock('biot_savart')
|
||||
|
||||
! allocate memory
|
||||
allocate(aux(nrxxs), j_of_g(1:ngm,3))
|
||||
|
||||
! transform current to reciprocal space
|
||||
j_of_g(:,:) = 0.d0
|
||||
do ispin = 1, nspin
|
||||
do ipol = 1, 3
|
||||
aux(1:nrxxs) = j_bare(1:nrxxs,ipol,jpol,ispin)
|
||||
call cft3s(aux, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -1)
|
||||
j_of_g(1:ngm,ipol) = j_of_g(1:ngm,ipol) + aux(nl(1:ngm))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! compute induced field in reciprocal space
|
||||
do ig = gstart, ngm
|
||||
fact = (0.d0,1.d0) * (fpi/c) / (gg(ig) * tpiba)
|
||||
b_ind(ig,1,jpol) = fact * (g(2,ig)*j_of_g(ig,3) - g(3,ig)*j_of_g(ig,2))
|
||||
b_ind(ig,2,jpol) = fact * (g(3,ig)*j_of_g(ig,1) - g(1,ig)*j_of_g(ig,3))
|
||||
b_ind(ig,3,jpol) = fact * (g(1,ig)*j_of_g(ig,2) - g(2,ig)*j_of_g(ig,1))
|
||||
enddo
|
||||
|
||||
! transform induced field in real space
|
||||
do ipol = 1, 3
|
||||
aux = (0.d0,0.d0)
|
||||
aux(nl(1:ngm)) = b_ind(1:ngm,ipol,jpol)
|
||||
call cft3s(aux, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 1)
|
||||
b_ind_r(1:nrxxs,ipol,jpol) = real(aux(1:nrxxs))
|
||||
enddo
|
||||
|
||||
deallocate(aux, j_of_g)
|
||||
call stop_clock('biot_savart')
|
||||
|
||||
END SUBROUTINE biot_savart
|
||||
|
||||
|
||||
|
||||
SUBROUTINE field_to_reciprocal_space
|
||||
USE kinds, ONLY : DP
|
||||
USE gvect, ONLY : ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, &
|
||||
nrx3, nrxx, nl, nlm, g, gg, ecutwfc, gcutm
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
IMPLICIT NONE
|
||||
complex(dp), allocatable :: aux(:)
|
||||
integer :: ipol, jpol
|
||||
|
||||
allocate(aux(nrxxs))
|
||||
b_ind(:,:,:) = 0.d0
|
||||
do ipol = 1, 3
|
||||
do jpol = 1, 3
|
||||
aux(1:nrxxs) = b_ind_r(1:nrxxs,ipol,jpol)
|
||||
call cft3s(aux, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -1)
|
||||
b_ind(1:ngm,ipol,jpol) = aux(nl(1:ngm))
|
||||
enddo
|
||||
enddo
|
||||
deallocate(aux)
|
||||
|
||||
END SUBROUTINE field_to_reciprocal_space
|
||||
|
|
@ -0,0 +1,107 @@
|
|||
!
|
||||
! 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 ch_psi_all (n, h, ah, e, ik, m)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! This routine applies the operator ( H - \epsilon S + alpha_pv P_v)
|
||||
! to a vector h. The result is given in Ah.
|
||||
!
|
||||
#include "f_defs.h"
|
||||
|
||||
use pwcom
|
||||
use becmod
|
||||
USE kinds, only : DP
|
||||
USE nmr_module
|
||||
implicit none
|
||||
|
||||
integer :: n, m, ik
|
||||
! input: the dimension of h
|
||||
! input: the number of bands
|
||||
! input: the k point
|
||||
|
||||
real(DP) :: e (m)
|
||||
! input: the eigenvalue
|
||||
|
||||
complex(DP) :: h (npwx, m), ah (npwx, m)
|
||||
! input: the vector
|
||||
! output: the operator applied to the vector
|
||||
!
|
||||
! local variables
|
||||
!
|
||||
integer :: ibnd, ikq, ig
|
||||
! counter on bands
|
||||
! the point k+q
|
||||
! counter on G vetors
|
||||
|
||||
complex(DP), allocatable :: ps (:,:), hpsi (:,:), spsi (:,:)
|
||||
! scalar products
|
||||
! the product of the Hamiltonian and h
|
||||
! the product of the S matrix and h
|
||||
|
||||
call start_clock ('ch_psi')
|
||||
allocate (ps ( nbnd , m))
|
||||
allocate (hpsi( npwx , m))
|
||||
allocate (spsi( npwx , m))
|
||||
hpsi (:,:) = (0.d0, 0.d0)
|
||||
spsi (:,:) = (0.d0, 0.d0)
|
||||
!
|
||||
! compute the product of the hamiltonian with the h vector
|
||||
!
|
||||
call h_psiq (npwx, n, m, h, hpsi, spsi)
|
||||
|
||||
call start_clock ('last')
|
||||
!
|
||||
! then we compute the operator H-epsilon S
|
||||
!
|
||||
do ibnd = 1, m
|
||||
do ig = 1, n
|
||||
ah (ig, ibnd) = hpsi (ig, ibnd) - e (ibnd) * spsi (ig, ibnd)
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! Here we compute the projector in the valence band
|
||||
!
|
||||
!!!if (lgamma) then
|
||||
ikq = ik
|
||||
!!!else
|
||||
!!! ikq = 2 * ik
|
||||
!!!endif
|
||||
ps (:,:) = (0.d0, 0.d0)
|
||||
|
||||
call ZGEMM ('C', 'N', nbnd_occ (ikq) , m, n, (1.d0, 0.d0) , evq, &
|
||||
npwx, spsi, npwx, (0.d0, 0.d0) , ps, nbnd)
|
||||
ps (:,:) = ps(:,:) * alpha_pv
|
||||
#ifdef __PARA
|
||||
call reduce (2 * nbnd * m, ps)
|
||||
#endif
|
||||
|
||||
hpsi (:,:) = (0.d0, 0.d0)
|
||||
call ZGEMM ('N', 'N', n, m, nbnd_occ (ikq) , (1.d0, 0.d0) , evq, &
|
||||
npwx, ps, nbnd, (1.d0, 0.d0) , hpsi, npwx)
|
||||
spsi(:,:) = hpsi(:,:)
|
||||
!
|
||||
! And apply S again
|
||||
!
|
||||
call ccalbec (nkb, npwx, n, m, becp, vkb, hpsi)
|
||||
|
||||
call s_psi (npwx, n, m, hpsi, spsi)
|
||||
do ibnd = 1, m
|
||||
do ig = 1, n
|
||||
ah (ig, ibnd) = ah (ig, ibnd) + spsi (ig, ibnd)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
deallocate (spsi)
|
||||
deallocate (hpsi)
|
||||
deallocate (ps)
|
||||
call stop_clock ('last')
|
||||
call stop_clock ('ch_psi')
|
||||
return
|
||||
end subroutine ch_psi_all
|
||||
|
|
@ -0,0 +1,178 @@
|
|||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 compute_sigma_bare(chi_bare, sigma_bare)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Compute the bare contribution to the chemical shift at the
|
||||
! ... position of the nuclei, given the induced field
|
||||
USE kinds, ONLY : DP
|
||||
USE klist, ONLY : xk
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk
|
||||
USE gvect, ONLY : ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, &
|
||||
nrx3, nrxx, nl, nlm, g, gg, ecutwfc, gcutm
|
||||
USE ions_base, ONLY : nat, tau, atm, ityp
|
||||
USE io_global, ONLY : stdout
|
||||
USE symme, ONLY : s, nsym, irt
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
! Arguments
|
||||
REAL(DP), INTENT(IN) :: chi_bare(3,3)
|
||||
real(dp), intent(out) :: sigma_bare(3,3,nat)
|
||||
|
||||
! Local
|
||||
integer :: na, ig
|
||||
real(dp) :: macroscopic_shape(3,3)
|
||||
real(dp) :: arg, tr_sigma
|
||||
complex(dp) :: tmp_sigma(3,3)
|
||||
|
||||
macroscopic_shape(:,:) = 2.0_dp/3.0_dp
|
||||
! like in paratec:
|
||||
macroscopic_shape(:,:) = 0.0_dp
|
||||
do na = 1, 3
|
||||
macroscopic_shape(na,na) = 2.0_dp/3.0_dp
|
||||
enddo
|
||||
|
||||
write(stdout,'(5X,''NMR chemical bare shifts in ppm:'')')
|
||||
write(stdout,*)
|
||||
|
||||
do na = 1, nat
|
||||
tmp_sigma(:,:) = 0.0_dp
|
||||
|
||||
do ig = gstart, ngm
|
||||
arg = (g(1,ig)*tau(1,na) + g(2,ig)*tau(2,na) + g(3,ig)*tau(3,na)) * tpi
|
||||
tmp_sigma(:,:) = tmp_sigma(:,:) &
|
||||
+ b_ind(ig,:,:) * cmplx(cos(arg),sin(arg))
|
||||
enddo
|
||||
|
||||
#if 1
|
||||
! this is the G = 0 term
|
||||
if (gstart == 2) then
|
||||
tmp_sigma(:,:) = tmp_sigma(:,:) &
|
||||
- (4.0_dp*pi) * macroscopic_shape(:,:) * chi_bare(:,:) !*TMPTMPTMP
|
||||
end if
|
||||
#endif
|
||||
|
||||
sigma_bare(:,:,na) = real(tmp_sigma(:,:))
|
||||
enddo
|
||||
#ifdef __PARA
|
||||
call reduce(9*na, sigma_bare)
|
||||
#endif
|
||||
|
||||
#if 0
|
||||
! symmetrize tensors ??
|
||||
do na = 1, nat
|
||||
call trntns (sigma_bare(1,1,na), at, bg, -1)
|
||||
enddo
|
||||
call symz(sigma_bare, nsym, s, nat, irt)
|
||||
do na = 1, nat
|
||||
call trntns (sigma_bare(1,1,na), at, bg, 1)
|
||||
enddo
|
||||
#endif
|
||||
|
||||
do na = 1, nat
|
||||
tr_sigma = (sigma_bare(1,1,na)+sigma_bare(2,2,na)+sigma_bare(3,3,na))/3.0_dp
|
||||
write(stdout,'(5X,''Atom'',I3,2X,A3,'' pos: ('',3(F10.6),&
|
||||
'') sigma: '',F14.4)') na, atm(ityp(na)), tau(:,na), tr_sigma*1e6_dp
|
||||
write(stdout, tens_fmt) sigma_bare(:,:,na) * 1e6_dp
|
||||
enddo
|
||||
|
||||
end subroutine compute_sigma_bare
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE compute_sigma_diamagnetic( sigma_diamagnetic )
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Compute the diamagnetic contribution to the chemical shift at the
|
||||
! ... position of the nuclei, given the induced field
|
||||
USE kinds, ONLY : DP
|
||||
USE klist, ONLY : xk
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk
|
||||
USE gvect, ONLY : ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, &
|
||||
nrx3, nrxx, nl, nlm, g, gg, ecutwfc, gcutm
|
||||
USE ions_base, ONLY : nat, tau, atm, ityp
|
||||
USE io_global, ONLY : stdout
|
||||
USE symme, ONLY : s, nsym, irt
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
! Arguments
|
||||
real(dp), intent(inout) :: sigma_diamagnetic(3,3,nat)
|
||||
|
||||
! Local
|
||||
integer :: na
|
||||
real(dp) :: tr_sigma
|
||||
|
||||
|
||||
write(stdout,'(5X,''NMR chemical diamagnetic shifts in ppm:'')')
|
||||
write(stdout,*)
|
||||
|
||||
! symmetrize tensors
|
||||
do na = 1, nat
|
||||
call trntns (sigma_diamagnetic(1,1,na), at, bg, -1)
|
||||
enddo
|
||||
call symz(sigma_diamagnetic, nsym, s, nat, irt)
|
||||
do na = 1, nat
|
||||
call trntns (sigma_diamagnetic(1,1,na), at, bg, 1)
|
||||
enddo
|
||||
|
||||
do na = 1, nat
|
||||
tr_sigma = (sigma_diamagnetic(1,1,na)+sigma_diamagnetic(2,2,na) &
|
||||
+sigma_diamagnetic(3,3,na))/3.0_dp
|
||||
write(stdout,'(5X,''Atom'',I3,2X,A3,'' pos: ('',3(F10.6),&
|
||||
'') sigma: '',F14.4)') na, atm(ityp(na)), tau(:,na), tr_sigma*1e6_dp
|
||||
write(stdout, tens_fmt) sigma_diamagnetic(:,:,na) * 1e6_dp
|
||||
enddo
|
||||
|
||||
end subroutine compute_sigma_diamagnetic
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE compute_sigma_paramagnetic( sigma_paramagnetic )
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Compute the paramagnetic contribution to the chemical shift at the
|
||||
! ... position of the nuclei, given the induced field
|
||||
USE kinds, ONLY : DP
|
||||
USE klist, ONLY : xk
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk
|
||||
USE gvect, ONLY : ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, &
|
||||
nrx3, nrxx, nl, nlm, g, gg, ecutwfc, gcutm
|
||||
USE ions_base, ONLY : nat, tau, atm, ityp
|
||||
USE io_global, ONLY : stdout
|
||||
USE symme, ONLY : s, nsym, irt
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
! Arguments
|
||||
real(dp), intent(inout) :: sigma_paramagnetic(3,3,nat)
|
||||
|
||||
! Local
|
||||
integer :: na
|
||||
real(dp) :: tr_sigma
|
||||
|
||||
write(stdout,'(5X,''NMR chemical paramagnetic shifts in ppm:'')')
|
||||
write(stdout,*)
|
||||
|
||||
! symmetrize tensors
|
||||
do na = 1, nat
|
||||
call trntns (sigma_paramagnetic(1,1,na), at, bg, -1)
|
||||
enddo
|
||||
call symz(sigma_paramagnetic, nsym, s, nat, irt)
|
||||
do na = 1, nat
|
||||
call trntns (sigma_paramagnetic(1,1,na), at, bg, 1)
|
||||
enddo
|
||||
|
||||
do na = 1, nat
|
||||
tr_sigma = (sigma_paramagnetic(1,1,na)+sigma_paramagnetic(2,2,na) &
|
||||
+sigma_paramagnetic(3,3,na))/3.0_dp
|
||||
write(stdout,'(5X,''Atom'',I3,2X,A3,'' pos: ('',3(F10.6),&
|
||||
'') sigma: '',F14.4)') na, atm(ityp(na)), tau(:,na), tr_sigma*1e6_dp
|
||||
write(stdout, tens_fmt) sigma_paramagnetic(:,:,na) * 1e6_dp
|
||||
enddo
|
||||
|
||||
end subroutine compute_sigma_paramagnetic
|
||||
|
|
@ -0,0 +1,186 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE compute_u_kq(ik, q)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE io_files, ONLY : nwordwfc, iunwfc
|
||||
USE cell_base, ONLY : at, bg, omega, tpiba, tpiba2
|
||||
USE wavefunctions_module, ONLY : evc
|
||||
USE klist, ONLY : nkstot, xk
|
||||
USE wvfct, ONLY : nbnd, nbndx, npwx, npw, igk, g2kin, current_k, btype
|
||||
USE nmr_module
|
||||
USE gvect, only: g,ngm,ecutwfc,ngl,nrxx, nr1, nr2, nr3, nrx1, nrx2, nrx3
|
||||
USE uspp, ONLY : nkb, vkb, okvan
|
||||
USE pwcom, ONLY : et
|
||||
USE becmod, ONLY : becp
|
||||
USE lsda_mod, ONLY : current_spin, lsda, isk
|
||||
USE g_psi_mod, ONLY : h_diag, s_diag
|
||||
USE scf, ONLY : vrs, vltot, vr
|
||||
USE complex_diis_module
|
||||
|
||||
IMPLICIT NONE
|
||||
INTEGER, INTENT(IN) :: ik
|
||||
REAL(DP), INTENT(IN) :: q(3)
|
||||
|
||||
real(dp), allocatable :: etq(:)
|
||||
real(dp) :: xk_plus_q(3)
|
||||
integer :: ig, ntry, iter, notconv
|
||||
real(dp) :: v_of_0
|
||||
real(dp) :: avg_iter, dav_iter, cg_iter, ethr
|
||||
integer :: diis_iter
|
||||
logical :: lrot
|
||||
|
||||
integer, parameter :: max_cg_iter = 200
|
||||
integer, parameter :: isolve = 0
|
||||
integer, parameter :: david = 4
|
||||
complex(dp) :: ZDOTC
|
||||
|
||||
CALL start_clock( 'u_kq' )
|
||||
v_of_0 = SUM( vltot(1:nrxx) ) / DBLE( nr1 * nr2 * nr3 )
|
||||
CALL reduce( 1, v_of_0 )
|
||||
|
||||
call gk_sort(xk(1,ik),ngm,g,ecutwfc/tpiba2,npw,igk,g2kin)
|
||||
current_k = ik
|
||||
if (lsda) current_spin = isk(ik)
|
||||
|
||||
xk_plus_q(:) = xk(:,ik) + q(:)
|
||||
call init_us_2(npw,igk,xk_plus_q,vkb)
|
||||
|
||||
! ... read in wavefunctions from the SCF calculation at k
|
||||
! ... and use it as a starting guess
|
||||
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
|
||||
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
|
||||
!g2kin(1:npw) = ( ( xk(1,ik) + g(1,igk(1:npw)) )**2 + &
|
||||
! ( xk(2,ik) + g(2,igk(1:npw)) )**2 + &
|
||||
! ( xk(3,ik) + g(3,igk(1:npw)) )**2 ) * tpiba2
|
||||
!call h_psi(npwx, npw, nbnd, evc, evq)
|
||||
!print*, ZDOTC(npw, evc(1,1), 1, evq(1,1), 1)*13.605692d0
|
||||
!print*, ZDOTC(npw, evc(1,2), 1, evq(1,2), 1)*13.605692d0
|
||||
! initial guess
|
||||
evq = evc
|
||||
|
||||
! ... sets the kinetic energy
|
||||
g2kin(1:npw) = ( ( xk_plus_q(1) + g(1,igk(1:npw)) )**2 + &
|
||||
( xk_plus_q(2) + g(2,igk(1:npw)) )**2 + &
|
||||
( xk_plus_q(3) + g(3,igk(1:npw)) )**2 ) * tpiba2
|
||||
|
||||
ALLOCATE( h_diag( npwx ) )
|
||||
ALLOCATE( s_diag( npwx ) )
|
||||
ALLOCATE(btype( nbnd, nkstot ) )
|
||||
|
||||
nbndx = nbnd
|
||||
IF ( isolve == 0 ) nbndx = david * nbnd
|
||||
btype(:,:) = 1
|
||||
iter = 1
|
||||
allocate(etq(1:nbndx))
|
||||
ethr = conv_threshold
|
||||
|
||||
!====================================================================
|
||||
! Diagonalize
|
||||
!====================================================================
|
||||
IF ( isolve == 1 ) THEN
|
||||
! ... Conjugate-Gradient diagonalization
|
||||
! ... h_diag is the precondition matrix
|
||||
!!!call errore('compute_u_kq','CG not working properly?',-1)
|
||||
h_diag = 1.D0
|
||||
FORALL( ig = 1 : npw )
|
||||
h_diag(ig) = 1.D0 + g2kin(ig) + SQRT( 1.D0 + ( g2kin(ig) - 1.D0 )**2 )
|
||||
END FORALL
|
||||
|
||||
ntry = 0
|
||||
CG_loop : DO
|
||||
IF ( ntry > 0 ) THEN
|
||||
CALL cinitcgg( npwx, npw, nbnd, nbnd, evq, evq, etq )
|
||||
avg_iter = avg_iter + 1.D0
|
||||
END IF
|
||||
|
||||
CALL ccgdiagg( npwx, npw, nbnd, evq, etq, btype(1,ik), &
|
||||
h_diag, ethr, max_cg_iter, .true., &
|
||||
notconv, cg_iter )
|
||||
|
||||
avg_iter = avg_iter + cg_iter
|
||||
ntry = ntry + 1
|
||||
print*, etq(:)*13.605692d0
|
||||
|
||||
! ... exit condition
|
||||
IF ( test_exit_cond() ) EXIT CG_loop
|
||||
END DO CG_loop
|
||||
|
||||
|
||||
ELSE IF ( isolve == 2 ) THEN
|
||||
! ... RMM-DIIS method
|
||||
call errore('compute_u_kq','DIIS not working properly?',-1)
|
||||
h_diag(1:npw) = g2kin(1:npw) + v_of_0
|
||||
CALL usnldiag( h_diag, s_diag )
|
||||
|
||||
ntry = 0
|
||||
RMMDIIS_loop: DO
|
||||
CALL cdiisg( npw, npwx, nbnd, evq, etq, &
|
||||
btype(1,ik), notconv, diis_iter, iter )
|
||||
|
||||
avg_iter = avg_iter + diis_iter
|
||||
ntry = ntry + 1
|
||||
|
||||
! ... exit condition
|
||||
IF ( test_exit_cond() ) EXIT RMMDIIS_loop
|
||||
END DO RMMDIIS_loop
|
||||
|
||||
|
||||
ELSE IF (isolve == 0) THEN
|
||||
! ... Davidson diagonalization
|
||||
! ... h_diag are the diagonal matrix elements of the
|
||||
! ... hamiltonian used in g_psi to evaluate the correction
|
||||
! ... to the trial eigenvectors
|
||||
h_diag(1:npw) = g2kin(1:npw) + v_of_0
|
||||
CALL usnldiag( h_diag, s_diag )
|
||||
|
||||
ntry = 0
|
||||
david_loop: DO
|
||||
lrot = ( iter == 1 )
|
||||
CALL cegterg( npw, npwx, nbnd, nbndx, evq, ethr, &
|
||||
okvan, etq, btype(1,ik), notconv, &
|
||||
lrot, dav_iter )
|
||||
|
||||
avg_iter = avg_iter + dav_iter
|
||||
ntry = ntry + 1
|
||||
|
||||
! ... exit condition
|
||||
IF ( test_exit_cond() ) EXIT david_loop
|
||||
END DO david_loop
|
||||
!!!print '(4F8.4)', etq(1:nbnd)*13.605692d0
|
||||
ELSE
|
||||
|
||||
call errore('compute_u_kq', 'diagonalization method?', -1)
|
||||
|
||||
END IF
|
||||
|
||||
IF ( notconv > MAX( 5, nbnd / 4 ) ) THEN
|
||||
CALL errore( 'compute_u_kq', 'too many bands are not converged', 1 )
|
||||
END IF
|
||||
|
||||
call stop_clock('u_kq')
|
||||
|
||||
DEALLOCATE( h_diag )
|
||||
DEALLOCATE( s_diag )
|
||||
DEALLOCATE(btype )
|
||||
|
||||
|
||||
CONTAINS
|
||||
|
||||
FUNCTION test_exit_cond()
|
||||
IMPLICIT NONE
|
||||
LOGICAL :: test_exit_cond
|
||||
test_exit_cond = ( ntry > 100 ) .OR. ( notconv <= 0 )
|
||||
END FUNCTION test_exit_cond
|
||||
|
||||
|
||||
END SUBROUTINE compute_u_kq
|
||||
|
|
@ -0,0 +1,924 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE g_tensor_crystal
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Compute the "bare" susceptibility as in Eq.(64-65) of
|
||||
! ... PRB 63, 245101 (2001)
|
||||
! ... add more comments
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE io_files, ONLY : nwordwfc, iunwfc
|
||||
USE cell_base, ONLY : at, bg, omega, tpiba, tpiba2
|
||||
USE wavefunctions_module, ONLY : evc
|
||||
USE klist, ONLY : nks, nkstot, wk, xk, nelec
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk, wg, g2kin, &
|
||||
current_k
|
||||
USE lsda_mod, ONLY : current_spin, lsda, isk, nspin
|
||||
USE becmod, ONLY : becp
|
||||
USE symme, ONLY : nsym, s, ftau
|
||||
USE scf, ONLY : vr, vltot, rho
|
||||
USE gvect, ONLY : ngm, nr1, nr2, nr3, nrx1, nrx2, &
|
||||
nrx3, nrxx, nlm, g
|
||||
USE mp_global, ONLY : my_pool_id
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!<apsi>
|
||||
USE paw, ONLY: paw_vkb, paw_becp, paw_nkb, aephi, &
|
||||
psphi, paw_nh, paw_nhtol, &
|
||||
paw_nhtom, paw_indv, paw_nbeta
|
||||
USE ions_base, ONLY : nat
|
||||
!</apsi>
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
IMPLICIT NONE
|
||||
complex(dp), allocatable, dimension(:,:,:) :: p_evc, vel_evc, g_vel_evc
|
||||
complex(dp), allocatable :: aux(:,:)
|
||||
|
||||
! Q tensor of eq. (65) (pGv => HH in Paratec, vGv => VV in Paratec)
|
||||
real(dp) :: q_pGv(3,3,-1:1), q_vGv(3,3,-1:1)
|
||||
|
||||
! F tensor of eq. (64)
|
||||
real(dp) :: f_pGv(3,3,-1:1), f_vGv(3,3,-1:1)
|
||||
|
||||
! chi_bare tensor of eq. (64)
|
||||
real(dp) :: chi_bare_pGv(3,3), chi_bare_vGv(3,3)
|
||||
|
||||
! f-sum rule
|
||||
real(dp) :: f_sum(3,3)
|
||||
|
||||
real(dp) :: q(3), braket, delta_rmc, gipaw_delta_rmc, rmc_gipaw
|
||||
integer :: ia, ib, ik, ipol, jpol, i, ibnd, isign, ispin
|
||||
complex(dp), external :: ZDOTC
|
||||
|
||||
integer :: s_maj, s_min
|
||||
real(dp) :: e_hartree, charge, s_weight, rho_diff, d_omega
|
||||
real(dp), allocatable :: grad_vr(:,:), v_local(:,:)
|
||||
real(dp), allocatable :: grad_vh(:,:), vh(:)
|
||||
real(dp), dimension ( 3, 3 ) :: delta_g_rmc, delta_g_bare, delta_g_soo, &
|
||||
delta_g_soo_2, delta_g_paramagn, delta_g_diamagn, delta_g_total, &
|
||||
delta_g_rmc_gipaw
|
||||
real(dp) :: g_e = 2.0023192778_DP, g_prime, units_Ry2Ha = 0.5_DP
|
||||
|
||||
logical :: tevaluate_chi, tcalculate_correct_delta_g_soo
|
||||
|
||||
real(dp) :: diamagnetic_corr_tensor(3,3)
|
||||
real(dp) :: paramagnetic_corr_tensor(3,3)
|
||||
real(dp) :: sigma_paramagnetic(3,3)
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
|
||||
!<apsi> TMPTMPTMP Until the reconstruction has been implemented
|
||||
delta_g_paramagn = 0.0_DP
|
||||
delta_g_diamagn = 0.0_DP
|
||||
!</apsi>
|
||||
|
||||
g_prime = 2 * ( g_e - 1 )
|
||||
|
||||
!<apsi> TMPTMPTMP Move to input or make it default?
|
||||
tcalculate_correct_delta_g_soo = .true.
|
||||
!</apsi>
|
||||
if ( tcalculate_correct_delta_g_soo ) then
|
||||
tevaluate_chi = .true.
|
||||
else
|
||||
tevaluate_chi = .false.
|
||||
end if
|
||||
|
||||
! Select majority and minority spin components
|
||||
rho_diff = SUM ( rho ( :, 1 ) - rho ( :, nspin ) )
|
||||
if ( rho_diff > +1.0e-3 ) then
|
||||
s_maj = 1
|
||||
s_min = nspin
|
||||
else if ( rho_diff < -1.0e-3 ) then
|
||||
s_maj = nspin
|
||||
s_min = 1
|
||||
else
|
||||
write ( stdout, * ) "WARNING: rho_diff zero!"
|
||||
end if
|
||||
|
||||
! allocate memory
|
||||
allocate ( p_evc(npwx,nbnd,3), vel_evc(npwx,nbnd,3), &
|
||||
aux(npwx,nbnd), g_vel_evc(npwx,nbnd,3) )
|
||||
|
||||
! zero the Q tensors
|
||||
q_pGv(:,:,:) = 0.0_dp
|
||||
q_vGv(:,:,:) = 0.0_dp
|
||||
|
||||
! zero the current and the field
|
||||
j_bare(:,:,:,:) = (0.0_dp,0.0_dp)
|
||||
b_ind(:,:,:) = (0.0_dp,0.0_dp)
|
||||
|
||||
delta_rmc = 0.0_DP
|
||||
gipaw_delta_rmc = 0.0_DP
|
||||
|
||||
write(stdout, '(5X,''Computing the magnetic susceptibility'')')
|
||||
|
||||
!====================================================================
|
||||
! loop over k-points
|
||||
!====================================================================
|
||||
do ik = 1, nks
|
||||
#ifdef __PARA
|
||||
if (my_pool_id == 0) &
|
||||
write(*, '(5X,''k-point #'',I5,'' of '',I5,'' (my_pool_id='',I3)') &
|
||||
ik, nks, my_pool_id
|
||||
#else
|
||||
write(*, '(5X,''k-point #'',I5,'' of '',I5)') ik, nks
|
||||
#endif
|
||||
current_k = ik
|
||||
current_spin = isk(ik)
|
||||
|
||||
if ( current_spin == s_maj ) then
|
||||
s_weight = +1
|
||||
else
|
||||
s_weight = -1
|
||||
end if
|
||||
|
||||
! initialize at k-point k
|
||||
call gk_sort(xk(1,ik), ngm, g, ecutwfc/tpiba2, npw, igk, g2kin)
|
||||
g2kin(:) = g2kin(:) * tpiba2
|
||||
call init_us_2(npw,igk,xk(1,ik),vkb)
|
||||
|
||||
! read wfcs from file and compute becp
|
||||
call davcio (evc, nwordwfc, iunwfc, ik, -1)
|
||||
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
|
||||
|
||||
!<apsi>
|
||||
call init_paw_2_no_phase (npw, igk, xk (1, ik), paw_vkb)
|
||||
call ccalbec (paw_nkb, npwx, npw, nbnd, paw_becp, paw_vkb, evc)
|
||||
diamagnetic_corr_tensor = 0.0
|
||||
call diamagnetic_correction ( diamagnetic_corr_tensor )
|
||||
delta_g_diamagn = delta_g_diamagn + s_weight * diamagnetic_corr_tensor
|
||||
!</apsi>
|
||||
|
||||
! this is the case q = 0 (like the case of the f-sum rule)
|
||||
q(:) = 0.0_dp
|
||||
!!!write(*,'(''q='',3(F12.4))') q
|
||||
|
||||
call compute_u_kq(ik, q)
|
||||
evc = evq
|
||||
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
delta_rmc = delta_rmc - s_weight &
|
||||
* wg(ibnd,ik) * SUM ( g2kin(:) * conjg(evc(:,ibnd)) * evc(:,ibnd) )
|
||||
end do
|
||||
|
||||
CALL relativistic_mass_correction ( rmc_gipaw )
|
||||
write(6,*) "CCC: ", rmc_gipaw
|
||||
gipaw_delta_rmc = gipaw_delta_rmc - s_weight * rmc_gipaw
|
||||
|
||||
! compute p_k|evc>, v_k|evc> and G_k v_{k,k}|evc>
|
||||
call apply_operators
|
||||
|
||||
!------------------------------------------------------------------
|
||||
! f-sum rule (pGv term only)
|
||||
!------------------------------------------------------------------
|
||||
do ia = 1, 3
|
||||
do ib = 1, 3
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
braket = 2.0_dp*real(ZDOTC(npw, p_evc(1,ibnd,ia), 1, &
|
||||
g_vel_evc(1,ibnd,ib), 1), DP)
|
||||
f_sum(ia,ib) = f_sum(ia,ib) + wg(ibnd,ik) * braket
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#ifdef __PARA
|
||||
call reduce(9, f_sum)
|
||||
#endif
|
||||
|
||||
!------------------------------------------------------------------
|
||||
! pGv and vGv contribution to chi_{bare}
|
||||
!------------------------------------------------------------------
|
||||
if ( tevaluate_chi ) then
|
||||
do i = 1, 3
|
||||
call add_to_tensor(q_pGv(:,:,0), p_evc, g_vel_evc)
|
||||
call add_to_tensor(q_vGv(:,:,0), vel_evc, g_vel_evc)
|
||||
enddo
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------
|
||||
! loop over -q and +q
|
||||
!------------------------------------------------------------------
|
||||
do isign = -1, 1, 2
|
||||
|
||||
! loop over cartesian directions
|
||||
do i = 1, 3
|
||||
if (iverbosity > 10) then
|
||||
write(stdout,*) " QQQ: ", i * isign
|
||||
end if
|
||||
! set the q vector
|
||||
q(:) = 0.0_dp
|
||||
q(i) = real(isign,dp) * q_nmr
|
||||
!!!write(*,'(''q='',3(F12.4))') q
|
||||
|
||||
! compute the wfcs at k+q
|
||||
call compute_u_kq(ik, q)
|
||||
|
||||
! compute p_k|evc>, v_k|evc> and G_{k+q} v_{k+q,k}|evc>
|
||||
call apply_operators
|
||||
|
||||
!<apsi>
|
||||
call init_paw_2_no_phase (npw, igk, xk(:,ik)+q(:), paw_vkb)
|
||||
call paramagnetic_correction ( paramagnetic_corr_tensor )
|
||||
paramagnetic_corr_tensor = paramagnetic_corr_tensor * s_weight
|
||||
call add_to_sigma_para ( paramagnetic_corr_tensor, delta_g_paramagn )
|
||||
!</apsi>
|
||||
|
||||
if ( tevaluate_chi ) then
|
||||
! pGv and vGv contribution to chi_bare
|
||||
call add_to_tensor(q_pGv(:,:,isign), p_evc, g_vel_evc)
|
||||
call add_to_tensor(q_vGv(:,:,isign), vel_evc, g_vel_evc)
|
||||
end if
|
||||
|
||||
! now the j_bare term
|
||||
call add_to_current(j_bare(:,:,:,current_spin), evc, g_vel_evc)
|
||||
enddo ! i
|
||||
|
||||
enddo ! isign
|
||||
enddo ! ik
|
||||
|
||||
! TODO: put a lot of poolreduce here
|
||||
#ifdef __PARA
|
||||
call poolreduce(9, f_sum)
|
||||
call poolreduce(9, q_pGv)
|
||||
call poolreduce(9, q_vGv)
|
||||
! TODO: non working yet in parallel!!!
|
||||
#endif
|
||||
|
||||
!====================================================================
|
||||
! print out results
|
||||
!====================================================================
|
||||
write(stdout,'(5X,''End of magnetic susceptibility calculation'',/)')
|
||||
|
||||
! f-sum rule
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(5X,''f-sum rule:'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') f_sum
|
||||
endif
|
||||
call sym_cart_tensor(f_sum)
|
||||
write(stdout, '(5X,''f-sum rule (symmetrized):'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') f_sum
|
||||
|
||||
!--------------------------------------------------------------------
|
||||
! now get the current, induced field and delta_g
|
||||
!--------------------------------------------------------------------
|
||||
chi_bare_pGv(:,:) = chi_bare_pGv(:,:) / omega
|
||||
j_bare(:,:,:,:) = j_bare(:,:,:,:) / (2.0_dp * q_nmr * tpiba * c * omega)
|
||||
|
||||
! either you symmetrize the current ...
|
||||
do ispin = 1, nspin
|
||||
call symmetrize_field(j_bare(:,:,:,ispin),1)
|
||||
end do
|
||||
|
||||
!
|
||||
! calculate the susceptibility
|
||||
!
|
||||
! F_{ij} = (2 - \delta_{ij}) Q_{ij}
|
||||
if ( tevaluate_chi ) then
|
||||
do ipol = 1, 3
|
||||
do jpol = 1, 3
|
||||
f_pGv(ipol,jpol,:) = 2.0_dp*q_pGv(ipol,jpol,:)
|
||||
if (ipol == jpol) f_pGv(ipol,jpol,:) = q_pGv(ipol,jpol,:)
|
||||
|
||||
f_vGv(ipol,jpol,:) = 2.0_dp*q_vGv(ipol,jpol,:)
|
||||
if (ipol == jpol) f_vGv(ipol,jpol,:) = q_vGv(ipol,jpol,:)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! compute chi_bare both pGv and vGv terms
|
||||
chi_bare_pGv(:,:) = f_pGv(:,:,1) - 2.0_dp*f_pGv(:,:,0) + f_pGv(:,:,-1)
|
||||
chi_bare_pGv(:,:) = -0.50_dp * chi_bare_pGv(:,:) / (c * q_nmr * tpiba)**2
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(5X,''chi_bare pGv (HH) in paratec units:'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_pGv(:,:) * c**2
|
||||
endif
|
||||
call sym_cart_tensor(chi_bare_pGv)
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_pGv(:,:) * c**2
|
||||
endif
|
||||
|
||||
chi_bare_vGv(:,:) = f_vGv(:,:,1) - 2.0_dp*f_vGv(:,:,0) + f_vGv(:,:,-1)
|
||||
chi_bare_vGv(:,:) = -0.50_dp * chi_bare_vGv(:,:) / (c * q_nmr * tpiba)**2
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(5X,''chi_bare vGv (VV) in paratec units:'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_vGv(:,:) * c**2
|
||||
endif
|
||||
call sym_cart_tensor(chi_bare_vGv)
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_vGv(:,:) * c**2
|
||||
endif
|
||||
end if
|
||||
|
||||
! compute induced field
|
||||
do ipol = 1, 3
|
||||
call biot_savart(ipol)
|
||||
enddo
|
||||
|
||||
! compute chemical shifts
|
||||
!call compute_sigma_bare(chi_bare_pGv)
|
||||
|
||||
|
||||
d_omega = omega / REAL ( nr1 * nr2 * nr3, DP )
|
||||
|
||||
!
|
||||
! ***************** spin-orbit-bare *******************
|
||||
!
|
||||
|
||||
! <apsi> TMPTMPTMP PLEASE CHECK FOR VANDERBILT/HARD GRIDS
|
||||
allocate ( grad_vr ( 3, nrxx ), v_local ( nrxx, nspin ) )
|
||||
do ispin = 1, nspin
|
||||
v_local(:,ispin) = vltot(:) + vr(:,ispin)
|
||||
end do
|
||||
call gradient ( nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, v_local, &
|
||||
ngm, g, nl, alat, grad_vr )
|
||||
grad_vr = grad_vr * units_Ry2Ha
|
||||
deallocate ( v_local )
|
||||
! </apsi>
|
||||
|
||||
do ipol = 1, 3
|
||||
delta_g_bare ( ipol, 1 ) = SUM ( &
|
||||
( j_bare(:,2,ipol,s_maj)-j_bare(:,2,ipol,s_min)) * grad_vr ( 3, : ) &
|
||||
-(j_bare(:,3,ipol,s_maj)-j_bare(:,3,ipol,s_min)) * grad_vr ( 2, : ) )
|
||||
delta_g_bare ( ipol, 2 ) = SUM ( &
|
||||
( j_bare(:,3,ipol,s_maj)-j_bare(:,3,ipol,s_min)) * grad_vr ( 1, : ) &
|
||||
-(j_bare(:,1,ipol,s_maj)-j_bare(:,1,ipol,s_min)) * grad_vr ( 3, : ) )
|
||||
delta_g_bare ( ipol, 3 ) = SUM ( &
|
||||
( j_bare(:,1,ipol,s_maj)-j_bare(:,1,ipol,s_min)) * grad_vr ( 2, : ) &
|
||||
-(j_bare(:,2,ipol,s_maj)-j_bare(:,2,ipol,s_min)) * grad_vr ( 1, : ) )
|
||||
end do
|
||||
deallocate ( grad_vr )
|
||||
|
||||
#ifdef __PARA
|
||||
call reduce(9, delta_g_bare)
|
||||
#endif
|
||||
|
||||
delta_g_bare = delta_g_bare * d_omega
|
||||
delta_g_bare = - delta_g_bare * g_prime / 2.0_dp / c * 1.0e6
|
||||
|
||||
!
|
||||
! ***************** spin-other-orbit *******************
|
||||
!
|
||||
|
||||
! This is the form used in the implementation in 'paratec'
|
||||
|
||||
! calculate the spin-other-orbit term a'la paratec:
|
||||
! int_r j_up(r) x v_h[n_unpaired] d^3r
|
||||
allocate ( grad_vh ( 3, nrxx ), vh ( nrxx ) )
|
||||
call v_h( rho(:,s_maj)-rho(:,s_min), nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, &
|
||||
nl, ngm, gg, gstart, 1, alat, omega, e_hartree, charge, vh )
|
||||
call gradient ( nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, vh, &
|
||||
ngm, g, nl, alat, grad_vh )
|
||||
grad_vh = grad_vh * units_Ry2Ha
|
||||
deallocate ( vh )
|
||||
|
||||
! -2*j_bare_dw(r) because of the self-interaction correction:
|
||||
! j_bare(r) - [j_bare_up(r)-j_bare_dw(r)] = -2*j_bare_dw(r)
|
||||
do ipol = 1, 3
|
||||
delta_g_soo ( ipol, 1 ) = 2 * SUM &
|
||||
( j_bare(:,2,ipol,s_min) * grad_vh ( 3, : ) &
|
||||
- j_bare(:,3,ipol,s_min) * grad_vh ( 2, : ) )
|
||||
delta_g_soo ( ipol, 2 ) = 2 * SUM &
|
||||
( j_bare(:,3,ipol,s_min) * grad_vh ( 1, : ) &
|
||||
- j_bare(:,1,ipol,s_min) * grad_vh ( 3, : ) )
|
||||
delta_g_soo ( ipol, 3 ) = 2 * SUM &
|
||||
( j_bare(:,1,ipol,s_min) * grad_vh ( 2, : ) &
|
||||
- j_bare(:,2,ipol,s_min) * grad_vh ( 1, : ) )
|
||||
end do
|
||||
deallocate ( grad_vh )
|
||||
|
||||
#ifdef __PARA
|
||||
call reduce ( 9, delta_g_soo )
|
||||
#endif
|
||||
|
||||
delta_g_soo = delta_g_soo * d_omega
|
||||
delta_g_soo = delta_g_soo * 2 / c * 1.0e6
|
||||
|
||||
! This is obtained using the equation (7) of Pickard et Mauri, PRL 88/086043
|
||||
if ( tcalculate_correct_delta_g_soo ) then
|
||||
!<apsi> The G=0 term is still not yet in
|
||||
!</apsi>
|
||||
do jpol = 1, 3
|
||||
do ipol = 1, 3
|
||||
delta_g_soo_2 ( ipol, jpol ) = SUM ( b_ind_r(:,ipol,jpol) &
|
||||
* ( rho(:,s_maj)-rho(:,s_min) ) )
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
||||
delta_g_soo_2 = delta_g_soo_2 * d_omega
|
||||
delta_g_soo_2 = delta_g_soo_2 * 2 * 1.0e6
|
||||
|
||||
!
|
||||
! ***************** relativistic-mass-correction *******************
|
||||
!
|
||||
|
||||
delta_rmc = delta_rmc / c ** 2 * g_e * units_Ry2Ha * 1.0e6
|
||||
delta_g_rmc = 0.0_DP
|
||||
do i = 1, 3
|
||||
delta_g_rmc(i,i) = delta_rmc
|
||||
end do
|
||||
|
||||
!
|
||||
! ***************** relativistic-mass-correction gipaw *******************
|
||||
!
|
||||
|
||||
write(6,*) "RMC: ", gipaw_delta_rmc, 1e6 / c ** 2
|
||||
gipaw_delta_rmc = gipaw_delta_rmc / c ** 2 * g_e * units_Ry2Ha * 1e6
|
||||
delta_g_rmc_gipaw = 0.0_DP
|
||||
do i = 1, 3
|
||||
delta_g_rmc_gipaw(i,i) = gipaw_delta_rmc
|
||||
end do
|
||||
|
||||
!
|
||||
! ***************** diamagnetic reconstruction *******************
|
||||
!
|
||||
|
||||
! symmetrize tensors
|
||||
!call trntns (delta_g_diamagn, at, bg, -1)
|
||||
!call symz(delta_g_diamagn, nsym, s, 1, irt)
|
||||
!call trntns (delta_g_diamagn, at, bg, 1)
|
||||
|
||||
delta_g_diamagn(:,:) = delta_g_diamagn(:,:) &
|
||||
* g_prime / 4 * units_Ry2Ha * 1e6
|
||||
|
||||
!
|
||||
! ***************** paramagnetic reconstruction *******************
|
||||
!
|
||||
|
||||
! symmetrize tensors
|
||||
!call trntns (delta_g_paramagn, at, bg, -1)
|
||||
!call symz(delta_g_paramagn, nsym, s, 1, irt)
|
||||
!call trntns (delta_g_paramagn, at, bg, 1)
|
||||
|
||||
delta_g_paramagn(:,:) = delta_g_paramagn(:,:) * g_prime / 2 * 1e6 &
|
||||
* units_Ry2Ha
|
||||
|
||||
!
|
||||
! ***************** total delta_g *******************
|
||||
!
|
||||
|
||||
delta_g_total = delta_g_rmc + delta_g_rmc_gipaw + delta_g_bare &
|
||||
+ delta_g_diamagn + delta_g_paramagn + delta_g_soo
|
||||
|
||||
write (stdout,*)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - relativistic-mass-correction'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_rmc(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - relativistic-mass-correction gipaw'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_rmc_gipaw(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - spin-orbit-bare'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_bare(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - spin-orbit diamagnetic correction (GIPAW)'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_diamagn(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - spin-orbit paramagnetic correction (GIPAW)'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_paramagn(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - spin-other-orbit'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_soo(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - spin-other-orbit, version 2'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_soo_2(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
write (stdout,*) 'Delta g - total'
|
||||
write (stdout,*)
|
||||
write ( stdout, '(3(5X,3(F12.6,2X)/))' ) delta_g_total(:,:)
|
||||
write (stdout,*) '**********************************************'
|
||||
write (stdout,*)
|
||||
|
||||
deallocate(p_evc, vel_evc, aux, g_vel_evc, j_bare, b_ind)
|
||||
|
||||
CONTAINS
|
||||
|
||||
!====================================================================
|
||||
! compute p_k|evc>, v_k|evc> and G_k v_{k+q,k}|evc>
|
||||
!====================================================================
|
||||
SUBROUTINE apply_operators
|
||||
implicit none
|
||||
integer ipol
|
||||
|
||||
do ipol = 1, 3
|
||||
call apply_p(evc, p_evc(1,1,ipol), ik, ipol, q)
|
||||
call apply_vel(evc, vel_evc(1,1,ipol), ik, ipol, q)
|
||||
! necessary because aux is overwritten by subroutine greenfunction
|
||||
aux(:,:) = vel_evc(:,:,ipol)
|
||||
call greenfunction(ik, aux, g_vel_evc(1,1,ipol), q)
|
||||
enddo
|
||||
END SUBROUTINE apply_operators
|
||||
|
||||
|
||||
!====================================================================
|
||||
! add contribution the Q tensors
|
||||
! Q_{\alpha,\beta} += <(e_i \times ul)_\alpha | (e_i \times ur)_\beta>
|
||||
!====================================================================
|
||||
SUBROUTINE add_to_tensor(qt, ul, ur)
|
||||
implicit none
|
||||
real(dp), intent(inout) :: qt(3,3)
|
||||
complex(dp), intent(in) :: ul(npwx,nbnd,3), ur(npwx,nbnd,3)
|
||||
real(dp) :: braket
|
||||
integer :: ibnd, ia, ib, comp_ia, comp_ib, ind(3,3), mult(3,3)
|
||||
|
||||
! index for the cross product
|
||||
ind(:,1) = (/ 1, 3, 2/); mult(:,1) = (/ 0,-1, 1 /)
|
||||
ind(:,2) = (/ 3, 2, 1/); mult(:,2) = (/ 1, 0,-1 /)
|
||||
ind(:,3) = (/ 2, 1, 3/); mult(:,3) = (/-1, 1, 0 /)
|
||||
|
||||
do ia = 1, 3 ! ia = alpha
|
||||
comp_ia = ind(ia,i)
|
||||
if (mult(ia,i) == 0) cycle
|
||||
|
||||
do ib = 1, 3 ! ib = beta
|
||||
comp_ib = ind(ib,i)
|
||||
if (mult(ib,i) == 0) cycle
|
||||
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
braket = real(ZDOTC(npw, ul(1,ibnd,comp_ia), 1, &
|
||||
ur(1,ibnd,comp_ib), 1), DP)
|
||||
qt(ia,ib) = qt(ia,ib) + wg(ibnd,ik) * &
|
||||
braket * mult(ia,i) * mult(ib,i)
|
||||
enddo ! ibnd
|
||||
|
||||
enddo ! ib
|
||||
enddo ! ia
|
||||
#ifdef __PARA
|
||||
call reduce(9, qt)
|
||||
#endif
|
||||
END SUBROUTINE add_to_tensor
|
||||
|
||||
|
||||
!====================================================================
|
||||
! add contribution the the current
|
||||
! j(r)_{\alpha,\beta} += <ul|J(r)|(B\times e_i \cdot ur)>
|
||||
!====================================================================
|
||||
SUBROUTINE add_to_current(j, ul, ur)
|
||||
implicit none
|
||||
real(dp), intent(inout) :: j(nrxxs,3,3)
|
||||
complex(dp), intent(in) :: ul(npwx,nbnd), ur(npwx,nbnd,3)
|
||||
real(dp) :: fact
|
||||
integer :: ibdir, icomp, ind(3,3), mult(3,3)
|
||||
|
||||
! index for the cross product
|
||||
ind(:,1) = (/ 1, 3, 2/); mult(:,1) = (/ 0,-1, 1 /)
|
||||
ind(:,2) = (/ 3, 2, 1/); mult(:,2) = (/ 1, 0,-1 /)
|
||||
ind(:,3) = (/ 2, 1, 3/); mult(:,3) = (/-1, 1, 0 /)
|
||||
|
||||
! loop over B direction
|
||||
do ibdir = 1, 3
|
||||
if (i == ibdir) cycle
|
||||
icomp = ind(ibdir, i)
|
||||
fact = real(mult(ibdir,i)*isign)
|
||||
call j_para(fact, ul(1,1), ur(1,1,icomp), ik, q, j(1,1,ibdir))
|
||||
enddo
|
||||
END SUBROUTINE add_to_current
|
||||
|
||||
!====================================================================
|
||||
! ...
|
||||
!====================================================================
|
||||
SUBROUTINE relativistic_mass_correction ( rmc_gipaw )
|
||||
|
||||
USE atom, ONLY : r, rab
|
||||
USE ions_base, ONLY : nat, ityp, ntyp => nsp
|
||||
USE nmr_module, ONLY : c
|
||||
|
||||
implicit none
|
||||
|
||||
! Arguments
|
||||
real(dp), intent(inout):: rmc_gipaw
|
||||
|
||||
integer :: l1, m1, lm1, l2, m2, lm2, ih, ikb, nbs1, jh, jkb, nbs2
|
||||
integer :: nt, ibnd, na, lm, nrc, ijkb0
|
||||
complex(dp) :: efg_corr
|
||||
complex(dp) :: bec_product
|
||||
|
||||
efg_corr = 0.0_dp
|
||||
|
||||
do ibnd = 1, nbnd
|
||||
ijkb0 = 0
|
||||
do nt = 1, ntyp
|
||||
do na = 1, nat
|
||||
if (ityp (na) .eq.nt) then
|
||||
do ih = 1, paw_nh (nt)
|
||||
ikb = ijkb0 + ih
|
||||
nbs1=paw_indv(ih,nt)
|
||||
l1=paw_nhtol(ih,nt)
|
||||
m1=paw_nhtom(ih,nt)
|
||||
lm1=m1+l1**2
|
||||
do jh = 1, paw_nh (nt)
|
||||
jkb = ijkb0 + jh
|
||||
nbs2=paw_indv(jh,nt)
|
||||
l2=paw_nhtol(jh,nt)
|
||||
m2=paw_nhtom(jh,nt)
|
||||
lm2=m2+l2**2
|
||||
|
||||
IF ( l1 /= l2 ) CYCLE
|
||||
IF ( m1 /= m2 ) CYCLE
|
||||
|
||||
bec_product = paw_becp(jkb,ibnd) &
|
||||
* CONJG(paw_becp(ikb,ibnd))
|
||||
|
||||
efg_corr = efg_corr &
|
||||
+ bec_product &
|
||||
* radial_integral_rmc(nbs1,nbs2,nt) &
|
||||
* wg(ibnd,ik)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
ijkb0 = ijkb0 + paw_nh (nt)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
rmc_gipaw = REAL ( efg_corr, dp )
|
||||
|
||||
END SUBROUTINE relativistic_mass_correction
|
||||
|
||||
|
||||
!====================================================================
|
||||
! ...
|
||||
!====================================================================
|
||||
SUBROUTINE diamagnetic_correction ( diamagnetic_tensor )
|
||||
|
||||
USE atom, ONLY : r, rab
|
||||
USE ions_base, ONLY : nat, ityp, ntyp => nsp
|
||||
USE nmr_module, ONLY : c
|
||||
|
||||
implicit none
|
||||
|
||||
! Arguments
|
||||
real(dp), intent(inout):: diamagnetic_tensor(3,3)
|
||||
|
||||
integer :: l1, m1, lm1, l2, m2, lm2, ih, ikb, nbs1, jh, jkb, nbs2
|
||||
integer :: nt, ibnd, na, lm, nrc, ijkb0
|
||||
complex(dp) , allocatable :: efg_corr(:)
|
||||
complex(dp) :: bec_product
|
||||
|
||||
allocate ( efg_corr ( lmaxx**2 ) )
|
||||
efg_corr = 0.0_dp
|
||||
|
||||
do ibnd = 1, nbnd
|
||||
ijkb0 = 0
|
||||
do nt = 1, ntyp
|
||||
do na = 1, nat
|
||||
if (ityp (na) .eq.nt) then
|
||||
do ih = 1, paw_nh (nt)
|
||||
ikb = ijkb0 + ih
|
||||
nbs1=paw_indv(ih,nt)
|
||||
l1=paw_nhtol(ih,nt)
|
||||
m1=paw_nhtom(ih,nt)
|
||||
lm1=m1+l1**2
|
||||
do jh = 1, paw_nh (nt)
|
||||
jkb = ijkb0 + jh
|
||||
nbs2=paw_indv(jh,nt)
|
||||
l2=paw_nhtol(jh,nt)
|
||||
m2=paw_nhtom(jh,nt)
|
||||
lm2=m2+l2**2
|
||||
|
||||
bec_product = paw_becp(jkb,ibnd) &
|
||||
* CONJG(paw_becp(ikb,ibnd))
|
||||
|
||||
!<apsi> s/non-trace-zero component
|
||||
! 2/3 to separate the non-trace vanishing component
|
||||
! 1/(2c^2) from the equation (59) in PM-PRB
|
||||
IF ( l1 == l2 .AND. m1 == m2 ) THEN
|
||||
diamagnetic_tensor(1,1) &
|
||||
= diamagnetic_tensor(1,1) &
|
||||
+ 2.0_dp / 3.0_dp * bec_product &
|
||||
* radial_integral_diamagnetic_so(nbs1,nbs2,nt) &
|
||||
* wg(ibnd,ik) / c**2
|
||||
diamagnetic_tensor(2,2) &
|
||||
= diamagnetic_tensor(2,2) &
|
||||
+ 2.0_dp / 3.0_dp * bec_product &
|
||||
* radial_integral_diamagnetic_so(nbs1,nbs2,nt) &
|
||||
* wg(ibnd,ik) / c**2
|
||||
diamagnetic_tensor(3,3) &
|
||||
= diamagnetic_tensor(3,3) &
|
||||
+ 2.0_dp / 3.0_dp * bec_product &
|
||||
* radial_integral_diamagnetic_so(nbs1,nbs2,nt) &
|
||||
* wg(ibnd,ik) / c**2
|
||||
END IF
|
||||
|
||||
! 2/3 to separate the non-trace vanishing component
|
||||
do lm = 5, 9
|
||||
efg_corr(lm) = efg_corr(lm) &
|
||||
+ bec_product / 3.0_dp &
|
||||
* radial_integral_diamagnetic_so(nbs1,nbs2,nt) &
|
||||
* ap(lm,lm1,lm2) * wg(ibnd,ik) / c**2
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
ijkb0 = ijkb0 + paw_nh (nt)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
write(6,'("CCC1",5F14.8)') efg_corr(5:9)
|
||||
|
||||
!
|
||||
! transform in cartesian coordinates
|
||||
!
|
||||
|
||||
efg_corr(5:9) = - sqrt(4.0_dp * pi/5.0_dp) * efg_corr(5:9)
|
||||
|
||||
diamagnetic_tensor(1,1) = diamagnetic_tensor(1,1) &
|
||||
+ sqrt(3.0_dp) * efg_corr(8) - efg_corr(5)
|
||||
diamagnetic_tensor(2,2) = diamagnetic_tensor(2,2) &
|
||||
- sqrt(3.0_dp) * efg_corr(8) - efg_corr(5)
|
||||
diamagnetic_tensor(3,3) = diamagnetic_tensor(3,3) &
|
||||
+ efg_corr(5) * 2.0_dp
|
||||
diamagnetic_tensor(1,2) = diamagnetic_tensor(1,2) &
|
||||
+ efg_corr(9) * sqrt(3.0_dp)
|
||||
diamagnetic_tensor(2,1) = diamagnetic_tensor(1,2)
|
||||
diamagnetic_tensor(1,3) = diamagnetic_tensor(1,3) &
|
||||
- efg_corr(6) * sqrt(3.0_dp)
|
||||
diamagnetic_tensor(3,1) = diamagnetic_tensor(1,3)
|
||||
diamagnetic_tensor(2,3) = diamagnetic_tensor(2,3) &
|
||||
- efg_corr(7) * sqrt(3.0_dp)
|
||||
diamagnetic_tensor(3,2) = diamagnetic_tensor(2,3)
|
||||
|
||||
! efg_corr(5) = 3z^2-1
|
||||
! efg_corr(6) = -xz
|
||||
! efg_corr(7) = -yz
|
||||
! efg_corr(8) = x^2-y^2
|
||||
! efg_corr(9) = xy
|
||||
|
||||
deallocate ( efg_corr )
|
||||
|
||||
END SUBROUTINE diamagnetic_correction
|
||||
|
||||
|
||||
!====================================================================
|
||||
! ...
|
||||
!====================================================================
|
||||
SUBROUTINE paramagnetic_correction ( paramagnetic_tensor )
|
||||
|
||||
USE ions_base, ONLY : nat, ityp, ntyp => nsp
|
||||
USE nmr_module, ONLY : c
|
||||
|
||||
implicit none
|
||||
|
||||
! Arguments
|
||||
real(dp), intent(inout):: paramagnetic_tensor(3,3)
|
||||
|
||||
integer :: l1, m1, lm1, l2, m2, lm2, ih, ikb, nbs1, jh, jkb, nbs2
|
||||
integer :: nt, ibnd, na, lm, j, ijkb0, ipol
|
||||
complex(dp) :: bec_product
|
||||
complex(dp) , allocatable :: efg_corr(:)
|
||||
|
||||
integer, parameter :: ng_ = 27, lmax2_ = 16
|
||||
integer :: mg, i1, i2, i3
|
||||
real(DP) :: g_ (3, ng_), gg_ (ng_)
|
||||
real(DP) :: ylm_ (ng_,lmax2_)
|
||||
|
||||
!--------------------------------------------------------------------------
|
||||
|
||||
allocate (efg_corr(3))
|
||||
|
||||
!
|
||||
! calculation of the reconstruction part
|
||||
!
|
||||
|
||||
do ipol = 1, 3
|
||||
|
||||
if ( ipol == i ) cycle !TESTTESTTEST
|
||||
|
||||
call ccalbec (paw_nkb, npwx, npw, nbnd, paw_becp2, paw_vkb, &
|
||||
g_vel_evc(1,1,ipol))
|
||||
|
||||
efg_corr = 0.0_dp
|
||||
|
||||
do ibnd = 1, nbnd
|
||||
ijkb0 = 0
|
||||
do nt = 1, ntyp
|
||||
do na = 1, nat
|
||||
|
||||
if (ityp (na) .eq.nt) then
|
||||
do ih = 1, paw_nh (nt)
|
||||
ikb = ijkb0 + ih
|
||||
nbs1=paw_indv(ih,nt)
|
||||
l1=paw_nhtol(ih,nt)
|
||||
m1=paw_nhtom(ih,nt)
|
||||
lm1=m1+l1**2
|
||||
|
||||
do jh = 1, paw_nh (nt)
|
||||
jkb = ijkb0 + jh
|
||||
nbs2=paw_indv(jh,nt)
|
||||
l2=paw_nhtol(jh,nt)
|
||||
m2=paw_nhtom(jh,nt)
|
||||
lm2=m2+l2**2
|
||||
|
||||
if ( l1 /= l2 ) cycle
|
||||
|
||||
bec_product = CONJG(paw_becp(ikb,ibnd)) &
|
||||
* paw_becp2(jkb,ibnd)
|
||||
|
||||
efg_corr(1) = efg_corr(1) &
|
||||
+ bec_product &
|
||||
* radial_integral_paramagnetic_so(nbs1,nbs2,nt) &
|
||||
* lx ( lm1, lm2 ) * wg(ibnd,ik) / c ** 2
|
||||
efg_corr(2) = efg_corr(2) &
|
||||
+ bec_product &
|
||||
* radial_integral_paramagnetic_so(nbs1,nbs2,nt) &
|
||||
* ly ( lm1, lm2 ) * wg(ibnd,ik) / c ** 2
|
||||
efg_corr(3) = efg_corr(3) &
|
||||
+ bec_product &
|
||||
* radial_integral_paramagnetic_so(nbs1,nbs2,nt) &
|
||||
* lz ( lm1, lm2 ) * wg(ibnd,ik) / c ** 2
|
||||
if (lz(lm1,lm2)/=0.and.ibnd==1.and.l1==-1) then
|
||||
write(6,*) "ZZZ3: ", &
|
||||
ibnd, lm1, lm2, &
|
||||
bec_product, &
|
||||
radial_integral_paramagnetic_so(nbs1,nbs2,nt)
|
||||
|
||||
end if
|
||||
enddo
|
||||
enddo
|
||||
ijkb0 = ijkb0 + paw_nh (nt)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
paramagnetic_tensor ( :, ipol ) = REAL ( efg_corr, dp )
|
||||
|
||||
write(6,'("DDD1",2I3,3(F16.7,2X)') &
|
||||
ipol, i*isign, REAL ( efg_corr(1:3) ) * 1e6
|
||||
! stop
|
||||
|
||||
end do
|
||||
|
||||
deallocate(efg_corr)
|
||||
|
||||
END SUBROUTINE paramagnetic_correction
|
||||
|
||||
!====================================================================
|
||||
! ...
|
||||
!====================================================================
|
||||
SUBROUTINE add_to_sigma_para( paramagnetic_correction, sigma_paramagnetic )
|
||||
implicit none
|
||||
real(dp), intent(in) :: paramagnetic_correction(3,3)
|
||||
real(dp), intent(inout) :: sigma_paramagnetic(3,3)
|
||||
real(dp) :: fact
|
||||
integer :: ibdir, icomp, ipol, ind(3,3), mult(3,3)
|
||||
|
||||
! index for the cross product
|
||||
ind(:,1) = (/ 1, 3, 2/); mult(:,1) = (/ 0,-1, 1 /)
|
||||
ind(:,2) = (/ 3, 2, 1/); mult(:,2) = (/ 1, 0,-1 /)
|
||||
ind(:,3) = (/ 2, 1, 3/); mult(:,3) = (/-1, 1, 0 /)
|
||||
|
||||
! loop over B direction
|
||||
do ibdir = 1, 3
|
||||
if (i == ibdir) cycle
|
||||
icomp = ind(ibdir, i)
|
||||
fact = real(mult(ibdir,i)*isign)
|
||||
|
||||
do ipol = 1, 3
|
||||
sigma_paramagnetic ( ipol, icomp ) &
|
||||
= sigma_paramagnetic ( ipol, icomp ) &
|
||||
+ fact * paramagnetic_correction ( ipol, ibdir ) &
|
||||
/ ( 2 * q_nmr * tpiba )
|
||||
|
||||
! Do iq=1, 3 ! loop over all q-points
|
||||
! ...
|
||||
! Do iperm0 =1,-1,-2
|
||||
! b0 = Mod(iq+iperm0+2,3) +1 ! gives different b0 for iperm0
|
||||
! p0 = Mod(iq-iperm0+2,3) +1 ! gives p0 = abs(q x b0)
|
||||
! ...
|
||||
! Call take_nonloc_deriv_kq_k(p0,k_gspace, rk(1),u_k(1,i),&
|
||||
! ...
|
||||
! para_corr_shift_tot(b0,:,:,:) = para_corr_shift_tot(b0,:,:,:)-&
|
||||
! Real(iperm0,dp)*Real(iqsign,dp)*para_corr_shift*&
|
||||
! kpoints%w(irk)/qmag/dtwo/Real(crys%nspin,dp)
|
||||
|
||||
end do
|
||||
enddo
|
||||
END SUBROUTINE add_to_sigma_para
|
||||
|
||||
END SUBROUTINE g_tensor_crystal
|
||||
|
|
@ -0,0 +1,144 @@
|
|||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE greenfunction(ik, psi, g_psi, q)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Apply the Green function operator
|
||||
! ... (formula here)
|
||||
! ... We use Hartree atomic units; since G is the inverse of an
|
||||
! ... energy: G => G / ryd_to_hartree
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE becmod, ONLY : becp
|
||||
USE wavefunctions_module, ONLY : evc
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- parameters ---------------------------------------------------------
|
||||
IMPLICIT none
|
||||
INTEGER, INTENT(IN) :: ik
|
||||
COMPLEX(DP), INTENT(INOUT) :: psi(npwx,nbnd) ! psi is changed on output!!!
|
||||
COMPLEX(DP), INTENT(OUT) :: g_psi(npwx,nbnd)
|
||||
REAL(DP) :: q(3)
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
real(dp), parameter :: ryd_to_hartree = 0.5d0
|
||||
complex(dp), allocatable :: ps(:,:), work (:,:)
|
||||
real(dp), allocatable :: h_diag (:,:), eprec (:)
|
||||
real(dp) :: anorm, thresh, gk(3), dxk(3)
|
||||
integer :: ibnd, ig, lter
|
||||
logical :: conv_root, q_is_zero
|
||||
complex(dp), external :: ZDOTC
|
||||
external ch_psi_all, cg_psi
|
||||
|
||||
|
||||
call start_clock ('greenf')
|
||||
|
||||
! allocate memory
|
||||
allocate (work(npwx, MAX(nkb,1)), ps(nbnd,nbnd), h_diag(npwx,nbnd), &
|
||||
eprec(nbnd))
|
||||
|
||||
! check if |q| is zero
|
||||
q_is_zero = .false.
|
||||
if (sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3)) < 1d-8) q_is_zero = .true.
|
||||
if (q_is_zero) evq = evc
|
||||
|
||||
|
||||
!====================================================================
|
||||
! apply -Q_{k+q} to the r.h.s.
|
||||
!====================================================================
|
||||
! project on <evq|: ps(i,j) = <evq(i)|psi(j)>
|
||||
ps = (0.d0,0.d0)
|
||||
CALL ZGEMM('C', 'N', nbnd_occ (ik), nbnd_occ (ik), npw, &
|
||||
(1.d0,0.d0), evq(1,1), npwx, psi(1,1), npwx, (0.d0,0.d0), &
|
||||
ps(1,1), nbnd)
|
||||
#ifdef __PARA
|
||||
call reduce (2 * nbnd * nbnd_occ (ik), ps(1,1))
|
||||
#endif
|
||||
|
||||
!! this is the case with overlap (ultrasoft)
|
||||
!! g_psi is used as work space to store S|evc>
|
||||
!!
|
||||
!!CALL ccalbec (nkb, npwx, npw, nbnd_occ(ik), becp, vkb, evc)
|
||||
!!CALL s_psi (npwx, npw, nbnd_occ(ik), evc, g_psi)
|
||||
!! |psi> = -(|psi> - S|evc><evc|psi>)
|
||||
!!
|
||||
!!CALL ZGEMM( 'N', 'N', npw, nbnd_occ(ik), nbnd_occ(ik), &
|
||||
!! (1.d0,0.d0), g_psi(1,1), npwx, ps(1,1), nbnd, (-1.d0,0.d0), &
|
||||
!! psi(1,1), npwx )
|
||||
|
||||
! |psi> = -(1 - |evq><evq|) |psi>
|
||||
CALL ZGEMM('N', 'N', npw, nbnd_occ(ik), nbnd_occ(ik), &
|
||||
(1.d0,0.d0), evq(1,1), npwx, ps(1,1), nbnd, (-1.d0,0.d0), &
|
||||
psi(1,1), npwx)
|
||||
|
||||
|
||||
!====================================================================
|
||||
! solve the linear system (apply G_{k+q})
|
||||
!====================================================================
|
||||
! convergence treshold
|
||||
!thresh = 1d-12
|
||||
thresh = sqrt(conv_threshold) ! sqrt(of that of PARATEC)
|
||||
|
||||
! preconditioning of the linear system
|
||||
do ibnd = 1, nbnd_occ (ik)
|
||||
conv_root = .true.
|
||||
do ig = 1, npw
|
||||
work (ig,1) = g2kin (ig) * evc (ig, ibnd)
|
||||
enddo
|
||||
eprec (ibnd) = 1.35d0 * ZDOTC (npw, evq (1, ibnd), 1, work, 1)
|
||||
enddo
|
||||
#ifdef __PARA
|
||||
call reduce (nbnd_occ (ik), eprec)
|
||||
#endif
|
||||
do ibnd = 1, nbnd_occ (ik)
|
||||
do ig = 1, npw
|
||||
h_diag (ig, ibnd) = 1.d0 / max (1.0d0, g2kin (ig) / eprec (ibnd) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! use the hamiltonian at k+q
|
||||
do ig = 1, npw
|
||||
gk(1) = (xk(1,ik) + g(1,igk(ig)) + q(1)) * tpiba
|
||||
gk(2) = (xk(2,ik) + g(2,igk(ig)) + q(2)) * tpiba
|
||||
gk(3) = (xk(3,ik) + g(3,igk(ig)) + q(3)) * tpiba
|
||||
g2kin (ig) = gk(1)**2 + gk(2)**2 + gk(3)**2
|
||||
enddo
|
||||
|
||||
if (.not. q_is_zero) then
|
||||
dxk = xk(:,ik) + q
|
||||
call init_us_2(npw, igk, dxk, vkb)
|
||||
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, psi)
|
||||
endif
|
||||
|
||||
! initial guess
|
||||
g_psi(:,:) = (0.d0, 0.d0)
|
||||
|
||||
! solve linear system
|
||||
call cgsolve_all (ch_psi_all, cg_psi, et(1,ik), psi, g_psi, &
|
||||
h_diag, npwx, npw, thresh, ik, lter, conv_root, anorm, &
|
||||
nbnd_occ(ik) )
|
||||
!!write(stdout, '(5X,''cgsolve_all converged in '',I3,'' iterations'')') &
|
||||
!! lter
|
||||
|
||||
if (.not.conv_root) WRITE( stdout, '(5x,"ik",i4," ibnd",i4, &
|
||||
& " linter: root not converged ",e10.3)') &
|
||||
ik, ibnd, anorm
|
||||
|
||||
! convert to Hartree
|
||||
g_psi(:,:) = g_psi(:,:) / ryd_to_hartree
|
||||
|
||||
call flush_unit( stdout )
|
||||
call stop_clock('greenf')
|
||||
|
||||
! free memory
|
||||
deallocate (work, h_diag, eprec, ps)
|
||||
|
||||
END SUBROUTINE greenfunction
|
||||
|
|
@ -0,0 +1,102 @@
|
|||
!
|
||||
! 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 h_psiq (lda, n, m, psi, hpsi, spsi)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! This routine computes the product of the Hamiltonian
|
||||
! and of the S matrix with a m wavefunctions contained
|
||||
! in psi. It first computes the bec matrix of these
|
||||
! wavefunctions and then with the routines hus_1psi and
|
||||
! s_psi computes for each band the required products
|
||||
!
|
||||
|
||||
use pwcom
|
||||
USE wavefunctions_module, ONLY: psic
|
||||
USE becmod, ONLY: becp
|
||||
USE kinds, only : DP
|
||||
use nmr_module
|
||||
implicit none
|
||||
!
|
||||
! Here the local variables
|
||||
!
|
||||
integer :: ibnd
|
||||
! counter on bands
|
||||
|
||||
integer :: lda, n, m
|
||||
! input: the leading dimension of the array psi
|
||||
! input: the real dimension of psi
|
||||
! input: the number of psi to compute
|
||||
integer :: j
|
||||
! do loop index
|
||||
|
||||
complex(DP) :: psi (lda, m), hpsi (lda, m), spsi (lda, m)
|
||||
! input: the functions where to apply H and S
|
||||
! output: H times psi
|
||||
! output: S times psi (Us PP's only)
|
||||
|
||||
|
||||
call start_clock ('h_psiq')
|
||||
call start_clock ('init')
|
||||
|
||||
call ccalbec (nkb, npwx, n, m, becp, vkb, psi)
|
||||
!
|
||||
! Here we apply the kinetic energy (k+G)^2 psi
|
||||
!
|
||||
do ibnd = 1, m
|
||||
do j = 1, n
|
||||
hpsi (j, ibnd) = g2kin (j) * psi (j, ibnd)
|
||||
enddo
|
||||
enddo
|
||||
call stop_clock ('init')
|
||||
!
|
||||
! the local potential V_Loc psi. First the psi in real space
|
||||
!
|
||||
|
||||
do ibnd = 1, m
|
||||
call start_clock ('firstfft')
|
||||
psic(:) = (0.d0, 0.d0)
|
||||
!do j = 1, n
|
||||
! psic (nls(igk(j))) = psi (j, ibnd)
|
||||
!enddo
|
||||
psic (nls(igk(1:n))) = psi(1:n, ibnd)
|
||||
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
|
||||
call stop_clock ('firstfft')
|
||||
!
|
||||
! and then the product with the potential vrs = (vltot+vr) on the smoo
|
||||
!
|
||||
!do j = 1, nrxxs
|
||||
! psic (j) = psic (j) * vrs (j, current_spin)
|
||||
!enddo
|
||||
psic (1:nrxxs) = psic (1:nrxxs) * vrs (1:nrxxs, current_spin)
|
||||
!
|
||||
! back to reciprocal space
|
||||
!
|
||||
call start_clock ('secondfft')
|
||||
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)
|
||||
!
|
||||
! addition to the total product
|
||||
!
|
||||
!do j = 1, n
|
||||
! hpsi (j, ibnd) = hpsi (j, ibnd) + psic (nls(igk(j)))
|
||||
!enddo
|
||||
hpsi (1:n, ibnd) = hpsi (1:n, ibnd) + psic (nls(igk(1:n)))
|
||||
call stop_clock ('secondfft')
|
||||
enddo
|
||||
!
|
||||
! Here the product with the non local potential V_NL psi
|
||||
!
|
||||
|
||||
call add_vuspsi (lda, n, m, psi, hpsi)
|
||||
|
||||
call s_psi (lda, n, m, psi, spsi)
|
||||
|
||||
call stop_clock ('h_psiq')
|
||||
return
|
||||
end subroutine h_psiq
|
||||
|
|
@ -0,0 +1,161 @@
|
|||
!
|
||||
! 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!
|
||||
!----------------------------------------------------------------------
|
||||
subroutine init_paw_2_no_phase (npw_, igk_, q_, vkb_)
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
! Calculates paw_beta functions (paw projectors), with
|
||||
! structure factor, for all atoms, in reciprocal space
|
||||
!
|
||||
USE kinds , ONLY : dp
|
||||
USE constants , ONLY : tpi
|
||||
USE wvfct , ONLY : npwx
|
||||
USE cell_base , ONLY : tpiba
|
||||
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau
|
||||
USE gvect , ONLY : eigts1, eigts2, eigts3, g, ig1, ig2, ig3
|
||||
USE us, ONLY : nqxq, dq
|
||||
USE paw, ONLY : paw_nkb, paw_lmaxkb, paw_nhm, paw_nh, paw_nhtol, &
|
||||
paw_nhtom, paw_indv, paw_tab, paw_nbeta, &
|
||||
paw_tab_d2y
|
||||
USE splinelib
|
||||
!
|
||||
implicit none
|
||||
!
|
||||
integer :: npw_, igk_ (npw_)
|
||||
! input: number of PW's
|
||||
! input: indices of q+G
|
||||
real(DP) :: q_(3)
|
||||
! input: q vector
|
||||
complex(DP) :: vkb_ (npwx, paw_nkb)
|
||||
! output: beta functions
|
||||
!
|
||||
! Local variables
|
||||
!
|
||||
integer :: i0,i1,i2,i3, ig, l, lm, na, nt, nb, ih, jkb
|
||||
|
||||
real(DP) :: px, ux, vx, wx, arg
|
||||
real(DP), allocatable :: gk (:,:), qg (:), vq (:), ylm (:,:), vkb1(:,:)
|
||||
|
||||
complex(DP) :: phase, pref
|
||||
complex(DP), allocatable :: sk(:)
|
||||
|
||||
!<apsi>
|
||||
integer :: startq, lastq, iq
|
||||
real(DP), allocatable :: xdata(:)
|
||||
!</apsi>
|
||||
!
|
||||
!
|
||||
if (paw_lmaxkb.lt.0) return
|
||||
call start_clock ('init_paw_2')
|
||||
allocate (vkb1( npw_,paw_nhm))
|
||||
allocate ( sk( npw_))
|
||||
allocate ( qg( npw_))
|
||||
allocate ( vq( npw_))
|
||||
allocate ( ylm( npw_, (paw_lmaxkb + 1) **2))
|
||||
allocate ( gk( 3, npw_))
|
||||
!
|
||||
|
||||
do ig = 1, npw_
|
||||
gk (1,ig) = q_(1) + g(1, igk_(ig) )
|
||||
gk (2,ig) = q_(2) + g(2, igk_(ig) )
|
||||
gk (3,ig) = q_(3) + g(3, igk_(ig) )
|
||||
qg (ig) = gk(1, ig)**2 + gk(2, ig)**2 + gk(3, ig)**2
|
||||
enddo
|
||||
!
|
||||
call ylmr2 ((paw_lmaxkb+1)**2, npw_, gk, qg, ylm)
|
||||
!
|
||||
! set now qg=|q+G| in atomic units
|
||||
!
|
||||
do ig = 1, npw_
|
||||
qg(ig) = sqrt(qg(ig))*tpiba
|
||||
enddo
|
||||
|
||||
!<ceres>
|
||||
call divide (nqxq, startq, lastq)
|
||||
startq = 1 !*TMPTMPTMP
|
||||
lastq = nqxq
|
||||
allocate(xdata(lastq-startq+1))
|
||||
do iq = startq, lastq
|
||||
xdata(iq) = (iq - 1) * dq
|
||||
enddo
|
||||
!</ceres>
|
||||
|
||||
jkb = 0
|
||||
do nt = 1, ntyp
|
||||
! calculate beta in G-space using an interpolation table
|
||||
do nb = 1, paw_nbeta (nt)
|
||||
do ig = 1, npw_
|
||||
!<ceres>
|
||||
!px = qg (ig) / dq - int (qg (ig) / dq)
|
||||
!ux = 1.d0 - px
|
||||
!vx = 2.d0 - px
|
||||
!wx = 3.d0 - px
|
||||
!i0 = qg (ig) / dq + 1
|
||||
!i1 = i0 + 1
|
||||
!i2 = i0 + 2
|
||||
!i3 = i0 + 3
|
||||
!vq (ig) = paw_tab (i0, nb, nt) * ux * vx * wx / 6.d0 + &
|
||||
! paw_tab (i1, nb, nt) * px * vx * wx / 2.d0 - &
|
||||
! paw_tab (i2, nb, nt) * px * ux * wx / 2.d0 + &
|
||||
! paw_tab (i3, nb, nt) * px * ux * vx / 6.d0
|
||||
vq(ig) = splint(xdata, paw_tab(:,nb,nt), paw_tab_d2y(:,nb,nt), qg(ig))
|
||||
!</ceres>
|
||||
enddo
|
||||
! add spherical harmonic part
|
||||
do ih = 1, paw_nh (nt)
|
||||
if (nb.eq.paw_indv (ih, nt) ) then
|
||||
l = paw_nhtol (ih, nt)
|
||||
lm = l * l + paw_nhtom (ih, nt)
|
||||
do ig = 1, npw_
|
||||
vkb1 (ig,ih) = ylm (ig, lm) * vq (ig)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! vkb1 contains all betas including angular part for type nt
|
||||
! now add the structure factor and factor (-i)^l
|
||||
!
|
||||
do na = 1, nat
|
||||
! ordering: first all betas for atoms of type 1
|
||||
! then all betas for atoms of type 2 and so on
|
||||
if (ityp (na) .eq.nt) then
|
||||
arg = (q_(1) * tau (1, na) + &
|
||||
q_(2) * tau (2, na) + &
|
||||
q_(3) * tau (3, na) ) * tpi
|
||||
phase = CMPLX (cos (arg), - sin (arg) )
|
||||
do ig = 1, npw_
|
||||
sk (ig) = eigts1 (ig1(igk_(ig)), na) * &
|
||||
eigts2 (ig2(igk_(ig)), na) * &
|
||||
eigts3 (ig3(igk_(ig)), na)
|
||||
enddo
|
||||
do ih = 1, paw_nh(nt)
|
||||
jkb = jkb + 1
|
||||
pref = (0.d0, -1.d0) ** paw_nhtol (ih, nt) ! * phase
|
||||
do ig = 1, npw_
|
||||
vkb_(ig, jkb) = vkb1 (ig,ih) * sk (ig) * pref
|
||||
enddo
|
||||
|
||||
enddo
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
deallocate (gk)
|
||||
deallocate (ylm)
|
||||
deallocate (vq)
|
||||
deallocate (qg)
|
||||
deallocate (sk)
|
||||
deallocate (vkb1)
|
||||
|
||||
call stop_clock ('init_paw_2')
|
||||
return
|
||||
end subroutine init_paw_2_no_phase
|
||||
|
||||
|
|
@ -0,0 +1,172 @@
|
|||
!
|
||||
! 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
|
||||
! for the time being, don't use splines, till it is implemented
|
||||
! also in PW
|
||||
#undef SPLINES
|
||||
#define SPLINES
|
||||
|
||||
!----------------------------------------------------------------------
|
||||
subroutine init_us_2_no_phase (npw_, igk_, q_, vkb_)
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
! Calculates beta functions (Kleinman-Bylander projectors), with
|
||||
! structure factor, for all atoms, in reciprocal space
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau
|
||||
USE cell_base, ONLY : tpiba
|
||||
USE constants, ONLY : tpi
|
||||
USE gvect, ONLY : eigts1, eigts2, eigts3, ig1, ig2, ig3, g
|
||||
USE wvfct, ONLY : npw, npwx, igk
|
||||
#ifdef SPLINES
|
||||
USE us, ONLY : nqxq, dq, tab, tab_d2y
|
||||
USE splinelib
|
||||
#else
|
||||
USE us, ONLY : nqxq, dq, tab
|
||||
#endif
|
||||
USE uspp, ONLY : nkb, vkb, nhtol, nhtolm, indv
|
||||
USE uspp_param, ONLY : lmaxkb, nbeta, nhm, nh
|
||||
!
|
||||
implicit none
|
||||
!
|
||||
integer :: npw_, igk_ (npw_)
|
||||
! input: number of PW's
|
||||
! input: indices of q+G
|
||||
real(DP) :: q_(3)
|
||||
! input: q vector
|
||||
complex(DP) :: vkb_ (npwx, nkb)
|
||||
! output: beta functions
|
||||
!
|
||||
! Local variables
|
||||
!
|
||||
integer :: i0,i1,i2,i3, ig, l, lm, na, nt, nb, ih, jkb, iq
|
||||
integer :: startq, lastq
|
||||
|
||||
real(DP) :: px, ux, vx, wx, arg
|
||||
real(DP), allocatable :: gk (:,:), qg (:), vq (:), ylm (:,:), vkb1(:,:)
|
||||
|
||||
complex(DP) :: phase, pref
|
||||
complex(DP), allocatable :: sk(:)
|
||||
real(DP), allocatable :: xdata(:)
|
||||
|
||||
!
|
||||
!
|
||||
if (lmaxkb.lt.0) return
|
||||
call start_clock ('init_us_2')
|
||||
allocate (vkb1( npw_,nhm))
|
||||
allocate ( sk( npw_))
|
||||
allocate ( qg( npw_))
|
||||
allocate ( vq( npw_))
|
||||
allocate ( ylm( npw_, (lmaxkb + 1) **2))
|
||||
allocate ( gk( 3, npw_))
|
||||
!
|
||||
do ig = 1, npw_
|
||||
gk (1,ig) = q_(1) + g(1, igk_(ig) )
|
||||
gk (2,ig) = q_(2) + g(2, igk_(ig) )
|
||||
gk (3,ig) = q_(3) + g(3, igk_(ig) )
|
||||
qg (ig) = gk(1, ig)**2 + gk(2, ig)**2 + gk(3, ig)**2
|
||||
enddo
|
||||
!
|
||||
call ylmr2 ((lmaxkb+1)**2, npw_, gk, qg, ylm)
|
||||
!
|
||||
! set now qg=|q+G| in atomic units
|
||||
!
|
||||
do ig = 1, npw_
|
||||
qg(ig) = sqrt(qg(ig))*tpiba
|
||||
enddo
|
||||
|
||||
#ifdef SPLINES
|
||||
!<ceres>
|
||||
startq = 1
|
||||
lastq = nqxq
|
||||
call divide (nqxq, startq, lastq)
|
||||
allocate(xdata(lastq-startq+1))
|
||||
do iq = startq, lastq
|
||||
xdata(iq) = (iq - 1) * dq
|
||||
enddo
|
||||
!</ceres>
|
||||
#endif
|
||||
|
||||
jkb = 0
|
||||
do nt = 1, ntyp
|
||||
! calculate beta in G-space using an interpolation table
|
||||
do nb = 1, nbeta (nt)
|
||||
do ig = 1, npw_
|
||||
#ifdef SPLINES
|
||||
vq(ig) = splint(xdata, tab(:,nb,nt), tab_d2y(:,nb,nt), qg(ig))
|
||||
#else
|
||||
px = qg (ig) / dq - int (qg (ig) / dq)
|
||||
ux = 1.d0 - px
|
||||
vx = 2.d0 - px
|
||||
wx = 3.d0 - px
|
||||
i0 = INT( qg (ig) / dq ) + 1
|
||||
i1 = i0 + 1
|
||||
i2 = i0 + 2
|
||||
i3 = i0 + 3
|
||||
!*apsi TMPTMPTMP
|
||||
if ( i3 > size(tab,1) ) then
|
||||
vq(ig) = 0.0_dp
|
||||
else
|
||||
vq (ig) = tab (i0, nb, nt) * ux * vx * wx / 6.d0 + &
|
||||
tab (i1, nb, nt) * px * vx * wx / 2.d0 - &
|
||||
tab (i2, nb, nt) * px * ux * wx / 2.d0 + &
|
||||
tab (i3, nb, nt) * px * ux * vx / 6.d0
|
||||
end if
|
||||
#endif
|
||||
enddo
|
||||
! add spherical harmonic part
|
||||
do ih = 1, nh (nt)
|
||||
if (nb.eq.indv (ih, nt) ) then
|
||||
l = nhtol (ih, nt)
|
||||
lm =nhtolm (ih, nt)
|
||||
do ig = 1, npw_
|
||||
vkb1 (ig,ih) = ylm (ig, lm) * vq (ig)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! vkb1 contains all betas including angular part for type nt
|
||||
! now add the structure factor and factor (-i)^l
|
||||
!
|
||||
do na = 1, nat
|
||||
! ordering: first all betas for atoms of type 1
|
||||
! then all betas for atoms of type 2 and so on
|
||||
if (ityp (na) .eq.nt) then
|
||||
!arg = (q_(1) * tau (1, na) + &
|
||||
! q_(2) * tau (2, na) + &
|
||||
! q_(3) * tau (3, na) ) * tpi
|
||||
!phase = CMPLX (cos (arg), - sin (arg) )
|
||||
do ig = 1, npw_
|
||||
sk (ig) = eigts1 (ig1(igk_(ig)), na) * &
|
||||
eigts2 (ig2(igk_(ig)), na) * &
|
||||
eigts3 (ig3(igk_(ig)), na)
|
||||
enddo
|
||||
do ih = 1, nh (nt)
|
||||
jkb = jkb + 1
|
||||
pref = (0.d0, -1.d0) **nhtol (ih, nt)! * phase
|
||||
do ig = 1, npw_
|
||||
vkb_(ig, jkb) = vkb1 (ig,ih) * sk (ig) * pref
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
deallocate (gk)
|
||||
deallocate (ylm)
|
||||
deallocate (vq)
|
||||
deallocate (qg)
|
||||
deallocate (sk)
|
||||
deallocate (vkb1)
|
||||
|
||||
call stop_clock ('init_us_2')
|
||||
return
|
||||
end subroutine init_us_2_no_phase
|
||||
|
||||
|
|
@ -0,0 +1,93 @@
|
|||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 j_para(fact, psi_n, psi_m, ik, q, j)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Compute the paramgnetic current between two states:
|
||||
! ... j_para(r') = (1/2) fact <psi_n| { p_k|r'><r'| + |r'><r'|p_{k+q} } |psi_m>
|
||||
! ... the result is added to j and is returned in real space
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE klist, ONLY : xk
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk, wg
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- parameters ---------------------------------------------------------
|
||||
IMPLICIT none
|
||||
INTEGER, INTENT(IN) :: ik ! k-point
|
||||
REAL(DP), INTENT(IN) :: fact ! multiplication factor
|
||||
REAL(DP), INTENT(IN) :: q(3)
|
||||
COMPLEX(DP), INTENT(IN) :: psi_n(npwx,nbnd), psi_m(npwx,nbnd)
|
||||
REAL(DP), INTENT(INOUT) :: j(nrxxs,3)
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
COMPLEX(DP), allocatable :: p_psic(:), psic(:), aux(:)
|
||||
REAL(DP) :: gk
|
||||
INTEGER :: ig, ipol, ibnd
|
||||
|
||||
call start_clock('j_para')
|
||||
|
||||
! allocate real space wavefunctions
|
||||
allocate(p_psic(nrxxs), psic(nrxxs), aux(npwx))
|
||||
|
||||
! loop over cartesian components
|
||||
do ipol = 1, 3
|
||||
|
||||
! loop over bands
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
|
||||
! apply p_k on the left
|
||||
do ig = 1, npw
|
||||
gk = xk(ipol,ik) + g(ipol,igk(ig))
|
||||
aux(ig) = gk * tpiba * psi_n(ig,ibnd)
|
||||
enddo
|
||||
|
||||
! transform to real space
|
||||
p_psic(:) = (0.d0,0.d0)
|
||||
p_psic(nls(igk(1:npw))) = aux(1:npw)
|
||||
call cft3s(p_psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
|
||||
|
||||
psic(:) = (0.d0,0.d0)
|
||||
psic(nls(igk(1:npw))) = psi_m(1:npw,ibnd)
|
||||
call cft3s(psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
|
||||
|
||||
! add to the current
|
||||
j(1:nrxxs,ipol) = j(1:nrxxs,ipol) + 0.5d0 * fact * wg(ibnd,ik) * &
|
||||
aimag(conjg(p_psic(1:nrxxs)) * psic(1:nrxxs))
|
||||
|
||||
! apply p_{k+q} on the right
|
||||
do ig = 1, npw
|
||||
gk = xk(ipol,ik) + g(ipol,igk(ig)) + q(ipol)
|
||||
aux(ig) = gk * tpiba * psi_m(ig,ibnd)
|
||||
enddo
|
||||
|
||||
! transform to real space
|
||||
p_psic(:) = (0.d0,0.d0)
|
||||
p_psic(nls(igk(1:npw))) = aux(1:npw)
|
||||
call cft3s(p_psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
|
||||
|
||||
psic(:) = (0.d0,0.d0)
|
||||
psic(nls(igk(1:npw))) = psi_n(1:npw,ibnd)
|
||||
call cft3s(psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
|
||||
|
||||
! add to the current
|
||||
j(1:nrxxs,ipol) = j(1:nrxxs,ipol) + 0.5d0 * fact * wg(ibnd,ik) * &
|
||||
aimag(conjg(psic(1:nrxxs)) * p_psic(1:nrxxs))
|
||||
|
||||
enddo ! ibnd
|
||||
enddo ! ipol
|
||||
|
||||
! free memory
|
||||
deallocate(p_psic, psic, aux)
|
||||
|
||||
call stop_clock('j_para')
|
||||
|
||||
END SUBROUTINE j_para
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,102 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
PROGRAM magn_main
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... This is the main driver of the magnetic resposne program.
|
||||
! ... It controls the initialization routines.
|
||||
! ... Features: NMR chemical shifts
|
||||
!... EPR g-tensor
|
||||
! ... Ported to Espresso by:
|
||||
! ... D. Ceresoli, A. P. Seitsonen, U. Gerstamnn and F. Mauri
|
||||
! ...
|
||||
! ... References (NMR):
|
||||
! ... F. Mauri and S. G. Louie Phys. Rev. Lett. 76, 4246 (1996)
|
||||
! ... F. Mauri, B. G. Pfrommer, S. G. Louie, Phys. Rev. Lett. 77, 5300 (1996)
|
||||
! ... T. Gregor, F. Mauri, and R. Car, J. Chem. Phys. 111, 1815 (1999)
|
||||
! ... C. J. Pickard and F. Mauri, Phys. Rev. B 63, 245101 (2001)
|
||||
! ... C. J. Pickard and F. Mauri, Phys. Rev. Lett. 91, 196401 (2003)
|
||||
! ...
|
||||
! ... References (g-tensor):
|
||||
! ... C. J. Pickard and F. Mauri, Phys. Rev. Lett. 88, 086403 (2002)
|
||||
! ...
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout, ionode, ionode_id
|
||||
USE wvfct, ONLY : gamma_only
|
||||
USE io_files, ONLY : prefix, tmp_dir, nd_nmbr
|
||||
USE klist, ONLY : nks
|
||||
USE mp, ONLY : mp_bcast
|
||||
USE cell_base, ONLY : tpiba
|
||||
USE nmr_module
|
||||
USE global_version, ONLY : version_number
|
||||
!<apsi>
|
||||
use paw, ONLY : read_recon, read_recon_paratec
|
||||
USE mp, ONLY : mp_barrier !*TMPTMPTMP
|
||||
!</apsi>
|
||||
!------------------------------------------------------------------------
|
||||
IMPLICIT NONE
|
||||
CHARACTER (LEN=9) :: code = 'MAGN'
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
CALL init_clocks( .TRUE. )
|
||||
|
||||
CALL start_clock( 'MAGN' )
|
||||
|
||||
! ... and begin with the initialization part
|
||||
CALL startup( nd_nmbr, code, version_number )
|
||||
CALL nmr_readin()
|
||||
|
||||
! read ground state wavefunctions
|
||||
call read_file
|
||||
call openfil
|
||||
|
||||
if (gamma_only) call errore('MAGN_main', 'gamma_only == .true.', 1)
|
||||
|
||||
!<apsi>
|
||||
! Read in qe format
|
||||
IF ( read_recon_in_paratec_fmt ) THEN
|
||||
! Read in paratec format
|
||||
CALL read_recon_paratec ( file_reconstruction )
|
||||
ELSE
|
||||
CALL read_recon ( file_reconstruction )
|
||||
END IF
|
||||
!</apsi>
|
||||
! write(0,*) "QQQ2: ", aephi(1,1)%psi(5), psphi(1,1)%psi(5)
|
||||
|
||||
CALL nmr_allocate()
|
||||
CALL nmr_setup()
|
||||
CALL nmr_summary()
|
||||
CALL nmr_openfil()
|
||||
CALL print_clock( 'MAGN' )
|
||||
|
||||
! convert q_nmr in units of tpiba
|
||||
q_nmr = q_nmr / tpiba
|
||||
|
||||
! calculation
|
||||
select case(trim(job))
|
||||
case('nmr')
|
||||
call suscept_crystal
|
||||
|
||||
case('g_tensor')
|
||||
call g_tensor_crystal
|
||||
|
||||
case('f-sum')
|
||||
call test_f_sum_rule
|
||||
|
||||
case default
|
||||
call errore('magn_main', 'wrong or undefined job in input', 1)
|
||||
end select
|
||||
|
||||
! print timings and stop the code
|
||||
CALL print_clock_nmr()
|
||||
CALL stop_code( .TRUE. )
|
||||
|
||||
STOP
|
||||
|
||||
END PROGRAM magn_main
|
||||
|
|
@ -0,0 +1,692 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
MODULE nmr_module
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... This module contains the variables used for NMR calculations
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE constants, ONLY : a0_to_cm => bohr_radius_cm
|
||||
USE parameters, ONLY: npk, ntypx, lmaxx
|
||||
IMPLICIT NONE
|
||||
SAVE
|
||||
|
||||
! speed of light in atomic units: c = 1/alpha
|
||||
REAL(DP), PARAMETER :: c = 137.03599911d0
|
||||
|
||||
! avogadro number
|
||||
REAL(DP), PARAMETER :: avogadro = 6.022142d23
|
||||
|
||||
! number of occupied bands at each k-point
|
||||
INTEGER :: nbnd_occ(npk)
|
||||
|
||||
! alpha shift of the projector on the valence wfcs
|
||||
REAL(DP) :: alpha_pv
|
||||
|
||||
! eigenvalues and eigenfunctions at k+q
|
||||
COMPLEX(DP), ALLOCATABLE :: evq(:,:)
|
||||
|
||||
! induced current (bare term) and induced magnetic field
|
||||
REAL(DP), ALLOCATABLE :: j_bare(:,:,:,:), b_ind_r(:,:,:)
|
||||
|
||||
! induced magnetic field in reciprocal space
|
||||
COMPLEX(DP), ALLOCATABLE :: b_ind(:,:,:)
|
||||
|
||||
! convergence threshold for diagonalizationa and greenfunction
|
||||
REAL(DP) :: conv_threshold
|
||||
|
||||
! q for the perturbation (in bohrradius^{-1})
|
||||
REAL(DP) :: q_nmr
|
||||
|
||||
! verbosity
|
||||
INTEGER :: iverbosity
|
||||
|
||||
! job: nmr, g_tensor
|
||||
CHARACTER(80) :: job
|
||||
|
||||
! format for a rank-2 tensor
|
||||
CHARACTER(*), PARAMETER :: tens_fmt = '(3(5X,3(F14.4,2X)/))'
|
||||
|
||||
! for plotting the induced current and induced field
|
||||
CHARACTER(80) :: filcurr, filfield
|
||||
|
||||
!<apsi>
|
||||
CHARACTER(256) :: file_reconstruction ( ntypx )
|
||||
LOGICAL :: read_recon_in_paratec_fmt
|
||||
REAL(dp) :: rc(ntypx,0:lmaxx)
|
||||
COMPLEX(dp), ALLOCATABLE :: paw_becp2 ( :, : )
|
||||
REAL(dp), ALLOCATABLE, DIMENSION ( :, : ) :: lx, ly, lz
|
||||
REAL(dp), ALLOCATABLE :: radial_integral_paramagnetic(:,:,:)
|
||||
REAL(dp), ALLOCATABLE :: radial_integral_diamagnetic(:,:,:)
|
||||
REAL(dp), ALLOCATABLE :: radial_integral_paramagnetic_so(:,:,:)
|
||||
REAL(dp), ALLOCATABLE :: radial_integral_diamagnetic_so(:,:,:)
|
||||
REAL(dp), ALLOCATABLE :: radial_integral_rmc(:,:,:)
|
||||
!<apsi>
|
||||
|
||||
CONTAINS
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! Read in the nmr input file.
|
||||
! Format: &inputmagn
|
||||
! prefix = '...' prefix of SCF calculation
|
||||
! tmp_dir = '...' scratch directory
|
||||
! job = 'nmr' or 'g_tensor'
|
||||
! conv_threshold = 1d-14
|
||||
! q_nmr = 0.01d0
|
||||
! filcurr = '...'
|
||||
! filfield = '...'
|
||||
! iverbosity = 0
|
||||
! /
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE nmr_readin()
|
||||
USE io_files, ONLY : nd_nmbr, prefix, tmp_dir
|
||||
USE io_global, ONLY : ionode
|
||||
IMPLICIT NONE
|
||||
INTEGER :: ios
|
||||
NAMELIST /inputmagn/ job, prefix, tmp_dir, conv_threshold, &
|
||||
q_nmr, iverbosity, filcurr, filfield, &
|
||||
read_recon_in_paratec_fmt, &
|
||||
file_reconstruction
|
||||
|
||||
if ( .not. ionode ) goto 400
|
||||
|
||||
job = ''
|
||||
prefix = 'pwscf'
|
||||
tmp_dir = './scratch/'
|
||||
conv_threshold = 1d-14
|
||||
q_nmr = 0.01d0
|
||||
iverbosity = 0
|
||||
filcurr = ''
|
||||
filfield = ''
|
||||
read_recon_in_paratec_fmt = .FALSE.
|
||||
file_reconstruction ( : ) = " "
|
||||
|
||||
read( 5, inputmagn, err = 200, iostat = ios )
|
||||
|
||||
200 call errore( 'nmr_readin', 'reading inputmagn namelist', abs( ios ) )
|
||||
|
||||
400 continue
|
||||
|
||||
call nmr_bcast_input
|
||||
|
||||
END SUBROUTINE nmr_readin
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! Broadcast input data to all processors
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE nmr_bcast_input
|
||||
#ifdef __PARA
|
||||
#include "f_defs.h"
|
||||
USE mp, ONLY : mp_bcast
|
||||
USE io_files, ONLY : prefix, tmp_dir
|
||||
implicit none
|
||||
integer :: root = 0
|
||||
call mp_bcast(job, root)
|
||||
call mp_bcast(prefix, root)
|
||||
call mp_bcast(tmp_dir, root)
|
||||
call mp_bcast(conv_threshold, root)
|
||||
call mp_bcast(q_nmr, root)
|
||||
call mp_bcast(iverbosity, root)
|
||||
call mp_bcast(filcurr, root)
|
||||
call mp_bcast(filfield, root)
|
||||
call mp_bcast(read_recon_in_paratec_fmt, root)
|
||||
call mp_bcast(file_reconstruction, root)
|
||||
#endif
|
||||
END SUBROUTINE nmr_bcast_input
|
||||
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! Allocate memory for NMR
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE nmr_allocate
|
||||
USE becmod, ONLY : becp
|
||||
USE lsda_mod, ONLY : nspin, lsda
|
||||
USE pwcom
|
||||
|
||||
!<apsi>
|
||||
USE paw, only: paw_vkb, paw_becp, paw_nkb
|
||||
!</apsi>
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
allocate(evq(npwx,nbnd))
|
||||
allocate(becp(nkb,nbnd))
|
||||
allocate(j_bare(nrxxs,3,3,nspin), b_ind_r(nrxxs,3,3), b_ind(ngm,3,3))
|
||||
|
||||
END SUBROUTINE nmr_allocate
|
||||
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! Print a short summary of the calculation
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE nmr_summary
|
||||
USE io_global, ONLY : stdout
|
||||
IMPLICIT NONE
|
||||
|
||||
CALL flush_unit( stdout )
|
||||
END SUBROUTINE nmr_summary
|
||||
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! Open files needed for NMR
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE nmr_openfil
|
||||
USE wvfct, ONLY : npwx
|
||||
USE uspp, ONLY : nkb
|
||||
IMPLICIT NONE
|
||||
|
||||
END SUBROUTINE nmr_openfil
|
||||
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! Print timings
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE print_clock_nmr
|
||||
USE io_global, ONLY : stdout
|
||||
IMPLICIT NONE
|
||||
|
||||
WRITE( stdout, * )
|
||||
call print_clock ('MAGN')
|
||||
WRITE( stdout, * ) ' INITIALIZATION: '
|
||||
call print_clock ('nmr_setup')
|
||||
WRITE( stdout, * )
|
||||
call print_clock ('greenf')
|
||||
call print_clock ('cgsolve')
|
||||
call print_clock ('ch_psi')
|
||||
call print_clock ('h_psiq')
|
||||
WRITE( stdout, * )
|
||||
call print_clock ('u_kq')
|
||||
call print_clock ('h_psi')
|
||||
WRITE( stdout, * )
|
||||
call print_clock ('apply_p')
|
||||
call print_clock ('apply_vel')
|
||||
WRITE( stdout, * )
|
||||
call print_clock ('j_para')
|
||||
call print_clock ('biot_savart')
|
||||
call print_clock ('c_sigma')
|
||||
WRITE( stdout, * )
|
||||
WRITE( stdout, * ) ' General routines'
|
||||
call print_clock ('ccalbec')
|
||||
call print_clock ('cft3')
|
||||
call print_clock ('cft3s')
|
||||
call print_clock ('cinterpolate')
|
||||
call print_clock ('davcio')
|
||||
call print_clock ('write_rec')
|
||||
WRITE( stdout, * )
|
||||
#ifdef __PARA
|
||||
WRITE( stdout, * ) ' Parallel routines'
|
||||
call print_clock ('reduce')
|
||||
call print_clock ('poolreduce')
|
||||
#endif
|
||||
END SUBROUTINE print_clock_nmr
|
||||
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! NMR setup
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE nmr_setup
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE ions_base, ONLY : tau, nat, ntyp => nsp
|
||||
USE atom, ONLY : nlcc
|
||||
USE wvfct, ONLY : nbnd, et, wg, npwx
|
||||
USE lsda_mod, ONLY : nspin, lsda
|
||||
USE scf, ONLY : vr, vrs, vltot, rho, rho_core
|
||||
USE gvect, ONLY : nrxx, ngm
|
||||
USE gsmooth, ONLY : doublegrid
|
||||
USE klist, ONLY : xk, degauss, ngauss, nks, nelec
|
||||
USE constants, ONLY : degspin, pi
|
||||
!<apsi>
|
||||
USE paw, ONLY : paw_nbeta, aephi, psphi, paw_vkb, &
|
||||
paw_becp, paw_nkb, vloc_present, &
|
||||
gipaw_ae_vloc, gipaw_ps_vloc
|
||||
USE atom, ONLY : r, rab
|
||||
!</apsi>
|
||||
|
||||
IMPLICIT none
|
||||
integer :: ik, nt, ibnd
|
||||
logical :: nlcc_any
|
||||
real(dp) :: emin, emax
|
||||
|
||||
!<apsi>
|
||||
integer :: il, lm, l, m, lm1, lm2, m1, m2, abs_m1, abs_m2
|
||||
integer :: sign_m1, sign_m2, il1, il2, l1, l2, j, kkpsi, nrc
|
||||
real(dp) :: alpha_lm, beta_lm
|
||||
integer, allocatable :: lm2l ( : ), lm2m ( : ), kkpsi_max
|
||||
real(dp), allocatable :: work(:), kinetic_aephi(:), kinetic_psphi(:)
|
||||
real(dp), allocatable :: aephi_dvloc_dr(:), psphi_dvloc_dr(:)
|
||||
!</apsi>
|
||||
|
||||
real(dp) :: mysum1 ( 3, lmaxx ) !TMPTMPTMP
|
||||
real(dp) :: mysum2 ( 3, 1:lmaxx ) !TMPTMPTMP
|
||||
|
||||
call start_clock ('nmr_setup')
|
||||
|
||||
! initialize pseudopotentials
|
||||
call init_us_1
|
||||
|
||||
!<apsi>
|
||||
! initialize paw
|
||||
do nt=1,ntyp
|
||||
do il = 1,paw_nbeta(nt)
|
||||
IF ( psphi(nt,il)%label%rc < -0.99_dp ) THEN
|
||||
rc(nt,psphi(nt,il)%label%l) = 1.6_dp
|
||||
rc(nt,aephi(nt,il)%label%l) = 1.6_dp
|
||||
psphi(nt,il)%label%rc = rc(nt,psphi(nt,il)%label%l)
|
||||
aephi(nt,il)%label%rc = rc(nt,aephi(nt,il)%label%l)
|
||||
ELSE
|
||||
rc(nt,psphi(nt,il)%label%l) = psphi(nt,il)%label%rc
|
||||
rc(nt,aephi(nt,il)%label%l) = aephi(nt,il)%label%rc
|
||||
END IF
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call init_paw_1
|
||||
|
||||
allocate (paw_vkb( npwx, paw_nkb))
|
||||
allocate (paw_becp(paw_nkb, nbnd))
|
||||
allocate (paw_becp2(paw_nkb, nbnd))
|
||||
|
||||
!<apsi>
|
||||
allocate ( radial_integral_diamagnetic(paw_nkb,paw_nkb,ntypx) )
|
||||
allocate ( radial_integral_paramagnetic(paw_nkb,paw_nkb,ntypx) )
|
||||
allocate ( radial_integral_diamagnetic_so(paw_nkb,paw_nkb,ntypx) )
|
||||
allocate ( radial_integral_paramagnetic_so(paw_nkb,paw_nkb,ntypx) )
|
||||
allocate ( radial_integral_rmc(paw_nkb,paw_nkb,ntypx) )
|
||||
radial_integral_diamagnetic = 0.0_dp
|
||||
radial_integral_paramagnetic = 0.0_dp
|
||||
radial_integral_diamagnetic_so = 0.0_dp
|
||||
radial_integral_paramagnetic_so = 0.0_dp
|
||||
radial_integral_rmc = 0.0_dp
|
||||
do nt=1, ntyp
|
||||
|
||||
do il1=1, paw_nbeta(nt)
|
||||
l1 = psphi(nt,il1)%label%l
|
||||
kkpsi = aephi(nt,il1)%kkpsi
|
||||
|
||||
nrc = psphi(nt,il1)%label%nrc
|
||||
write(0,*) "ZZZtmptmptmp: ", kkpsi, nrc
|
||||
|
||||
allocate ( work(kkpsi) )
|
||||
|
||||
do il2 = 1, paw_nbeta(nt)
|
||||
l2 = psphi(nt,il2)%label%l
|
||||
|
||||
IF ( l1 /= l2 ) CYCLE
|
||||
|
||||
!
|
||||
! NMR shielding, diamagnetic
|
||||
!
|
||||
!work = 0.0_dp
|
||||
do j = 1, nrc
|
||||
work(j) = (aephi(nt,il1)%psi(j)*aephi(nt,il2)%psi(j)-&
|
||||
psphi(nt,il1)%psi(j)*psphi(nt,il2)%psi(j))/r(j,nt)
|
||||
enddo
|
||||
|
||||
!work(1) = 0.0_dp
|
||||
CALL simpson( nrc, work, rab(:nrc,nt), &
|
||||
radial_integral_diamagnetic(il1,il2,nt) )
|
||||
|
||||
!
|
||||
! NMR shielding, paramagnetic
|
||||
!
|
||||
! calculate radial integration on atom site
|
||||
! <aephi|1/r^3|aephi>-<psphi|1/r^3|psphi>
|
||||
!
|
||||
!work(1) = 0.0_dp
|
||||
do j = 1, nrc
|
||||
work(j) = &
|
||||
( aephi(nt,il1)%psi(j) * aephi(nt,il2)%psi(j) &
|
||||
- psphi(nt,il1)%psi(j) * psphi(nt,il2)%psi(j) ) &
|
||||
/ r(j,nt) ** 3
|
||||
end do
|
||||
|
||||
!work(1) = 0.0_dp
|
||||
call simpson( nrc, work, rab(:,nt), &
|
||||
radial_integral_paramagnetic(il1,il2,nt) )
|
||||
write(stdout,*) "WWW1: ", l2, il1, il2, &
|
||||
radial_integral_paramagnetic(il1,il2,nt) / c**2 * 1e6 * 4
|
||||
|
||||
! Calculate the radial integral only if the radial potential
|
||||
! is present
|
||||
IF ( .NOT. vloc_present ( nt ) ) CYCLE
|
||||
|
||||
!
|
||||
! g tensor, relativistic mass correction
|
||||
!
|
||||
ALLOCATE ( kinetic_aephi ( kkpsi ), kinetic_psphi ( kkpsi ) )
|
||||
CALL radial_kinetic_energy ( l2, r(:nrc,nt), &
|
||||
aephi(nt,il2)%psi(:nrc), kinetic_aephi(:nrc) )
|
||||
CALL radial_kinetic_energy ( l2, r(:nrc,nt), &
|
||||
psphi(nt,il2)%psi(:nrc), kinetic_psphi(:nrc) )
|
||||
|
||||
do j = 1, nrc
|
||||
work(j) = ( aephi(nt,il1)%psi(j) * kinetic_aephi(j) &
|
||||
- psphi(nt,il1)%psi(j) * kinetic_psphi(j) )
|
||||
end do
|
||||
DEALLOCATE ( kinetic_aephi, kinetic_psphi )
|
||||
|
||||
!work(1) = 0.0_dp
|
||||
CALL simpson ( nrc, work, rab(:,nt), &
|
||||
radial_integral_rmc(il1,il2,nt) )
|
||||
write(stdout,*) "WWW2: ", l2, il1, il2, &
|
||||
radial_integral_rmc(il1,il2,nt)
|
||||
|
||||
ALLOCATE ( aephi_dvloc_dr ( nrc ), psphi_dvloc_dr ( nrc ) )
|
||||
|
||||
CALL radial_derivative ( r(:nrc,nt), &
|
||||
gipaw_ae_vloc ( :nrc, nt ), &
|
||||
aephi_dvloc_dr ( :nrc ) )
|
||||
CALL radial_derivative ( r(:nrc,nt), &
|
||||
gipaw_ps_vloc ( :nrc, nt ), &
|
||||
psphi_dvloc_dr ( :nrc ) )
|
||||
|
||||
!
|
||||
! g tensor, diamagnetic
|
||||
!
|
||||
do j = 1, nrc
|
||||
work(j) = ( aephi(nt,il1)%psi(j) * aephi_dvloc_dr(j) &
|
||||
* aephi(nt,il2)%psi(j) - psphi(nt,il1)%psi(j) &
|
||||
* psphi_dvloc_dr(j) * psphi(nt,il2)%psi(j) ) &
|
||||
* r(j,nt)
|
||||
end do
|
||||
|
||||
!work(1) = 0.0_dp
|
||||
call simpson( nrc, work, rab(:,nt), &
|
||||
radial_integral_diamagnetic_so(il1,il2,nt) )
|
||||
write(stdout,*) "WWW3: ", l2, il1, il2, &
|
||||
radial_integral_diamagnetic_so(il1,il2,nt) / c
|
||||
|
||||
!
|
||||
! g tensor, paramagnetic
|
||||
!
|
||||
!work(1) = 0.0_dp
|
||||
do j = 1, nrc
|
||||
work(j) = ( aephi(nt,il1)%psi(j) * aephi_dvloc_dr(j) &
|
||||
* aephi(nt,il2)%psi(j) - psphi(nt,il1)%psi(j) &
|
||||
* psphi_dvloc_dr(j) * psphi(nt,il2)%psi(j) ) &
|
||||
/ r(j,nt)
|
||||
end do
|
||||
if ( l1 == 0 ) then
|
||||
do j = 1, nrc
|
||||
write(90,*) r(j,nt), work(j)*r(j,nt)**2
|
||||
end do
|
||||
write(90,*) ""
|
||||
end if
|
||||
|
||||
!work(1) = 0.0_dp
|
||||
call simpson( nrc,work,rab(:,nt), &
|
||||
radial_integral_paramagnetic_so(il1,il2,nt) )
|
||||
write(stdout,*) "WWW4: ", l2, il1, il2, &
|
||||
radial_integral_paramagnetic_so(il1,il2,nt) / c
|
||||
|
||||
DEALLOCATE ( aephi_dvloc_dr, psphi_dvloc_dr )
|
||||
|
||||
enddo
|
||||
|
||||
DEALLOCATE ( work )
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!</apsi>
|
||||
|
||||
! terms for paramagnetic
|
||||
|
||||
!<apsi>
|
||||
allocate ( lx ( lmaxx**2, lmaxx**2 ) )
|
||||
allocate ( ly ( lmaxx**2, lmaxx**2 ) )
|
||||
allocate ( lz ( lmaxx**2, lmaxx**2 ) )
|
||||
|
||||
allocate ( lm2l ( lmaxx**2 ), lm2m ( lmaxx**2 ) )
|
||||
|
||||
lm = 0
|
||||
do l = 0, lmaxx - 1
|
||||
do m = 0, l
|
||||
lm = lm + 1
|
||||
lm2l ( lm ) = l
|
||||
lm2m ( lm ) = m
|
||||
if ( m /= 0 ) then
|
||||
lm = lm + 1
|
||||
lm2l ( lm ) = l
|
||||
lm2m ( lm ) = - m
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
lx = 0.0_dp
|
||||
ly = 0.0_dp
|
||||
lz = 0.0_dp
|
||||
do lm2 = 1, lmaxx**2
|
||||
do lm1 = 1, lmaxx**2
|
||||
if ( lm2l ( lm1 ) /= lm2l ( lm2 ) ) cycle
|
||||
|
||||
l = lm2l ( lm1 )
|
||||
|
||||
m1 = lm2m ( lm1 )
|
||||
m2 = lm2m ( lm2 )
|
||||
|
||||
! L_x, L_y
|
||||
if ( m2 == 0 ) then
|
||||
if ( m1 == -1 ) then
|
||||
lx ( lm1, lm2 ) = - sqrt(real(l*(l+1),dp)) / sqrt(2.0_dp)
|
||||
else if ( m1 == +1 ) then
|
||||
ly ( lm1, lm2 ) = + sqrt(real(l*(l+1),dp)) / sqrt(2.0_dp)
|
||||
end if
|
||||
else if ( m1 == 0 ) then
|
||||
if ( m2 == -1 ) then
|
||||
lx ( lm1, lm2 ) = + sqrt(real(l*(l+1),dp)) / sqrt(2.0_dp)
|
||||
else if ( m2 == +1 ) then
|
||||
ly ( lm1, lm2 ) = - sqrt(real(l*(l+1),dp)) / sqrt(2.0_dp)
|
||||
end if
|
||||
else
|
||||
abs_m1 = abs ( m1 )
|
||||
abs_m2 = abs ( m2 )
|
||||
sign_m1 = sign ( 1, m1 )
|
||||
sign_m2 = sign ( 1, m2 )
|
||||
alpha_lm = sqrt(real(l*(l+1)-abs_m2*(abs_m2+1),dp))
|
||||
beta_lm = sqrt(real(l*(l+1)-abs_m2*(abs_m2-1),dp))
|
||||
if ( abs_m1 == abs_m2 + 1 ) then
|
||||
lx ( lm1, lm2 ) =-( sign_m2 - sign_m1 ) * alpha_lm / 4.0_dp
|
||||
ly ( lm1, lm2 ) = ( sign_m2 + sign_m1 ) * alpha_lm / 4.0_dp &
|
||||
/ sign_m2
|
||||
else if ( abs_m1 == abs_m2 - 1 ) then
|
||||
lx ( lm1, lm2 ) =-( sign_m2 - sign_m1 ) * beta_lm / 4.0_dp
|
||||
ly ( lm1, lm2 ) =-( sign_m2 + sign_m1 ) * beta_lm / 4.0_dp &
|
||||
/ sign_m2
|
||||
end if
|
||||
end if
|
||||
|
||||
! L_z
|
||||
if ( m1 == - m2 ) then
|
||||
lz ( lm1, lm2 ) = - m2
|
||||
end if
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
write(stdout,'(A)') "lx:"
|
||||
write(stdout,'(9F8.5)') lx
|
||||
|
||||
write(stdout,'(A)') "ly:"
|
||||
write(stdout,'(9F8.5)') ly
|
||||
|
||||
write(stdout,'(A)') "lz:"
|
||||
write(stdout,'(9F8.5)') lz
|
||||
|
||||
! Checks
|
||||
mysum1 = 0
|
||||
mysum2 = 0
|
||||
do lm2 = 1, lmaxx**2
|
||||
do lm1 = 1, lmaxx**2
|
||||
if ( lm2l ( lm1 ) /= lm2l ( lm2 ) ) cycle
|
||||
l = lm2l ( lm2 )
|
||||
|
||||
mysum1(1,l+1) = mysum1(1,l+1) + lx(lm1,lm2)
|
||||
mysum2(1,l+1) = mysum2(1,l+1) + lx(lm1,lm2)**2
|
||||
mysum1(2,l+1) = mysum1(2,l+1) + ly(lm1,lm2)
|
||||
mysum2(2,l+1) = mysum2(2,l+1) + ly(lm1,lm2)**2
|
||||
mysum1(3,l+1) = mysum1(3,l+1) + lz(lm1,lm2)
|
||||
mysum2(3,l+1) = mysum2(3,l+1) + lz(lm1,lm2)**2
|
||||
end do
|
||||
end do
|
||||
write(stdout,'(A,9F8.4)') "DDD1x: ", mysum1(1,:)
|
||||
write(stdout,'(A,9F8.4)') "DDD1y: ", mysum1(2,:)
|
||||
write(stdout,'(A,9F8.4)') "DDD1z: ", mysum1(3,:)
|
||||
write(stdout,'(A,9F8.4)') "DDD2x: ", mysum2(1,:)
|
||||
write(stdout,'(A,9F8.4)') "DDD2y: ", mysum2(2,:)
|
||||
write(stdout,'(A,9F8.4)') "DDD2z: ", mysum2(3,:)
|
||||
|
||||
deallocate ( lm2l, lm2m )
|
||||
|
||||
!
|
||||
! calculate radial integration on atom site
|
||||
! <aephi|1/r^3|aephi>-<psphi|1/r^3|psphi>
|
||||
!
|
||||
! allocate ( radial_integral_paramagnetic(paw_nkb,paw_nkb,ntypx) )
|
||||
! radial_integral_paramagnetic = 0.0_dp
|
||||
! do nt=1,ntyp
|
||||
! do il1=1,paw_nbeta(nt)
|
||||
! l1 = psphi(nt,il1)%label%l
|
||||
! kkpsi = aephi(nt,il1)%kkpsi
|
||||
! allocate ( work(kkpsi) )
|
||||
! nrc = psphi(nt,il1)%label%nrc
|
||||
! do il2=1,paw_nbeta(nt)
|
||||
! l2 = psphi(nt,il2)%label%l
|
||||
! IF ( l1 /= l2 ) CYCLE
|
||||
! work=0.0_dp
|
||||
! do j = 1, nrc
|
||||
! work(j) = (aephi(nt,il1)%psi(j)*aephi(nt,il2)%psi(j)-&
|
||||
! psphi(nt,il1)%psi(j)*psphi(nt,il2)%psi(j))/r(j,nt)**3
|
||||
! enddo
|
||||
!
|
||||
! !work(1) = 0.0_dp
|
||||
! call simpson( nrc,work,rab(:,nt), &
|
||||
! radial_integral_paramagnetic(il1,il2,nt) )
|
||||
!
|
||||
! if ( il1 == 1 .and. il2 == 1 ) then
|
||||
! write(stdout,*) "WWW: ", il1, il2, &
|
||||
! radial_integral_paramagnetic(il1,il2,nt) / c**2 * 1e6 * 4
|
||||
! end if
|
||||
!
|
||||
! enddo
|
||||
!
|
||||
! deallocate(work)
|
||||
! enddo
|
||||
! enddo
|
||||
!</apsi>
|
||||
! stop "DONE"
|
||||
|
||||
! computes the total local potential (external+scf) on the smooth grid
|
||||
call set_vrs (vrs, vltot, vr, nrxx, nspin, doublegrid)
|
||||
|
||||
! compute the D for the pseudopotentials
|
||||
call newd
|
||||
|
||||
! set non linear core correction stuff
|
||||
nlcc_any = .false.
|
||||
do nt = 1, ntyp
|
||||
nlcc_any = nlcc_any.or.nlcc (nt)
|
||||
enddo
|
||||
!!if (nlcc_any) allocate (drc( ngm, ntyp))
|
||||
|
||||
!! setup all gradient correction stuff
|
||||
!!call setup_dgc
|
||||
|
||||
! computes the number of occupied bands for each k point
|
||||
nbnd_occ (:) = 0
|
||||
do ik = 1, nks
|
||||
do ibnd = 1, nbnd
|
||||
if ( wg(ibnd,ik) > 1e-6 ) then
|
||||
nbnd_occ(ik) = ibnd
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
! computes alpha_pv
|
||||
emin = et (1, 1)
|
||||
do ik = 1, nks
|
||||
do ibnd = 1, nbnd
|
||||
emin = min (emin, et (ibnd, ik) )
|
||||
enddo
|
||||
enddo
|
||||
#ifdef __PARA
|
||||
! find the minimum across pools
|
||||
call poolextreme (emin, -1)
|
||||
#endif
|
||||
if (degauss.ne.0.d0) then
|
||||
call errore('nmr_setup', 'implemented only for insulators', -1)
|
||||
else
|
||||
emax = et (1, 1)
|
||||
do ik = 1, nks
|
||||
do ibnd = 1, nbnd
|
||||
emax = max (emax, et (ibnd, ik) )
|
||||
enddo
|
||||
enddo
|
||||
#ifdef __PARA
|
||||
! find the maximum across pools
|
||||
call poolextreme (emax, + 1)
|
||||
#endif
|
||||
alpha_pv = 2.d0 * (emax - emin)
|
||||
endif
|
||||
! avoid zero value for alpha_pv
|
||||
alpha_pv = max (alpha_pv, 1.0d-2)
|
||||
|
||||
call stop_clock('nmr_setup')
|
||||
|
||||
CONTAINS
|
||||
|
||||
SUBROUTINE radial_kinetic_energy ( l, rdata, ydata, kin_ydata )
|
||||
USE splinelib
|
||||
|
||||
INTEGER, INTENT ( IN ) :: l
|
||||
REAL(dp), INTENT ( IN ) :: rdata ( : ), ydata ( : )
|
||||
REAL(dp), INTENT ( OUT ) :: kin_ydata ( : )
|
||||
|
||||
REAL(dp) :: d1
|
||||
|
||||
d1 = ( ydata(2) - ydata(1) ) / ( rdata(2) - rdata(1) )
|
||||
CALL spline ( rdata, ydata, 0.d0, d1, kin_ydata )
|
||||
|
||||
kin_ydata = - kin_ydata + l*(l+1) * ydata / rdata ** 2
|
||||
|
||||
END SUBROUTINE radial_kinetic_energy
|
||||
|
||||
SUBROUTINE radial_derivative ( rdata, ydata, dydata_dr )
|
||||
USE splinelib
|
||||
|
||||
! ydata passed as y * r
|
||||
REAL(dp), INTENT ( IN ) :: rdata ( : ), ydata ( : )
|
||||
REAL(dp), INTENT ( OUT ) :: dydata_dr ( : )
|
||||
|
||||
INTEGER :: j
|
||||
REAL(dp) :: d1, tab_d2y ( SIZE ( ydata ) )
|
||||
|
||||
d1 = ( ydata(2) - ydata(1) ) / ( rdata(2) - rdata(1) )
|
||||
CALL spline ( rdata, ydata, 0.d0, d1, tab_d2y )
|
||||
|
||||
DO j = 1, SIZE ( ydata )
|
||||
dydata_dr ( j ) = &
|
||||
( splint_deriv ( rdata, ydata, tab_d2y, rdata ( j ) ) &
|
||||
- ydata ( j ) / rdata ( j ) ) / rdata ( j )
|
||||
END DO
|
||||
|
||||
END SUBROUTINE radial_derivative
|
||||
|
||||
END SUBROUTINE nmr_setup
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
END MODULE nmr_module
|
||||
!-----------------------------------------------------------------------
|
||||
|
|
@ -0,0 +1,46 @@
|
|||
!
|
||||
! Copyright (C) 2001-2004 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 stop_code( flag )
|
||||
!----------------------------------------------------------------------------
|
||||
!
|
||||
! ... Synchronize processes before stopping.
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE mp, ONLY : mp_end, mp_barrier
|
||||
!
|
||||
USE parallel_include
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
LOGICAL :: flag
|
||||
!
|
||||
!
|
||||
CALL mp_barrier()
|
||||
!
|
||||
CALL mp_end()
|
||||
!
|
||||
#if defined (__T3E)
|
||||
!
|
||||
! ... set streambuffers off
|
||||
!
|
||||
CALL set_d_stream( 0 )
|
||||
!
|
||||
#endif
|
||||
!
|
||||
IF ( flag ) THEN
|
||||
!
|
||||
STOP
|
||||
!
|
||||
ELSE
|
||||
!
|
||||
STOP 1
|
||||
!
|
||||
ENDIF
|
||||
!
|
||||
END SUBROUTINE stop_code
|
||||
|
|
@ -0,0 +1,649 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE suscept_crystal
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Compute the "bare" susceptibility as in Eq.(64-65) of
|
||||
! ... PRB 63, 245101 (2001)
|
||||
! ... add more comments
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE io_files, ONLY : nwordwfc, iunwfc
|
||||
USE cell_base, ONLY : at, bg, omega, tpiba, tpiba2
|
||||
USE wavefunctions_module, ONLY : evc
|
||||
USE klist, ONLY : nks, nkstot, wk, xk, nelec
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk, wg, g2kin, &
|
||||
current_k
|
||||
USE lsda_mod, ONLY : current_spin, lsda, isk
|
||||
USE becmod, ONLY : becp
|
||||
USE symme, ONLY : nsym, s, ftau
|
||||
USE mp_global, ONLY : my_pool_id, npool, me_pool, root_pool
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!<apsi>
|
||||
USE paw, ONLY: paw_vkb, paw_becp, paw_nkb, aephi, &
|
||||
psphi, paw_nh, paw_nhtol, &
|
||||
paw_nhtom, paw_indv, paw_nbeta
|
||||
USE ions_base, ONLY : nat
|
||||
!</apsi>
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
IMPLICIT NONE
|
||||
complex(dp), allocatable, dimension(:,:,:) :: p_evc, vel_evc, g_vel_evc
|
||||
complex(dp), allocatable :: aux(:,:)
|
||||
|
||||
! Q tensor of eq. (65) (pGv => HH in Paratec, vGv => VV in Paratec)
|
||||
real(dp) :: q_pGv(3,3,-1:1), q_vGv(3,3,-1:1)
|
||||
|
||||
! F tensor of eq. (64)
|
||||
real(dp) :: f_pGv(3,3,-1:1), f_vGv(3,3,-1:1)
|
||||
|
||||
! chi_bare tensor of eq. (64)
|
||||
real(dp) :: chi_bare_pGv(3,3), chi_bare_vGv(3,3)
|
||||
|
||||
! f-sum rule
|
||||
real(dp) :: f_sum(3,3)
|
||||
|
||||
integer :: ia, ib, ik, ipol, jpol, i, ibnd, isign
|
||||
real(dp) :: tmp(3,3), q(3), braket, sigma_bare(3,3,nat)
|
||||
real(dp) :: diamagnetic_corr_tensor(3,3,nat)
|
||||
real(dp) :: paramagnetic_corr_tensor(3,3,nat)
|
||||
real(dp) :: sigma_diamagnetic(3,3,nat)
|
||||
real(dp) :: sigma_paramagnetic(3,3,nat)
|
||||
complex(dp), external :: ZDOTC
|
||||
!-----------------------------------------------------------------------
|
||||
|
||||
! allocate memory
|
||||
allocate ( p_evc(npwx,nbnd,3), vel_evc(npwx,nbnd,3), &
|
||||
aux(npwx,nbnd), g_vel_evc(npwx,nbnd,3) )
|
||||
|
||||
! zero the Q tensors
|
||||
q_pGv(:,:,:) = 0.0_dp
|
||||
q_vGv(:,:,:) = 0.0_dp
|
||||
|
||||
! zero the current and the field
|
||||
j_bare(:,:,:,:) = (0.0_dp,0.0_dp)
|
||||
b_ind(:,:,:) = (0.0_dp,0.0_dp)
|
||||
|
||||
sigma_diamagnetic = 0.0_dp
|
||||
sigma_paramagnetic = 0.0_dp
|
||||
|
||||
write(stdout, '(5X,''Computing the magnetic susceptibility'')')
|
||||
!====================================================================
|
||||
! loop over k-points
|
||||
!====================================================================
|
||||
do ik = 1, nks
|
||||
#ifdef __PARA
|
||||
if (me_pool == root_pool) &
|
||||
write(*, '(5X,''k-point #'',I5,'' of '',I5,6X,''pool #'',I3)') &
|
||||
ik, nks, my_pool_id+1
|
||||
#else
|
||||
write(stdout, '(5X,''k-point #'',I5,'' of '',I5)') ik, nks
|
||||
#endif
|
||||
current_k = ik
|
||||
current_spin = isk(ik)
|
||||
|
||||
! initialize at k-point k
|
||||
call gk_sort(xk(1,ik), ngm, g, ecutwfc/tpiba2, npw, igk, g2kin)
|
||||
g2kin(:) = g2kin(:) * tpiba2
|
||||
call init_us_2(npw,igk,xk(1,ik),vkb)
|
||||
|
||||
! read wfcs from file and compute becp
|
||||
call davcio (evc, nwordwfc, iunwfc, ik, -1)
|
||||
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
|
||||
|
||||
!<apsi>
|
||||
call init_paw_2_no_phase (npw, igk, xk (1, ik), paw_vkb)
|
||||
call ccalbec (paw_nkb, npwx, npw, nbnd, paw_becp, paw_vkb, evc)
|
||||
diamagnetic_corr_tensor = 0.0
|
||||
call diamagnetic_correction ( diamagnetic_corr_tensor )
|
||||
sigma_diamagnetic = sigma_diamagnetic + diamagnetic_corr_tensor
|
||||
!</apsi>
|
||||
|
||||
! this is the case q = 0 (like the case of the f-sum rule)
|
||||
q(:) = 0.0_dp
|
||||
!!!write(*,'(''q='',3(F12.4))') q
|
||||
|
||||
call compute_u_kq(ik, q)
|
||||
evc = evq
|
||||
|
||||
! compute p_k|evc>, v_k|evc> and G_k v_{k,k}|evc>
|
||||
call apply_operators
|
||||
|
||||
!------------------------------------------------------------------
|
||||
! f-sum rule (pGv term only)
|
||||
!------------------------------------------------------------------
|
||||
do ia = 1, 3
|
||||
do ib = 1, 3
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
braket = 2.0_dp*real(ZDOTC(npw, p_evc(1,ibnd,ia), 1, &
|
||||
g_vel_evc(1,ibnd,ib), 1), DP)
|
||||
f_sum(ia,ib) = f_sum(ia,ib) + wg(ibnd,ik) * braket
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!------------------------------------------------------------------
|
||||
! pGv and vGv contribution to chi_{bare}
|
||||
!------------------------------------------------------------------
|
||||
do i = 1, 3
|
||||
call add_to_tensor(q_pGv(:,:,0), p_evc, g_vel_evc)
|
||||
call add_to_tensor(q_vGv(:,:,0), vel_evc, g_vel_evc)
|
||||
enddo
|
||||
|
||||
!------------------------------------------------------------------
|
||||
! loop over -q and +q
|
||||
!------------------------------------------------------------------
|
||||
do isign = -1, 1, 2
|
||||
|
||||
! loop over cartesian directions
|
||||
do i = 1, 3
|
||||
! set the q vector
|
||||
q(:) = 0.0_dp
|
||||
q(i) = dble(isign) * q_nmr
|
||||
!!!write(*,'(''q='',3(F12.4))') q
|
||||
|
||||
! compute the wfcs at k+q
|
||||
call compute_u_kq(ik, q)
|
||||
|
||||
! compute p_k|evc>, v_k|evc> and G_{k+q} v_{k+q,k}|evc>
|
||||
call apply_operators
|
||||
|
||||
!<apsi>
|
||||
call init_paw_2_no_phase (npw, igk, xk(:,ik)+q(:), paw_vkb)
|
||||
call paramagnetic_correction ( paramagnetic_corr_tensor )
|
||||
call add_to_sigma_para ( paramagnetic_corr_tensor, sigma_paramagnetic )
|
||||
!</apsi>
|
||||
|
||||
! pGv and vGv contribution to chi_bare
|
||||
call add_to_tensor(q_pGv(:,:,isign), p_evc, g_vel_evc)
|
||||
call add_to_tensor(q_vGv(:,:,isign), vel_evc, g_vel_evc)
|
||||
|
||||
! now the j_bare term
|
||||
call add_to_current(j_bare(:,:,:,current_spin), evc, g_vel_evc)
|
||||
enddo ! i
|
||||
|
||||
enddo ! isign
|
||||
enddo ! ik
|
||||
|
||||
#ifdef __PARA
|
||||
! reduce over G-vectors
|
||||
call reduce(9, f_sum)
|
||||
call reduce(3*9, q_pGv)
|
||||
call reduce(3*9, q_vGv)
|
||||
#endif
|
||||
|
||||
#ifdef __PARA
|
||||
! reduce over k-points
|
||||
call poolreduce(9, f_sum)
|
||||
call poolreduce(3*9, q_pGv)
|
||||
call poolreduce(3*9, q_vGv)
|
||||
call poolreduce(nrxxs*nspin*9, j_bare)
|
||||
call poolreduce(9*nat, sigma_diamagnetic)
|
||||
call poolreduce(9*nat, sigma_paramagnetic)
|
||||
#endif
|
||||
|
||||
!====================================================================
|
||||
! print out results
|
||||
!====================================================================
|
||||
write(stdout,'(5X,''End of magnetic susceptibility calculation'')')
|
||||
write(stdout,*)
|
||||
|
||||
! f-sum rule
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(5X,''f-sum rule:'')')
|
||||
write(stdout, tens_fmt) f_sum
|
||||
endif
|
||||
call sym_cart_tensor(f_sum)
|
||||
write(stdout, '(5X,''f-sum rule (symmetrized):'')')
|
||||
write(stdout, tens_fmt) f_sum
|
||||
|
||||
! F_{ij} = (2 - \delta_{ij}) Q_{ij}
|
||||
do ipol = 1, 3
|
||||
do jpol = 1, 3
|
||||
f_pGv(ipol,jpol,:) = 2.0_dp*q_pGv(ipol,jpol,:)
|
||||
if (ipol == jpol) f_pGv(ipol,jpol,:) = q_pGv(ipol,jpol,:)
|
||||
|
||||
f_vGv(ipol,jpol,:) = 2.0_dp*q_vGv(ipol,jpol,:)
|
||||
if (ipol == jpol) f_vGv(ipol,jpol,:) = q_vGv(ipol,jpol,:)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! compute chi_bare both pGv and vGv terms
|
||||
chi_bare_pGv(:,:) = f_pGv(:,:,1) - 2.0_dp*f_pGv(:,:,0) + f_pGv(:,:,-1)
|
||||
chi_bare_pGv(:,:) = -0.5_dp * chi_bare_pGv(:,:) / (c * q_nmr * tpiba)**2
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(5X,''chi_bare pGv (HH) in paratec units:'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_pGv(:,:) * c**2
|
||||
endif
|
||||
call sym_cart_tensor(chi_bare_pGv)
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_pGv(:,:) * c**2
|
||||
endif
|
||||
|
||||
chi_bare_vGv(:,:) = f_vGv(:,:,1) - 2.0_dp*f_vGv(:,:,0) + f_vGv(:,:,-1)
|
||||
chi_bare_vGv(:,:) = -0.5_dp * chi_bare_vGv(:,:) / (c * q_nmr * tpiba)**2
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(5X,''chi_bare vGv (VV) in paratec units:'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_vGv(:,:) * c**2
|
||||
endif
|
||||
call sym_cart_tensor(chi_bare_vGv)
|
||||
if (iverbosity > 0) then
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') chi_bare_vGv(:,:) * c**2
|
||||
endif
|
||||
|
||||
! convert from atomic units to 10^{-6} cm^3 / mol
|
||||
tmp(:,:) = chi_bare_pGv(:,:) * 1e6_dp * a0_to_cm**3.0_dp * avogadro
|
||||
write(stdout, '(5X,''chi_bare pGv (HH) in 10^{-6} cm^3/mol:'')')
|
||||
write(stdout, tens_fmt) tmp(:,:)
|
||||
|
||||
tmp(:,:) = chi_bare_vGv(:,:) * 1e6_dp * a0_to_cm**3.0_dp * avogadro
|
||||
write(stdout, '(5X,''chi_bare vGv (VV) in 10^{-6} cm^3/mol:'')')
|
||||
write(stdout, tens_fmt) tmp(:,:)
|
||||
|
||||
!--------------------------------------------------------------------
|
||||
! now get the current, induced field and chemical shifts
|
||||
!--------------------------------------------------------------------
|
||||
chi_bare_pGv(:,:) = chi_bare_pGv(:,:) / omega
|
||||
j_bare(:,:,:,:) = j_bare(:,:,:,:) / (2.0_dp * q_nmr * tpiba * c * omega)
|
||||
|
||||
!nsym = 1
|
||||
! either you symmetrize the current ...
|
||||
do i = 1, nspin
|
||||
#ifdef __PARA
|
||||
call psymmetrize_field(j_bare(:,:,:,i),1)
|
||||
#else
|
||||
call symmetrize_field(j_bare(:,:,:,i),1)
|
||||
#endif
|
||||
enddo
|
||||
|
||||
! compute induced field
|
||||
do ipol = 1, 3
|
||||
call biot_savart(ipol)
|
||||
enddo
|
||||
|
||||
do i = 1, nspin
|
||||
if (trim(filcurr) /= '') &
|
||||
call write_tensor_field(filcurr, i, j_bare(1,1,1,i))
|
||||
enddo
|
||||
if (trim(filfield) /= '') &
|
||||
call write_tensor_field(filfield, 0, b_ind_r)
|
||||
|
||||
! ... or you symmetrize the induced field
|
||||
!call symmetrize_field(b_ind_r,0)
|
||||
!call field_to_reciprocal_space
|
||||
|
||||
! compute chemical shifts
|
||||
call compute_sigma_bare(chi_bare_pGv, sigma_bare)
|
||||
|
||||
call compute_sigma_diamagnetic( sigma_diamagnetic )
|
||||
|
||||
call compute_sigma_paramagnetic( sigma_paramagnetic )
|
||||
|
||||
deallocate(p_evc, vel_evc, aux, g_vel_evc, j_bare, b_ind)
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
||||
!====================================================================
|
||||
! compute p_k|evc>, v_k|evc> and G_k v_{k+q,k}|evc>
|
||||
!====================================================================
|
||||
SUBROUTINE apply_operators
|
||||
implicit none
|
||||
integer ipol
|
||||
|
||||
do ipol = 1, 3
|
||||
call apply_p(evc, p_evc(1,1,ipol), ik, ipol, q)
|
||||
call apply_vel(evc, vel_evc(1,1,ipol), ik, ipol, q)
|
||||
! necessary because aux is overwritten by subroutine greenfunction
|
||||
aux(:,:) = vel_evc(:,:,ipol)
|
||||
call greenfunction(ik, aux, g_vel_evc(1,1,ipol), q)
|
||||
enddo
|
||||
END SUBROUTINE apply_operators
|
||||
|
||||
|
||||
|
||||
!====================================================================
|
||||
! add contribution the Q tensors
|
||||
! Q_{\alpha,\beta} += <(e_i \times ul)_\alpha | (e_i \times ur)_\beta>
|
||||
!====================================================================
|
||||
SUBROUTINE add_to_tensor(qt, ul, ur)
|
||||
implicit none
|
||||
real(dp), intent(inout) :: qt(3,3)
|
||||
complex(dp), intent(in) :: ul(npwx,nbnd,3), ur(npwx,nbnd,3)
|
||||
real(dp) :: braket
|
||||
integer :: ibnd, ia, ib, comp_ia, comp_ib, ind(3,3), mult(3,3)
|
||||
|
||||
! index for the cross product
|
||||
ind(:,1) = (/ 1, 3, 2/); mult(:,1) = (/ 0,-1, 1 /)
|
||||
ind(:,2) = (/ 3, 2, 1/); mult(:,2) = (/ 1, 0,-1 /)
|
||||
ind(:,3) = (/ 2, 1, 3/); mult(:,3) = (/-1, 1, 0 /)
|
||||
|
||||
do ia = 1, 3 ! ia = alpha
|
||||
comp_ia = ind(ia,i)
|
||||
if (mult(ia,i) == 0) cycle
|
||||
|
||||
do ib = 1, 3 ! ib = beta
|
||||
comp_ib = ind(ib,i)
|
||||
if (mult(ib,i) == 0) cycle
|
||||
|
||||
do ibnd = 1, nbnd_occ(ik)
|
||||
braket = real(ZDOTC(npw, ul(1,ibnd,comp_ia), 1, &
|
||||
ur(1,ibnd,comp_ib), 1), DP)
|
||||
qt(ia,ib) = qt(ia,ib) + wg(ibnd,ik) * &
|
||||
braket * mult(ia,i) * mult(ib,i)
|
||||
enddo ! ibnd
|
||||
|
||||
enddo ! ib
|
||||
enddo ! ia
|
||||
!#ifdef __PARA
|
||||
! call reduce(9, qt)
|
||||
!#endif
|
||||
END SUBROUTINE add_to_tensor
|
||||
|
||||
|
||||
|
||||
!====================================================================
|
||||
! add contribution the the current
|
||||
! j(r)_{\alpha,\beta} += <ul|J(r)|(B\times e_i \cdot ur)>
|
||||
!====================================================================
|
||||
SUBROUTINE add_to_current(j, ul, ur)
|
||||
implicit none
|
||||
real(dp), intent(inout) :: j(nrxxs,3,3)
|
||||
complex(dp), intent(in) :: ul(npwx,nbnd), ur(npwx,nbnd,3)
|
||||
real(dp) :: braket, fact
|
||||
integer :: ibdir, icomp, ind(3,3), mult(3,3)
|
||||
|
||||
! index for the cross product
|
||||
ind(:,1) = (/ 1, 3, 2/); mult(:,1) = (/ 0,-1, 1 /)
|
||||
ind(:,2) = (/ 3, 2, 1/); mult(:,2) = (/ 1, 0,-1 /)
|
||||
ind(:,3) = (/ 2, 1, 3/); mult(:,3) = (/-1, 1, 0 /)
|
||||
|
||||
! loop over B direction
|
||||
do ibdir = 1, 3
|
||||
if (i == ibdir) cycle
|
||||
icomp = ind(ibdir, i)
|
||||
fact = real(mult(ibdir,i)*isign)
|
||||
call j_para(fact, ul(1,1), ur(1,1,icomp), ik, q, j(1,1,ibdir))
|
||||
enddo
|
||||
END SUBROUTINE add_to_current
|
||||
|
||||
|
||||
!====================================================================
|
||||
! ...
|
||||
!====================================================================
|
||||
SUBROUTINE diamagnetic_correction ( diamagnetic_tensor )
|
||||
|
||||
USE atom, ONLY : r, rab
|
||||
USE ions_base, ONLY : nat, ityp, ntyp => nsp
|
||||
USE nmr_module, ONLY : c
|
||||
|
||||
implicit none
|
||||
|
||||
! Arguments
|
||||
real(dp), intent(inout):: diamagnetic_tensor(3,3,nat)
|
||||
|
||||
integer :: l1, m1, lm1, l2, m2, lm2, ih, ikb, nbs1, jh, jkb, nbs2
|
||||
integer :: nt, ibnd, na, lm, nrc, ijkb0
|
||||
complex(dp) , allocatable :: efg_corr(:,:)
|
||||
complex(dp) :: bec_product
|
||||
|
||||
allocate (efg_corr(lmaxx**2,nat))
|
||||
efg_corr = 0.0_dp
|
||||
|
||||
! rc=1.60_dp
|
||||
! nrc=count(r(1:msh(1),1).le.rc)
|
||||
!
|
||||
! calculate radial integration on atom site
|
||||
! <aephi|1/r^3|aephi>-<psphi|1/r^3|psphi>
|
||||
!
|
||||
|
||||
!
|
||||
! calculation of the reconstruction part
|
||||
!
|
||||
|
||||
do ibnd = 1, nbnd
|
||||
ijkb0 = 0
|
||||
do nt = 1, ntyp
|
||||
do na = 1, nat
|
||||
if (ityp (na) .eq.nt) then
|
||||
do ih = 1, paw_nh (nt)
|
||||
ikb = ijkb0 + ih
|
||||
nbs1=paw_indv(ih,nt)
|
||||
l1=paw_nhtol(ih,nt)
|
||||
m1=paw_nhtom(ih,nt)
|
||||
lm1=m1+l1**2
|
||||
do jh = 1, paw_nh (nt)
|
||||
jkb = ijkb0 + jh
|
||||
nbs2=paw_indv(jh,nt)
|
||||
l2=paw_nhtol(jh,nt)
|
||||
m2=paw_nhtom(jh,nt)
|
||||
lm2=m2+l2**2
|
||||
|
||||
bec_product = paw_becp(jkb,ibnd) &
|
||||
* CONJG(paw_becp(ikb,ibnd))
|
||||
|
||||
!<apsi> s/non-trace-zero component
|
||||
! 2/3 to separate the non-trace vanishing component
|
||||
! 1/(2c^2) from the equation (59) in PM-PRB
|
||||
IF ( l1 == l2 .AND. m1 == m2 ) THEN
|
||||
diamagnetic_tensor(1,1,na) &
|
||||
= diamagnetic_tensor(1,1,na) &
|
||||
+ 2.0_dp / 3.0_dp * bec_product &
|
||||
* radial_integral_diamagnetic(nbs1,nbs2,nt) &
|
||||
* wg(ibnd,ik) / (2.0_dp*c**2)
|
||||
diamagnetic_tensor(2,2,na) &
|
||||
= diamagnetic_tensor(2,2,na) &
|
||||
+ 2.0_dp / 3.0_dp * bec_product &
|
||||
* radial_integral_diamagnetic(nbs1,nbs2,nt) &
|
||||
* wg(ibnd,ik) / (2.0_dp*c**2)
|
||||
diamagnetic_tensor(3,3,na) &
|
||||
= diamagnetic_tensor(3,3,na) &
|
||||
+ 2.0_dp / 3.0_dp * bec_product &
|
||||
* radial_integral_diamagnetic(nbs1,nbs2,nt) &
|
||||
* wg(ibnd,ik) / (2.0_dp*c**2)
|
||||
END IF
|
||||
|
||||
! 2/3 to separate the non-trace vanishing component
|
||||
do lm = 5, 9
|
||||
efg_corr(lm,na) = efg_corr(lm,na) &
|
||||
+ bec_product / 3.0_dp &
|
||||
* radial_integral_diamagnetic(nbs1,nbs2,nt) &
|
||||
* ap(lm,lm1,lm2) * wg(ibnd,ik) / (2.0_dp*c**2)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
ijkb0 = ijkb0 + paw_nh (nt)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! write(6,'("DDD1",5F14.8)') efg_corr(5:9,1)
|
||||
! write(6,'("DDD2",5F14.8)') efg_corr(5:9,2)
|
||||
|
||||
!
|
||||
! transform in cartesian coordinates
|
||||
!
|
||||
|
||||
efg_corr(5:9,:nat) = - sqrt(4.0_dp * pi/5.0_dp) * efg_corr(5:9,:nat)
|
||||
|
||||
diamagnetic_tensor(1,1,:) = diamagnetic_tensor(1,1,:) &
|
||||
+ sqrt(3.0_dp) * efg_corr(8,:) - efg_corr(5,:)
|
||||
diamagnetic_tensor(2,2,:) = diamagnetic_tensor(2,2,:) &
|
||||
- sqrt(3.0_dp) * efg_corr(8,:) - efg_corr(5,:)
|
||||
diamagnetic_tensor(3,3,:) = diamagnetic_tensor(3,3,:) &
|
||||
+ efg_corr(5,:) * 2.0_dp
|
||||
diamagnetic_tensor(1,2,:) = diamagnetic_tensor(1,2,:) &
|
||||
+ efg_corr(9,:) * sqrt(3.0_dp)
|
||||
diamagnetic_tensor(2,1,:) = diamagnetic_tensor(1,2,:)
|
||||
diamagnetic_tensor(1,3,:) = diamagnetic_tensor(1,3,:) &
|
||||
- efg_corr(6,:) * sqrt(3.0_dp)
|
||||
diamagnetic_tensor(3,1,:) = diamagnetic_tensor(1,3,:)
|
||||
diamagnetic_tensor(2,3,:) = diamagnetic_tensor(2,3,:) &
|
||||
- efg_corr(7,:) * sqrt(3.0_dp)
|
||||
diamagnetic_tensor(3,2,:) = diamagnetic_tensor(2,3,:)
|
||||
|
||||
! efg_corr(5,:) = 3z^2-1
|
||||
! efg_corr(6,:) = -xz
|
||||
! efg_corr(7,:) = -yz
|
||||
! efg_corr(8,:) = x^2-y^2
|
||||
! efg_corr(9,:) = xy
|
||||
|
||||
deallocate ( efg_corr )
|
||||
|
||||
END SUBROUTINE diamagnetic_correction
|
||||
|
||||
|
||||
!====================================================================
|
||||
! ...
|
||||
!====================================================================
|
||||
SUBROUTINE paramagnetic_correction ( paramagnetic_tensor )
|
||||
|
||||
USE ions_base, ONLY : nat, ityp, ntyp => nsp
|
||||
USE nmr_module, ONLY : c
|
||||
|
||||
implicit none
|
||||
|
||||
! Arguments
|
||||
real(dp), intent(inout):: paramagnetic_tensor(3,3,nat)
|
||||
|
||||
integer :: l1, m1, lm1, l2, m2, lm2, ih, ikb, nbs1, jh, jkb, nbs2
|
||||
integer :: nt, ibnd, na, lm, j, ijkb0, ipol
|
||||
complex(dp) :: bec_product
|
||||
complex(dp) , allocatable :: efg_corr(:,:)
|
||||
|
||||
integer, parameter :: ng_ = 27, lmax2_ = 16
|
||||
integer :: mg, i1, i2, i3
|
||||
real(DP) :: g_ (3, ng_), gg_ (ng_)
|
||||
real(DP) :: ylm_ (ng_,lmax2_)
|
||||
|
||||
!--------------------------------------------------------------------------
|
||||
|
||||
allocate (efg_corr(3,nat))
|
||||
|
||||
!
|
||||
! calculation of the reconstruction part
|
||||
!
|
||||
|
||||
do ipol = 1, 3
|
||||
|
||||
if ( ipol == i ) cycle !TESTTESTTEST
|
||||
|
||||
call ccalbec (paw_nkb, npwx, npw, nbnd, paw_becp2, paw_vkb, &
|
||||
g_vel_evc(1,1,ipol))
|
||||
|
||||
efg_corr = 0.0_dp
|
||||
|
||||
do ibnd = 1, nbnd
|
||||
ijkb0 = 0
|
||||
do nt = 1, ntyp
|
||||
do na = 1, nat
|
||||
|
||||
if (ityp (na) .eq.nt) then
|
||||
do ih = 1, paw_nh (nt)
|
||||
ikb = ijkb0 + ih
|
||||
nbs1=paw_indv(ih,nt)
|
||||
l1=paw_nhtol(ih,nt)
|
||||
m1=paw_nhtom(ih,nt)
|
||||
lm1=m1+l1**2
|
||||
|
||||
do jh = 1, paw_nh (nt)
|
||||
jkb = ijkb0 + jh
|
||||
nbs2=paw_indv(jh,nt)
|
||||
l2=paw_nhtol(jh,nt)
|
||||
m2=paw_nhtom(jh,nt)
|
||||
lm2=m2+l2**2
|
||||
|
||||
if ( l1 /= l2 ) cycle
|
||||
|
||||
bec_product = CONJG(paw_becp(ikb,ibnd)) &
|
||||
* paw_becp2(jkb,ibnd)
|
||||
|
||||
efg_corr(1,na) = efg_corr(1,na) &
|
||||
+ bec_product &
|
||||
* radial_integral_paramagnetic(nbs1,nbs2,nt) &
|
||||
* lx ( lm1, lm2 ) * wg(ibnd,ik) / c ** 2
|
||||
efg_corr(2,na) = efg_corr(2,na) &
|
||||
+ bec_product &
|
||||
* radial_integral_paramagnetic(nbs1,nbs2,nt) &
|
||||
* ly ( lm1, lm2 ) * wg(ibnd,ik) / c ** 2
|
||||
efg_corr(3,na) = efg_corr(3,na) &
|
||||
+ bec_product &
|
||||
* radial_integral_paramagnetic(nbs1,nbs2,nt) &
|
||||
* lz ( lm1, lm2 ) * wg(ibnd,ik) / c ** 2
|
||||
|
||||
enddo
|
||||
enddo
|
||||
ijkb0 = ijkb0 + paw_nh (nt)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
paramagnetic_tensor ( :, ipol, : ) = REAL ( efg_corr, dp )
|
||||
|
||||
write(6,'("DDD1",2I3,3(F16.7,2X)') &
|
||||
ipol, i*isign, REAL ( efg_corr(1:3,1) ) * 1e6
|
||||
!write(6,'("DDD2",2I3,3(F16.7,2X)') &
|
||||
! ipol, i*isign, REAL ( efg_corr(1:3,2) ) * 1e6
|
||||
|
||||
end do
|
||||
|
||||
deallocate(efg_corr)
|
||||
|
||||
END SUBROUTINE paramagnetic_correction
|
||||
|
||||
!====================================================================
|
||||
! ...
|
||||
!====================================================================
|
||||
SUBROUTINE add_to_sigma_para( paramagnetic_correction, sigma_paramagnetic )
|
||||
implicit none
|
||||
real(dp), intent(in) :: paramagnetic_correction(3,3,nat)
|
||||
real(dp), intent(inout) :: sigma_paramagnetic(3,3,nat)
|
||||
real(dp) :: fact
|
||||
integer :: ibdir, icomp, ipol, ind(3,3), mult(3,3)
|
||||
|
||||
! index for the cross product
|
||||
ind(:,1) = (/ 1, 3, 2/); mult(:,1) = (/ 0,-1, 1 /)
|
||||
ind(:,2) = (/ 3, 2, 1/); mult(:,2) = (/ 1, 0,-1 /)
|
||||
ind(:,3) = (/ 2, 1, 3/); mult(:,3) = (/-1, 1, 0 /)
|
||||
|
||||
! loop over B direction
|
||||
do ibdir = 1, 3
|
||||
if (i == ibdir) cycle
|
||||
icomp = ind(ibdir, i)
|
||||
fact = real(mult(ibdir,i)*isign)
|
||||
|
||||
do ipol = 1, 3
|
||||
sigma_paramagnetic ( ipol, icomp, : ) &
|
||||
= sigma_paramagnetic ( ipol, icomp, : ) &
|
||||
+ fact * paramagnetic_correction ( ipol, ibdir, : ) &
|
||||
/ ( 2 * q_nmr * tpiba )
|
||||
|
||||
! Do iq=1, 3 ! loop over all q-points
|
||||
! ...
|
||||
! Do iperm0 =1,-1,-2
|
||||
! b0 = Mod(iq+iperm0+2,3) +1 ! gives different b0 for iperm0
|
||||
! p0 = Mod(iq-iperm0+2,3) +1 ! gives p0 = abs(q x b0)
|
||||
! ...
|
||||
! Call take_nonloc_deriv_kq_k(p0,k_gspace, rk(1),u_k(1,i),&
|
||||
! ...
|
||||
! para_corr_shift_tot(b0,:,:,:) = para_corr_shift_tot(b0,:,:,:)-&
|
||||
! Real(iperm0,dp)*Real(iqsign,dp)*para_corr_shift*&
|
||||
! kpoints%w(irk)/qmag/dtwo/Real(crys%nspin,dp)
|
||||
|
||||
end do
|
||||
enddo
|
||||
END SUBROUTINE add_to_sigma_para
|
||||
|
||||
END SUBROUTINE suscept_crystal
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 sym_cart_tensor(tens)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... symmetrize a rank-2 tensor in cartesian coordinates
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE cell_base, ONLY : at, bg
|
||||
USE symme, ONLY : s, nsym
|
||||
|
||||
IMPLICIT NONE
|
||||
REAL(DP), INTENT(INOUT) :: tens(3,3)
|
||||
|
||||
! cartesian to crystal
|
||||
call trntns (tens, at, bg, -1)
|
||||
|
||||
! symmetrize
|
||||
call symtns (tens, nsym, s)
|
||||
|
||||
! crystal to cartesian
|
||||
call trntns (tens, at, bg, 1)
|
||||
|
||||
END SUBROUTINE sym_cart_tensor
|
||||
|
|
@ -0,0 +1,195 @@
|
|||
!
|
||||
! 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 .
|
||||
!
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE symmetrize_field(field, iflag)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! symmetrize a tensor field (e.g. current or induced magnetic field)
|
||||
! iflag = 0 => tensor (e.g. induced B field)
|
||||
! iflag = 1 => pseudo-tensor (e.g. induced current)
|
||||
!
|
||||
! don't use nrxxs: in the parallel case nrx1s*nrx2s*nrx3s /= nrxxs
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE symme, ONLY : s, nsym
|
||||
USE cell_base, ONLY : at, bg
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- parameters ------------------------------------------------------
|
||||
IMPLICIT NONE
|
||||
REAL(DP), INTENT(INOUT) :: field(nrx1s*nrx2s*nrx3s,3,3)
|
||||
INTEGER :: iflag
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
complex(dp), allocatable :: aux(:)
|
||||
real(dp) :: tmp(3,3)
|
||||
integer :: i, ipol, jpol
|
||||
|
||||
! if no symmetries, return
|
||||
if (nsym.eq.1) return
|
||||
|
||||
! cartesian to crystal
|
||||
do i = 1, nrx1s*nrx2s*nrx3s
|
||||
tmp(:,:) = field(i,:,:)
|
||||
call trntns (tmp, at, bg, -1)
|
||||
field(i,:,:) = tmp(:,:)
|
||||
enddo
|
||||
|
||||
! symmetrize
|
||||
call syme2(field, iflag)
|
||||
|
||||
! crystal to cartesian
|
||||
do i = 1, nrx1s*nrx2s*nrx3s
|
||||
tmp(:,:) = field(i,:,:)
|
||||
call trntns (tmp, at, bg, 1)
|
||||
field(i,:,:) = tmp(:,:)
|
||||
enddo
|
||||
|
||||
END SUBROUTINE symmetrize_field
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE psymmetrize_field(field, iflag)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! symmetrize a tensor field (e.g. current or induced magnetic field)
|
||||
! (parallel version)
|
||||
! iflag = 0 => tensor (e.g. induced B field)
|
||||
! iflag = 1 => pseudo-tensor (e.g. induced current)
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE pfft, ONLY : nxx
|
||||
USE mp_global, ONLY : me_pool
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- parameters ------------------------------------------------------
|
||||
IMPLICIT NONE
|
||||
REAL(DP), INTENT(INOUT) :: field(nxx,3,3)
|
||||
INTEGER :: iflag
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
real(dp), allocatable :: aux(:,:,:)
|
||||
integer :: i, j
|
||||
|
||||
! if no symmetries, return
|
||||
if (nsym.eq.1) return
|
||||
|
||||
allocate( aux(nrx1s*nrx2s*nrx3s,3,3) )
|
||||
do i = 1, 3
|
||||
do j = 1, 3
|
||||
call gather(field(1,i,j), aux(1,i,j))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if ( me_pool == 0 ) call symmetrize_field(aux, iflag)
|
||||
|
||||
do i = 1, 3
|
||||
do j = 1, 3
|
||||
call scatter(aux(1,i,j), field(1,i,j))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
deallocate(aux)
|
||||
END SUBROUTINE psymmetrize_field
|
||||
|
||||
|
||||
!
|
||||
! 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 .
|
||||
!
|
||||
!
|
||||
!#include "f_defs.h"
|
||||
!---------------------------------------------------------------------
|
||||
subroutine syme2 (dvsym, iflag)
|
||||
!-------------------------------------------------------------------
|
||||
!
|
||||
!
|
||||
use kinds, only : DP
|
||||
use pwcom
|
||||
implicit none
|
||||
|
||||
real(DP) :: dvsym (nrx1, nrx2, nrx3, 3, 3)
|
||||
real(DP), allocatable :: aux (:,:,:,:,:)
|
||||
! the function to symmetrize
|
||||
! auxiliary space
|
||||
|
||||
integer :: ix, jx, kx, ri, rj, rk, irot, ip, jp, lp, mp, iflag
|
||||
! define a real-space point on the grid
|
||||
! the rotated points
|
||||
! counter on symmetries
|
||||
! counter on polarizations
|
||||
real(dp) :: det(48), sc(3,3), d
|
||||
|
||||
if (nsym.eq.1) return
|
||||
allocate (aux(nrx1 , nrx2 , nrx3 , 3, 3))
|
||||
|
||||
call DCOPY (nrx1 * nrx2 * nrx3 * 9, dvsym, 1, aux, 1)
|
||||
|
||||
! compute determinants of transformation matrixes
|
||||
do irot = 1, nsym
|
||||
if (iflag == 1) then ! pseudo-tensor
|
||||
sc(:,:) = dble(s(:,:,irot))
|
||||
! crystal to cartesian
|
||||
call trntns (sc, at, bg, 1)
|
||||
d = sc(1,1)*sc(2,2)*sc(3,3) + &
|
||||
sc(1,2)*sc(2,3)*sc(3,1) + &
|
||||
sc(1,3)*sc(2,1)*sc(3,2) - &
|
||||
sc(1,3)*sc(2,2)*sc(3,1) - &
|
||||
sc(1,2)*sc(2,1)*sc(3,3) - &
|
||||
sc(1,1)*sc(2,3)*sc(3,2)
|
||||
det(irot) = sign(1.d0,d)
|
||||
else ! tensor
|
||||
det(irot) = 1.d0
|
||||
endif
|
||||
enddo
|
||||
|
||||
dvsym (:,:,:,:,:) = 0.d0
|
||||
!
|
||||
! symmmetrize
|
||||
!
|
||||
do kx = 1, nr3
|
||||
do jx = 1, nr2
|
||||
do ix = 1, nr1
|
||||
do irot = 1, nsym
|
||||
call ruotaijk(s (1, 1, irot), ftau (1, irot), ix, jx, kx, &
|
||||
nr1, nr2, nr3, ri, rj, rk)
|
||||
!
|
||||
! ruotaijk finds the rotated of ix,jx,kx with the inverse of S
|
||||
!
|
||||
do ip = 1, 3
|
||||
do jp = 1, 3
|
||||
do lp = 1, 3
|
||||
do mp = 1, 3
|
||||
dvsym (ix, jx, kx, ip, jp) = &
|
||||
dvsym (ix, jx, kx, ip, jp) + det(irot)*&
|
||||
DBLE (s (ip, lp, irot))* &
|
||||
DBLE (s (jp, mp, irot))* &
|
||||
aux (ri, rj, rk, lp, mp)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call DSCAL (nrx1 * nrx2 * nrx3 * 9, 1.d0 / DBLE (nsym), dvsym , 1)
|
||||
|
||||
deallocate (aux)
|
||||
return
|
||||
end subroutine syme2
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,113 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE test_f_sum_rule
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... Test the f-sum rule
|
||||
! ... 2/N_k sum_k sum_n^{occ}
|
||||
! ... <u_{nk}| p_{k,\alpha} G_k v_{k,\beta} | u_{nk}> =
|
||||
! ... = -N_{el} \delta_{\alpha,\beta}
|
||||
! ...
|
||||
! ... where: p_k = -i\nabla + k and v_k = -i [r, H_k]
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE io_files, ONLY : nwordwfc, iunwfc
|
||||
USE cell_base, ONLY : at, bg, omega, tpiba, tpiba2
|
||||
USE wavefunctions_module, ONLY : evc
|
||||
USE klist, ONLY : nks, nkstot, wk, xk, nelec
|
||||
USE wvfct, ONLY : nbnd, npwx, npw, igk, wg, g2kin, &
|
||||
current_k
|
||||
USE lsda_mod, ONLY : current_spin, lsda, isk
|
||||
USE becmod, ONLY : becp
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
|
||||
!-- local variables ----------------------------------------------------
|
||||
IMPLICIT NONE
|
||||
complex(dp), allocatable, dimension(:,:,:) :: p_evc, vel_evc, g_vel_evc
|
||||
real(dp) :: f_sum(3,3), f_sum_k(3,3), q(3)
|
||||
integer :: ik, ipol, jpol, ibnd, ig
|
||||
complex(dp), external :: ZDOTC
|
||||
|
||||
! allocate memory
|
||||
allocate ( p_evc(npwx,nbnd,3), &
|
||||
vel_evc(npwx,nbnd,3), &
|
||||
g_vel_evc(npwx,nbnd,3) )
|
||||
|
||||
! zero the f-sum
|
||||
f_sum(:,:) = 0.d0
|
||||
q(:) = 0.d0
|
||||
|
||||
write(stdout, '(5X,''Computing the f-sum rule'')')
|
||||
|
||||
!====================================================================
|
||||
! loop over k-points
|
||||
!====================================================================
|
||||
do ik = 1, nks
|
||||
current_k = ik
|
||||
current_spin = isk(ik)
|
||||
|
||||
! initialize at k-point k
|
||||
call gk_sort(xk(1,ik), ngm, g, ecutwfc/tpiba2, npw, igk, g2kin)
|
||||
g2kin(:) = g2kin(:) * tpiba2
|
||||
call init_us_2(npw,igk,xk(1,ik),vkb)
|
||||
|
||||
! read wfcs from file and compute becp
|
||||
CALL davcio (evc, nwordwfc, iunwfc, ik, -1)
|
||||
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
|
||||
|
||||
!q = 0.d0; q(1) = 0.0d0
|
||||
!write(stdout, '(5X,''computing wfcs at k + q'')')
|
||||
!call compute_u_kq(ik, q)
|
||||
!evc = evq
|
||||
|
||||
! compute p_k|evc>, v_k|evc> and G_k v_k|evc>
|
||||
do ipol = 1, 3
|
||||
call apply_p(evc, p_evc(1,1,ipol), ik, ipol, q)
|
||||
call apply_vel(evc, vel_evc(1,1,ipol), ik, ipol, q)
|
||||
call greenfunction(ik, vel_evc(1,1,ipol), g_vel_evc(1,1,ipol), q)
|
||||
enddo
|
||||
|
||||
! k-point contribution to the f-sum rule
|
||||
f_sum_k = 0.0
|
||||
|
||||
! loop over cartesian directions
|
||||
do jpol = 1, 3
|
||||
do ipol = 1, 3
|
||||
do ibnd = 1, nbnd_occ (ik)
|
||||
f_sum_k(ipol,jpol) = f_sum_k(ipol,jpol) + wg(ibnd,ik) * &
|
||||
2.d0 * real(ZDOTC(npw, p_evc(1,ibnd,ipol), 1, &
|
||||
g_vel_evc(1,ibnd,jpol), 1))
|
||||
enddo
|
||||
enddo ! ipol
|
||||
enddo ! jpol
|
||||
|
||||
write(stdout, '(5X,''f-sum rule (ik='',I5,''):'')') ik
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') f_sum_k
|
||||
|
||||
f_sum(:,:) = f_sum(:,:) + f_sum_k(:,:)
|
||||
enddo ! ik
|
||||
#ifdef __PARA
|
||||
call reduce(9, f_sum)
|
||||
call poolreduce(9, f_sum)
|
||||
#endif
|
||||
write(stdout, '(5X,''f-sum rule:'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') f_sum(:,:)
|
||||
|
||||
call sym_cart_tensor(f_sum)
|
||||
write(stdout, '(5X,''f-sum rule (symmetrized):'')')
|
||||
write(stdout, '(3(5X,3(F12.6,2X)/))') f_sum
|
||||
|
||||
deallocate(p_evc, vel_evc, g_vel_evc)
|
||||
|
||||
END SUBROUTINE test_f_sum_rule
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,130 @@
|
|||
!
|
||||
! Copyright (C) 2001-2005 Quantum-ESPRESSO 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 .
|
||||
!
|
||||
#include "f_defs.h"
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE write_tensor_field(name, ispin, field)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! ... write the tensor field in xcrysden format
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE mp_global, ONLY : me_pool
|
||||
USE cell_base, ONLY : at, bg, alat
|
||||
USE ions_base, ONLY : nat, tau, atm, ityp
|
||||
USE pwcom
|
||||
USE nmr_module
|
||||
!--------------------------------------------------------------------
|
||||
character*(*) name
|
||||
integer :: ispin
|
||||
real(dp) :: field(nrx1,nrx2,nrx3,3)
|
||||
!--------------------------------------------------------------------
|
||||
integer, parameter :: ounit = 48
|
||||
character*80 :: fname
|
||||
integer :: ios, ipol
|
||||
|
||||
if (me_pool /= 0) return
|
||||
|
||||
do ipol = 1, 3
|
||||
! form the name of the output file
|
||||
if (ispin == 0) fname = trim(name)//'_'
|
||||
if (ispin == 1) fname = trim(name)//'_UP_'
|
||||
if (ispin == 2) fname = trim(name)//'_DW_'
|
||||
|
||||
if (ipol == 1) fname = trim(fname)//'X.xsf'
|
||||
if (ipol == 2) fname = trim(fname)//'Y.xsf'
|
||||
if (ipol == 3) fname = trim(fname)//'Z.xsf'
|
||||
write(stdout, '(5X,''write_tensor_field: '',A40)') fname
|
||||
|
||||
open(unit=ounit, file=fname, iostat=ios, form='formatted', &
|
||||
status='unknown')
|
||||
if (ios /= 0) &
|
||||
call errore('write_tensor_field', 'error opening '//fname, ounit)
|
||||
|
||||
call xsf_struct (alat, at, nat, tau, atm, ityp, nr1*nr2*nr3, ounit)
|
||||
call xsf_vector_3d(field(1,1,1,ipol), &
|
||||
nr1, nr2, nr3, nrx1, nrx2, nrx3, at, bg, alat, ounit)
|
||||
close(unit=48)
|
||||
enddo
|
||||
end subroutine write_tensor_field
|
||||
|
||||
|
||||
|
||||
!
|
||||
! Copyright (C) 2003 Tone Kokalj
|
||||
! 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 .
|
||||
!
|
||||
! This file holds XSF (=Xcrysden Structure File) utilities.
|
||||
! Routines written by Tone Kokalj on Mon Jan 27 18:51:17 CET 2003
|
||||
!
|
||||
! -------------------------------------------------------------------
|
||||
! this routine writes the crystal structure in XSF format
|
||||
! -------------------------------------------------------------------
|
||||
subroutine xsf_struct (alat, at, nat, tau, atm, ityp, nr, ounit)
|
||||
USE kinds, only : DP
|
||||
implicit none
|
||||
integer :: nat, ityp (nat), nr, ounit
|
||||
character(len=3) :: atm(*)
|
||||
real(DP) :: alat, tau (3, nat), at (3, 3)
|
||||
! --
|
||||
integer :: i, j, n
|
||||
real(DP) :: at1 (3, 3)
|
||||
! convert lattice vectors to ANGSTROM units ...
|
||||
do i=1,3
|
||||
do j=1,3
|
||||
at1(j,i) = at(j,i)*alat*0.529177d0
|
||||
enddo
|
||||
enddo
|
||||
|
||||
write(ounit,*) 'CRYSTAL'
|
||||
write(ounit,*) 'PRIMVEC'
|
||||
write(ounit,'(2(3F15.9/),3f15.9)') at1
|
||||
write(ounit,*) 'PRIMCOORD'
|
||||
write(ounit,*) nat + nr, 1
|
||||
|
||||
do n=1,nat
|
||||
! positions are in Angstroms
|
||||
write(ounit,'(a3,3x,3f15.9)') atm(ityp(n)), &
|
||||
tau(1,n)*alat*0.529177d0, &
|
||||
tau(2,n)*alat*0.529177d0, &
|
||||
tau(3,n)*alat*0.529177d0
|
||||
enddo
|
||||
return
|
||||
end subroutine xsf_struct
|
||||
|
||||
|
||||
|
||||
! -------------------------------------------------------------------
|
||||
! this routine writes a 3D vector field
|
||||
! -------------------------------------------------------------------
|
||||
subroutine xsf_vector_3d(v, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
|
||||
at, bg, alat, ounit)
|
||||
USE kinds, only : DP
|
||||
implicit none
|
||||
integer :: nrx1, nrx2, nrx3, nr1, nr2, nr3, ounit
|
||||
real(DP) :: at(3,3), bg(3,3), x(3), alat, v(nrx1,nrx2,nrx3,3)
|
||||
integer :: i1, i2, i3
|
||||
|
||||
do i1 = 1, nr1
|
||||
do i2 = 1, nr2
|
||||
do i3 = 1, nr3
|
||||
! coordinate in angstrom
|
||||
x(1) = dble(i1)/dble(nr1)
|
||||
x(2) = dble(i2)/dble(nr2)
|
||||
x(3) = dble(i3)/dble(nr3)
|
||||
! crystal to cartesian
|
||||
call trnvect (x, at, bg, 1)
|
||||
x = x * alat * 0.529177d0
|
||||
write(ounit,'(''X '',3x,3f15.9,2x,3e12.4)') x, v(i1,i2,i3,1:3)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end subroutine xsf_vector_3d
|
||||
Loading…
Reference in New Issue