mirror of https://gitlab.com/QEF/q-e.git
misc. cleanup, postprocessing adapted to last changes
git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@323 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz
parent
1d34b5c0a8
commit
1361c7ba39
11
PP/bands.f90
11
PP/bands.f90
|
|
@ -75,7 +75,16 @@ subroutine punch_band (filband)
|
|||
#ifdef __PARA
|
||||
use para, only: me
|
||||
#endif
|
||||
use pwcom
|
||||
use atom
|
||||
use basis
|
||||
use brilz
|
||||
use constants, only: rytoev
|
||||
use gvect
|
||||
use klist
|
||||
use units
|
||||
use wvfct
|
||||
use wavefunctions
|
||||
use us
|
||||
use becmod
|
||||
|
||||
implicit none
|
||||
|
|
|
|||
|
|
@ -16,6 +16,7 @@ program projwfc
|
|||
! outdir temporary directory where files resides
|
||||
!
|
||||
use pwcom
|
||||
use parameters, only : DP
|
||||
use io_files, only: nd_nmbr, prefix, tmp_dir
|
||||
#ifdef __PARA
|
||||
use para, only: me
|
||||
|
|
@ -84,7 +85,20 @@ subroutine projwave (io_choice,Emin, Emax, DeltaE, smoothing)
|
|||
!-----------------------------------------------------------------------
|
||||
!
|
||||
#include "machine.h"
|
||||
use pwcom
|
||||
use atom
|
||||
use basis
|
||||
use brilz
|
||||
use constants, only: rytoev
|
||||
use gvect
|
||||
use klist
|
||||
use ldaU
|
||||
use lsda_mod
|
||||
use symme, only: nsym, irt
|
||||
use varie, only: newpseudo
|
||||
use wvfct
|
||||
use units
|
||||
use us
|
||||
use wavefunctions
|
||||
use becmod
|
||||
use io_files, only: nd_nmbr, prefix, tmp_dir
|
||||
#ifdef __PARA
|
||||
|
|
|
|||
|
|
@ -6,115 +6,115 @@
|
|||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
subroutine compute_dip(dip, dipion, z0)
|
||||
!
|
||||
! This routine computes the integral 1/Omega \int \rho(r) r d^3r
|
||||
! and gives as output the projection of the dipole in the direction of
|
||||
! the electric field. (This routine is called only if tefield is true)
|
||||
! The direction is the reciprocal lattice vector bg(.,edir)
|
||||
!
|
||||
use pwcom
|
||||
!
|
||||
! This routine computes the integral 1/Omega \int \rho(r) r d^3r
|
||||
! and gives as output the projection of the dipole in the direction of
|
||||
! the electric field. (This routine is called only if tefield is true)
|
||||
! The direction is the reciprocal lattice vector bg(.,edir)
|
||||
!
|
||||
use pwcom
|
||||
#ifdef __PARA
|
||||
use para
|
||||
use para
|
||||
#endif
|
||||
implicit none
|
||||
real(kind=dp) :: dip, dipion
|
||||
real(kind=dp), allocatable :: rrho (:), aux(:), rws(:,:)
|
||||
real(kind=dp) :: dipol_ion(3), dipol(3)
|
||||
implicit none
|
||||
real(kind=dp) :: dip, dipion
|
||||
real(kind=dp), allocatable :: rrho (:), aux(:), rws(:,:)
|
||||
real(kind=dp) :: dipol_ion(3), dipol(3)
|
||||
|
||||
integer:: ipol, i, j, k, i1, j1, k11, na, is, ir
|
||||
integer:: nrws, nrwsx
|
||||
real(kind=dp) :: deltax, deltay, deltaz, rijk(3), bmod, proj, x0(3), z0
|
||||
real(kind=dp) :: weight, wsweight
|
||||
!
|
||||
! calculate ionic dipole
|
||||
!
|
||||
x0=0.d0
|
||||
x0(3)=-z0
|
||||
dipol_ion=0.d0
|
||||
do na=1,nat
|
||||
do ipol=1,3
|
||||
dipol_ion(ipol)=dipol_ion(ipol)+zv(ityp(na))*(tau(ipol,na)+x0(ipol))*alat
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! collect the charge density: sum over spin and collect in parallel case
|
||||
!
|
||||
allocate (rrho(nrx1*nrx2*nrx3))
|
||||
integer:: ipol, i, j, k, i1, j1, k11, na, is, ir
|
||||
integer:: nrws, nrwsx
|
||||
real(kind=dp) :: deltax, deltay, deltaz, rijk(3), bmod, proj, x0(3), z0
|
||||
real(kind=dp) :: weight, wsweight
|
||||
!
|
||||
! calculate ionic dipole
|
||||
!
|
||||
x0=0.d0
|
||||
x0(3)=-z0
|
||||
dipol_ion=0.d0
|
||||
do na=1,nat
|
||||
do ipol=1,3
|
||||
dipol_ion(ipol)=dipol_ion(ipol)+zv(ityp(na))*(tau(ipol,na)+x0(ipol))*alat
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! collect the charge density: sum over spin and collect in parallel case
|
||||
!
|
||||
allocate (rrho(nrx1*nrx2*nrx3))
|
||||
rrho(:) = 0.d0
|
||||
#ifdef __PARA
|
||||
allocate(aux(nrx1*nrx2*nrx3))
|
||||
rrho=0.d0
|
||||
do is=1,nspin
|
||||
call gather (rho(1,is), aux)
|
||||
rrho=rrho+aux
|
||||
enddo
|
||||
deallocate(aux)
|
||||
if (me.eq.1) then
|
||||
allocate(aux(nrxx))
|
||||
aux(:) =0.d0
|
||||
do is=1,nspin
|
||||
aux(:) = aux(:) + rho(:,is)
|
||||
enddo
|
||||
call gather (aux(1), rrho(1))
|
||||
deallocate(aux)
|
||||
if (me.eq.1) then
|
||||
#else
|
||||
rrho=0.d0
|
||||
do is=1,nspin
|
||||
rrho=rrho+rho(:,is)
|
||||
enddo
|
||||
do is=1,nspin
|
||||
rrho=rrho+rho(:,is)
|
||||
enddo
|
||||
#endif
|
||||
|
||||
nrwsx=125
|
||||
allocate(rws(0:3,nrwsx))
|
||||
call wsinit(rws,nrwsx,nrws,at)
|
||||
nrwsx=125
|
||||
allocate(rws(0:3,nrwsx))
|
||||
call wsinit(rws,nrwsx,nrws,at)
|
||||
|
||||
deltax=1.d0/real(nr1,dp)
|
||||
deltay=1.d0/real(nr2,dp)
|
||||
deltaz=1.d0/real(nr3,dp)
|
||||
dipol=0.d0
|
||||
do i1 = -nr1/2, nr1/2
|
||||
i=i1+1
|
||||
if (i.lt.1) i=i+nr1
|
||||
do j1 = -nr2/2,nr2/2
|
||||
j=j1+1
|
||||
if (j.lt.1) j=j+nr2
|
||||
do k11 = -nr3/2, nr3/2
|
||||
k=k11+1
|
||||
if (k.lt.1) k=k+nr3
|
||||
ir=i + (j-1)*nrx1 + (k-1)*nrx1*nrx2
|
||||
do ipol=1,3
|
||||
rijk(ipol) = real(i1,dp)*at(ipol,1)*deltax + &
|
||||
real(j1,dp)*at(ipol,2)*deltay + &
|
||||
real(k11,dp)*at(ipol,3)*deltaz
|
||||
enddo
|
||||
weight=wsweight(rijk,rws,nrws)
|
||||
do ipol=1,3
|
||||
dipol(ipol)=dipol(ipol)+weight*(rijk(ipol)+x0(ipol))*rrho(ir)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
deltax=1.d0/real(nr1,dp)
|
||||
deltay=1.d0/real(nr2,dp)
|
||||
deltaz=1.d0/real(nr3,dp)
|
||||
dipol=0.d0
|
||||
do i1 = -nr1/2, nr1/2
|
||||
i=i1+1
|
||||
if (i.lt.1) i=i+nr1
|
||||
do j1 = -nr2/2,nr2/2
|
||||
j=j1+1
|
||||
if (j.lt.1) j=j+nr2
|
||||
do k11 = -nr3/2, nr3/2
|
||||
k=k11+1
|
||||
if (k.lt.1) k=k+nr3
|
||||
ir=i + (j-1)*nrx1 + (k-1)*nrx1*nrx2
|
||||
do ipol=1,3
|
||||
rijk(ipol) = real(i1,dp)*at(ipol,1)*deltax + &
|
||||
real(j1,dp)*at(ipol,2)*deltay + &
|
||||
real(k11,dp)*at(ipol,3)*deltaz
|
||||
enddo
|
||||
weight=wsweight(rijk,rws,nrws)
|
||||
do ipol=1,3
|
||||
dipol(ipol)=dipol(ipol)+weight*(rijk(ipol)+x0(ipol))*rrho(ir)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
dipol=dipol*alat*omega/nr1/nr2/nr3
|
||||
write(6,'(5x,"electron", 3f15.5)') dipol(1), dipol(2), dipol(3)
|
||||
write(6,'(5x,"ion ", 3f15.5)') dipol_ion(1), dipol_ion(2), dipol_ion(3)
|
||||
write(6,'(5x,"total ", 3f15.5)') dipol_ion(1)-dipol(1), &
|
||||
dipol_ion(2)-dipol(2), &
|
||||
dipol_ion(3)-dipol(3)
|
||||
dipol=dipol*alat*omega/nr1/nr2/nr3
|
||||
write(6,'(5x,"electron", 3f15.5)') dipol(1), dipol(2), dipol(3)
|
||||
write(6,'(5x,"ion ", 3f15.5)') dipol_ion(1), dipol_ion(2), dipol_ion(3)
|
||||
write(6,'(5x,"total ", 3f15.5)') dipol_ion(1)-dipol(1), &
|
||||
dipol_ion(2)-dipol(2), &
|
||||
dipol_ion(3)-dipol(3)
|
||||
|
||||
bmod=sqrt(bg(1,edir)**2+bg(2,edir)**2+bg(3,edir)**2)
|
||||
proj=0.d0
|
||||
do ipol=1,3
|
||||
proj=proj+(dipol_ion(ipol)-dipol(ipol))*bg(ipol,edir)
|
||||
enddo
|
||||
proj=proj/bmod
|
||||
dip= fpi*proj/omega
|
||||
bmod=sqrt(bg(1,edir)**2+bg(2,edir)**2+bg(3,edir)**2)
|
||||
proj=0.d0
|
||||
do ipol=1,3
|
||||
proj=proj+(dipol_ion(ipol)-dipol(ipol))*bg(ipol,edir)
|
||||
enddo
|
||||
proj=proj/bmod
|
||||
dip= fpi*proj/omega
|
||||
|
||||
proj=0.d0
|
||||
do ipol=1,3
|
||||
proj=proj+dipol_ion(ipol)*bg(ipol,edir)
|
||||
enddo
|
||||
proj=proj/bmod
|
||||
dipion= fpi*proj/omega
|
||||
deallocate(rws)
|
||||
proj=0.d0
|
||||
do ipol=1,3
|
||||
proj=proj+dipol_ion(ipol)*bg(ipol,edir)
|
||||
enddo
|
||||
proj=proj/bmod
|
||||
dipion= fpi*proj/omega
|
||||
deallocate(rws)
|
||||
|
||||
#ifdef __PARA
|
||||
endif
|
||||
endif
|
||||
#endif
|
||||
|
||||
deallocate (rrho)
|
||||
deallocate (rrho)
|
||||
|
||||
return
|
||||
end
|
||||
return
|
||||
end subroutine compute_dip
|
||||
|
|
|
|||
|
|
@ -7,54 +7,53 @@
|
|||
!
|
||||
|
||||
subroutine delta_e (nr1, nr2, nr3, nrxx, rho, v_in, v_out, omega, &
|
||||
de, deband, nspin)
|
||||
use parameters
|
||||
implicit none
|
||||
de, deband, nspin)
|
||||
use parameters
|
||||
implicit none
|
||||
|
||||
integer :: nspin, nr1, nr2, nr3, nrxx, i
|
||||
! nspin=1 if LDA, nspin=2 if LSDA
|
||||
! fft grid dimensions
|
||||
! counter
|
||||
integer :: nspin, nr1, nr2, nr3, nrxx, i
|
||||
! nspin=1 if LDA, nspin=2 if LSDA
|
||||
! fft grid dimensions
|
||||
! counter
|
||||
|
||||
real(kind=DP) :: rho (nrxx, nspin), v_in (nrxx, nspin), v_out (nrxx, &
|
||||
nspin), omega, de, deband
|
||||
! charge density
|
||||
! in and ...
|
||||
! ... out potentials from potential mixing
|
||||
! cell volume
|
||||
! total energy and band energy corrections
|
||||
|
||||
real(kind=DP) :: rho (nrxx, nspin), v_in (nrxx, nspin), v_out (nrxx, &
|
||||
nspin), omega, de, deband
|
||||
! charge density
|
||||
! in and ...
|
||||
! ... out potentials from potential mixing
|
||||
! cell volume
|
||||
! total energy and band energy corrections
|
||||
de = 0.d0
|
||||
de = 0.d0
|
||||
deband = 0.d0
|
||||
if (nspin.eq.1) then
|
||||
!
|
||||
! spin-unpolarized case
|
||||
!
|
||||
do i = 1, nrxx
|
||||
de = de + rho(i,nspin) * ( v_out(i,1) - v_in(i,1) )
|
||||
deband = deband - rho(i,nspin) * v_in(i,1)
|
||||
enddo
|
||||
else
|
||||
!
|
||||
! spin-polarized case
|
||||
!
|
||||
do i = 1, nrxx
|
||||
de = de + rho(i,1) * ( v_out(i,1) - v_in(i,1) ) + &
|
||||
rho(i,2) * ( v_out(i,2) - v_in(i,2) )
|
||||
deband = deband - rho(i,1) * v_in(i,1) - rho(i,2) * v_in(i,2)
|
||||
enddo
|
||||
|
||||
deband = 0.d0
|
||||
if (nspin.eq.1) then
|
||||
!
|
||||
! LDA case
|
||||
!
|
||||
do i = 1, nrxx
|
||||
de = de + rho(i,nspin) * ( v_out(i,1) - v_in(i,1) )
|
||||
deband = deband - rho(i,nspin) * v_in(i,1)
|
||||
enddo
|
||||
else
|
||||
!
|
||||
! LSDA case
|
||||
!
|
||||
do i = 1, nrxx
|
||||
de = de + rho(i,1) * ( v_out(i,1) - v_in(i,1) ) + &
|
||||
rho(i,2) * ( v_out(i,2) - v_in(i,2) )
|
||||
deband = deband - rho(i,1) * v_in(i,1) - rho(i,2) * v_in(i,2)
|
||||
enddo
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
de = omega * de / (nr1 * nr2 * nr3)
|
||||
deband = omega * deband / (nr1 * nr2 * nr3)
|
||||
de = omega * de / (nr1 * nr2 * nr3)
|
||||
deband = omega * deband / (nr1 * nr2 * nr3)
|
||||
|
||||
#ifdef __PARA
|
||||
call reduce (1, de)
|
||||
call reduce (1, deband)
|
||||
call reduce (1, de)
|
||||
call reduce (1, deband)
|
||||
#endif
|
||||
|
||||
return
|
||||
return
|
||||
|
||||
end subroutine delta_e
|
||||
|
|
|
|||
|
|
@ -20,12 +20,11 @@ subroutine gather (f_in, f_out)
|
|||
use parameters, only : DP
|
||||
implicit none
|
||||
|
||||
real (8) :: f_in (nxx), f_out ( * )
|
||||
real (kind=DP) :: f_in (nxx), f_out ( * )
|
||||
include 'mpif.h'
|
||||
|
||||
|
||||
integer :: root, proc, info, displs (nprocp), recvcount (nprocp)
|
||||
root = 0
|
||||
integer :: root = 0
|
||||
integer :: proc, info, displs (nprocp), recvcount (nprocp)
|
||||
|
||||
call start_clock ('gather')
|
||||
do proc = 1, nprocp
|
||||
recvcount (proc) = ncplane * npp (proc)
|
||||
|
|
@ -34,10 +33,9 @@ subroutine gather (f_in, f_out)
|
|||
else
|
||||
displs (proc) = displs (proc - 1) + recvcount (proc - 1)
|
||||
endif
|
||||
|
||||
enddo
|
||||
call mpi_barrier (MPI_COMM_POOL, info)
|
||||
|
||||
call mpi_barrier (MPI_COMM_POOL, info)
|
||||
call mpi_gatherv (f_in, recvcount (me), MPI_REAL8, f_out, &
|
||||
recvcount, displs, MPI_REAL8, root, MPI_COMM_POOL, info)
|
||||
call errore ('gather', 'info<>0', info)
|
||||
|
|
|
|||
110
PW/move_ions.f90
110
PW/move_ions.f90
|
|
@ -91,64 +91,64 @@ subroutine move_ions
|
|||
end subroutine move_ions
|
||||
!-------------------------------------------------------------------
|
||||
subroutine new_force (dg, dg2)
|
||||
!-------------------------------------------------------------------
|
||||
!
|
||||
! find the lagrange multiplier lambda for the problem with one const
|
||||
!
|
||||
! force*dg
|
||||
! lambda = - --------,
|
||||
! |dg|^2
|
||||
!
|
||||
! and redefine the forces:
|
||||
!
|
||||
! force = force + lambda*dg
|
||||
!
|
||||
! where dg is the gradient of the constraint function
|
||||
!
|
||||
use pwcom
|
||||
integer :: na, i, ipol
|
||||
!-------------------------------------------------------------------
|
||||
!
|
||||
! find the lagrange multiplier lambda for the problem with one const
|
||||
!
|
||||
! force*dg
|
||||
! lambda = - --------,
|
||||
! |dg|^2
|
||||
!
|
||||
! and redefine the forces:
|
||||
!
|
||||
! force = force + lambda*dg
|
||||
!
|
||||
! where dg is the gradient of the constraint function
|
||||
!
|
||||
use pwcom
|
||||
integer :: na, i, ipol
|
||||
|
||||
real(kind=DP) :: dg (3, nat), lambda, dg2, DDOT, sum
|
||||
real(kind=DP) :: dg (3, nat), lambda, dg2, DDOT, sum
|
||||
|
||||
lambda = 0.d0
|
||||
if (dg2.ne.0.d0) then
|
||||
lambda = - DDOT (3 * nat, force, 1, dg, 1) / dg2
|
||||
call DAXPY (3 * nat, lambda, dg, 1, force, 1)
|
||||
if (DDOT (3 * nat, force, 1, dg, 1) **2.gt.1.d-30) then
|
||||
call errore ('new_force', 'force is not orthogonal to constrain', - 1)
|
||||
print *, DDOT (3 * nat, force, 1, dg, 1) **2
|
||||
endif
|
||||
do ipol = 1, 3
|
||||
sum = 0.d0
|
||||
do na = 1, nat
|
||||
sum = sum + force (ipol, na)
|
||||
enddo
|
||||
!
|
||||
! impose total force = 0
|
||||
!
|
||||
do na = 1, nat
|
||||
force (ipol, na) = force (ipol, na) - sum / nat
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! resymmetrize (should not be needed, but...)
|
||||
!
|
||||
if (nsym.gt.1) then
|
||||
do na = 1, nat
|
||||
call trnvect (force (1, na), at, bg, - 1)
|
||||
enddo
|
||||
call symvect (nat, force, nsym, s, irt)
|
||||
do na = 1, nat
|
||||
call trnvect (force (1, na), at, bg, 1)
|
||||
enddo
|
||||
endif
|
||||
write (6, '(/5x,"Constrained forces")')
|
||||
do na = 1, nat
|
||||
write (6, '(3f14.8)') (force (i, na) , i = 1, 3)
|
||||
enddo
|
||||
lambda = 0.d0
|
||||
if (dg2.ne.0.d0) then
|
||||
lambda = - DDOT (3 * nat, force, 1, dg, 1) / dg2
|
||||
call DAXPY (3 * nat, lambda, dg, 1, force, 1)
|
||||
if (DDOT (3 * nat, force, 1, dg, 1) **2.gt.1.d-30) then
|
||||
call errore ('new_force', 'force is not orthogonal to constrain', - 1)
|
||||
print *, DDOT (3 * nat, force, 1, dg, 1) **2
|
||||
endif
|
||||
do ipol = 1, 3
|
||||
sum = 0.d0
|
||||
do na = 1, nat
|
||||
sum = sum + force (ipol, na)
|
||||
enddo
|
||||
!
|
||||
! impose total force = 0
|
||||
!
|
||||
do na = 1, nat
|
||||
force (ipol, na) = force (ipol, na) - sum / nat
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! resymmetrize (should not be needed, but...)
|
||||
!
|
||||
if (nsym.gt.1) then
|
||||
do na = 1, nat
|
||||
call trnvect (force (1, na), at, bg, - 1)
|
||||
enddo
|
||||
call symvect (nat, force, nsym, s, irt)
|
||||
do na = 1, nat
|
||||
call trnvect (force (1, na), at, bg, 1)
|
||||
enddo
|
||||
endif
|
||||
write (6, '(/5x,"Constrained forces")')
|
||||
do na = 1, nat
|
||||
write (6, '(3f14.8)') (force (i, na) , i = 1, 3)
|
||||
enddo
|
||||
|
||||
endif
|
||||
return
|
||||
endif
|
||||
return
|
||||
|
||||
end subroutine new_force
|
||||
!---------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -10,7 +10,7 @@
|
|||
subroutine print_clock_pw
|
||||
!
|
||||
! this routine prints out the clocks at the end of the run
|
||||
! it tries to construct the calling tree of the program..
|
||||
! it tries to construct the calling tree of the program.
|
||||
|
||||
use pwcom
|
||||
implicit none
|
||||
|
|
|
|||
|
|
@ -35,7 +35,6 @@ program pwscf
|
|||
call ions
|
||||
if (conv_ions) goto 10
|
||||
call hinit1
|
||||
|
||||
enddo
|
||||
|
||||
10 call punch
|
||||
|
|
|
|||
|
|
@ -47,8 +47,8 @@ subroutine read_file
|
|||
kunittmp = 0
|
||||
|
||||
! here we read the variables that dimension the system
|
||||
! in parallel execution, only root proc read the file
|
||||
! and then broadcast the values to all ather procs
|
||||
! in parallel execution, only root proc reads the file
|
||||
! and then broadcasts the values to all other procs
|
||||
!
|
||||
call readfile_new( 'dim', iunpun, rdum, rdum, kunittmp, 0, 0, ierr )
|
||||
IF( ierr /= 0 ) THEN
|
||||
|
|
|
|||
|
|
@ -17,9 +17,10 @@ subroutine scatter (f_in, f_out)
|
|||
#ifdef __PARA
|
||||
#include "machine.h"
|
||||
use para
|
||||
use parameters, only : DP
|
||||
implicit none
|
||||
|
||||
real (8) :: f_in ( * ), f_out (nxx)
|
||||
real (kind=DP) :: f_in ( * ), f_out (nxx)
|
||||
include 'mpif.h'
|
||||
|
||||
|
||||
|
|
@ -38,7 +39,7 @@ subroutine scatter (f_in, f_out)
|
|||
call mpi_barrier (MPI_COMM_POOL, info)
|
||||
call mpi_scatterv (f_in, sendcount, displs, MPI_REAL8, f_out, &
|
||||
sendcount (me), MPI_REAL8, root, MPI_COMM_POOL, info)
|
||||
call errore ('gather_pot', 'info<>0', info)
|
||||
call errore ('scatter', 'info<>0', info)
|
||||
if (sendcount(me).ne.nxx) f_out(sendcount(me)+1:nxx) = 0.d0
|
||||
|
||||
call stop_clock ('scatter')
|
||||
|
|
|
|||
|
|
@ -177,12 +177,12 @@ subroutine setup
|
|||
|
||||
call recips ( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) )
|
||||
!
|
||||
! If lxkcry = .true. , the input k-point components in crystallogra
|
||||
! units are transformed in cartesian coordinates
|
||||
! If lxkcry = .true. , the input k-point components in crystal
|
||||
! axis are transformed in cartesian coordinates
|
||||
!
|
||||
if (lxkcry) call cryst_to_cart (nks, xk, bg, 1)
|
||||
!
|
||||
! Test that atomic coordinates are different
|
||||
! Test that atoms do not overlap
|
||||
!
|
||||
if (.not. (lchk_tauxk (nat, tau, bg) ) ) call errore ('setup', &
|
||||
'Wrong atomic coordinates ', 1)
|
||||
|
|
@ -290,7 +290,7 @@ subroutine setup
|
|||
.and. .not. ( calc.eq.'mm' .or. calc.eq.'nm' ) ) &
|
||||
call errore ('setup', 'Dynamics, you should have no symmetries', -1)
|
||||
!
|
||||
! Automatic generation of k-points
|
||||
! Calculate quantities used in tetrahedra method
|
||||
!
|
||||
|
||||
if (ltetra) then
|
||||
|
|
|
|||
|
|
@ -152,7 +152,7 @@ subroutine sgam_at (nrot, s, nat, tau, ityp, at, bg, nr1, nr2, &
|
|||
ft1 = ft (1) * nr1
|
||||
ft2 = ft (2) * nr2
|
||||
ft3 = ft (3) * nr3
|
||||
! check if the fractional transaltions are commensurate
|
||||
! check if the fractional translations are commensurate
|
||||
! with the FFT grid, discard sym.op. if not
|
||||
if (abs (ft1 - nint (ft1) ) / nr1.gt.1.0d-5 .or. &
|
||||
abs (ft2 - nint (ft2) ) / nr2.gt.1.0d-5 .or. &
|
||||
|
|
|
|||
128
PW/symrho.f90
128
PW/symrho.f90
|
|
@ -8,81 +8,67 @@
|
|||
!
|
||||
!-----------------------------------------------------------------------
|
||||
subroutine symrho (rho, nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, s, &
|
||||
ftau)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! symmetrize the charge density.
|
||||
!
|
||||
ftau)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! symmetrize the charge density.
|
||||
!
|
||||
#include "machine.h"
|
||||
use parameters
|
||||
implicit none
|
||||
!
|
||||
! first the dummy variables
|
||||
!
|
||||
integer :: nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, s (3, 3, 48), &
|
||||
ftau (3, 48)
|
||||
!
|
||||
!
|
||||
! input: dimensions of the FFT mesh
|
||||
!
|
||||
!
|
||||
!
|
||||
! input: the number of symmetry
|
||||
! input: the symmetry matrices
|
||||
! input: the fractionary translation
|
||||
real(kind=DP) :: rho (nrx1, nrx2, nrx3)
|
||||
! inp/out: the charge density
|
||||
integer , allocatable :: symflag (:,:,:)
|
||||
integer :: ri (48), rj (48), rk (48), &
|
||||
i, j, k, isym
|
||||
! tells if the point has been se
|
||||
!
|
||||
! the rotated of each mesh point
|
||||
!
|
||||
!
|
||||
! counters on mesh points
|
||||
!
|
||||
! counter on symmetries
|
||||
use parameters
|
||||
implicit none
|
||||
!
|
||||
! first the dummy variables
|
||||
!
|
||||
integer :: nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, s (3, 3, 48), &
|
||||
ftau (3, 48)
|
||||
!
|
||||
! input: dimensions of the FFT mesh
|
||||
! input: the number of symmetries
|
||||
! input: the symmetry matrices
|
||||
! input: the fractionary translations
|
||||
!
|
||||
real(kind=DP) :: rho (nrx1, nrx2, nrx3)
|
||||
! inp/out: the charge density
|
||||
integer , allocatable :: symflag (:,:,:)
|
||||
integer :: ri (48), rj (48), rk (48), i, j, k, isym
|
||||
real(kind=DP) :: sum
|
||||
|
||||
real(kind=DP) :: sum
|
||||
! auxiliary variable
|
||||
if (nsym.eq.1) return
|
||||
|
||||
if (nsym.eq.1) return
|
||||
allocate (symflag(nrx1, nrx2, nrx3))
|
||||
do k = 1, nr3
|
||||
do j = 1, nr2
|
||||
do i = 1, nr1
|
||||
symflag (i, j, k) = 0
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do k = 1, nr3
|
||||
do j = 1, nr2
|
||||
do i = 1, nr1
|
||||
if (symflag (i, j, k) .eq.0) then
|
||||
sum = 0.d0
|
||||
do isym = 1, nsym
|
||||
call ruotaijk (s (1, 1, isym), ftau (1, isym), i, j, k, nr1, &
|
||||
nr2, nr3, ri (isym), rj (isym), rk (isym) )
|
||||
sum = sum + rho (ri (isym), rj (isym), rk (isym) )
|
||||
enddo
|
||||
sum = sum / nsym
|
||||
!
|
||||
! sum contains the symmetrized charge density at point r.
|
||||
! now fill the star of r with this sum.
|
||||
!
|
||||
do isym = 1, nsym
|
||||
rho (ri (isym), rj (isym), rk (isym) ) = sum
|
||||
symflag (ri (isym), rj (isym), rk (isym) ) = 1
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
allocate (symflag(nrx1, nrx2, nrx3))
|
||||
do k = 1, nr3
|
||||
do j = 1, nr2
|
||||
do i = 1, nr1
|
||||
symflag (i, j, k) = 0
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do k = 1, nr3
|
||||
do j = 1, nr2
|
||||
do i = 1, nr1
|
||||
if (symflag (i, j, k) .eq.0) then
|
||||
sum = 0.d0
|
||||
do isym = 1, nsym
|
||||
call ruotaijk (s (1, 1, isym), ftau (1, isym), i, j, k, nr1, &
|
||||
nr2, nr3, ri (isym), rj (isym), rk (isym) )
|
||||
sum = sum + rho (ri (isym), rj (isym), rk (isym) )
|
||||
enddo
|
||||
sum = sum / nsym
|
||||
!
|
||||
! sum contains the symmetrized charge density at point r.
|
||||
! now fill the star of r with this sum.
|
||||
!
|
||||
do isym = 1, nsym
|
||||
rho (ri (isym), rj (isym), rk (isym) ) = sum
|
||||
symflag (ri (isym), rj (isym), rk (isym) ) = 1
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
deallocate(symflag)
|
||||
return
|
||||
deallocate(symflag)
|
||||
return
|
||||
end subroutine symrho
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue