mirror of https://gitlab.com/QEF/q-e.git
Introdotto il calcolo del dipolo di cariche localizzate al centro della
cella unitaria in chdens (utile solo per molecole isolate). Introdotta la possibilita' di aggiungere un potenziale a forma di dente di sega al potenziale degli ioni per simulare un campo elettrico finito. (by J. Tobik) git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@141 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
dalcorso
parent
139c026205
commit
513f48608f
|
|
@ -108,6 +108,7 @@ MODULES = ../Modules/*.o
|
|||
|
||||
PWOBJS = ../PW/pwcom.o \
|
||||
../PW/aainit.o \
|
||||
../PW/add_efield.o \
|
||||
../PW/addusdens.o \
|
||||
../PW/addusforce.o \
|
||||
../PW/addusstress.o \
|
||||
|
|
|
|||
|
|
@ -27,12 +27,14 @@ stm.o \
|
|||
stop_pp.o \
|
||||
subtract_vxc.o \
|
||||
voronoy.o \
|
||||
work_function.o
|
||||
work_function.o \
|
||||
wsweight.o
|
||||
|
||||
MODULES = ../Modules/*.o
|
||||
|
||||
PWOBJS = ../PW/pwcom.o \
|
||||
../PW/aainit.o \
|
||||
../PW/add_efield.o \
|
||||
../PW/addusdens.o \
|
||||
../PW/addusforce.o \
|
||||
../PW/addusstress.o \
|
||||
|
|
|
|||
260
PP/chdens.f90
260
PP/chdens.f90
|
|
@ -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,
|
||||
|
|
@ -70,6 +70,18 @@ subroutine do_chdens
|
|||
! is written (in postproc format) for further processing
|
||||
! (macroscopic average)
|
||||
!
|
||||
!----IF iflag=3,30 and plot_out=1
|
||||
!
|
||||
! idpol 1 the ionic and electronic dipole moment of the charge
|
||||
! is computed.
|
||||
! 2 only the electronic dipole moment is computed
|
||||
!
|
||||
! NB: This option is to be used for an isolated molecule in
|
||||
! a box and the molecule must be at the center of the box.
|
||||
! The code computes the dipole on the Wigner-Seitz cell of
|
||||
! the Bravais lattice. The 3d box must contain this cell
|
||||
! otherwise meaningless numbers are printed.
|
||||
!
|
||||
!----END_IF
|
||||
!
|
||||
!-/ END of namelist &input
|
||||
|
|
@ -146,16 +158,17 @@ subroutine do_chdens
|
|||
real(kind=DP), allocatable :: taus (:,:)
|
||||
integer :: ibravs, nrx1sa, nrx2sa, nrx3sa, nr1sa, nr2sa, nr3sa, &
|
||||
ntyps, nats
|
||||
integer :: idpol ! dipol moment flag
|
||||
integer, allocatable :: ityps (:)
|
||||
character (len=3) :: atms(ntypx)
|
||||
character (len=80) :: filepp(nfilemax)
|
||||
real(kind=DP) :: rhodum
|
||||
real(kind=DP) :: rhodum, dipol(0:3)
|
||||
complex(kind=DP), allocatable:: vgc (:)
|
||||
! rho or polarization in G space
|
||||
logical :: fast3d
|
||||
|
||||
namelist /input/ &
|
||||
nfile, filepp, weight, iflag, &
|
||||
nfile, filepp, weight, iflag, idpol, &
|
||||
plot_out, output_format, fileout, epsilon, filepol
|
||||
|
||||
!
|
||||
|
|
@ -170,8 +183,9 @@ subroutine do_chdens
|
|||
output_format = 0
|
||||
fileout = ' '
|
||||
epsilon = 1.0d0
|
||||
idpol = 0
|
||||
filepol = ' '
|
||||
fast3d = .false.
|
||||
fast3d = .false.
|
||||
|
||||
!
|
||||
! read and check input data
|
||||
|
|
@ -188,6 +202,7 @@ subroutine do_chdens
|
|||
if (nfile.le.0.or.nfile.gt.nfilemax) &
|
||||
call errore ('chdens ', 'nfile is wrong ', 1)
|
||||
|
||||
|
||||
! check for iflag
|
||||
if (iflag.eq.30) then
|
||||
iflag = 3
|
||||
|
|
@ -196,6 +211,17 @@ subroutine do_chdens
|
|||
if ((iflag.lt.1.or.iflag.gt.4) .and. iflag.ne.31) &
|
||||
call errore ('chdens', 'iflag not implemented', 1)
|
||||
|
||||
! check for idpol
|
||||
if (idpol.eq.1.or.idpol.eq.2) then
|
||||
if (iflag.ne.3) then
|
||||
idpol=0
|
||||
call errore("chdens","dipole computed only if iflag=3 or 30",-1)
|
||||
endif
|
||||
if (plot_out.ne.1) then
|
||||
idpol=0
|
||||
call errore("chdens","dipole computed only if plot_out=1",-1)
|
||||
endif
|
||||
endif
|
||||
|
||||
! reading the rest of input (spanning vectors, origin, number-of points)
|
||||
if (iflag.lt.4) then
|
||||
|
|
@ -480,12 +506,15 @@ subroutine do_chdens
|
|||
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, + 1)
|
||||
call plot_fast (celldm (1), at, nat, tau, atm, ityp, &
|
||||
nrx1, nrx2, nrx3, nr1, nr2, nr3, psic, &
|
||||
bg, m1, m2, m3, x0, e, output_format, ounit)
|
||||
bg, m1, m2, m3, x0, e, output_format, ounit, dipol(0))
|
||||
else
|
||||
call plot_3d (celldm (1), at, nat, tau, atm, ityp, ngm, g, vgc, &
|
||||
nx, ny, nz, m1, m2, m3, x0, e, output_format, ounit)
|
||||
nx, ny, nz, m1, m2, m3, x0, e, output_format, ounit, dipol(0))
|
||||
end if
|
||||
|
||||
if (idpol.eq.1.or.idpol.eq.2) &
|
||||
call write_dipol(dipol(0),tau,nat,alat,zv,ntyp,ityp,idpol)
|
||||
|
||||
elseif (iflag.eq.4) then
|
||||
radius = radius / alat
|
||||
|
||||
|
|
@ -591,8 +620,8 @@ subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
|
|||
|
||||
rhoim = 0.d0
|
||||
do i = 1, nx
|
||||
rhomin = min (rhomin, dreal (carica (i) ) )
|
||||
rhomax = max (rhomax, dreal (carica (i) ) )
|
||||
rhomin = min (rhomin, DREAL (carica (i) ) )
|
||||
rhomax = max (rhomax, DREAL (carica (i) ) )
|
||||
rhoim = rhoim + abs (DIMAG (carica (i) ) )
|
||||
|
||||
enddo
|
||||
|
|
@ -702,8 +731,8 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
|
|||
rhoim = 0.d0
|
||||
do i = 1, nx
|
||||
do j = 1, ny
|
||||
rhomin = min (rhomin, dreal (carica (i, j) ) )
|
||||
rhomax = max (rhomax, dreal (carica (i, j) ) )
|
||||
rhomin = min (rhomin, DREAL (carica (i, j) ) )
|
||||
rhomax = max (rhomax, DREAL (carica (i, j) ) )
|
||||
rhoim = rhoim + abs (dimag (carica (i, j) ) )
|
||||
enddo
|
||||
|
||||
|
|
@ -723,7 +752,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
|
|||
!
|
||||
! write(ounit,'(2i6)') nx,ny
|
||||
do i = 1, nx
|
||||
write (ounit, '(e25.14)') (dreal (carica (i, j) ) , j = 1, ny)
|
||||
write (ounit, '(e25.14)') (DREAL (carica (i, j) ) , j = 1, ny)
|
||||
write (ounit, * )
|
||||
enddo
|
||||
elseif (output_format.eq.1) then
|
||||
|
|
@ -731,7 +760,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
|
|||
! contour.x format
|
||||
!
|
||||
write (ounit, '(3i5,2e25.14)') nx, ny, 1, deltax, deltay
|
||||
write (ounit, '(4e25.14)') ( (dreal (carica (i, j) ) , j = 1, &
|
||||
write (ounit, '(4e25.14)') ( (DREAL (carica (i, j) ) , j = 1, &
|
||||
ny) , i = 1, nx)
|
||||
elseif (output_format.eq.2) then
|
||||
!
|
||||
|
|
@ -740,7 +769,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
|
|||
write (ounit, '(2i4)') nx - 1, ny - 1
|
||||
write (ounit, '(8f8.4)') (deltax * (i - 1) , i = 1, nx)
|
||||
write (ounit, '(8f8.4)') (deltay * (j - 1) , j = 1, ny)
|
||||
write (ounit, '(6e12.4)') ( (dreal (carica (i, j) ) , i = 1, &
|
||||
write (ounit, '(6e12.4)') ( (DREAL (carica (i, j) ) , i = 1, &
|
||||
nx) , j = 1, ny)
|
||||
write (ounit, '(3f8.4)') x0
|
||||
write (ounit, '(3f8.4)') (m1 * e (i, 1) , i = 1, 3)
|
||||
|
|
@ -838,8 +867,8 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
|
|||
rhoim = 0.d0
|
||||
do i = 1, nx
|
||||
do j = 1, ny
|
||||
rhomin = min (rhomin, dreal (carica (i, j) ) )
|
||||
rhomax = max (rhomax, dreal (carica (i, j) ) )
|
||||
rhomin = min (rhomin, DREAL (carica (i, j) ) )
|
||||
rhomax = max (rhomax, DREAL (carica (i, j) ) )
|
||||
rhoim = rhoim + abs (dimag (carica (i, j) ) )
|
||||
enddo
|
||||
|
||||
|
|
@ -857,14 +886,14 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
|
|||
!
|
||||
write (ounit, '(2i8)') nx, ny
|
||||
do i = 1, nx
|
||||
write (ounit, '(e25.14)') (dreal (carica (i, j) ) , j = 1, ny)
|
||||
write (ounit, '(e25.14)') (DREAL (carica (i, j) ) , j = 1, ny)
|
||||
enddo
|
||||
elseif (output_format.eq.1) then
|
||||
!
|
||||
! contour.x format
|
||||
!
|
||||
write (ounit, '(3i5,2e25.14)') nx, ny, 1, deltax, deltay
|
||||
write (ounit, '(4e25.14)') ( (dreal (carica (i, j) ) , j = 1, &
|
||||
write (ounit, '(4e25.14)') ( (DREAL (carica (i, j) ) , j = 1, &
|
||||
ny) , i = 1, nx)
|
||||
else
|
||||
call errore ('plot_2ds', 'not implemented plot', 1)
|
||||
|
|
@ -878,7 +907,7 @@ end subroutine plot_2ds
|
|||
!
|
||||
!-----------------------------------------------------------------------
|
||||
subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
||||
nx, ny, nz, m1, m2, m3, x0, e, output_format, ounit)
|
||||
nx, ny, nz, m1, m2, m3, x0, e, output_format, ounit, dipol)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
use parameters, only : DP
|
||||
|
|
@ -893,7 +922,7 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
|||
character(len=3) :: atm(*)
|
||||
|
||||
real(kind=DP) :: alat, tau (3, nat), at (3, 3), g (3, ngm), e (3, 3), &
|
||||
x0 (3), m1, m2, m3
|
||||
x0 (3), m1, m2, m3, dipol(0:3)
|
||||
! lattice parameter
|
||||
! atomic positions
|
||||
! lattice vectors
|
||||
|
|
@ -911,12 +940,15 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
|||
! steps along e1, e2, e3
|
||||
real(kind=DP), parameter :: pi = 3.14159265358979d0
|
||||
complex(kind=DP), allocatable :: eigx (:), eigy (:), eigz (:)
|
||||
real(kind=DP), allocatable :: carica (:,:,:)
|
||||
real(kind=DP), allocatable :: carica (:,:,:), fact(:,:,:), rws(:,:)
|
||||
real(kind=dp) :: wsweight, r(3), rijk, omega, suma
|
||||
integer :: ipol, na, nrwsx, nrws
|
||||
|
||||
allocate (eigx( nx))
|
||||
allocate (eigy( ny))
|
||||
allocate (eigz( nz))
|
||||
allocate (carica( nx , ny , nz))
|
||||
allocate (fact( nx , ny , nz))
|
||||
|
||||
deltax = m1 / (nx - 1)
|
||||
deltay = m2 / (ny - 1)
|
||||
|
|
@ -929,23 +961,23 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
|||
! These factors are calculated and stored in order to save CPU time
|
||||
!
|
||||
do i = 1, nx
|
||||
eigx (i) = exp ( (0.d0, 1.d0) * 2.d0 * pi * ( (i - 1) * deltax * &
|
||||
(e (1, 1) * g (1, ig) + e (2, 1) * g (2, ig) + e (3, 1) * g (3, ig)) &
|
||||
+ (x0 (1) * g (1, ig) + x0 (2) * g (2, ig) + x0 (3) * g (3, ig) ) ) )
|
||||
eigx (i) = exp((0.d0, 1.d0)*2.d0*pi*((i - 1) * deltax * &
|
||||
(e(1,1)*g(1,ig)+e(2,1)*g(2,ig)+e(3,1)*g(3,ig)) &
|
||||
+(x0(1)*g(1,ig)+x0(2)*g(2,ig)+x0(3)*g(3,ig) ) ) )
|
||||
enddo
|
||||
do j = 1, ny
|
||||
eigy (j) = exp ( (0.d0, 1.d0) * 2.d0 * pi * (j - 1) * deltay * &
|
||||
(e (1, 2) * g (1, ig) + e (2, 2) * g (2, ig) + e (3, 2) * g (3, ig) ) )
|
||||
eigy (j) = exp ( (0.d0,1.d0)*2.d0*pi*(j-1)*deltay* &
|
||||
(e(1,2)*g(1,ig)+e(2,2)*g(2,ig)+e(3,2)*g(3,ig) ) )
|
||||
enddo
|
||||
do k = 1, nz
|
||||
eigz (k) = exp ( (0.d0, 1.d0) * 2.d0 * pi * (k - 1) * deltaz * &
|
||||
(e (1, 3) * g (1, ig) + e (2, 3) * g (2, ig) + e (3, 3) * g (3, ig) ) )
|
||||
eigz (k) = exp ( (0.d0,1.d0)*2.d0*pi*(k-1)*deltaz* &
|
||||
(e(1,3)*g(1,ig)+e(2,3)*g(2,ig)+e(3,3)*g(3,ig)) )
|
||||
enddo
|
||||
do k = 1, nz
|
||||
do j = 1, ny
|
||||
do i = 1, nx
|
||||
carica (i, j, k) = carica (i, j, k) + &
|
||||
dreal (vgc (ig) * eigz (k) * eigy (j) * eigx (i) )
|
||||
DREAL (vgc (ig) * eigz (k) * eigy (j) * eigx (i) )
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
|
@ -953,11 +985,35 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
|||
enddo
|
||||
!
|
||||
! Here we check the value of the resulting charge
|
||||
! and compute the dipole of the charge
|
||||
!
|
||||
nrwsx=125
|
||||
allocate(rws(0:3,nrwsx))
|
||||
call wsinit(rws,nrwsx,nrws,at)
|
||||
|
||||
fact=0.d0
|
||||
suma=0.d0
|
||||
do k = 1, nz
|
||||
do j = 1, ny
|
||||
do i = 1, nx
|
||||
do ipol=1,3
|
||||
r(ipol)=x0(ipol)+(i-1)*e(ipol,1)*deltax + &
|
||||
(j-1)*e(ipol,2)*deltay + &
|
||||
(k-1)*e(ipol,3)*deltaz
|
||||
enddo
|
||||
fact(i,j,k)=wsweight(r,rws,nrws)
|
||||
suma=suma+fact(i,j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call volume(alat,at(1,1),at(1,2),at(1,3),omega)
|
||||
|
||||
rhomin = 1.0d10
|
||||
rhomax =-1.0d10
|
||||
rhotot = 0.d0
|
||||
rhoabs = 0.d0
|
||||
dipol = 0.d0
|
||||
do k = 1, nz
|
||||
do j = 1, ny
|
||||
do i = 1, nx
|
||||
|
|
@ -965,12 +1021,22 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
|||
rhomax = max (rhomax, carica (i, j, k) )
|
||||
rhotot = rhotot + carica (i, j, k)
|
||||
rhoabs = rhoabs + abs (carica (i, j, k) )
|
||||
dipol(0) = dipol(0) + fact(i,j,k)*carica (i, j, k)
|
||||
do ipol=1,3
|
||||
rijk=x0(ipol)+(i-1)*e(ipol,1)*deltax + &
|
||||
(j-1)*e(ipol,2)*deltay + &
|
||||
(k-1)*e(ipol,3)*deltaz
|
||||
dipol(ipol)=dipol(ipol)+fact(i,j,k)*rijk*carica(i,j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
rhotot = rhotot / nx / ny / nz * m1 * m2 * m3 * alat**3
|
||||
rhoabs = rhoabs / nx / ny / nz * m1 * m2 * m3 * alat**3
|
||||
do ipol=1,3
|
||||
dipol(ipol)=dipol(ipol) / suma * omega * alat
|
||||
enddo
|
||||
print '(/5x,"Min, Max, Total, Abs charge: ",4f10.6)', rhomin, &
|
||||
rhomax, rhotot, rhoabs
|
||||
|
||||
|
|
@ -993,6 +1059,8 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
|||
endif
|
||||
|
||||
deallocate (carica)
|
||||
deallocate (fact)
|
||||
deallocate (rws)
|
||||
deallocate (eigz)
|
||||
deallocate (eigy)
|
||||
deallocate (eigx)
|
||||
|
|
@ -1000,9 +1068,9 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
|
|||
end subroutine plot_3d
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
subroutine plot_fast (alat, at, nat, tau, atm, ityp, &
|
||||
subroutine plot_fast (alat, at, nat, tau, atm, ityp,&
|
||||
nrx1, nrx2, nrx3, nr1, nr2, nr3, rho, &
|
||||
bg, m1, m2, m3, x0, e, output_format, ounit)
|
||||
bg, m1, m2, m3, x0, e, output_format, ounit, dipol)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
use parameters, only : DP
|
||||
|
|
@ -1012,11 +1080,14 @@ subroutine plot_fast (alat, at, nat, tau, atm, ityp, &
|
|||
character(len=3) :: atm(*)
|
||||
|
||||
real(kind=DP) :: alat, tau (3, nat), at (3, 3), rho(2, nrx1,nrx2,nrx3), &
|
||||
bg (3, 3), e (3, 3), x0 (3), m1, m2, m3
|
||||
bg (3, 3), e (3, 3), x0 (3), m1, m2, m3, dipol(0:3)
|
||||
|
||||
integer :: nx, ny, nz, nx0, ny0, nz0, nx1, ny1, nz1, i, j, k, i1, j1, k1
|
||||
real(kind=DP) :: rhomin, rhomax, rhotot, rhoabs
|
||||
real(kind=DP), allocatable :: carica (:,:,:)
|
||||
real(kind=DP), allocatable :: carica (:,:,:), fact(:,:,:), rws(:,:)
|
||||
real(kind=DP) :: rijk, deltax, deltay, deltaz, debye
|
||||
real(kind=dp) :: wsweight, r(3), omega, suma
|
||||
integer :: ipol, na, nrwsx, nrws
|
||||
|
||||
! find FFT grid point closer to X0 (origin of the parallelepiped)
|
||||
! (add 1 because r=0 correspond to n=1)
|
||||
|
|
@ -1030,18 +1101,19 @@ subroutine plot_fast (alat, at, nat, tau, atm, ityp, &
|
|||
e(1,3) .ne. 0.d0 .or. e(2,3) .ne. 0.d0 ) &
|
||||
call errore ('plot_fast','need vectors along x,y,z',1)
|
||||
|
||||
! find FFT grid points closer to X0 + e1, X0 + e2, X(0 + e3
|
||||
! find FFT grid points closer to X0 + e1, X0 + e2, X0 + e3
|
||||
! (the opposite vertex of the parallelepiped)
|
||||
|
||||
nx1 = nint ((x0(1)+m1)*bg(1,1)*nr1 + x0(2)*bg(2,1)*nr1 + x0(3)*bg(3,1)*nr1 ) + 1
|
||||
ny1 = nint ( x0(1)*bg(1,2)*nr2 +(x0(2)+m2)*bg(2,2)*nr2 + x0(3)*bg(3,2)*nr2 ) + 1
|
||||
nz1 = nint ( x0(1)*bg(1,3)*nr3 + x0(2)*bg(2,3)*nr3 +(x0(3)+m3)*bg(3,3)*nr3 ) + 1
|
||||
nx1 = nint ((x0(1)+m1)*bg(1,1)*nr1+x0(2)*bg(2,1)*nr1+x0(3)*bg(3,1)*nr1)+1
|
||||
ny1 = nint (x0(1)*bg(1,2)*nr2+(x0(2)+m2)*bg(2,2)*nr2+x0(3)*bg(3,2)*nr2)+1
|
||||
nz1 = nint (x0(1)*bg(1,3)*nr3+x0(2)*bg(2,3)*nr3+(x0(3)+m3)*bg(3,3)*nr3)+1
|
||||
|
||||
nx = nx1 - nx0 + 1
|
||||
ny = ny1 - ny0 + 1
|
||||
nz = nz1 - nz0 + 1
|
||||
|
||||
allocate (carica( nx, ny, nz))
|
||||
allocate (fact( nx, ny, nz))
|
||||
|
||||
carica = 0.d0
|
||||
do k = nz0, nz1
|
||||
|
|
@ -1071,17 +1143,45 @@ subroutine plot_fast (alat, at, nat, tau, atm, ityp, &
|
|||
! consistent with the FFT grid
|
||||
!
|
||||
write(6,'(5x,"Requested parallelepiped origin: ",3f8.4)') x0
|
||||
x0(1) = (nx0-1) * at(1,1)/nr1 + (ny0-1) * at(1,2)/nr2 + (nz0-1) * at(1,3)/nr3
|
||||
x0(2) = (nx0-1) * at(2,1)/nr1 + (ny0-1) * at(2,2)/nr2 + (nz0-1) * at(2,3)/nr3
|
||||
x0(3) = (nx0-1) * at(3,1)/nr1 + (ny0-1) * at(3,2)/nr2 + (nz0-1) * at(3,3)/nr3
|
||||
x0(1) = (nx0-1)*at(1,1)/nr1+(ny0-1)*at(1,2)/nr2+(nz0-1)*at(1,3)/nr3
|
||||
x0(2) = (nx0-1)*at(2,1)/nr1+(ny0-1)*at(2,2)/nr2+(nz0-1)*at(2,3)/nr3
|
||||
x0(3) = (nx0-1)*at(3,1)/nr1+(ny0-1)*at(3,2)/nr2+(nz0-1)*at(3,3)/nr3
|
||||
write(6,'(5x,"Redefined parallelepiped origin: ",3f8.4)') x0
|
||||
|
||||
deltax = m1/(nx - 1)
|
||||
deltay = m2/(ny - 1)
|
||||
deltaz = m3/(nz - 1)
|
||||
!
|
||||
! Here we check the value of the resulting charge
|
||||
! and compute the dipole
|
||||
!
|
||||
nrwsx=125
|
||||
allocate(rws(0:3,nrwsx))
|
||||
call wsinit(rws,nrwsx,nrws,at)
|
||||
|
||||
fact=0.d0
|
||||
suma=0.d0
|
||||
do k = 1, nz
|
||||
do j = 1, ny
|
||||
do i = 1, nx
|
||||
do ipol=1,3
|
||||
r(ipol)=x0(ipol)+(i-1)*e(ipol,1)*deltax + &
|
||||
(j-1)*e(ipol,2)*deltay + &
|
||||
(k-1)*e(ipol,3)*deltaz
|
||||
enddo
|
||||
fact(i,j,k)=wsweight(r,rws,nrws)
|
||||
suma=suma+fact(i,j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call volume(alat,at(1,1),at(1,2),at(1,3),omega)
|
||||
|
||||
rhomin = 1.0d10
|
||||
rhomax =-1.0d10
|
||||
rhotot = 0.d0
|
||||
rhoabs = 0.d0
|
||||
dipol=0.d0
|
||||
do k = 1, nz
|
||||
do j = 1, ny
|
||||
do i = 1, nx
|
||||
|
|
@ -1089,20 +1189,35 @@ subroutine plot_fast (alat, at, nat, tau, atm, ityp, &
|
|||
rhomax = max (rhomax, carica (i, j, k) )
|
||||
rhotot = rhotot + carica (i, j, k)
|
||||
rhoabs = rhoabs + abs (carica (i, j, k) )
|
||||
dipol(0) = dipol(0) + fact(i,j,k)*carica (i, j, k)
|
||||
do ipol=1,3
|
||||
rijk=x0(ipol)+(i-1)*e(ipol,1)*deltax + &
|
||||
(j-1)*e(ipol,2)*deltay + &
|
||||
(k-1)*e(ipol,3)*deltaz
|
||||
dipol(ipol)=dipol(ipol)+fact(i,j,k)*rijk*carica(i,j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
rhotot = rhotot / nx / ny / nz * m1 * m2 * m3 * alat**3
|
||||
rhoabs = rhoabs / nx / ny / nz * m1 * m2 * m3 * alat**3
|
||||
rhotot = rhotot / (nx-1) / (ny-1) / (nz-1) * m1 * m2 * m3 * alat**3
|
||||
rhoabs = rhoabs / (nx-1) / (ny-1) / (nz-1) * m1 * m2 * m3 * alat**3
|
||||
dipol(0) = dipol(0) / suma * omega
|
||||
|
||||
if (omega.gt.m1*m2*m3*alat**3) &
|
||||
write(6,*) 'Warning: the box is too small to calculate dipole'
|
||||
|
||||
print '(/5x,"Min, Max, Total, Abs charge: ",4f10.6)', rhomin, &
|
||||
rhomax, rhotot, rhoabs
|
||||
|
||||
do ipol=1,3
|
||||
dipol(ipol)=dipol(ipol) / suma *omega*alat
|
||||
enddo
|
||||
|
||||
if (output_format.eq.4) then
|
||||
!
|
||||
! "gopenmol" file
|
||||
!
|
||||
|
||||
call write_openmol_file (alat, at, nat, tau, atm, ityp, x0, &
|
||||
m1, m2, m3, nx, ny, nz, rhomax, carica, ounit)
|
||||
else
|
||||
|
|
@ -1112,9 +1227,10 @@ subroutine plot_fast (alat, at, nat, tau, atm, ityp, &
|
|||
call xsf_struct (alat, at, nat, tau, atm, ityp, ounit)
|
||||
call xsf_datagrid_3d (carica, nx, ny, nz, m1, m2, m3, x0, e, alat, ounit)
|
||||
endif
|
||||
|
||||
!
|
||||
deallocate (carica)
|
||||
deallocate (fact)
|
||||
deallocate (rws)
|
||||
return
|
||||
|
||||
end subroutine plot_fast
|
||||
|
|
@ -1202,3 +1318,61 @@ subroutine write_openmol_file (alat, at, nat, tau, atm, ityp, x0, &
|
|||
!
|
||||
return
|
||||
end subroutine write_openmol_file
|
||||
|
||||
subroutine write_dipol(dipol,tau,nat,alat,zv,ntyp,ityp,idpol)
|
||||
use parameters, only : dp
|
||||
implicit none
|
||||
|
||||
integer :: nat, ntyp, ityp(nat), idpol
|
||||
real(kind=dp) :: dipol(0:3), tau(3,nat), zv(ntyp), alat
|
||||
|
||||
real(kind=dp) :: debye, dipol_ion(3)
|
||||
|
||||
integer :: na, ipol
|
||||
!
|
||||
! compute ion dipole moments
|
||||
!
|
||||
if (idpol.eq.1) then
|
||||
dipol_ion=0.d0
|
||||
do na=1,nat
|
||||
do ipol=1,3
|
||||
dipol_ion(ipol)=dipol_ion(ipol)+zv(ityp(na))*tau(ipol,na)*alat
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
!
|
||||
! Charge inside the Wigner-Seitz cell
|
||||
!
|
||||
write(6, '(/4x," Charge density inside the Wigner-Seitz cell:",3f14.8," el.")') &
|
||||
dipol(0)
|
||||
|
||||
!
|
||||
! print the electron dipole moment calculated by the plotting 3d routines
|
||||
! A positive dipole goes from the - charge to the + charge.
|
||||
!
|
||||
write(6, '(/4x,"Electrons dipole moments",3f14.8," a.u.")') &
|
||||
(-dipol(ipol),ipol=1,3)
|
||||
!
|
||||
! print the ionic and total dipole moment
|
||||
!
|
||||
if (idpol.eq.1) then
|
||||
write(6, '(4x," Ions dipole moments",3f14.8," a.u.")') &
|
||||
(dipol_ion(ipol),ipol=1,3)
|
||||
write(6,'(4x," Total dipole moments",3f14.8," a.u.")') &
|
||||
((-dipol(ipol)+dipol_ion(ipol)),ipol=1,3)
|
||||
endif
|
||||
!
|
||||
! Print the same information in Debye
|
||||
!
|
||||
debye=2.54176d0
|
||||
|
||||
write(6,'(/4x,"Electrons dipole moments",3f14.8," Debye")') &
|
||||
(-dipol(ipol)*debye,ipol=1,3)
|
||||
if (idpol.eq.1) then
|
||||
write(6,'(4x," Ions dipole moments",3f14.8," Debye")') &
|
||||
(dipol_ion(ipol)*debye,ipol=1,3)
|
||||
write(6,'(4x," Total dipole moments",3f14.8," Debye")') &
|
||||
((-dipol(ipol)+dipol_ion(ipol))*debye,ipol=1,3)
|
||||
endif
|
||||
|
||||
end subroutine write_dipol
|
||||
|
|
|
|||
|
|
@ -0,0 +1,50 @@
|
|||
!
|
||||
!-----------------------------------------------------------------------
|
||||
subroutine wsinit(rws,nrwsx,nrws,atw)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
use parameters, only : DP
|
||||
implicit none
|
||||
integer i, ii, ir, jr, kr, nrws, nrwsx, nx
|
||||
real(kind=dp) rt, eps, rws(0:3,nrwsx), atw(3,3)
|
||||
parameter (eps=1.0d-6,nx=2)
|
||||
ii = 1
|
||||
do ir=-nx,nx
|
||||
do jr=-nx,nx
|
||||
do kr=-nx,nx
|
||||
do i=1,3
|
||||
rws(i,ii) = atw(i,1)*ir + atw(i,2)*jr + atw(i,3)*kr
|
||||
end do
|
||||
rws(0,ii)=rws(1,ii)*rws(1,ii)+rws(2,ii)*rws(2,ii)+ &
|
||||
rws(3,ii)*rws(3,ii)
|
||||
rws(0,ii)=0.5d0*rws(0,ii)
|
||||
if (rws(0,ii).gt.eps) ii = ii + 1
|
||||
if (ii.gt.nrwsx) call errore('wsinit', 'ii.gt.nrwsx',1)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
nrws = ii - 1
|
||||
return
|
||||
end subroutine wsinit
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
function wsweight(r,rws,nrws)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
use parameters, only : dp
|
||||
implicit none
|
||||
integer ir, nreq, nrws
|
||||
real(kind=dp) r(3), rrt, ck, eps, rws(0:3,nrws), wsweight
|
||||
parameter (eps=1.0d-6)
|
||||
!
|
||||
wsweight = 0.d0
|
||||
nreq = 1
|
||||
do ir =1,nrws
|
||||
rrt = r(1)*rws(1,ir) + r(2)*rws(2,ir) + r(3)*rws(3,ir)
|
||||
ck = rrt-rws(0,ir)
|
||||
if ( ck .gt. eps ) return
|
||||
if ( abs(ck) .lt. eps ) nreq = nreq + 1
|
||||
end do
|
||||
wsweight = 1.d0/float(nreq)
|
||||
return
|
||||
end function wsweight
|
||||
|
|
@ -7,6 +7,7 @@ include ../make.sys
|
|||
PWOBJS = pwcom.o \
|
||||
para.o \
|
||||
aainit.o \
|
||||
add_efield.o \
|
||||
addusdens.o \
|
||||
addusforce.o \
|
||||
addusstress.o \
|
||||
|
|
|
|||
|
|
@ -0,0 +1,163 @@
|
|||
! Copyright (C) 2003 J. Tobik
|
||||
! 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 add_efield
|
||||
!--------------------------------------------------------------------------
|
||||
!
|
||||
! This routine adds an electric field to the ionic potential. The
|
||||
! field is made artificially periodic by introducing a saw-tooth
|
||||
! potential. The field is parallel to a reciprocal lattice vector bg,
|
||||
! according to the index edir.
|
||||
!
|
||||
!
|
||||
#include "machine.h"
|
||||
use pwcom
|
||||
#ifdef __PARA
|
||||
use para
|
||||
#endif
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: npoints, nmax, ndesc
|
||||
integer :: ii, ij, ik, itmp, ir, izlb, izub
|
||||
real(kind=dp) :: length, vamp, value
|
||||
real(kind=dp), parameter :: eps=1.d-8
|
||||
|
||||
#ifndef __PARA
|
||||
integer me, npp(1)
|
||||
me=1
|
||||
npp(1)=nr3
|
||||
#endif
|
||||
if(edir.eq.1) then
|
||||
npoints=nr1
|
||||
length=alat/sqrt(bg(1,1)**2+bg(2,1)**2+bg(3,1)**2)
|
||||
else if (edir.eq.2) then
|
||||
npoints=nr2
|
||||
length=alat/sqrt(bg(1,2)**2+bg(2,2)**2+bg(3,2)**2)
|
||||
elseif (edir.eq.3) then
|
||||
npoints=nr3
|
||||
length=alat/sqrt(bg(1,3)**2+bg(2,3)**2+bg(3,3)**2)
|
||||
else
|
||||
call errore('setlocal',' wrong edir',1)
|
||||
endif
|
||||
|
||||
nmax =int(real(npoints,dp)*(emaxpos-eps))+1
|
||||
if (nmax.lt.1.or.nmax.gt.npoints) &
|
||||
call errore('setlocal','nmax out of range',1)
|
||||
|
||||
ndesc=int(real(npoints,dp)*(eopreg-eps))+1
|
||||
if (ndesc.lt.1.or.ndesc.gt.npoints) &
|
||||
call errore('setlocal','ndesc out of range',1)
|
||||
!
|
||||
! The electric field is assumed in a.u.( 1 a.u. of field change the
|
||||
! potential energy of an electron of 1 Hartree in a distance of 1 Bohr.
|
||||
! The factor 2 converts potential energy to Ry.
|
||||
!
|
||||
vamp=2.0d0*eamp*length*real(npoints-ndesc,dp)/real(npoints,dp)
|
||||
|
||||
write(6,*)
|
||||
write(6,*) 'Adding an external electric field'
|
||||
write(6,*) 'Intensity [a.u.]: ', eamp
|
||||
write(6,*) 'Amplitude vamp [Ry]', vamp
|
||||
write(6,*) 'Total length [points] ', npoints
|
||||
write(6,*) 'Total length [bohr rad] ', length
|
||||
write(6,*) 'Field is reversed between points ', nmax, nmax+ndesc
|
||||
!
|
||||
! in this case in x direction
|
||||
!
|
||||
if(edir.eq.1) then
|
||||
do ij=1,nr2
|
||||
do ik=1,npp(me)
|
||||
do ii=nmax,nmax+ndesc-1
|
||||
value=vamp*(real(nmax+ndesc-ii,dp)/real(ndesc,dp)-0.5d0)
|
||||
itmp=ii
|
||||
if (itmp.gt.nr1) itmp=itmp-nr1
|
||||
ir=itmp+(ij-1)*nrx1+(ik-1)*nrx1*nrx2
|
||||
vltot(ir)=vltot(ir)+value
|
||||
end do
|
||||
do ii=nmax+ndesc,nmax+nr1-1
|
||||
value=vamp*(real(ii-nmax-ndesc,dp)/real(nr1-ndesc,dp)-0.5d0)
|
||||
itmp=ii
|
||||
if (itmp.gt.nr1) itmp=itmp-nr1
|
||||
ir=itmp+(ij-1)*nrx1+(ik-1)*nrx1*nrx2
|
||||
vltot(ir)=vltot(ir)+value
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!
|
||||
! in this case in y direction
|
||||
!
|
||||
else if (edir.eq.2) then
|
||||
do ii=1,nr1
|
||||
do ik=1,npp(me)
|
||||
do ij=nmax,nmax+ndesc-1
|
||||
value=vamp*(real(nmax+ndesc-ij,dp)/real(ndesc,dp)-0.5d0)
|
||||
itmp=ij
|
||||
if (itmp.gt.nr2) itmp=itmp-nr2
|
||||
ir=ii+(itmp-1)*nrx1+(ik-1)*nrx1*nrx2
|
||||
vltot(ir)=vltot(ir)+value
|
||||
end do
|
||||
do ij=nmax+ndesc,nmax+nr2-1
|
||||
value=vamp*(real(ij-nmax-ndesc,dp)/real(nr2-ndesc,dp)-0.5d0)
|
||||
itmp=ij
|
||||
if (itmp.gt.nr2) itmp=itmp-nr2
|
||||
ir=ii+(itmp-1)*nrx1+(ik-1)*nrx1*nrx2
|
||||
vltot(ir)=vltot(ir)+value
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!
|
||||
! and in other cases in z direction
|
||||
!
|
||||
elseif (edir.eq.3) then
|
||||
#ifdef __PARA
|
||||
izub=0
|
||||
do itmp=1,me
|
||||
izlb=izub+1
|
||||
izub=izub+npp(itmp)
|
||||
end do
|
||||
#else
|
||||
izlb=1
|
||||
izub=nr3
|
||||
#endif
|
||||
!
|
||||
! now we have set up boundaries - let's calculate potential
|
||||
!
|
||||
do ii=1,nr1
|
||||
do ij=1,nr2
|
||||
do ik=nmax,nmax+ndesc-1
|
||||
value=vamp*(real(nmax+ndesc-ik,dp)/real(ndesc,dp)-0.5d0)
|
||||
itmp=ik
|
||||
if (itmp.gt.nr3) itmp=itmp-nr3
|
||||
if((itmp.ge.izlb).and.(itmp.le.izub)) then
|
||||
!
|
||||
! Yes - this point belongs to me
|
||||
!
|
||||
itmp=itmp-izlb+1
|
||||
ir=ii+(ij-1)*nrx1+(itmp-1)*nrx1*nrx2
|
||||
vltot(ir)=vltot(ir)+value
|
||||
end if
|
||||
end do
|
||||
do ik=nmax+ndesc,nmax+nr3-1
|
||||
value=vamp*(real(ik-nmax-ndesc,dp)/real(nr3-ndesc,dp)-0.5d0)
|
||||
itmp=ik
|
||||
if (itmp.gt.nr3) itmp=itmp-nr3
|
||||
if((itmp.ge.izlb).and.(itmp.le.izub)) then
|
||||
itmp=itmp-izlb+1
|
||||
ir=ii+(ij-1)*nrx1+(itmp-1)*nrx1*nrx2
|
||||
vltot(ir)=vltot(ir)+value
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
else
|
||||
call errore('setlocal', 'wrong edir', 1)
|
||||
endif
|
||||
return
|
||||
end subroutine add_efield
|
||||
|
||||
31
PW/input.f90
31
PW/input.f90
|
|
@ -26,7 +26,8 @@ subroutine iosys
|
|||
press, cmass, calc, cell_factor, xqq, alat, ntypx, &
|
||||
lmovecell, imix, at, omega, ityp, tau, nks, xk, wk, uakbar, amconv, &
|
||||
force, at_old, omega_old, starting_scf_threshold, title, crystal, &
|
||||
atm, nk1, nk2, nk3, k1, k2, k3
|
||||
atm, nk1, nk2, nk3, k1, k2, k3, &
|
||||
tefield, edir, emaxpos, eopreg, eamp
|
||||
use io, only : tmp_dir, prefix, pseudo_dir, pseudop
|
||||
use constants, only: pi
|
||||
#ifdef __PARA
|
||||
|
|
@ -53,7 +54,7 @@ subroutine iosys
|
|||
NAMELIST / control / title, calculation, verbosity, &
|
||||
restart_mode, nstep, iprint, isave, tstress, tprnfor, &
|
||||
dt, ndr, ndw, outdir, prefix, max_seconds, ekin_conv_thr,&
|
||||
etot_conv_thr, forc_conv_thr, pseudo_dir, disk_io
|
||||
etot_conv_thr, forc_conv_thr, pseudo_dir, disk_io, tefield
|
||||
|
||||
! SYSTEM namelist
|
||||
|
||||
|
|
@ -65,7 +66,8 @@ subroutine iosys
|
|||
nr1b, nr2b, nr3b, nosym, starting_magnetization, &
|
||||
occupations, degauss, smearing, &
|
||||
nelup, neldw, nspin, ecfixed, qcutz, q2sigma, xc_type, &
|
||||
lda_plus_U, Hubbard_U, Hubbard_alpha
|
||||
lda_plus_U, Hubbard_U, Hubbard_alpha, &
|
||||
edir, emaxpos, eopreg, eamp
|
||||
|
||||
! ELECTRONS namelist
|
||||
|
||||
|
|
@ -155,6 +157,7 @@ subroutine iosys
|
|||
etot_conv_thr = 1.d-4
|
||||
forc_conv_thr = 1.d-3
|
||||
disk_io = 'default'
|
||||
tefield=.false.
|
||||
noinv = .false. ! not actually used
|
||||
!
|
||||
#ifdef __T3E
|
||||
|
|
@ -198,6 +201,10 @@ subroutine iosys
|
|||
Hubbard_U(:) = 0.d0
|
||||
Hubbard_alpha(:) = 0.d0
|
||||
starting_magnetization(:) = 0.d0
|
||||
edir=1
|
||||
emaxpos=0.5d0
|
||||
eopreg=0.1d0
|
||||
eamp=1.d-3
|
||||
|
||||
! ... Variables initialization for ELECTRONS
|
||||
!
|
||||
|
|
@ -325,6 +332,7 @@ subroutine iosys
|
|||
startingpot = 'atomic'
|
||||
end if
|
||||
|
||||
|
||||
READ (unit, system, iostat = ios )
|
||||
if (ios /= 0) call errore ('reading','namelist &system',2)
|
||||
READ (unit, electrons, iostat = ios )
|
||||
|
|
@ -344,6 +352,16 @@ subroutine iosys
|
|||
READ (unit, phonon, iostat = ios )
|
||||
if (ios /= 0) call errore ('reading','namelist &phonon',5)
|
||||
end if
|
||||
|
||||
if (tefield.and..not.nosym) then
|
||||
nosym=.true.
|
||||
write(6,'(5x,"Presently no symmetry can be used with electric field",/)')
|
||||
endif
|
||||
if (tefield.and.(tstress.or.tprnfor)) then
|
||||
tstress=.false.
|
||||
tprnfor=.false.
|
||||
write(6,'(5x,"Presently stress and forces not available with electric field",/)')
|
||||
endif
|
||||
#ifdef __PARA
|
||||
end if
|
||||
|
||||
|
|
@ -370,6 +388,7 @@ subroutine iosys
|
|||
CALL mp_bcast( forc_conv_thr, ionode_id )
|
||||
CALL mp_bcast( pseudo_dir, ionode_id )
|
||||
CALL mp_bcast( disk_io, ionode_id )
|
||||
CALL mp_bcast( tefield, ionode_id )
|
||||
!
|
||||
! ... SYSTEM Variables Broadcast
|
||||
!
|
||||
|
|
@ -405,7 +424,11 @@ subroutine iosys
|
|||
CALL mp_bcast( Hubbard_U, ionode_id )
|
||||
CALL mp_bcast( Hubbard_alpha, ionode_id )
|
||||
CALL mp_bcast( starting_magnetization, ionode_id )
|
||||
|
||||
CALL mp_bcast( edir, ionode_id )
|
||||
CALL mp_bcast( emaxpos, ionode_id )
|
||||
CALL mp_bcast( eopreg, ionode_id )
|
||||
CALL mp_bcast( eamp, ionode_id )
|
||||
!
|
||||
! ... ELECTRONS Variables Broadcast
|
||||
!
|
||||
CALL mp_bcast( emass, ionode_id )
|
||||
|
|
|
|||
14
PW/pwcom.f90
14
PW/pwcom.f90
|
|
@ -583,6 +583,19 @@ module ldaU
|
|||
conv_ns ! .true. if ns are converged
|
||||
end module ldaU
|
||||
|
||||
module efield_0
|
||||
use parameters
|
||||
logical :: tefield ! if true a finite electric field is added to the
|
||||
! local potential
|
||||
integer :: edir ! direction of the field
|
||||
real(kind=DP) :: &
|
||||
emaxpos, & ! position of the maximum of the field (0<emaxpos<1)
|
||||
eopreg, & ! amplitude of the inverse region (0<eopreg<1)
|
||||
eamp ! field amplitude (in a.u.) (1 a.u.=51.44 10^11 V/m)
|
||||
|
||||
end module efield_0
|
||||
|
||||
|
||||
module pwcom
|
||||
use constants, only: e2, degspin, rytoev, amconv, uakbar, pi, tpi, fpi
|
||||
use brilz
|
||||
|
|
@ -611,5 +624,6 @@ module pwcom
|
|||
use filnam
|
||||
use us
|
||||
use ldaU
|
||||
use efield_0
|
||||
end module pwcom
|
||||
!
|
||||
|
|
|
|||
|
|
@ -14,6 +14,7 @@ subroutine setlocal
|
|||
!
|
||||
#include "machine.h"
|
||||
use pwcom
|
||||
|
||||
implicit none
|
||||
complex(kind=DP), allocatable :: aux (:)
|
||||
! auxiliary variable
|
||||
|
|
@ -22,12 +23,13 @@ subroutine setlocal
|
|||
! counter on g vectors
|
||||
! counter on r vectors
|
||||
!
|
||||
|
||||
allocate (aux( nrxx))
|
||||
!
|
||||
call setv (2 * nrxx, 0.d0, aux, 1)
|
||||
do nt = 1, ntyp
|
||||
do ng = 1, ngm
|
||||
aux (nl (ng) ) = aux (nl (ng) ) + vloc (igtongl (ng), nt) * strf (ng, nt)
|
||||
aux (nl(ng))=aux(nl(ng)) + vloc (igtongl (ng), nt) * strf (ng, nt)
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
|
|
@ -38,7 +40,12 @@ subroutine setlocal
|
|||
do ir = 1, nrxx
|
||||
vltot (ir) = DREAL (aux (ir) )
|
||||
enddo
|
||||
!
|
||||
!
|
||||
! If required add an electric field to the local potential
|
||||
!
|
||||
|
||||
if (tefield) call add_efield
|
||||
|
||||
deallocate(aux)
|
||||
return
|
||||
end subroutine setlocal
|
||||
|
|
|
|||
Loading…
Reference in New Issue