scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Trying to stay up-to-date with the recent CVS changes. Gipaw compiles 

cleanly but doesn't not working at all!
I give up: I'm going in vacation for the next few days. (D.C.)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3814 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
ceresoli 2007-02-23 15:32:08 +00:00
parent 542c94701f
commit 50d7de6418
8 changed files with 150 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
GIPAW/apply_vel.f90
View File

@ -63,6 +63,7 @@ SUBROUTINE apply_vel(psi, vel_psi, ik, ipol, q)
! compute (1/2|dk|) ( V^{NL}_{k+dk+q,k+dk} |psi> -
! V^{NL}_{k-dk+q,k-dk} |psi> )
!====================================================================
allocate(becp(nkb,nbnd))
do isign = -1,1,2
dxk(:) = xk(:,ik)
dxk(ipol) = dxk(ipol) + isign * dk ! k + dk
@ -82,6 +83,7 @@ SUBROUTINE apply_vel(psi, vel_psi, ik, ipol, q)
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
deallocate(becp)
#else
do isign = -1,1,2

6
GIPAW/ch_psi_all2.f90
View File

@ -77,11 +77,7 @@ subroutine ch_psi_all2 (n, h, ah, e, ik, m)
!
! Here we compute the projector in the valence band
!
if (lgamma) then
ikq = ik
else
ikq = 2 * ik
endif
ikq = ik
ps (:,:) = (0.d0, 0.d0)
IF (noncolin) THEN

116
GIPAW/compute_u_kq.f90 Normal file
View File

@ -0,0 +1,116 @@
! Copyright (C) 2001-2007 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)
!----------------------------------------------------------------------------
!
! ... diagonalize at k+q
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE io_files, ONLY : iunigk, nwordatwfc, iunsat, iunwfc, &
nwordwfc
USE klist, ONLY : nkstot, nks, xk, ngk
USE uspp, ONLY : vkb, nkb
USE gvect, ONLY : g, nrxx, nr1, nr2, nr3
USE wvfct, ONLY : et, nbnd, npwx, igk, npw, g2kin, &
current_k, nbndx, btype
USE control_flags, ONLY : ethr, isolve, io_level
USE ldaU, ONLY : lda_plus_u, swfcatom
USE lsda_mod, ONLY : current_spin, lsda, isk
USE noncollin_module, ONLY : noncolin, npol
USE wavefunctions_module, ONLY : evc
USE gvect, ONLY : g, ngm, ecutwfc, ngl, nrxx, &
nr1, nr2, nr3, nrx1, nrx2, nrx3
USE cell_base, ONLY : at, bg, omega, tpiba, tpiba2
USE bp, ONLY : lelfield
USE control_flags, ONLY : iverbosity
USE becmod, ONLY : becp
USE control_flags, ONLY : ethr, isolve, max_cg_iter
USE gipaw_module
IMPLICIT NONE
INTEGER :: ik, iter ! k-point, current iterations
REAL(DP) :: q(3) ! q-vector
REAL(DP) :: avg_iter
INTEGER :: ig
REAL(DP) :: xkold(3)
REAL(DP), allocatable :: et_old(:,:)
CALL start_clock( 'c_bands' )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
iter = 1
max_cg_iter = 200
isolve = 1 ! CG
ethr = conv_threshold
if (allocated(btype)) deallocate(btype)
allocate(btype(nkstot,nbndx))
btype(:,:) = 1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! save eigenvalues
allocate( et_old(nbnd,nkstot) )
et_old = et
! debug
!!WRITE(stdout, '(5X,"compute_u_kq: q = (",F10.4,",",F10.4,",",F10.4,")")') q
!!WRITE(stdout, '(5X," isolve = ", I2)') isolve
avg_iter = 0.D0
current_k = ik
IF ( lsda ) current_spin = isk(ik)
! same sorting of G-vector at k+q
call gk_sort(xk(1,ik),ngm,g,ecutwfc/tpiba2,npw,igk,g2kin)
! set the k-point
xkold(:) = xk(:,ik)
xk(:,ik) = xk(:,ik) + q(:)
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
! various initializations
IF ( nkb > 0 ) CALL init_us_2( npw, igk, xk(1,ik), vkb )
! read in wavefunctions from the previous iteration
!!IF ( nks > 1 .OR. (io_level > 1) .OR. lelfield ) &
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
! Needed for LDA+U
IF ( lda_plus_u ) CALL davcio( swfcatom, nwordatwfc, iunsat, ik, -1 )
! diagonalization of bands for k-point ik
call diag_bands ( iter, ik, avg_iter )
call poolreduce( 1, avg_iter )
avg_iter = avg_iter / nkstot
! debug
!!WRITE( stdout, &
!! '( 5X,"ethr = ",1PE9.2,", avg # of iterations =",0PF5.1 )' ) &
!! ethr, avg_iter
! restore the k-point and eigenvalues
print*, et(1:nbnd,ik)*13.6056982d0
xk(:,ik) = xkold(:)
et = et_old
deallocate(et_old)
! restore wavefunctions
!!evq = evc
!!IF ( nks > 1 .OR. (io_level > 1) .OR. lelfield ) &
!! CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
CALL stop_clock( 'c_bands' )
RETURN
END SUBROUTINE compute_u_kq

4
GIPAW/gipaw_module.f90
View File

@ -174,7 +174,9 @@ CONTAINS
IMPLICIT NONE
allocate(evq(npwx,nbnd))
allocate(becp(nkb,nbnd))
! DO NOT not allocate becp here! it's allocated and deallocated
! by a lot of routines (crazy!)
! allocate(becp(nkb,nbnd))
allocate(j_bare(nrxxs,3,3,nspin), b_ind_r(nrxxs,3,3), b_ind(ngm,3,3))
allocate ( paw_recon(ntyp) )

30
GIPAW/greenfunction.f90
View File

@ -44,14 +44,13 @@ SUBROUTINE greenfunction(ik, psi, g_psi, q)
! allocate memory
allocate (work(npwx), ps(nbnd,nbnd), h_diag(npwx,nbnd), &
eprec(nbnd))
eprec(nbnd), becp(nkb,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.
!====================================================================
@ -85,9 +84,16 @@ SUBROUTINE greenfunction(ik, psi, g_psi, q)
! solve the linear system (apply G_{k+q})
!====================================================================
! convergence treshold
!thresh = 1d-12
thresh = sqrt(conv_threshold) ! sqrt(of that of PARATEC)
! 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
! preconditioning of the linear system
do ibnd = 1, nbnd_occ (ik)
do ig = 1, npw
@ -104,28 +110,26 @@ SUBROUTINE greenfunction(ik, psi, g_psi, q)
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)
!!call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, psi)
else
call init_us_2(npw, igk, xk(1,ik), vkb)
endif
!!call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, psi)
! initial guess
g_psi(:,:) = (0.d0, 0.d0)
! solve linear system
conv_root = .true.
!!PRINT*, et(1:nbnd,ik)*13.6056982d0
call cgsolve_all (ch_psi_all2, cg_psi, et(1,ik), psi, g_psi, &
h_diag, npwx, npw, thresh, ik, lter, conv_root, anorm, &
nbnd_occ(ik), 1 )
! debug
!!write(stdout, '(5X,''cgsolve_all converged in '',I3,'' iterations'')') &
!! lter
@ -140,6 +144,6 @@ SUBROUTINE greenfunction(ik, psi, g_psi, q)
call stop_clock('greenf')
! free memory
deallocate (work, h_diag, eprec, ps)
deallocate (work, h_diag, eprec, ps, becp)
END SUBROUTINE greenfunction

4
GIPAW/make.depend
View File

@ -21,14 +21,13 @@ compute_sigma.o : ../Modules/kind.o
compute_sigma.o : ../PW/pwcom.o
compute_sigma.o : gipaw_module.o
compute_u_kq.o : ../Modules/cell_base.o
compute_u_kq.o : ../Modules/control_flags.o
compute_u_kq.o : ../Modules/io_files.o
compute_u_kq.o : ../Modules/io_global.o
compute_u_kq.o : ../Modules/kind.o
compute_u_kq.o : ../Modules/uspp.o
compute_u_kq.o : ../Modules/wavefunctions.o
compute_u_kq.o : ../PW/becmod.o
compute_u_kq.o : ../PW/diis_base.o
compute_u_kq.o : ../PW/g_psi_mod.o
compute_u_kq.o : ../PW/noncol.o
compute_u_kq.o : ../PW/pwcom.o
compute_u_kq.o : gipaw_module.o
@ -58,6 +57,7 @@ g_tensor_crystal.o : ../PW/paw.o
g_tensor_crystal.o : ../PW/pwcom.o
g_tensor_crystal.o : gipaw_module.o
gipaw_main.o : ../Modules/cell_base.o
gipaw_main.o : ../Modules/control_flags.o
gipaw_main.o : ../Modules/io_files.o
gipaw_main.o : ../Modules/io_global.o
gipaw_main.o : ../Modules/kind.o

4
GIPAW/suscept_crystal.f90
View File

@ -104,7 +104,7 @@ SUBROUTINE suscept_crystal
! read wfcs from file and compute becp
call davcio (evc, nwordwfc, iunwfc, ik, -1)
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
!!call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
!<apsi>
call init_paw_2_no_phase (npw, igk, xk (1, ik), paw_vkb)
@ -119,7 +119,7 @@ SUBROUTINE suscept_crystal
!!!write(*,'(''q='',3(F12.4))') q
call compute_u_kq(ik, q)
evc = evq
!!evc = evq
! compute p_k|evc>, v_k|evc> and G_k v_{k,k}|evc>
call apply_operators

17
GIPAW/test_sum_rule.f90
View File

@ -26,7 +26,6 @@ SUBROUTINE test_f_sum_rule
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 gipaw_module
@ -60,15 +59,9 @@ SUBROUTINE test_f_sum_rule
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)
! read wfcs from file
call davcio (evc, nwordwfc, iunwfc, ik, -1)
!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)
@ -84,8 +77,12 @@ SUBROUTINE test_f_sum_rule
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, &
2.d0 * real(ZDOTC(npw, evc(1,ibnd), 1, &
!! 2.d0 * real(ZDOTC(npw, p_evc(1,ibnd,ipol), 1, &
g_vel_evc(1,ibnd,jpol), 1))
PRINT*, ibnd,ipol,jpol, 2.d0 * real(ZDOTC(npw, evc(1,ibnd), 1, &
g_vel_evc(1,ibnd,jpol), 1))
enddo
enddo ! ipol
enddo ! jpol