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

D3/Makefile : missing add_efield added 

cegterg, regterg: use GEMM instead of GEMV, do not require
any longer evc, et to be dimensioned nbndx (et is now
dimensioned (nbnd,nkstot)) - misc. cleanup

Il calcolo di (H-eS)*psi ('update') nella diagonalizzazione iterativa
prendeva un tempo esagerato. L'ho modificata in modo da usare prodotti
matrice-matrice su tutti gli psi invece che matrice-vettore su ogni
psi, se piu' di 1/4 dei vettori non e' a convergenza. La cosa e' fatta
a naso e richiede ulteriori di prove, ma mi sembra che apporti dei
miglioramenti.

In TODO ho messo una lista di cose da fare.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@145 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2003-04-07 12:55:29 +00:00
parent 6e3853d809
commit ab5b6af5e5
27 changed files with 402 additions and 362 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
D3/Makefile
View File

@ -165,6 +165,7 @@ PWOBJS = ../PW/pwcom.o \
../PW/addusdens.o \
../PW/addusforce.o \
../PW/addusstress.o \
../PW/add_efield.o \
../PW/add_vuspsi.o \
../PW/allocate_fft.o \
../PW/allocate_locpot.o \

13
Gamma/allocate_wfc.f90
View File

@ -19,15 +19,16 @@ subroutine allocate_wfc
!
! Allocate memory
!
allocate (et( nbndx, nkstot))
allocate (wg( nbnd, nkstot))
allocate (evc( npwx, nbndx))
allocate (et( nbnd, nkstot))
allocate (wg( nbnd, nkstot))
allocate (evc(npwx, nbndx))
allocate(becp(nkb, nbndx))
! Needed with LDA+U
!
! Needed for LDA+U
!
if (lda_plus_u) allocate (swfcatom( npwx, natomwfc))
call setv (nbndx * nkstot, 0.d0, et, 1)
et(:,:) = 0.d0
write (6, 100) nbndx, nbnd, natomwfc, npwx, nelec, nkb, ngl
100 format (/5x,'nbndx = ',i5,' nbnd = ',i5,' natomwfc = ',i5, &

4
Gamma/electrons.f90
View File

@ -126,7 +126,7 @@ implicit none
if (.not.lscf) then
conv_elec=.true.
#ifdef __PARA
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
#endif
do ik = 1, nkstot
@ -248,7 +248,7 @@ implicit none
enddo
call ireduce (nks, ngkp)
call ipoolrecover (ngkp, 1, nkstot, nks)
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
#endif
do ik = 1, nkstot

158
Gamma/regterg.f90
View File

@ -13,10 +13,10 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
!
! iterative solution of the eigenvalue problem:
!
! ( H - e S ) * psi = 0
! ( H - e S ) * evc = 0
!
! where H is an hermitean operator, e is a real scalar,
! S is an overlap matrix, psi is a complex vector
! S is an overlap matrix, evc is a complex vector
! (real wavefunctions with only half plane waves stored)
!
#include "machine.h"
@ -25,7 +25,7 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
! on INPUT
integer :: ndim, ndmx, nvec, nvecx, gstart
! dimension of the matrix to be diagonalized
! leading dimension of matrix psi, as declared in the calling pgm unit
! leading dimension of matrix evc, as declared in the calling pgm unit
! integer number of searched low-lying roots
! maximum dimension of the reduced basis set
! (the basis set is refreshed when its dimension would exceed nvecx)
@ -58,16 +58,16 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
! S matrix on the reduced basis
! eigenvectors of the Hamiltonian
! eigenvalues of the reduced hamiltonian
real(kind=DP) :: DDOT
real(kind=DP), external :: DDOT
complex(kind=DP), allocatable :: psi(:,:), hpsi (:,:),spsi (:,:)
! work space, contains psi
! the product of H and psi
! the product of S and psi
integer, allocatable :: conv (:)
! if 1 the root is converged if 0 no
logical, allocatable :: conv (:)
! if .true. the root is converged
!
! Called routines:
external h_psi, g_psi, DDOT
external h_psi, g_psi
! h_psi(ndmx,ndim,nvec,psi,hpsi,spsi)
! calculates nvec H|psi> and S|psi> (if needed) products.
! Vectors psi,hpsi,spsi are dimensioned (ndmx,nvec)
@ -88,13 +88,13 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
allocate (ew( nvecx))
allocate (conv( nvec))
if (nvec.gt.nvecx / 2) call errore ('regter', 'nvecx is too small',1)
if (nvec > nvecx / 2) call errore ('regter', 'nvecx is too small',1)
!
! prepare the hamiltonian for the first iteration
!
notcnv = nvec
nbase = nvec
call DCOPY (2*ndmx*nvec, evc, 1, psi, 1)
psi(:, 1:nvec) = evc(:, 1:nvec)
!
! hpsi contains h times the basis vectors
!
@ -103,8 +103,8 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
! hr contains the projection of the hamiltonian onto the reduced space
! vr contains the eigenvectors of hr
!
call setv (nvecx * nvecx, 0.d0, hr, 1)
call setv (nvecx * nvecx, 0.d0, vr, 1)
hr (:,:) = 0.d0
vr (:,:) = 0.d0
call DGEMM ('t', 'n', nbase, nbase, 2*ndim, 2.d0 , psi, &
2*ndmx, hpsi, 2*ndmx, 0.d0, hr, nvecx)
if (gstart.eq.2) call DGER (nbase, nbase, -1.d0, psi, 2*ndmx, &
@ -112,7 +112,7 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
#ifdef __PARA
call reduce (nbase * nvecx, hr)
#endif
call setv (nvecx * nvecx, 0.d0, sr, 1)
sr(:,:) = 0.d0
call DGEMM ('t', 'n', nbase, nbase, 2*ndim, 2.d0, psi, &
2*ndmx, spsi, 2*ndmx, 0.d0, sr, nvecx)
if (gstart.eq.2) call DGER (nbase, nbase, -1.d0, psi, 2*ndmx, &
@ -123,31 +123,47 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
do n = 1, nbase
e (n) = hr (n, n)
conv (n) = 0
conv (n) = .false.
vr (n, n) = 1.d0
enddo
!
! iterate
!
do kter = 1, maxter
iter = kter+1
np = nbase
iter = kter
call start_clock ('update')
do n = 1, nvec
if (conv (n) .eq.0) then
!
! this root not yet converged ... set a new basis vector
! (position np) to (h-es)psi ...
!
np = np + 1
ew(np) = e (n)
call setv (2 * ndim, 0.d0, psi (1, np), 1)
call DGEMV ('n', 2*ndim, nbase, 1.d0, hpsi, 2*ndmx, &
vr (1, n) , 1, 0.d0, psi (1, np) , 1)
call DGEMV ('n', 2*ndim, nbase, -e(n),spsi, 2*ndmx, &
vr (1, n) , 1, 1.d0, psi (1, np) , 1)
endif
enddo
if (notcnv < nvec) then
np = nbase
do n = 1, nvec
if ( .not.conv (n) ) then
!
! this root not yet converged ... set a new basis vector
! (position np) to (h-es)psi ...
!
np = np + 1
! for use in g_psi
ew(np) = e (n)
call DGEMV ('n', 2*ndim, nbase, 1.d0, hpsi, 2*ndmx, &
vr (1, n) , 1, 0.d0, psi (1, np) , 1)
call DGEMV ('n', 2*ndim, nbase, -e(n),spsi, 2*ndmx, &
vr (1, n) , 1, 1.d0, psi (1, np) , 1)
endif
enddo
else
!
! expand the basis set with new basis vectors (h-es)psi ...
!
call DGEMM ('n', 'n', 2*ndim, nvec, nbase, 1.d0, spsi, &
2*ndmx, vr, nvecx, 0.d0, psi (1, nbase+1), 2*ndmx)
do n = 1, nvec
! for use in g_psi
ew (nbase+n) = e (n)
psi (:,nbase+n) = - e(n) * psi(:,nbase+n)
end do
call DGEMM ('n', 'n', 2*ndim, nvec, nbase, 1.d0, hpsi, &
2*ndmx, vr, nvecx, 1.d0, psi (1, nbase+1), 2*ndmx)
end if
call stop_clock ('update')
!
! approximate inverse iteration
@ -212,62 +228,80 @@ subroutine regterg (ndim, ndmx, nvec, nvecx, evc, ethr, gstart, &
!
notcnv = 0
do n = 1, nvec
if (conv (n) .eq.0) then
if (.not.abs (ew (n) - e (n) ) .le.ethr) then
notcnv = notcnv + 1
else
! root converged
conv (n) = 1
endif
endif
!!! conv (n) = conv(n) .or. ( abs (ew (n) - e (n) ) <= ethr )
conv (n) = ( abs (ew (n) - e (n) ) <= ethr )
if ( .not. conv(n) ) notcnv = notcnv + 1
e (n) = ew (n)
end do
enddo
!!! TEST : update all eigenvectors if more than 1/4 are converged
if (notcnv > nvec/4) then
notcnv = nvec
conv(:) = .false.
end if
!!! END OF TEST
!
! if overall convergence has been achieved or the dimension of the
! reduced basis set is becoming too large, refresh the basis set
! i.e. replace the first nvec elements with the current estimate
! of the eigenvectors and set the basis dimension to nvec.
! if overall convergence has been achieved, OR
! the dimension of the reduced basis set is becoming too large, OR
! in any case if we are at the last iteration
! refresh the basis set. i.e. replace the first nvec elements
! with the current estimate of the eigenvectors;
! set the basis dimension to nvec.
!
if (notcnv.eq.0 .or. iter.gt.maxter .or. nbase+notcnv.gt.nvecx) then
if (notcnv == 0 .or. nbase+notcnv > nvecx .or. iter == maxter) then
call start_clock ('last')
call DGEMM ('n', 'n', 2*ndim, nvec, nbase, 1.d0, psi, &
2*ndmx, vr, nvecx, 0.d0, evc, 2*ndmx)
if (notcnv.eq.0.or. iter.gt.maxter) then
if (notcnv == 0) then
!
! all roots converged: return
!
call stop_clock ('last')
goto 10
else if (iter == maxter) then
!
! last iteration, some roots not converged: return
!
#ifdef DEBUG_DAVIDSON
do n = 1, nvec
if ( .not.conv (n) ) write (6, '(" WARNING: e(",i3,") =",&
f10.5," is not converged to within ",1pe8.1)') n, e(n), ethr
enddo
#else
write (6, '(" WARNING: ",i5," eigenvalues not converged")') &
notcnv
#endif
call stop_clock ('last')
goto 10
end if
call DCOPY (2*ndmx*nvec, evc, 1, psi, 1)
!
! refresh psi, H*psi and S*psi
!
psi(:, 1:nvec) = evc(:, 1:nvec)
call DGEMM ('n', 'n', 2*ndim, nvec, nbase, 1.d0, spsi, &
2*ndmx, vr, nvecx, 0.d0, psi(1,nvec+1), 2*ndmx)
call DCOPY (2 * ndmx * nvec, psi(1,nvec+1), 1, spsi, 1)
2*ndmx, vr, nvecx, 0.d0, psi(1, nvec + 1), 2*ndmx)
spsi(:, 1:nvec) = psi(:, nvec+1:2*nvec)
call DGEMM ('n', 'n', 2*ndim, nvec, nbase, 1.d0, hpsi, &
2*ndmx, vr, nvecx, 0.d0, psi(1,nvec+1), 2*ndmx)
call DCOPY (2 * ndmx * nvec, psi(1,nvec+1), 1, hpsi, 1)
2*ndmx, vr, nvecx, 0.d0, psi(1, nvec + 1), 2*ndmx)
hpsi(:, 1:nvec) = psi(:, nvec+1:2*nvec)
!
! modify the reduced hamiltonian accordingly
! refresh the reduced hamiltonian
!
call setv (nvecx * nvecx, 0.d0, hr, 1)
call DCOPY (nbase, e, 1, hr, nvecx + 1 )
call setv (nvecx * nvecx, 0.d0, sr, 1)
call setv (nvecx, 1.d0, sr, nvecx + 1 )
call stop_clock ('last')
nbase = nvec
call setv (nvecx * nbase, 0.d0, vr, 1)
hr (:, 1:nbase) = 0.d0
sr (:, 1:nbase) = 0.d0
vr (:, 1:nbase) = 0.d0
do n = 1, nbase
hr (n, n) = e(n)
sr (n, n) = 1.d0
vr (n, n) = 1.d0
enddo
end do
call stop_clock ('last')
endif
enddo
10 continue
do n = 1, nvec
if (conv (n) .eq.0) write (6, '(" WARNING: e(",i3,") =",&
& f10.5," is not converged to within ",1pe8.1)') n, e(n), ethr
enddo
deallocate (conv)
deallocate (ew)

6
Gamma/sum_band.f90
View File

@ -50,10 +50,10 @@ subroutine sum_band
!
elseif (ltetra) then
#ifdef __PARA
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
if (me.eq.1.and.mypool.eq.1) then
#endif
call tweights (nkstot, nspin, nbndx, nbnd, nelec, ntetra, &
call tweights (nkstot, nspin, nbnd, nelec, ntetra, &
tetra, et, ef, wg)
#ifdef __PARA
endif
@ -62,7 +62,7 @@ subroutine sum_band
call broadcast (1, ef)
#endif
elseif (lgauss) then
call gweights (nks, wk, nbndx, nbnd, nelec, degauss, ngauss, &
call gweights (nks, wk, nbnd, nelec, degauss, ngauss, &
et, ef, demet, wg)
endif
!

2
Gamma/wfcinit.f90
View File

@ -129,7 +129,7 @@ subroutine wfcinit
enddo
if (iprint.eq.1) then
#ifdef __PARA
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
#endif
do ik = 1, nkstot
write (6, 9010) (xk (ipol, ik), ipol = 1, 3)

12
PH/elphon.f90
View File

@ -331,9 +331,8 @@ subroutine elphsum
! Note that the weights of k+q points must be set to zero for the
! following call to yield correct results
!
call efermig (et, nbndx, nbnd, nks, nelec, wk, degauss1, ngauss1, &
ef1)
dosef = dos_ef (ngauss1, degauss1, ef1, et, nbndx, wk, nks, nbnd)
call efermig (et, nbnd, nks, nelec, wk, degauss1, ngauss1, ef1)
dosef = dos_ef (ngauss1, degauss1, ef1, et, wk, nks, nbnd)
! N(Ef) is the DOS per spin, not summed over spin
dosef = dosef / 2.d0
!
@ -442,15 +441,14 @@ subroutine elphsum
end subroutine elphsum
!
!-----------------------------------------------------------------------
function dos_ef (ngauss, degauss, ef, et, nbndx, wk, nks, &
nbnd)
function dos_ef (ngauss, degauss, ef, et, wk, nks, nbnd)
!-----------------------------------------------------------------------
!
use parameters, only : DP
implicit none
real(kind=DP) :: dos_ef
integer :: ngauss, nbndx, nbnd, nks
real(kind=DP) :: et (nbndx, nks), wk (nks), ef, degauss
integer :: ngauss, nbnd, nks
real(kind=DP) :: et (nbnd, nks), wk (nks), ef, degauss
!
integer :: ik, ibnd
real(kind=DP) :: w0gauss

4
PP/dos.f90
View File

@ -116,9 +116,9 @@ subroutine dos (nodenumber)
do n= 1, ndos
E = Emin + (n - 1) * DeltaE
if (ltetra) then
call dos_t(et,nspin,nbndx,nbnd, nks,ntetra,tetra, E, DOSofE)
call dos_t(et,nspin,nbnd, nks,ntetra,tetra, E, DOSofE)
else
call dos_g(et,nspin,nbndx,nbnd, nks,wk,degauss,ngauss, E, DOSofE)
call dos_g(et,nspin,nbnd, nks,wk,degauss,ngauss, E, DOSofE)
endif
if (nspin.eq.1) then
DOSint = DOSint + DOSofE (1) * DeltaE

9
PP/dosg.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,15 +7,14 @@
!
!
!--------------------------------------------------------------------
subroutine dos_g (et, nspin, nbndx, nbnd, nks, wk, Degauss, &
ngauss, E, dosg)
subroutine dos_g (et, nspin, nbnd, nks, wk, Degauss, ngauss, E, dosg)
!--------------------------------------------------------------------
!
use parameters, only : DP
implicit none
integer :: nspin, nks, nbndx, nbnd, ngauss
integer :: nspin, nks, nbnd, ngauss
real(kind=DP) :: wk (nks), et (nbndx, nks), Degauss, E, dosg (2)
real(kind=DP) :: wk (nks), et (nbnd, nks), Degauss, E, dosg (2)
real(kind=DP) :: w0gauss
integer :: n, ns, nk0, nk, ik
external w0gauss

8
PP/dost.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,14 +7,14 @@
!
!
!--------------------------------------------------------------------
subroutine dos_t (et, nspin, nbndx, nbnd, nks, ntetra, tetra, e, dost)
subroutine dos_t (et, nspin, nbnd, nks, ntetra, tetra, e, dost)
!------------------------------------------------------------------
!
use parameters, only : DP
implicit none
integer :: nspin, nbndx, nbnd, nks, ntetra, tetra (4, ntetra)
integer :: nspin, nbnd, nks, ntetra, tetra (4, ntetra)
real(kind=DP) :: et (nbndx, nks), e, dost (2)
real(kind=DP) :: et (nbnd, nks), e, dost (2)
integer :: itetra (4), nk, ns, nt, ibnd, i
real(kind=DP) :: etetra (4), e1, e2, e3, e4

2
PP/projwave.f90
View File

@ -274,7 +274,7 @@ subroutine projwave (io_choice,Emin, Emax, DeltaE, smoothing)
!
! recover the vector proj over the pools
!
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
call poolrecover (proj, nbnd * natomwfc, nkstot, nks)
!
if (me.eq.1.and.mypool.eq.1) then

50
PP/stm.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -12,15 +12,15 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
!
! This routine calculates an stm image defined as the local density
! of states at the fermi energy.
! To this end uses a matched exponential behaved wavefunction, in th
! To this end uses a matched exponential behaved wavefunction, in the
! spirit of the Tersoff and Hamann approximation (PRB 31, 805 (1985)
! The matching with the true wavefunction is decided by the variable
! in celldm(1) units, then stm images are calculateted every dz.
! z in celldm(1) units, then stm images are calculateted every dz.
! The bias of the sample is decided by sample_bias, states between
! ef and ef + sample_bias are taken into account.
! It needs the workfunction wf. On output wf contains the number of
! states used to compute the image.
! The slab must be orientated with the main axis along celldm(3).
! The slab must be oriented with the main axis along celldm(3).
! It may not properly work if the slab has two symmetric surfaces.
!
#include "machine.h"
@ -38,7 +38,7 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
integer :: istates, igs, npws, ir, ir1, irx, iry, irz, ig, ibnd, &
ik, nbnd_ocp
! the number of states to compute the im
! the number of states to compute the image
! counter on surface G vectors
! number of surfac g-vectors
! counters on 3D r points
@ -67,11 +67,10 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
write ( * , * ) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
!
! if matching is .true. then matches the wfc's and uses their
! exponential behaviour, otherwise uses the true wfc's on the fft gr
! exponential behaviour, otherwise uses the true wfc's on fft grid
!
call setv (nrx1 * nrx2 * nrx3, 0.d0, stmdos, 1)
if (.not.stm_wfc_matching) call setv (nrxx, 0.d0, rho, 1)
stmdos(:) = 0.d0
if (.not.stm_wfc_matching) rho(:,:) = 0.d0
if (stm_wfc_matching) then
z = z * alat
@ -122,7 +121,7 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
endif
!
! take only the states in the energy window above or below the fermi
! as determined by the bias of the sample
! energy as determined by the bias of the sample
!
if (sample_bias.gt.0) then
up = ef + sample_bias
@ -170,14 +169,15 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
uguale = .false.
do igs = 1, npws
!
! if the surface part of G is equal to at least one of the surface v
! already found then uguale = .true.
! if the surface part of G is equal to at least one of the
! surface vectors already found then uguale = .true.
!
uguale = uguale.or. (g (1, igk (ig) ) .eq.gs (1, igs) .and.g ( &
2, igk (ig) ) .eq.gs (2, igs) )
enddo
!
! if G is not equal to any surface vector then G is a new surface ve
! if G is not equal to any surface vector then G is a new
! surface vector
!
if (.not.uguale) then
npws = npws + 1
@ -196,12 +196,11 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
!
istates = istates + 1
!
! for this state the work function is modified accordingly to its en
!
!
! find the coefficients of the matching wfcs
!
if (stm_wfc_matching) then
! for this state the work function is modified accordingly
! to its energy
wf1 = wf - (et (ibnd, ik) - ef - sample_bias)
do igs = 1, npws
a (igs) = (0.d0, 0.d0)
@ -219,18 +218,19 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
enddo
enddo
!
! reconstruct the wfc for the z in interest for this k-point and thi
! uses nr3/2 planes only, the other nr3/2 are empty -> only the uppe
! surface is used. N.B. it may not properly work if the upper surfac
! is connected to the lower surface by a symmetry operation (one
! should take the average of the two surfaces...).
! reconstruct the wfc for the z of interest for this k-point
! and this band. Uses nr3/2 planes only, the other nr3/2 are
! empty -> only the upper surface is used.
! N.B. it may not properly work if the upper surfacw
! is connected to the lower surface by a symmetry operation
! (one should take the average of the two surfaces...).
!
do irz = 2, nr3 / 2
!
! zz is the new z
!
zz = z + dz * (irz - 2)
call setv (2 * nrx1 * nrx2, 0.d0, psi, 1)
psi(:,:) = (0.d0, 0.d0)
do igs = 1, npws
fac = exp ( - sqrt (wf1 + ( (xk (1, ik) + gs (1, igs) ) **2 + &
(xk (2, ik) + gs (2, igs) ) **2) * tpiba2) * zz)
@ -256,8 +256,8 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
call reduce (2 * nrx1 * nrx2, psi)
#endif
!
! now sum for each k-point and for each band the square modulus of t
! wfc times the weighting factor
! now sum for each k-point and for each band the square
! modulus of the wfc times the weighting factor
!
do iry = 1, nr2
do irx = 1, nr1
@ -271,7 +271,7 @@ subroutine stm (wf, sample_bias, z, dz, stm_wfc_matching, stmdos)
!
! do not match
!
call setv (2 * nrxx, 0.d0, psic, 1)
psic(:) = (0.d0, 0.d0)
do ig = 1, npw
psic (nl (igk (ig) ) ) = evc (ig, ibnd)

15
PW/allocate_wfc.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -19,15 +19,16 @@ subroutine allocate_wfc
!
! Allocate memory
!
allocate (et( nbndx, nkstot))
allocate (wg( nbnd, nkstot))
allocate (evc( npwx, nbndx))
allocate (et( nbnd, nkstot))
allocate (wg( nbnd, nkstot))
allocate (evc(npwx, nbndx))
allocate(becp(nkb, nbndx))
! Needed with LDA+U
!
! Needed for LDA+U
!
if (lda_plus_u) allocate (swfcatom( npwx, natomwfc))
call setv (nbndx * nkstot, 0.d0, et, 1)
et(:,:) = 0.d0
write (6, 100) nbndx, nbnd, natomwfc, npwx, nelec, nkb, ngl
100 format (/5x,'nbndx = ',i5,' nbnd = ',i5,' natomwfc = ',i5, &

236
PW/cegterg.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,16 +7,16 @@
!
!
!----------------------------------------------------------------------
subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
subroutine cegterg (ndim, ndmx, nvec, nvecx, evc, ethr, overlap, &
e, notcnv, iter)
!----------------------------------------------------------------------
!
! iterative solution of the eigenvalue problem:
!
! ( H - e S ) * psi = 0
! ( H - e S ) * evc = 0
!
! where H is an hermitean operator, e is a real scalar,
! S is an overlap matrix, psi is a complex vector
! S is an overlap matrix, evc is a complex vector
!
#include "machine.h"
use parameters, only : DP
@ -24,11 +24,11 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
! on INPUT
integer :: ndim, ndmx, nvec, nvecx
! dimension of the matrix to be diagonalized
! leading dimension of matrix psi, as declared in the calling pgm unit
! leading dimension of matrix evc, as declared in the calling pgm unit
! integer number of searched low-lying roots
! maximum dimension of the reduced basis set
! (the basis set is refreshed when its dimension would exceed nvecx)
complex(kind=DP) :: psi (ndmx, nvecx)
complex(kind=DP) :: evc (ndmx, nvec)
real(kind=DP) :: ethr
! energy threshold for convergence
! root improvement is stopped, when two consecutive estimates of the root
@ -37,9 +37,8 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
! if .true. : use overlap matrix
! if .false. : use orthogonalization
! on OUTPUT
! complex*16 psi the first nvec columns contain the
! refined estimates of the eigenvectors
real(kind=DP) :: e (nvecx)
! evc contains the refined estimates of the eigenvectors
real(kind=DP) :: e (nvec)
! contains the estimated roots.
integer :: iter, notcnv
@ -58,21 +57,22 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
! dimension of the reduced basis
! counter on the reduced basis vectors
! do-loop counters
complex(kind=DP), allocatable :: hc (:,:), smat (:,:), vc (:,:), aux (:)
complex(kind=DP), allocatable :: hc (:,:), sc (:,:), vc (:,:)
! Hamiltonian on the reduced basis
! S matrix on the reduced basis
! the eigenvectors of the Hamiltonian
! work space
complex(kind=DP), allocatable :: hpsi (:,:), spsi (:,:)
complex(kind=DP), allocatable :: psi(:,:), hpsi (:,:), spsi (:,:)
! work space, contains psi
! the product of H and psi
! the product of S and psi
complex(kind=DP) :: ZDOTC, eau
complex(kind=DP), external :: ZDOTC
! scalar product routine
complex(kind=DP) :: eau
! auxiliary complex variable
real(kind=DP), allocatable :: ew (:)
! eigenvalues of the reduced hamiltonian
integer, allocatable :: conv (:)
! if 1 the root is converged if 0 no
logical, allocatable :: conv (:)
! true if the root is converged
!
! Called routines:
external h_psi, g_psi
@ -82,28 +82,28 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
! g_psi(ndmx,ndim,notcnv,psi,e)
! calculates (diag(h)-e)^-1 * psi, diagonal approx. to (h-e)^-1*psi
! the first nvec columns contain the trial eigenvectors
external ZDOTC
!
! allocate the work arrays
!
call start_clock ('cegterg')
if (overlap) allocate(smat(nvecx, nvecx))
allocate(hc (nvecx, nvecx))
allocate(vc (nvecx, nvecx))
allocate( psi (ndmx, nvecx))
allocate(hpsi (ndmx, nvecx))
allocate(spsi (ndmx, nvecx))
allocate(aux (ndmx * nvec))
if (overlap) allocate(sc(nvecx, nvecx))
allocate(hc (nvecx, nvecx))
allocate(vc (nvecx, nvecx))
allocate(ew (nvecx))
allocate(conv (nvec))
! write(6,*) 'eneter cegter',hc,vc,hpsi
if (nvec.gt.nvecx / 2) call errore ('cegter', 'nvecx is too small',1)
if (nvec > nvecx / 2) call errore ('cegter', 'nvecx is too small',1)
!
! prepare the hamiltonian for the first iteration
!
notcnv = nvec
nbase = nvec
psi(:, 1:nvec) = evc(:, 1:nvec)
!
! hpsi contains h times the basis vectors
!
@ -120,17 +120,17 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
call reduce (2 * nbase * nvecx, hc)
#endif
if (overlap) then
smat(:,:) = (0.d0, 0.d0)
sc(:,:) = (0.d0, 0.d0)
call ZGEMM ('c', 'n', nbase, nbase, ndim, (1.d0, 0.d0) , psi, &
ndmx, spsi, ndmx, (0.d0, 0.d0) , smat, nvecx)
ndmx, spsi, ndmx, (0.d0, 0.d0) , sc, nvecx)
#ifdef __PARA
call reduce (2 * nbase * nvecx, smat)
call reduce (2 * nbase * nvecx, sc)
#endif
endif
do n = 1, nbase
e (n) = hc (n, n)
conv (n) = 0
conv (n) = .false.
vc (n, n) = (1.d0, 0.d0)
enddo
!
@ -138,34 +138,49 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
!
do kter = 1, maxter
iter = kter
np = nbase
call start_clock ('update')
do n = 1, nvec
if (conv (n) .eq.0) then
!
! this root not yet converged ... set a new basis vector
! (position np) to (h-es)psi ...
!
np = np + 1
e (np) = e (n)
psi(:,np) = (0.d0,0.d0)
call ZGEMV ('n', ndim, nbase, (1.d0, 0.d0) , hpsi, ndmx, vc (1, &
n) , 1, (0.d0, 0.d0) , psi (1, np) , 1)
eau = DCMPLX ( - e (n), 0.d0)
call ZGEMV ('n', ndim, nbase, eau, spsi, ndmx, vc (1, n) , 1, &
(1.d0, 0.d0) , psi (1, np) , 1)
endif
enddo
if (notcnv < nvec) then
np = nbase
do n = 1, nvec
if ( .not.conv (n) ) then
!
! this root not yet converged ... set a new basis vector
! (position np) to (h-es)psi ...
!
np = np + 1
! for use in g_psi
ew (np) = e (n)
call ZGEMV ('n', ndim, nbase, (1.d0, 0.d0) , hpsi, ndmx, &
vc (1, n) , 1, (0.d0, 0.d0) , psi (1, np) , 1)
eau = DCMPLX ( - e (n), 0.d0)
call ZGEMV ('n', ndim, nbase, eau, spsi, ndmx, vc (1, n) , 1, &
(1.d0, 0.d0) , psi (1, np) , 1)
endif
enddo
else
!
! expand the basis set with new basis vectors (h-es)psi ...
!
call ZGEMM ('n', 'n', ndim, nvec, nbase, (1.d0, 0.d0), spsi, &
ndmx, vc, nvecx, (0.d0, 0.d0), psi (1, nbase+1), ndmx)
do n = 1, nvec
! for use in g_psi
ew (nbase+n) = e (n)
psi (:,nbase+n) = - e(n) * psi(:,nbase+n)
end do
call ZGEMM ('n', 'n', ndim, nvec, nbase, (1.d0, 0d0), hpsi, &
ndmx, vc, nvecx, (1.d0, 0.d0), psi (1, nbase+1), ndmx)
end if
call stop_clock ('update')
!
! approximate inverse iteration
!
call g_psi (ndmx, ndim, notcnv, psi (1, nbase+1), e (nbase+1) )
call g_psi (ndmx, ndim, notcnv, psi (1, nbase+1), ew (nbase+1) )
!
! "normalize" correction vectors psi(*,nbase+1:nbase+notcnv) in order to
! the numerical stability of cdiaghg (use ew as work array for norm of p
! "normalize" correction vectors psi(*,nbase+1:nbase+notcnv) in order
! to improve numerical stability of subspace diagonalization rdiaghg
! ew is used as work array : ew = <psi_i|psi_i>, i=nbase+1,nbase+notcnv
!
do n = 1, notcnv
ew (n) = ZDOTC (ndim, psi (1, nbase+n), 1, psi (1, nbase+n), 1)
@ -200,17 +215,18 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
call reduce (2 * nvecx * notcnv, hc (1, nbase+1) )
#endif
if (overlap) then
call ZGEMM ('c', 'n', nbase+notcnv, notcnv, ndim, (1.d0, 0.d0) &
, psi, ndmx, spsi (1, nbase+1) , ndmx, (0.d0, 0.d0) , smat (1, &
nbase+1) , nvecx)
call ZGEMM ('c', 'n', nbase+notcnv, notcnv, ndim, (1.d0, 0.d0), &
psi, ndmx, spsi (1, nbase+1) , ndmx, (0.d0, 0.d0) , &
sc (1, nbase+1) , nvecx)
#ifdef __PARA
call reduce (2 * nvecx * notcnv, smat (1, nbase+1) )
call reduce (2 * nvecx * notcnv, sc (1, nbase+1) )
#endif
endif
call stop_clock ('overlap')
nbase = nbase+notcnv
do n = 1, nbase
! the diagonal of hc must be strictly real
hc (n, n) = DCMPLX (DREAL (hc (n, n) ), 0.d0)
do m = n + 1, nbase
hc (m, n) = conjg (hc (n, m) )
@ -221,12 +237,12 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
!
if (overlap) then
do n = 1, nbase
smat (n, n) = DCMPLX (DREAL (smat (n, n) ), 0.d0)
sc (n, n) = DCMPLX (DREAL (sc (n, n) ), 0.d0)
do m = n + 1, nbase
smat (m, n) = conjg (smat (n, m) )
sc (m, n) = conjg (sc (n, m) )
enddo
enddo
call cdiaghg (nbase, nvec, hc, smat, nvecx, ew, vc)
call cdiaghg (nbase, nvec, hc, sc, nvecx, ew, vc)
#ifdef DEBUG_DAVIDSON
write (6,'(a,18f10.6)') 'EIG=',(e(n),n=1,nvec)
write (6,'(a,18f10.6)') 'EIG=',(ew(n),n=1,nvec)
@ -240,87 +256,93 @@ subroutine cegterg (ndim, ndmx, nvec, nvecx, psi, ethr, overlap, &
!
notcnv = 0
do n = 1, nvec
! if (conv (n) .eq.0) then
if (.not.abs (ew (n) - e (n) ) .le.ethr) then
e (n) = ew (n)
conv (n) = 0
notcnv = notcnv + 1
else
e (n) = ew (n)
! root converged
conv (n) = 1
endif
! endif
!!! conv (n) = conv(n) .or. ( abs (ew (n) - e (n) ) <= ethr )
conv (n) = ( abs (ew (n) - e (n) ) <= ethr )
if ( .not. conv(n) ) notcnv = notcnv + 1
e (n) = ew (n)
enddo
! if overall convergence has been achieved or the dimension of the
! reduced basis set is becoming too large, then refresh the basis
! set. i.e. replace the first nvec elements with the current estimat
! of the eigenvectors and set the basis dimension to nvec.
if (notcnv.eq.0.or.nbase+notcnv.gt.nvecx) then
!!! TEST : update all eigenvectors if more than 1/4 are converged
if (notcnv > nvec/4) then
notcnv = nvec
conv(:) = .false.
end if
!!! END OF TEST
!
! if overall convergence has been achieved, OR
! the dimension of the reduced basis set is becoming too large, OR
! in any case if we are at the last iteration
! refresh the basis set. i.e. replace the first nvec elements
! with the current estimate of the eigenvectors;
! set the basis dimension to nvec.
!
if ( notcnv == 0 .or. nbase+notcnv > nvecx .or. iter == maxter) then
call start_clock ('last')
aux(:) = (0.d0,0.d0)
call ZGEMM ('n', 'n', ndim, nvec, nbase, (1.d0, 0.d0) , psi, &
ndmx, vc, nvecx, (0.d0, 0.d0) , aux, ndmx)
call ZCOPY (ndmx * nvec, aux, 1, psi, 1)
if (notcnv.eq.0) then
ndmx, vc, nvecx, (0.d0, 0.d0) , evc, ndmx)
if (notcnv == 0) then
!
! all roots converged: return
!
call stop_clock ('last')
goto 10
else if (iter == maxter) then
!
! last iteration, some roots not converged: return
!
#ifdef DEBUG_DAVIDSON
do n = 1, nvec
if ( .not.conv (n) ) write (6, '(" WARNING: e(",i3,") =",&
f10.5," is not converged to within ",1pe8.1)') n, e(n), ethr
enddo
#else
write (6, '(" WARNING: ",i5," eigenvalues not converged")') &
notcnv
#endif
call stop_clock ('last')
goto 10
end if
!
! refresh psi, H*psi and S*psi
!
psi(:, 1:nvec) = evc(:, 1:nvec)
call ZGEMM ('n', 'n', ndim, nvec, nbase, (1.d0, 0.d0) , spsi, &
ndmx, vc, nvecx, (0.d0, 0.d0) , aux, ndmx)
call ZCOPY (ndmx * nvec, aux, 1, spsi, 1)
ndmx, vc, nvecx, (0.d0, 0.d0) , psi(1, nvec + 1), ndmx)
spsi(:, 1:nvec) = psi(:, nvec+1:2*nvec)
call ZGEMM ('n', 'n', ndim, nvec, nbase, (1.d0, 0.d0) , hpsi, &
ndmx, vc, nvecx, (0.d0, 0.d0) , psi (1, nvec + 1) , ndmx)
call ZCOPY (ndmx * nvec, psi (1, nvec + 1), 1, hpsi, 1)
hpsi(:, 1:nvec) = psi(:, nvec+1:2*nvec)
!
! modify the reduced hamiltonian accordingly
! refresh the reduced hamiltonian
!
hc(:,:) = (0.d0, 0.d0)
call DCOPY (nbase, e, 1, hc, 2 * (nvecx + 1) )
nbase = nvec
hc (:, 1:nbase) = (0.d0, 0.d0)
vc (:, 1:nbase) = (0.d0, 0.d0)
do n = 1, nbase
hc (n, n) = e(n)
vc (n, n) = (1.d0, 0.d0)
enddo
if (overlap) then
smat = (0.d0, 0.d0)
do n = 1, nvecx
smat (n, n) = (1.d0, 0.d0)
sc (:, 1:nbase) = (0.d0, 0.d0)
do n = 1, nbase
sc (n, n) = (1.d0, 0.d0)
enddo
endif
call stop_clock ('last')
nbase = nvec
vc(:,1:nbase) = (0.d0, 0.d0)
do n = 1, nbase
vc (n, n) = (1.d0, 0.d0)
enddo
endif
enddo
!
! replace the first nvec elements with the current estimate of
! the eigenvectors before quiting also when overall convergence
! has not been achieved
!
aux(:) = (0.d0,0.d0)
call ZGEMM ('n', 'n', ndim, nvec, nbase, (1.d0, 0.d0) , psi, &
ndmx, vc, nvecx, (0.d0, 0.d0) , aux, ndmx)
call ZCOPY (ndmx * nvec, aux, 1, psi, 1)
do n = 1, nvec
if (conv (n) .eq.0) write (6, '(" WARNING: e(",i3,") =",&
& f10.5," is not converged to within ",1pe8.1)') n, e(n), ethr
endif
enddo
10 continue
deallocate (conv)
deallocate (ew)
deallocate (aux)
deallocate (spsi)
deallocate (hpsi)
deallocate (vc)
deallocate ( psi)
if (overlap) deallocate (sc)
deallocate (hc)
if (overlap) deallocate (smat)
deallocate (vc)
deallocate (ew)
deallocate (conv)
call stop_clock ('cegterg')
return

19
PW/efermig.f90
View File

@ -1,23 +1,22 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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 efermig (et, nbndx, nbnd, nks, nelec, wk, Degauss, &
Ngauss, Ef)
subroutine efermig (et, nbnd, nks, nelec, wk, Degauss, Ngauss, Ef)
!--------------------------------------------------------------------
!
! Finds the Fermi energy - Gaussian Broadening (Methfessel-Paxton)
!
use parameters
implicit none
integer :: nks, nbndx, nbnd, i, kpoint, Ngauss
real(kind=DP) :: wk (nks), et (nbndx, nks), Degauss, Ef, Eup, Elw, &
sumkg, sumkup, sumklw, sumkmid, nelec
external sumkg
integer :: nks, nbnd, i, kpoint, Ngauss
real(kind=DP) :: wk (nks), et (nbnd, nks), Degauss, Ef, Eup, Elw, &
sumkup, sumklw, sumkmid, nelec
real(kind=DP), external:: sumkg
!
! find bounds for the Fermi energy. Very safe choice!
!
@ -39,13 +38,13 @@ subroutine efermig (et, nbndx, nbnd, nks, nelec, wk, Degauss, &
!
! Bisection method
!
sumkup = sumkg (et, nbndx, nbnd, nks, wk, Degauss, Ngauss, Eup)
sumklw = sumkg (et, nbndx, nbnd, nks, wk, Degauss, Ngauss, Elw)
sumkup = sumkg (et, nbnd, nks, wk, Degauss, Ngauss, Eup)
sumklw = sumkg (et, nbnd, nks, wk, Degauss, Ngauss, Elw)
if ( (sumkup - nelec) .lt. - 1.0e-10.or. (sumklw - nelec) &
.gt.1.0e-10) call errore ('Efermi', 'unexpected error', 1)
do i = 1, 50
Ef = (Eup + Elw) / 2.0
sumkmid = sumkg (et, nbndx, nbnd, nks, wk, Degauss, Ngauss, Ef)
sumkmid = sumkg (et, nbnd, nks, wk, Degauss, Ngauss, Ef)
if (abs (sumkmid-nelec) .lt.1.0e-10) then
return
elseif ( (sumkmid-nelec) .lt. - 1.0e-10) then

18
PW/efermit.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,7 +7,7 @@
!
!
!--------------------------------------------------------------------
subroutine efermit (et, nbndx, nbnd, nks, nelec, nspin, ntetra, &
subroutine efermit (et, nbnd, nks, nelec, nspin, ntetra, &
tetra, ef)
!--------------------------------------------------------------------
!
@ -15,14 +15,13 @@ subroutine efermit (et, nbndx, nbnd, nks, nelec, nspin, ntetra, &
!
use parameters
implicit none
integer :: nks, nbndx, nbnd, nspin, ntetra, tetra (4, ntetra)
integer :: nks, nbnd, nspin, ntetra, tetra (4, ntetra)
! input: the number of k points
! input: the maximum number of bands
! input: the number of bands
! input: the number of spin components
! input: the number of tetrahedra
! input: the vertices of a tetrahedron
real(kind=DP) :: et (nbndx, nks), nelec, ef
real(kind=DP) :: et (nbnd, nks), nelec, ef
! input: the eigenvalues
! input: the number of electrons
! output: the fermi energy
@ -67,14 +66,14 @@ subroutine efermit (et, nbndx, nbnd, nks, nelec, nspin, ntetra, &
!
! Bisection method
!
sumkup = sumkt (et, nbndx, nbnd, nks, nspin, ntetra, tetra, eup)
sumklw = sumkt (et, nbndx, nbnd, nks, nspin, ntetra, tetra, elw)
sumkup = sumkt (et, nbnd, nks, nspin, ntetra, tetra, eup)
sumklw = sumkt (et, nbnd, nks, nspin, ntetra, tetra, elw)
if ( (sumkup - nelec) .lt. - eps.or. (sumklw - nelec) .gt.eps) &
call errore ('efermit', 'unexpected error', 1)
better = 1.0e+10
do iter = 1, maxiter
ef = (eup + elw) / 2.0
sumkmid = sumkt (et, nbndx, nbnd, nks, nspin, ntetra, tetra, ef)
sumkmid = sumkt (et, nbnd, nks, nspin, ntetra, tetra, ef)
if (abs (sumkmid-nelec) .lt.better) then
better = abs (sumkmid-nelec)
efbetter = ef
@ -93,8 +92,7 @@ subroutine efermit (et, nbndx, nbnd, nks, nelec, nspin, ntetra, &
! unconverged exit:
! the best available ef is used . Needed in some difficult cases
ef = efbetter
sumkmid = sumkt (et, nbndx, nbnd, nks, nspin, ntetra, tetra, ef)
sumkmid = sumkt (et, nbnd, nks, nspin, ntetra, tetra, ef)
write (6, 9010) ef * rydtoev, sumkmid
! converged exit:

6
PW/electrons.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -129,7 +129,7 @@ implicit none
if (.not.lscf) then
conv_elec=.true.
#ifdef __PARA
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
#endif
do ik = 1, nkstot
@ -259,7 +259,7 @@ implicit none
enddo
call ireduce (nks, ngkp)
call ipoolrecover (ngkp, 1, nkstot, nks)
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
#endif
do ik = 1, nkstot

22
PW/gweights.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,34 +7,28 @@
!
!
!--------------------------------------------------------------------
subroutine gweights (nks, wk, nbndx, nbnd, nelec, degauss, ngauss, &
subroutine gweights (nks, wk, nbnd, nelec, degauss, ngauss, &
et, ef, demet, wg)
!--------------------------------------------------------------------
! calculates weights with the gaussian spreading technique
use parameters
implicit none
!
integer :: nks, nbndx, nbnd, ngauss
real(kind=DP) :: wk (nks), et (nbndx, nks), nelec, degauss
integer :: nks, nbnd, ngauss
real(kind=DP) :: wk (nks), et (nbnd, nks), nelec, degauss
real(kind=DP) :: wg (nbnd, nks), ef, demet
real(kind=DP) :: wgauss, w1gauss
!
integer :: kpoint, ibnd
real(kind=DP) , external :: wgauss, w1gauss
external efermig, wgauss, w1gauss
! Calculate the Fermi energy ef
call efermig (et, nbndx, nbnd, nks, nelec, wk, degauss, ngauss, &
ef)
call efermig (et, nbnd, nks, nelec, wk, degauss, ngauss, ef)
demet = 0.d0
do kpoint = 1, nks
do ibnd = 1, nbnd
! Calculate the gaussian weights
! Calculate the gaussian weights
wg (ibnd, kpoint) = wk (kpoint) * wgauss ( (ef - et (ibnd, kpoint) &
) / degauss, ngauss)

4
PW/iweights.f90
View File

@ -17,11 +17,9 @@ subroutine iweights (nks, wk, nbnd, nelec, wg)
integer :: nks, nbnd
real(kind=DP) :: wk (nks), nelec
real(kind=DP) :: wg (nbnd, nks)
real(kind=DP) :: degspin
real(kind=DP), parameter :: degspin = 2.d0
integer :: kpoint, ibnd
parameter (degspin = 2.d0)
do kpoint = 1, nks
do ibnd = 1, nbnd
if (ibnd.le.nint (nelec) / degspin) then

8
PW/punch.f90
View File

@ -18,7 +18,7 @@ subroutine punch
!
!
use pwcom, only: nks, filpun, reduce_io, evc, nwordwfc, iunwfc, lscf, &
rho, nspin, iunpun, et, wg, nbnd, nkstot, nbndx
rho, nspin, iunpun, et, wg, nbnd, nkstot
use io, only: prefix
#ifdef __PARA
use para
@ -57,8 +57,8 @@ subroutine punch
! while eigenvalues et and weights wg must be
! explicitely collected to the first node
!
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (wg, nbnd , nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
call poolrecover (wg, nbnd, nkstot, nks)
!
! In parallel execution, only the first node writes this file
!
@ -127,7 +127,7 @@ subroutine punch
! while eigenvalues et and weights wg must be
! explicitely collected to the first node
!
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
call poolrecover (wg, nbnd , nkstot, nks)
!
! In parallel execution, only the first node writes this file

10
PW/sum_band.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -51,11 +51,10 @@ subroutine sum_band
!
elseif (ltetra) then
#ifdef __PARA
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
if (me.eq.1.and.mypool.eq.1) then
#endif
call tweights (nkstot, nspin, nbndx, nbnd, nelec, ntetra, &
tetra, et, ef, wg)
call tweights (nkstot, nspin, nbnd, nelec, ntetra, tetra, et, ef, wg)
#ifdef __PARA
endif
call poolscatter (nbnd, nkstot, wg, nks, wg)
@ -63,8 +62,7 @@ subroutine sum_band
call broadcast (1, ef)
#endif
elseif (lgauss) then
call gweights (nks, wk, nbndx, nbnd, nelec, degauss, ngauss, &
et, ef, demet, wg)
call gweights (nks, wk, nbnd, nelec, degauss, ngauss, et, ef, demet, wg)
endif
!
! Needed for LDA+U

84
PW/sumkg.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,49 +7,47 @@
!
!
!-----------------------------------------------------------------------
function sumkg (et, nbndx, nbnd, nks, wk, degauss, ngauss, e)
!-----------------------------------------------------------------------
!
! This function computes the number of states under a given energy e
!
!
use parameters
implicit none
!
! First the dummy variables.
!
integer :: nks, nbndx, nbnd, ngauss
! input: the total number of K points
! input: the maximum number of bands
! input: the number of bands
! input: the type of smearing
real(kind=DP) :: wk (nks), et (nbndx, nks), degauss, sumkg, e
! input: the weight of the k points
!input: the energy eigenvalues
! input: the smearing function
! output: the sum of the eigenvalues
! input: the energy to check
!
! the local variables
!
real(kind=DP) :: wgauss, sum1
! function which compute the smearing
! auxiliary function
integer :: ik, ibnd
! counter on k points
! counter on the band energy
!
sumkg = 0.d0
do ik = 1, nks
sum1 = 0.d0
do ibnd = 1, nbnd
sum1 = sum1 + wgauss ( (e-et (ibnd, ik) ) / degauss, ngauss)
enddo
sumkg = sumkg + wk (ik) * sum1
enddo
function sumkg (et, nbnd, nks, wk, degauss, ngauss, e)
!-----------------------------------------------------------------------
!
! This function computes the number of states under a given energy e
!
!
use parameters
implicit none
!
real(kind=DP) :: sumkg
!
integer :: nks, nbnd, ngauss
! input: the total number of K points
! input: the number of bands
! input: the type of smearing
real(kind=DP) :: wk (nks), et (nbnd, nks), degauss, e
! input: the weight of the k points
! input: the energy eigenvalues
! input: gaussian broadening
! input: the energy to check
!
! the local variables
!
real(kind=DP), external :: wgauss
! function which compute the smearing
real(kind=DP) ::sum1
integer :: ik, ibnd
! counter on k points
! counter on the band energy
!
sumkg = 0.d0
do ik = 1, nks
sum1 = 0.d0
do ibnd = 1, nbnd
sum1 = sum1 + wgauss ( (e-et (ibnd, ik) ) / degauss, ngauss)
enddo
sumkg = sumkg + wk (ik) * sum1
enddo
#ifdef __PARA
call poolreduce (1, sumkg)
call poolreduce (1, sumkg)
#endif
return
return
end function sumkg

16
PW/sumkt.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,19 +7,17 @@
!
!
!--------------------------------------------------------------------
real(8) function sumkt (et, nbndx, nbnd, nks, nspin, ntetra, &
tetra, e)
function sumkt (et, nbnd, nks, nspin, ntetra, tetra, e)
!--------------------------------------------------------------------
!
use parameters
implicit none
integer :: nbndx, nbnd, nks, nspin, ntetra, tetra (4, ntetra)
real(kind=DP) :: et (nbndx, nks), e
integer :: nt, nk, ns, ibnd, i
real(kind=DP) :: sumkt
integer :: nbnd, nks, nspin, ntetra, tetra (4, ntetra)
real(kind=DP) :: et (nbnd, nks), e
!
real(kind=DP) :: etetra (4), e1, e2, e3, e4
integer :: nt, nk, ns, ibnd, i
sumkt = 0.0
do ns = 1, nspin

19
PW/tweights.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -7,34 +7,29 @@
!
!
!--------------------------------------------------------------------
subroutine tweights (nks, nspin, nbndx, nbnd, nelec, ntetra, &
subroutine tweights (nks, nspin, nbnd, nelec, ntetra, &
tetra, et, ef, wg)
!--------------------------------------------------------------------
! calculates weights with the tetrahedron method (Bloechl version)
use parameters
implicit none
!
integer :: nks, nspin, nbndx, nbnd, ntetra, tetra (4, ntetra)
real(kind=DP) :: et (nbndx, nks), nelec
integer :: nks, nspin, nbnd, ntetra, tetra (4, ntetra)
real(kind=DP) :: et (nbnd, nks), nelec
real(kind=DP) :: wg (nbnd, nks), ef
real(kind=DP) :: e1, e2, e3, e4, c1, c2, c3, c4, etetra (4), dosef
integer :: ik, ibnd, nt, nk, ns, i, kp1, kp2, kp3, kp4, itetra (4)
! Calculate the Fermi energy ef
call efermit (et, nbndx, nbnd, nks, nelec, nspin, ntetra, tetra, &
ef)
call efermit (et, nbnd, nks, nelec, nspin, ntetra, tetra, ef)
do ik = 1, nks
do ibnd = 1, nbnd
wg (ibnd, ik) = 0.d0
enddo
enddo
do ns = 1, nspin
!
! nk is used to select k-points with up (ns=1) or down (ns=2) spin

9
PW/wfcinit.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -24,10 +24,9 @@ subroutine wfcinit
! " " polarization
! number of starting wavefunctions
real(kind=DP) :: rndm, rr, arg
real(kind=DP) ::rr, arg
real(kind=DP), external :: rndm
! random function generation
external rndm
!
! state what is going to happen
!
@ -128,7 +127,7 @@ subroutine wfcinit
enddo
if (iprint.eq.1) then
#ifdef __PARA
call poolrecover (et, nbndx, nkstot, nks)
call poolrecover (et, nbnd, nkstot, nks)
#endif
do ik = 1, nkstot
write (6, 9010) (xk (ipol, ik), ipol = 1, 3)

23
TODO
View File

@ -1,8 +1,8 @@
TODO LIST - 2 Apr 2003
TODO LIST - 7 Apr 2003
INSTALLATION
- compile fftw by default on request for PW and CP as well
- compile fftw by default for PW and CP as well
COMMON
@ -25,15 +25,16 @@ COMMON
PW
- wavefunctions should be dimensioned nbnd, not nbndx
- bfgs, md : atomic positions should be written in the same
format as they are read (but this should not spoil scripts
that extract coordinates from output file)
- remove analytical pps and related variables
- remove analytical PPs and related variables
- remove residual direct calls to MPI routines
use routines in mp.f90 instead
merge with existing wrapper if needed
- remove residual direct calls to MPI routines,
use (or merge with) existing routines in mp.f90 instead
- remove potential mixing, save and start from rho instead of V
@ -55,6 +56,8 @@ PW
- memory estimator tool, or memory report
- ultrasoft PP: newd and addusdens are too slow
POSTPROCESSING
- bands.x must either be finished or removed
@ -64,7 +67,7 @@ POSTPROCESSING
- stm in non-scf calculations to be verified
- add more scripts that process output file
- add more scripts that process output files
- dos, projected dos, etc: input data should be more uniform
@ -72,12 +75,15 @@ PH
- Tone: ntyp in input needed for phonon GUI ?
- better algorithm for electron-phonon (Malgorzata?)
- better algorithm for electron-phonon (Malgorzata)
DOCUMENTATION
- expand, update, verify, etc
- examples for many features are missing
examples should be quicker and easier to verify
- add a list of FAQ, or AAQ (Already Answered Questions)
CPV:
@ -86,6 +92,7 @@ CPV:
with ultrasoft pps and box grids. If it is a bug, it is serious.
- check on input ipp and pp type. Even better: remove ipp
(leave either UPF or old-format Vanderbilt)
- remove eispack routines, replace with lapack or parallel diag

6
pwtools/dynmat.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 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,
@ -180,7 +180,7 @@
end if
if (asr) then
!
! ASR sulle cariche efficaci
! ASR on effective charges
!
do i=1,3
do j=1,3
@ -194,7 +194,7 @@
end do
end do
!
! ASR sulla matrice dinamica
! ASR on dynamical matrix
!
do i=1,3
do j=1,3