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

- common xc and gga functionals subroutines for all code PW/CP/FPMD 

- Few more functionals has been added to file more_functionals.f90
  since in PW BLYP/LSDA were not implemented .
  In the same file, temporary, are stored old CP90 subroutines
  for testing purpose ( in CP functionals were vectorized!! ).
  The small program PP/xctest.f90 can be used to make a comparative
  test between new and old routines, in case someone suspect a problem
  there.
- buon anno!


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1533 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
cavazzon 2004-12-31 11:14:32 +00:00
parent a925df0b01
commit e32ce9b074
29 changed files with 2521 additions and 4203 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CPV/Makefile
View File

@ -49,7 +49,6 @@ fields_type.o \
fnl.o \
forces.o \
fromscra.o \
gga.o \
greenf.o \
grid.o \
gsmesh.o \
@ -61,10 +60,9 @@ interfaces_main.o \
interfaces.o \
ions.o \
ksstates.o \
lda.o \
lsd.o \
macdep.o \
main.o \
mainvar.o \
main_loops.o \
modules.o \
nl_base.o \

1557
CPV/cplib.f90
View File

File diff suppressed because it is too large Load Diff

129
CPV/cpr.f90
View File

@ -172,6 +172,8 @@
USE electrons_module, ONLY: cp_eigs
USE print_out_module, ONLY: cp_print_rho
USE cp_main_variables
!
!
implicit none
@ -181,64 +183,16 @@
real(kind=8) :: tau(3,*)
real(kind=8) :: fion_out(3,*)
real(kind=8) :: etot_out
!
!
! control variables
!
logical tbump
logical tfirst, tlast
logical tstop
logical tconv
real(kind=8) :: delta_etot
!
! structure factors e^{-ig*R}
!
complex(kind=8), allocatable:: ei1(:,:,:), ei2(:,:,:), ei3(:,:,:)
complex(kind=8), allocatable:: eigr(:,:,:)
!
! structure factors (summed over atoms of the same kind)
!
complex(kind=8), allocatable:: sfac(:,:)
!
! indexes, positions, and structure factors for the box grid
!
integer irb(3,natx,nsx)
real(kind=8) taub(3,natx)
complex(kind=8), allocatable:: eigrb(:,:,:)
!
! charge densities and potentials
! rhog = charge density in g space
! rhor = charge density in r space (dense grid)
! rhos = charge density in r space (smooth grid)
! rhoc = core charge density in real space (dense grid)
!
complex(kind=8), allocatable:: rhog(:,:)
real(kind=8), allocatable:: rhor(:,:), rhos(:,:), rhoc(:)
!
! nonlocal projectors:
! bec = scalar product of projectors and wave functions
! betae = nonlocal projectors in g space = beta x e^(-ig.R)
! becdr = <betae|g|psi> used in force calculation
! rhovan= \sum_i f(i) <psi(i)|beta_l><beta_m|psi(i)>
! deeq = \int V_eff(r) q_lm(r) dr
!
real(kind=8), allocatable:: bec(:,:), becdr(:,:,:)
real(kind=8), allocatable:: bephi(:,:), becp(:,:)
!
! mass preconditioning
!
real(kind=8), allocatable:: ema0bg(:)
real(kind=8), allocatable:: emadt2(:)
real(kind=8), allocatable:: emaver(:)
!
! constraints (lambda at t, lambdam at t-dt, lambdap at t+dt)
!
real(kind=8), allocatable:: lambda(:,:), lambdam(:,:), lambdap(:,:)
logical tbump, tfirst, tlast, tstop, tconv
logical :: ttprint ! logical variable used to control printout
!
! ionic positions, center of mass position
!
real(kind=8) cdm0(3)
real(kind=8) cdm(3)
!
! forces on ions
!
@ -246,8 +200,6 @@
!
! work variables
!
real(kind=8) acc(nacx)
complex(kind=8), allocatable:: c2(:), c3(:)
complex(kind=8) speed
real(kind=8) &
& tempp, verl1, verl2, verl3, &
@ -256,27 +208,23 @@
& ettt, ccc, bigr, dt2, dt2by2, twodel, dt2bye, dt2hbe
real(kind=8) ekinc0, ekinp, ekinpr, ekincm, ekinc
real(kind=8) temps(nsx)
integer is, nacc, ia, j, iter, nfi, i, isa, ipos
real(kind=8) ekinh, temphc, temp1, temp2, randy
real(kind=8) :: delta_etot
real(kind=8) ftmp, enb, enbi
integer :: is, nacc, ia, j, iter, nfi, i, isa, ipos
integer :: k, ii, l, m, ibeg
!
! work variables, 2
!
real(kind=8) hnew(3,3),velh(3,3),hgamma(3,3), temphh(3,3)
real(kind=8) fcell(3,3)
real(kind=8) cdm(3)
!
integer :: k, ii, l, m, ibeg
real(kind=8) ekinh, temphc, temp1, temp2, randy
real(kind=8) ftmp, enb, enbi
character(len=256) :: filename
character(len=256) :: dirname
integer :: strlen, dirlen
real(kind=8) :: b1(3), b2(3), b3(3)
real(kind=8) :: stress_gpa(3,3), thstress(3,3)
logical :: ttprint ! logical variable used to control printout
real(kind=8), allocatable :: tauw( :, : ) ! temporary array used
! to printout positions
! temporary array used to printout positions
real(kind=8), allocatable :: tauw( :, : )
!
! ==================================================================
@ -364,35 +312,18 @@
! allocation of all arrays not already allocated in init and nlinit
! ==================================================================
!
CALL allocate_mainvar &
( ngw, ngb, ngs, ng, nr1, nr2, nr3, nnr, nnrsx, nax, nsp, nspin, n, nx, nhsa, nlcc_any )
allocate(c0(ngw,nx,1,1))
allocate(cm(ngw,nx,1,1))
allocate(phi(ngw,nx,1,1))
allocate(wrk2(ngw,max(nax,n)))
allocate(eigr(ngw,nax,nsp))
allocate(eigrb(ngb,nax,nsp))
allocate(sfac(ngs,nsp))
allocate(rhops(ngs,nsp))
allocate(vps(ngs,nsp))
allocate(rhor(nnr,nspin))
allocate(rhos(nnrsx,nspin))
allocate(rhog(ng,nspin))
if ( nlcc_any ) allocate(rhoc(nnr))
allocate(wrk1(nnr))
allocate(qv(nnrb))
allocate(c2(ngw))
allocate(c3(ngw))
allocate(ema0bg(ngw))
allocate(lambda(nx,nx))
allocate(lambdam(nx,nx))
allocate(lambdap(nx,nx))
allocate(ei1(-nr1:nr1,nax,nsp))
allocate(ei2(-nr2:nr2,nax,nsp))
allocate(ei3(-nr3:nr3,nax,nsp))
allocate(betae(ngw,nhsa))
allocate(becdr(nhsa,n,3))
allocate(bec (nhsa,n))
allocate(bephi(nhsa,n))
allocate(becp (nhsa,n))
allocate(deeq(nhm,nhm,nat,nspin))
allocate(rhovan(nhm*(nhm+1)/2,nat,nspin))
allocate(dbec (nhsa,n,3,3))
@ -492,10 +423,10 @@
!
call from_restart &
( sfac, eigr, ei1, ei2, ei3, bec, becdr, tfirst, eself, fion, &
taub, irb, eigrb, b1, b2, b3, nfi, rhog, rhor, rhos, rhoc, enl, ekin, stress, &
detot, enthal, etot, lambda, lambdam, lambdap, ema0bg, dbec, delt, &
bephi, becp, velh, dt2bye, iforce, fionm, nbeg, xnhe0, xnhem, vnhe, ekincm )
( sfac, eigr, ei1, ei2, ei3, bec, becdr, tfirst, fion, &
taub, irb, eigrb, b1, b2, b3, nfi, rhog, rhor, rhos, rhoc, stress, &
detot, enthal, lambda, lambdam, lambdap, ema0bg, dbec, &
bephi, becp, velh, dt2bye, fionm, ekincm )
!
end if
!==============================================end of if(nbeg.lt.0)====
@ -1226,27 +1157,7 @@
!
1977 format(5x,//'====================== end cprvan ======================',//)
IF( ALLOCATED( ei1 ) ) DEALLOCATE( ei1 )
IF( ALLOCATED( ei2 ) ) DEALLOCATE( ei2 )
IF( ALLOCATED( ei3 ) ) DEALLOCATE( ei3 )
IF( ALLOCATED( eigr ) ) DEALLOCATE( eigr )
IF( ALLOCATED( sfac ) ) DEALLOCATE( sfac )
IF( ALLOCATED( eigrb ) ) DEALLOCATE( eigrb )
IF( ALLOCATED( rhor ) ) DEALLOCATE( rhor )
IF( ALLOCATED( rhos ) ) DEALLOCATE( rhos )
IF( ALLOCATED( rhog ) ) DEALLOCATE( rhog )
IF( ALLOCATED( rhoc ) ) DEALLOCATE( rhoc )
IF( ALLOCATED( bec ) ) DEALLOCATE( bec )
IF( ALLOCATED( becdr ) ) DEALLOCATE( becdr )
IF( ALLOCATED( bephi ) ) DEALLOCATE( bephi )
IF( ALLOCATED( becp ) ) DEALLOCATE( becp )
IF( ALLOCATED( ema0bg ) ) DEALLOCATE( ema0bg )
IF( ALLOCATED( lambda ) ) DEALLOCATE( lambda )
IF( ALLOCATED( lambdam ) ) DEALLOCATE( lambdam )
IF( ALLOCATED( lambdap ) ) DEALLOCATE( lambdap )
IF( ALLOCATED( c2 ) ) DEALLOCATE( c2 )
IF( ALLOCATED( c3 ) ) DEALLOCATE( c3 )
CALL deallocate_mainvar()
CALL deallocate_cvan()
CALL deallocate_efield( )
CALL deallocate_ensemble_dft()

40
CPV/cprsub.f90
View File

@ -108,10 +108,10 @@
use gvecp, only: ng => ngm
use grid_dimensions, only: nr1, nr2, nr3, nnr => nnrx
use cell_base, only: ainv
!use cell_module
use cell_base, only: omega
use control_flags, only: tpre
use derho
use mp, only: mp_sum
!
implicit none
! input
@ -134,27 +134,15 @@
call fillgrad(nspin,rhog,gradr)
end if
!
exc=0.0
if (iexch==1.and.icorr==1.and.igcx==0.and.igcc==0) then
! LDA (Perdew-Zunger)
call expxc(nnr,nspin,rhor,exc)
else if (iexch==1.and.icorr==4.and.igcx==2.and.igcc==2) then
! PW91
call ggapwold(nspin,rhog,gradr,rhor,exc)
else if (iexch==1.and.icorr==3.and.igcx==1.and.igcc==3) then
! BLYP
call ggablyp4(nspin,rhog,gradr,rhor,exc)
else if (iexch==1.and.icorr==4.and.igcx==3.and.igcc==4) then
! PBE
call ggapbe(nspin,rhog,gradr,rhor,exc)
else
call errore('exc-cor','no such exch-corr',1)
end if
CALL exch_corr_cp(nnr,nspin,gradr,rhor,exc)
call mp_sum( exc )
exc=exc*omega/dble(nr1*nr2*nr3)
!
! exchange-correlation contribution to pressure
!
dxc(:,:) = 0.0
dxc = 0.0
if (tpre) then
if (nspin.ne.1) call errore('exc-cor','spin not implemented',1)
!
@ -166,9 +154,7 @@
dxc(i,j)=omega/(nr1*nr2*nr3)*dxc(i,j)
end do
end do
#ifdef __PARA
call reduce (9,dxc)
#endif
call mp_sum ( dxc )
do j=1,3
do i=1,3
dxc(i,j) = dxc(i,j) + exc*ainv(j,i)
@ -179,16 +165,10 @@
! second part of the xc-potential
!
if (igcx /= 0 .or. igcc /= 0) then
call gradh(nspin,gradr,rhog,rhor,dexc)
call gradh( nspin, gradr, rhog, rhor, dexc)
if (tpre) then
#ifdef __PARA
call reduce (9,dexc)
#endif
do j=1,3
do i=1,3
dxc(i,j) = dxc(i,j) + dexc(i,j)
end do
end do
call mp_sum ( dexc )
dxc = dxc + dexc
end if
deallocate(gradr)
end if

286
CPV/exch_corr.f90
View File

@ -6,83 +6,267 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
MODULE exchange_correlation
!=----------------------------------------------------------------------------=!
MODULE exchange_correlation
!=----------------------------------------------------------------------------=!
USE kinds
USE local_density_approximation
USE local_spin_density
USE kinds, ONLY: dbl
IMPLICIT NONE
SAVE
PRIVATE
LOGICAL :: use_table
! ... Gradient Correction & exchange and correlation
PUBLIC :: v2gc, vofxc_lsd, exch_corr_energy, vofxc_lda
LOGICAL :: tgc
REAL(dbl), PARAMETER :: small_rho = 1.0d-10
PUBLIC :: v2gc, exch_corr_energy
PUBLIC :: exch_corr_setup, exch_corr_print_info
PUBLIC :: deallocate_lda
PUBLIC :: tgc, narray
PUBLIC :: tgc
CONTAINS
!=----------------------------------------------------------------------------=!
CONTAINS
!=----------------------------------------------------------------------------=!
SUBROUTINE exch_corr_setup(xc_type, narray_i, rmxxc_i)
CHARACTER(LEN = 20) xc_type
INTEGER, INTENT(IN) :: narray_i
REAL(dbl), INTENT(IN) :: rmxxc_i
CALL lsd_setup(xc_type)
CALL lda_setup(narray_i, rmxxc_i)
IF( xc_type == 'LDA' ) THEN
use_table = .TRUE.
ELSE
use_table = .FALSE.
END IF
RETURN
END SUBROUTINE exch_corr_setup
SUBROUTINE exch_corr_setup( xc_type )
USE funct, ONLY: dft, which_dft, igcx, igcc
CHARACTER(LEN = 20) :: xc_type
dft = TRIM( xc_type )
CALL which_dft( dft )
tgc = .FALSE.
if( ( igcx > 0 ) .OR. ( igcc > 0 ) ) tgc = .TRUE.
RETURN
END SUBROUTINE exch_corr_setup
SUBROUTINE exch_corr_print_info(iunit)
INTEGER, INTENT(IN) :: iunit
WRITE(iunit,800)
IF( use_table ) THEN
CALL lda_print_info(iunit)
ELSE
CALL lsd_print_info(iunit)
!=----------------------------------------------------------------------------=!
SUBROUTINE exch_corr_print_info(iunit)
USE funct, ONLY: dft, iexch, icorr, igcx, igcc
INTEGER, INTENT(IN) :: iunit
CHARACTER(LEN = 20) :: xcgc_type
CHARACTER(LEN = 60) :: exch_info
CHARACTER(LEN = 60) :: corr_info
CHARACTER(LEN = 60) :: exgc_info
CHARACTER(LEN = 60) :: cogc_info
WRITE(iunit,800)
! ... iexch => Exchange functional form
! ... icorr => Correlation functional form
! ... igcx => Gradient Correction to the Exchange potential
! ... igcc => Gradient Correction to the Correlation potential
SELECT CASE ( iexch )
CASE (0)
exch_info = 'NONE'
CASE (1)
exch_info = 'SLATER'
CASE (2)
exch_info = 'SLATER (alpha=1)'
CASE DEFAULT
exch_info = 'UNKNOWN'
END SELECT
SELECT CASE ( icorr )
CASE (0)
corr_info = 'NONE'
CASE (1)
corr_info = 'PERDEW AND ZUNGER'
CASE (2)
corr_info = 'VOSKO, WILK AND NUSAIR'
CASE (3)
corr_info = 'LEE, YANG, AND PARR'
CASE (4)
corr_info = 'PERDEW AND WANG'
CASE (9)
corr_info = 'PADE APPROXIMATION'
CASE DEFAULT
corr_info = 'UNKNOWN'
END SELECT
SELECT CASE ( igcx )
CASE (0)
exgc_info = 'NONE'
CASE (1)
exgc_info = 'BECKE'
CASE (2)
exgc_info = 'PERDEW'
CASE (3)
exgc_info = 'PERDEW BURKE ERNZERHOF'
CASE DEFAULT
exgc_info = 'UNKNOWN'
END SELECT
SELECT CASE ( igcc )
CASE (0)
cogc_info = 'NONE'
CASE (1)
cogc_info = 'PERDEW'
CASE (2)
cogc_info = 'LEE, YANG AND PARR'
CASE (3)
cogc_info = 'PERDEW AND WANG'
CASE (4)
cogc_info = 'PERDEW BURKE ERNZERHOF'
CASE DEFAULT
cogc_info = 'UNKNOWN'
END SELECT
WRITE(iunit,910)
WRITE(iunit,fmt='(5X,"Exchange functional: ",A)') exch_info
WRITE(iunit,fmt='(5X,"Correlation functional: ",A)') corr_info
IF(tgc) THEN
WRITE(iunit,810)
WRITE(iunit,fmt='(5X,"Exchange functional: ",A)') exgc_info
WRITE(iunit,fmt='(5X,"Correlation functional: ",A)') cogc_info
END IF
800 FORMAT(//,3X,'Exchange and correlations functionals',/ &
,3X,'-------------------------------------')
WRITE( iunit, '(5x,"Exchange-correlation = ",a, &
& " (",4i1,")")') trim(dft) , iexch, icorr, igcx, igcc
800 FORMAT(//,3X,'Exchange and correlations functionals',/ &
,3X,'-------------------------------------')
810 FORMAT( 3X,'Using Generalized Gradient Corrections with')
910 FORMAT( 3X,'Using Local Density Approximation with')
RETURN
END SUBROUTINE exch_corr_print_info
!=----------------------------------------------------------------------------=!
!=----------------------------------------------------------------------------=!
SUBROUTINE v2gc(v2xc, grho, rhoer, vpot, gv)
USE fft
USE stick, ONLY: dfftp
USE cell_base, ONLY: tpiba
USE cp_types
USE mp_global
!
implicit none
!
REAL(dbl) :: vpot(:,:,:,:)
REAL(dbl), intent(in) :: v2xc(:,:,:,:,:)
REAL(dbl), intent(in) :: grho(:,:,:,:,:)
REAL(dbl), intent(in) :: rhoer(:,:,:,:)
type (recvecs), intent(in) :: GV
!
integer ig, ipol, nxl, nyl, nzl, i, j, k, is, js, nspin
COMPLEX(dbl), allocatable :: psi(:,:,:)
COMPLEX(dbl), allocatable :: vtemp(:)
COMPLEX(dbl), allocatable :: vtemp_pol(:)
REAL(dbl), ALLOCATABLE :: v(:,:,:)
REAL(dbl) :: fac
INTEGER :: kk(2,2), nn(2,2)
! ...
nxl = dfftp%nr1
nyl = dfftp%nr2
nzl = dfftp%npl
nspin = SIZE(rhoer,4)
kk(1,1)=1
kk(1,2)=3
kk(2,1)=3
kk(2,2)=2
nn(1,1)=1
nn(1,2)=2
nn(2,1)=1
nn(2,2)=2
fac = REAL(nspin) / 2.0d0
DO js = 1, nspin
allocate(vtemp(gv%ng_l))
allocate(vtemp_pol(gv%ng_l))
allocate(psi(nxl,nyl,nzl))
vtemp = CMPLX(0.0d0,0.0d0)
!
! DO is = 1, nspin
! DO ipol = 1, 3
! DO k = 1, nzl
! DO j = 1, nyl
! DO i = 1, nxl
! psi(i,j,k) = fac * v2xc(i,j,k,kk(js,is)) * grho(i,j,k,ipol,nn(js,is))
! END DO
! END DO
! END DO
! CALL pfwfft(vtemp_pol, psi)
! DO ig = gv%gstart, gv%ng_l
! vtemp(ig) = vtemp(ig) + vtemp_pol(ig) * CMPLX( 0.d0, tpiba * gv%gx_l( ipol, ig ) )
! END DO
! END DO
! END DO
DO ipol = 1, 3
DO is = 1, nspin
DO k = 1, nzl
DO j = 1, nyl
DO i = 1, nxl
psi(i,j,k) = fac * v2xc(i,j,k,js,is) * grho(i,j,k,ipol,is)
END DO
END DO
END DO
CALL pfwfft(vtemp_pol, psi)
DO ig = gv%gstart, gv%ng_l
vtemp(ig) = vtemp(ig) + vtemp_pol(ig) * CMPLX( 0.d0, tpiba * gv%gx_l( ipol, ig ) )
END DO
END DO
END DO
DEALLOCATE(psi)
DEALLOCATE(vtemp_pol)
ALLOCATE(v(nxl,nyl,nzl))
CALL pinvfft(v, vtemp)
DO k = 1, nzl
DO j = 1, nyl
DO i = 1, nxl
vpot(i,j,k,js) = vpot(i,j,k,js) - v(i,j,k)
END DO
END DO
END DO
DEALLOCATE(vtemp)
DEALLOCATE(v)
END DO
RETURN
END SUBROUTINE exch_corr_print_info
END SUBROUTINE
!=----------------------------------------------------------------------------=!
SUBROUTINE exch_corr_energy(rhoetr, rhoetg, grho, vpot, &
sxc, vxc, v2xc, v2xc2, gv)
SUBROUTINE exch_corr_energy(rhoetr, rhoetg, grho, vpot, sxc, vxc, v2xc, gv)
USE kinds, ONLY: dbl
USE cp_types, ONLY: recvecs
REAL (dbl) :: rhoetr(:,:,:,:)
COMPLEX(dbl) :: rhoetg(:,:)
REAL (dbl) :: grho(:,:,:,:,:)
REAL (dbl) :: vpot(:,:,:,:)
REAL (dbl) :: sxc, vxc
REAL (dbl) :: v2xc(:,:,:,:)
REAL (dbl) :: v2xc2(:,:,:)
REAL(dbl) :: sxc ! E_xc energy
REAL(dbl) :: vxc ! SUM ( v(r) * rho(r) )
REAL (dbl) :: v2xc(:,:,:,:,:)
TYPE (recvecs), INTENT(IN) :: gv
INTEGER :: nspin, ispin
nspin = SIZE(rhoetr, 4)
IF(.NOT.tgc) THEN
! ... no gradient correction
! ... vpot = vxc(rhoetr); vpot(r) <-- u(r)
IF( (nspin .EQ. 1) .AND. use_table ) THEN
CALL vofxc_lda(rhoetr(:,:,:,1), vpot(:,:,:,1), sxc, vxc)
ELSE
CALL vofxc_lsd(rhoetr, vpot, sxc, vxc)
END IF
ELSE
! ... gradient correction
CALL vofxc_lsd(rhoetr, vpot, grho, v2xc, v2xc2, sxc, vxc)
CALL v2gc(v2xc, v2xc2, grho, rhoetr, vpot, gv)
INTEGER :: nxl, nyl, nzl, nspin, nnr
! ... vpot = vxc(rhoetr); vpot(r) <-- u(r)
nxl = SIZE(rhoetr,1)
nyl = SIZE(rhoetr,2)
nzl = SIZE(rhoetr,3)
nspin = SIZE(rhoetr,4)
nnr = nxl * nyl * nzl
!
CALL exch_corr_wrapper( nnr, nspin, grho, rhoetr, sxc, vpot, v2xc )
!
vxc = SUM( vpot * rhoetr )
!
IF(tgc) THEN
! ... vpot additional term for gradient correction
CALL v2gc( v2xc, grho, rhoetr, vpot, gv )
END If
RETURN
END SUBROUTINE
!=----------------------------------------------------------------------------=!
END MODULE exchange_correlation
!=----------------------------------------------------------------------------=!

706
CPV/gga.f90
View File

@ -1,706 +0,0 @@
!
! Copyright (C) 2002 FPMD group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
module gradient_corrections
USE kinds
IMPLICIT NONE
SAVE
PRIVATE
! ... Gradient Correction & exchange and correlation
logical :: tgc
integer :: mfxcx,mfxcc,mgcx,mgcc
REAL(dbl) :: salpha,bbeta,betapp
! PUBLIC :: gc_setup
contains
subroutine gc_setup(tgc_inp,type_inp)
logical, intent(in) :: tgc_inp
character*20, intent(in) :: type_inp
tgc = tgc_inp
IF(tgc) THEN
IF(type_inp.EQ.'BLYP') THEN
mfxcx = 1
mgcx = 1
mfxcc = 3
mgcc = 2
bbeta = 0.0042d0
betapp= 0.0d0
salpha= 0.6666666666d0
ELSE IF(type_inp.EQ.'BP') THEN
mfxcx = 1
mgcx = 1
mfxcc = 1
mgcc = 1
bbeta = 0.0042d0
betapp= 0.0d0
salpha= 0.6666666666d0
ELSE
call errore(' module gradient corrections ', &
' tgc is true but type is wrong ',0)
ENDIF
ENDIF
return
end subroutine gc_setup
! ==================================================================
SUBROUTINE XC(RHO,EX,EC,VX,VC)
! ==--------------------------------------------------------------==
! == LDA EXCHANGE AND CORRELATION FUNCTIONALS ==
! == ==
! == EXCHANGE : SLATER alpha ==
! == CORRELATION : CEPERLEY & ALDER (PERDEW-ZUNGER PARAMETERS) ==
! == VOSKO, WILK & NUSSAIR ==
! == LEE, YANG & PARR ==
! == PERDEW & WANG ==
! == WIGNER ==
! == HEDIN & LUNDQVIST ==
! == ORTIZ & BALLONE (PERDEW-ZUNGER FORMULA) ==
! == ORTIZ & BALLONE (PERDEW-WANG FORMULA) ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
! ... Gradient Correction & exchange and correlation
REAL(dbl) :: SMALL, PI34, THIRD
PARAMETER (SMALL=1.D-10)
PARAMETER (PI34= 0.75D0 / 3.141592653589793D+00 ,THIRD=1.D0/3.D0)
! ==--------------------------------------------------------------==
!..Exchange
IF(MFXCX.EQ.1) THEN
CALL SLATERX(RHO,EX,VX,SALPHA)
ELSE
EX=0.0D0
VX=0.0D0
ENDIF
IF(RHO.LE.SMALL) THEN
EC = 0.0D0
VC = 0.0D0
EX = 0.0D0
VX = 0.0D0
ELSE IF(MFXCC.EQ.1) THEN
RS=(PI34/RHO)**THIRD
IFLG=2
IF(RS.LT.1.0D0) IFLG=1
CALL PZ(RS,EC,VC,IFLG)
ELSEIF(MFXCC.EQ.2) THEN
RS = (PI34/RHO)**THIRD
CALL VWN(RS,EC,VC)
ELSEIF(MFXCC.EQ.3) THEN
CALL LYP(RHO,EC,VC)
ELSEIF(MFXCC.EQ.4) THEN
RS=(PI34/RHO)**THIRD
IFLG=2
IF(RS.LT.0.5D0) IFLG=1
IF(RS.GT.100.D0) IFLG=3
CALL PW(RS,EC,VC,IFLG)
ELSEIF(MFXCC.EQ.5) THEN
CALL WIGNER(RHO,EC,VC)
ELSEIF(MFXCC.EQ.6) THEN
CALL HEDIN(RHO,EC,VC)
ELSEIF(MFXCC.EQ.7) THEN
RS=(PI34/RHO)**THIRD
IFLG=2
IF(RS.LT.1.0D0) IFLG=1
CALL OBPZ(RS,EC,VC,IFLG)
ELSEIF(MFXCC.EQ.8) THEN
RS=(PI34/RHO)**THIRD
IFLG=2
IF(RS.LT.0.5D0) IFLG=1
IF(RS.GT.100.D0) IFLG=3
CALL OBPW(RS,EC,VC,IFLG)
ELSEIF(MFXCC.EQ.9) THEN
RS=(PI34/RHO)**THIRD
CALL PADE(RS,EC,VC)
ELSE
EC=0.0D0
VC=0.0D0
ENDIF
RETURN
END SUBROUTINE XC
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE GCXC(RHO,GRHO,SX,SC,V1X,V2X,V1C,V2C)
! ==--------------------------------------------------------------==
! == GRADIENT CORRECTIONS FOR EXCHANGE AND CORRELATION ==
! == ==
! == EXCHANGE : BECKE88 ==
! == GGAX ==
! == PBEX ==
! == CORRELATION : PERDEW86 ==
! == LEE, YANG & PARR ==
! == GGAC ==
! == PBEC ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: SMALL
PARAMETER (SMALL=1.D-10)
! ... Gradient Correction & exchange and correlation
! ==--------------------------------------------------------------==
!..Exchange
IF(RHO.LE.SMALL) THEN
SX = 0.0D0
V1X = 0.0D0
V2X = 0.0D0
ELSEIF(MGCX.EQ.1) THEN
CALL BECKE88(bbeta,RHO,GRHO,SX,V1X,V2X)
ELSEIF(MGCX.EQ.2) THEN
CALL GGAX(RHO,GRHO,SX,V1X,V2X)
ELSEIF(MGCX.EQ.3) THEN
CALL PBEX(RHO,GRHO,SX,V1X,V2X)
ELSE
SX=0.0D0
V1X=0.0D0
V2X=0.0D0
ENDIF
!..Correlation
IF(RHO.LE.SMALL) THEN
SC = 0.0D0
V1C = 0.0D0
V2C = 0.0D0
ELSEIF(MGCC.EQ.1) THEN
CALL PERDEW86(RHO,GRHO,SC,V1C,V2C)
ELSEIF(MGCC.EQ.2) THEN
CALL GLYP(RHO,GRHO,SC,V1C,V2C)
ELSEIF(MGCC.EQ.3) THEN
CALL GGAC(RHO,GRHO,SC,V1C,V2C)
ELSEIF(MGCC.EQ.4) THEN
CALL PBEC(RHO,GRHO,SC,V1C,V2C)
ELSE
SC=0.0D0
V1C=0.0D0
V2C=0.0D0
ENDIF
RETURN
END SUBROUTINE GCXC
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE SLATERX(RHO,EX,VX,ALPHA)
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: SMALL, F1, THIRD, F43
PARAMETER (SMALL=1.D-10)
PARAMETER (F1 = -1.10783814957303361D0)
PARAMETER (THIRD=1.D0/3.D0,F43=4.D0/3.D0)
! ==--------------------------------------------------------------==
IF(RHO.LE.SMALL) THEN
EX = 0.0D0
VX = 0.0D0
ELSE
RS = RHO**THIRD
EX = F1*ALPHA*RS
VX = F43*F1*ALPHA*RS
ENDIF
RETURN
END SUBROUTINE SLATERX
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE PZ(RS,EPZ,VPZ,IFLG)
! ==--------------------------------------------------------------==
! == J.P. PERDEW AND ALEX ZUNGER PRB 23, 5048 (1981) ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: A, B, C, D, GC, B1, B2
PARAMETER (A=0.0311D0,B=-0.048D0,C=0.0020D0,D=-0.0116D0, &
GC=-0.1423D0,B1=1.0529D0,B2=0.3334D0)
! ==--------------------------------------------------------------==
EPWC=0.0D0
VPWC=0.0D0
IF(IFLG.EQ.1) THEN
!..High density formula
XLN=LOG(RS)
EPZ=A*XLN+B+C*RS*XLN+D*RS
VPZ=A*XLN+(B-A/3.D0)+2.D0/3.D0*C*RS*XLN+ (2.D0*D-C)/3.D0*RS
ELSEIF(IFLG.EQ.2) THEN
!..Interpolation formula
RS1=SQRT(RS)
RS2=RS
OX=1.D0+B1*RS1+B2*RS2
DOX=1.D0+7./6.*B1*RS1+4./3.*B2*RS2
EPZ=GC/OX
VPZ=EPZ*DOX/OX
ENDIF
RETURN
END SUBROUTINE PZ
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE VWN(RS,EVWN,VVWN)
! ==--------------------------------------------------------------==
! == S.H VOSKO, L.WILK, AND M. NUSAIR, ==
! == CAN. J. PHYS. 58 1200 (1980) ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: a, b, c, x0, two
PARAMETER (A=0.0310907,B=3.72744,C=12.9352,X0=-0.10498)
PARAMETER (TWO=2.0D0)
! ==--------------------------------------------------------------==
Q = SQRT(4.D0*C-B*B)
F1 = TWO*B/Q
F2 = B*X0/(X0*X0+B*X0+C)
F3 = TWO*(TWO*X0+B)/Q
X = SQRT(RS)
FX = X*X+B*X+C
QX = ATAN(Q/(TWO*X+B))
EVWN=A*(LOG(RS/FX)+F1*QX-F2*(LOG((X-X0)**2/FX)+F3*QX))
TXPB=TWO*X+B
TTQQ=TXPB*TXPB+Q*Q
VVWN=EVWN - X*A/6.*(TWO/X-TXPB/FX-4.*B/TTQQ-F2*(TWO/(X-X0) &
-TXPB/FX-4.*(TWO*X0+B)/TTQQ))
RETURN
END SUBROUTINE VWN
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE LYP(RHO,ELYP,VLYP)
! ==--------------------------------------------------------------==
! == C. LEE, W. YANG, AND R.G. PARR, PRB 37, 785 (1988) ==
! == THIS IS ONLY THE LDA PART ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: a, b, c, d, cf
PARAMETER (A=0.04918,B=0.132,C=0.2533,D=0.349)
PARAMETER (CF=2.87123400018819108)
! ==--------------------------------------------------------------==
RS=RHO**(-1./3.)
ECRS=B*CF*EXP(-C*RS)
OX=1.D0/(1.D0+D*RS)
ELYP=-A*OX*(1.D0+ECRS)
VLYP=ELYP-RS/3.*A*OX*(D*OX+ECRS*(D*OX+C))
RETURN
END SUBROUTINE LYP
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE PW(RS,EPWC,VPWC,IFLG)
! ==--------------------------------------------------------------==
! == J.P. PERDEW AND YUE WANG PRB 45, 13244 (1992) ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: a, a1, b1, b2, b3, b4, c0, c1, c2, c3, d0, d1
PARAMETER (A=0.031091,A1=0.21370,B1=7.5957,B2=3.5876,B3=1.6382, &
B4=0.49294,C0=A,C1=0.046644,C2=0.00664,C3=0.01043, &
D0=0.4335,D1=1.4408)
! ==--------------------------------------------------------------==
EPWC=0.0D0
VPWC=0.0D0
IF(IFLG.EQ.1) THEN
!..High density formula
XLN=LOG(RS)
EPWC=C0*XLN-C1+C2*RS*XLN-C3*RS
VPWC=C0*XLN-(C1+C0/3.D0)+2.D0/3.D0*C2*RS*XLN- (2.D0*C3+C2)/3.D0*RS
ELSEIF(IFLG.EQ.2) THEN
!..Interpolation formula
RS1=SQRT(RS)
RS2=RS
RS3=RS2*RS1
RS4=RS2*RS2
OM=2.D0*A*(B1*RS1+B2*RS2+B3*RS3+B4*RS4)
DOM=2.D0*A*(0.5*B1*RS1+B2*RS2+1.5D0*B3*RS3+2.D0*B4*RS4)
OLOG=LOG(1.D0+1.0D0/OM)
EPWC=-2.D0*A*(1+A1*RS)*OLOG
VPWC=-2.*A*(1.D0+2./3.*A1*RS)*OLOG-2./3.*A*(1.+A1*RS)*DOM/(OM*(OM+1.))
ELSEIF(IFLG.EQ.3) THEN
!..Low density formula
EPWC=-D0/RS+D1/RS**1.5D0
VPWC=-4.D0/3.D0*D0/RS+1.5D0*D1/RS**1.5
ENDIF
RETURN
END SUBROUTINE PW
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE WIGNER(RHO,EXC,FXC)
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
RH=RHO
X=RH**0.33333333333333333D0
FXC=-X*((0.943656+8.8963*X)/(1.0+12.57*X)**2)
EXC=-0.738*X*(0.959/(1.0+12.57*X))
RETURN
END SUBROUTINE WIGNER
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE HEDIN(RHO,ECP,FCP)
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
SAVE RH
IF(RH .EQ. 0.0D0) RETURN
RH=RHO
RSM1=0.62035049D0*RH**(0.3333333333333333D0)
ALN=DLOG(1.0D0 + 21.0D0*RSM1)
X=21.0D0/RSM1
ECP = ALN+(X**3*ALN-X*X)+X/2-1.0D0/3.0D0
ECP = -0.0225*ECP
FCP = -0.0225*ALN
RETURN
END SUBROUTINE HEDIN
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE OBPZ(RS,EPZ,VPZ,IFLG)
! ==--------------------------------------------------------------==
! == G.ORTIZ AND P. BALLONE PRB 50, 1391 (1994) ==
! == PERDEW-ZUNGER FORMULA ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER(I-N)
REAL(dbl) :: a, b, c, d, gc, b1, b2
PARAMETER (A=0.031091D0,B=-0.046644D0,C=0.00419D0,D=-0.00983D0, &
GC=-0.103756D0,B1=0.56371D0,B2=0.27358D0)
! ==--------------------------------------------------------------==
EPWC=0.0D0
VPWC=0.0D0
IF(IFLG.EQ.1) THEN
!..High density formula
XLN=LOG(RS)
EPZ=A*XLN+B+C*RS*XLN+D*RS
VPZ=A*XLN+(B-A/3.D0)+2.D0/3.D0*C*RS*XLN+ (2.D0*D-C)/3.D0*RS
ELSEIF(IFLG.EQ.2) THEN
!..Interpolation formula
RS1=SQRT(RS)
RS2=RS
OX=1.D0+B1*RS1+B2*RS2
DOX=1.D0+7./6.*B1*RS1+4./3.*B2*RS2
EPZ=GC/OX
VPZ=EPZ*DOX/OX
ENDIF
RETURN
END SUBROUTINE OBPZ
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE OBPW(RS,EPWC,VPWC,IFLG)
! ==--------------------------------------------------------------==
! == G.ORTIZ AND P. BALLONE PRB 50, 1391 (1994) ==
! == PERDEW-WANG FORMULA ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER(I-N)
REAL(dbl) :: a, a1, b1, b2, b3, b4, c0, c1, c2, c3, d0, d1
PARAMETER (A=0.031091,A1=0.026481,B1=7.5957,B2=3.5876,B3=-0.46647, &
B4=0.13354,C0=A,C1=0.046644,C2=0.00664,C3=0.01043, &
D0=0.4335,D1=1.4408)
! ==--------------------------------------------------------------==
EPWC=0.0D0
VPWC=0.0D0
IF(IFLG.EQ.1) THEN
!..High density formula
XLN=LOG(RS)
EPWC=C0*XLN-C1+C2*RS*XLN-C3*RS
VPWC=C0*XLN-(C1+C0/3.D0)+2.D0/3.D0*C2*RS*XLN- (2.D0*C3+C2)/3.D0*RS
ELSEIF(IFLG.EQ.2) THEN
!..Interpolation formula
RS1=SQRT(RS)
RS2=RS
RS3=RS2*RS1
RS4=RS2*RS2
OM=2.D0*A*(B1*RS1+B2*RS2+B3*RS3+B4*RS4)
DOM=2.D0*A*(0.5*B1*RS1+B2*RS2+1.5D0*B3*RS3+2.D0*B4*RS4)
OLOG=LOG(1.D0+1.0D0/OM)
EPWC=-2.D0*A*(1+A1*RS)*OLOG
VPWC=-2.*A*(1.D0+2./3.*A1*RS)*OLOG &
-2./3.*A*(1.+A1*RS)*DOM/(OM*(OM+1.))
ELSEIF(IFLG.EQ.3) THEN
!..Low density formula
EPWC=-D0/RS+D1/RS**1.5D0
VPWC=-4.D0/3.D0*D0/RS+1.5D0*D1/RS**1.5
ENDIF
RETURN
END SUBROUTINE OBPW
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE PADE(RS,EC,VC)
! ==--------------------------------------------------------------==
! == PADE APPROXIMATION ==
! == S. GOEDECKER, M. TETER, J. HUTTER, PRB in press ==
! ==--------------------------------------------------------------==
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: a0, a1, a2, a3, b1, b2, b3, b4, O3
PARAMETER (A0=0.4581652932831429D0,A1=2.217058676663745D0, &
A2=0.7405551735357053D0,A3=0.01968227878617998D0)
PARAMETER (B1=1.0000000000000000D0,B2=4.504130959426697D0, &
B3=1.110667363742916D0,B4=0.02359291751427506D0)
PARAMETER (O3=1.D0/3.D0)
! ==--------------------------------------------------------------==
TOP=A0+RS*(A1+RS*(A2+RS*A3))
DTOP=A1+RS*(2.D0*A2+3.D0*A3*RS)
BOT=RS*(B1+RS*(B2+RS*(B3+RS*B4)))
DBOT=B1+RS*(2.D0*B2+RS*(3.D0*B3+RS*4.D0*B4))
EC=-TOP/BOT
VC=EC+RS*O3*(DTOP/BOT-TOP*DBOT/(BOT*BOT))
RETURN
END SUBROUTINE PADE
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE BECKE88(B1,RHO,GRHO,SX,V1X,V2X)
! ==--------------------------------------------------------------==
! BECKE EXCHANGE: PRA 38, 3098 (1988)
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: OB3
PARAMETER(OB3=1.D0/3.D0)
! ==--------------------------------------------------------------==
TWO13 = 2.0D0**(1.D0/3.D0)
AA = GRHO
A = SQRT(AA)
BR1 = RHO**OB3
BR2 = BR1*BR1
BR4 = BR2*BR2
XS = TWO13*A/BR4
XS2 = XS*XS
SA2B8 = SQRT(1.0D0+XS2)
SHM1 = LOG(XS+SA2B8)
DD = 1.0D0 + 6.0D0*B1*XS*SHM1
DD2 = DD*DD
EE = 6.0D0*B1*XS2/SA2B8 - 1.D0
SX = TWO13*AA/BR4*(-B1/DD)
V1X = -(4./3.)/TWO13*XS2*B1*BR1*EE/DD2
V2X = TWO13*B1*(EE-DD)/(BR4*DD2)
RETURN
END SUBROUTINE BECKE88
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE GGAX(RHO,GRHO,SX,V1X,V2X)
! ==--------------------------------------------------------------==
! J.P.PERDEW ET AL. PRB 46 6671 (1992)
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: f1, f2, f3, f4, f5, pi
PARAMETER(F1=0.19645,F2=7.7956,F3=0.2743,F4=0.1508,F5=0.004)
PARAMETER(PI=3.141592653589793D+00)
! ==--------------------------------------------------------------==
FP1 = -3./(16.*PI)*(3.*PI*PI)**(-1./3.)
FP2 = 0.5*(3.*PI*PI)**(-1./3.)
AA = GRHO
A = SQRT(AA)
RR = RHO**(-4./3.)
S = FP2*A*RR
S2 = S*S
S3 = S2*S
S4 = S2*S2
EXPS = F4*EXP(-100.*S2)
AS = F3-EXPS-F5*S2
SA2B8 = SQRT(1.0D0+F2*F2*S2)
SHM1 = LOG(F2*S+SA2B8)
BS = 1.D0+F1*S*SHM1+F5*S4
DAS = 200.*S*EXPS-2.*S*F5
DBS = F1*(SHM1+F2*S/SA2B8)+4.*F5*S3
DLS = (DAS/AS-DBS/BS)
SX = FP1*AA*RR*AS/BS
V1X = -4./3.*SX/RHO*(1.+S*DLS)
V2X = FP1*RR*AS/BS*(2.+S*DLS)
RETURN
END SUBROUTINE GGAX
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE PBEX(RHO,GRHO,SX,V1X,V2X)
! ==--------------------------------------------------------------==
! J.P.PERDEW ET AL. PRL XX XXXX (1996)
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: us, ax, um, uk, ul
PARAMETER(US=0.161620459673995492D0,AX=-0.738558766382022406D0, &
UM=0.2195149727645171D0,UK=0.8040D0,UL=UM/UK)
! ==--------------------------------------------------------------==
AA = GRHO
RR = RHO**(-4./3.)
EX = AX/RR
S2 = AA*RR*RR*US*US
PO = 1.D0/(1.D0 + UL*S2)
FX = UK-UK*PO
SX = EX*FX
DFX = 2.D0*UK*UL*PO*PO
V1X = 1.33333333333333D0*AX*RHO**0.333333333333D0*(FX-S2*DFX)
V2X = EX*DFX*(US*RR)**2
RETURN
END SUBROUTINE PBEX
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE PERDEW86(RHO,GRHO,SC,V1C,V2C)
! ==--------------------------------------------------------------==
! PERDEW CORRELATION: PRB 33, 8822 (1986)
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: p1, p2, p3, p4, pc1, pc2, pci, ob3, fpi
PARAMETER(P1=0.023266D0,P2=7.389D-6,P3=8.723D0,P4=0.472D0)
PARAMETER(PC1=0.001667D0,PC2=0.002568,PCI=PC1+PC2)
PARAMETER(OB3=1.D0/3.D0, FPI=4.0D0*3.141592653589793D0)
! ==--------------------------------------------------------------==
AA = GRHO
A = SQRT(AA)
BR1 = RHO**OB3
BR2 = BR1*BR1
BR4 = BR2*BR2
RS = (3./(FPI*RHO))**OB3
RS2 = RS*RS
RS3 = RS*RS2
CNA = PC2+P1*RS+P2*RS2
CNB = 1.+P3*RS+P4*RS2+1.D4*P2*RS3
CN = PC1 + CNA/CNB
DRS = -OB3*(3./FPI)**OB3 / BR4
DCNA = (P1+2.*P2*RS)*DRS
DCNB = (P3+2.*P4*RS+3.D4*P2*RS2)*DRS
DCN = DCNA/CNB - CNA/(CNB*CNB)*DCNB
PHI = 0.192*PCI/CN*A*RHO**(-7./6.)
EPHI = EXP(-PHI)
SC = AA/BR4*CN*EPHI
V1C = SC*((1.+PHI)*DCN/CN -((4./3.)-(7./6.)*PHI)/RHO)
V2C = CN*EPHI/BR4*(2.-PHI)
RETURN
END SUBROUTINE PERDEW86
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE GLYP(RHO,GRHO,SC,V1C,V2C)
! ==--------------------------------------------------------------==
! LEE, YANG PARR: GRADIENT CORRECTION PART
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: a, b, c, d
PARAMETER(A=0.04918,B=0.132,C=0.2533,D=0.349)
! ==--------------------------------------------------------------==
AA = GRHO
R = RHO**(-1.d0/3.d0)
OM = EXP(-C*R)/(1.d0+D*R)
R5 = R**5
XL = 1.d0+(7.d0/3.d0)*(C*R + D*R/(1.d0+D*R))
FF = A*B*AA/24.d0
SC = FF*R5*OM*XL
DR5 = 5.d0*R*R*R*R
DOM = -OM*(C+D+C*D*R)/(1.d0+D*R)
DXL = (7.d0/3.d0)*(C+D+2.d0*C*D*R+C*D*D*R*R)/(1.d0+D*R)**2
V1C = -FF*(R*R*R*R)/3.d0*( DR5*OM*XL + R5*DOM*XL + R5*OM*DXL)
V2C = A*B*R5*OM*XL/12.d0
RETURN
END SUBROUTINE GLYP
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE GGAC(RHO,GRHO,SC,V1C,V2C)
! ==--------------------------------------------------------------==
! PERDEW & WANG GGA CORRELATION PART
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: al, pa, pb, pc, pd, cx, cxc0, cc0, ob3, pi
PARAMETER(AL=0.09,PA=0.023266D0,PB=7.389D-6,pc=8.723D0,PD=0.472D0)
PARAMETER(CX=-0.001667D0,CXC0=0.002568,CC0=-CX+CXC0)
PARAMETER(OB3=1.D0/3.D0, PI=3.141592653589793D0)
! ==--------------------------------------------------------------==
XNU = 16./PI*(3.*PI*PI)**OB3
BE = XNU*CC0
CALL XC(RHO,EX,EC,VX,VC)
AA = GRHO
A = SQRT(AA)
RS = (3./(4.*PI*RHO))**OB3
RS2 = RS*RS
RS3 = RS*RS2
XKF = (9.*PI/4.)**OB3/RS
XKS = SQRT(4.*XKF/PI)
T = A/(2.*XKS*RHO)
EXPE = EXP(-2.*AL*EC/(BE*BE))
AF = 2.*AL/BE * (1./(EXPE-1.))
BF = EXPE*(VC-EC)
Y = AF*T*T
XY = (1.+Y)/(1.+Y+Y*Y)
QY = Y*Y*(2.+Y)/(1.+Y+Y*Y)**2
S1 = 1.+2.*AL/BE*T*T*XY
H0 = BE*BE/(2.*AL) * LOG(S1)
DH0 = BE*T*T/S1*(-7./3.*XY-QY*(AF*BF/BE-7./3.))
DDH0 = BE/(2.*XKS*XKS*RHO)*(XY-QY)/S1
EE = -100.*(XKS/XKF*T)**2
CNA = CXC0+PA*RS+PB*RS2
DCNA = -(PA*RS+2.*PB*RS2)/3.
CNB = 1.+pc*RS+PD*RS2+1.D4*PB*RS3
DCNB = -(pc*RS+2.*PD*RS2+3.D4*PB*RS3)/3.
CN = CNA/CNB - CX
DCN = DCNA/CNB - CNA*DCNB/(CNB*CNB)
H1 = XNU*(CN-CC0-3./7.*CX)*T*T*EXP(EE)
DH1 = -OB3*(H1*(7.+8.*EE)+XNU*T*T*EXP(EE)*DCN)
DDH1 = 2.*H1*(1.+EE)*RHO/AA
SC = RHO*(H0+H1)
V1C = H0+H1+DH0+DH1
V2C = DDH0+DDH1
RETURN
END SUBROUTINE GGAC
! ==--------------------------------------------------------------==
! ==================================================================
SUBROUTINE PBEC(RHO,GRHO,SC,V1C,V2C)
! ==--------------------------------------------------------------==
! PBE Correlation functional
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER(I-N)
REAL(dbl) :: be, ga, cx, cxc0, cc0, ob3, pi
PARAMETER(BE=0.06672455060314922D0,GA=0.031090690869654895D0)
PARAMETER(CX=-0.001667D0,CXC0=0.002568,CC0=-CX+CXC0)
PARAMETER(OB3=1.D0/3.D0, PI=3.141592653589793D0)
! ==--------------------------------------------------------------==
CALL XC(RHO,EX,EC,VX,VC)
AA = GRHO
A = SQRT(AA)
RS = (3./(4.*PI*RHO))**OB3
XKF = (9.*PI/4.)**OB3/RS
XKS = SQRT(4.*XKF/PI)
T = A/(2.*XKS*RHO)
EXPE = EXP(-EC/GA)
AF = BE/GA * (1./(EXPE-1.))
Y = AF*T*T
XY = (1.+Y)/(1.+Y+Y*Y)
S1 = 1.+BE/GA*T*T*XY
H0 = GA * LOG(S1)
DTDR = -T*7./(6.*RHO)
DADR = AF*AF*EXPE/BE*(VC-EC)/RHO
DSDA = -BE/GA * AF * T**6 * (2.+Y) / (1.+Y+Y*Y)**2
DSDT = 2.*BE/GA * T * (1.+2.*Y) / (1.+Y+Y*Y)**2
DSDR = DSDA*DADR + DSDT*DTDR
DHDT = GA/S1*DSDT
DHDR = GA/S1*DSDR
SC = RHO*H0
V1C = H0+DHDR*RHO
V2C = RHO*DHDT*T/AA
RETURN
END SUBROUTINE PBEC
! ==--------------------------------------------------------------==
! ==================================================================
end module gradient_corrections

4
CPV/input.f90
View File

@ -1318,7 +1318,7 @@ MODULE input_fpmd
tturbo_inp, nturbo_inp, diis_achmix, electron_damping, vhrmax_inp, &
vhasse_inp, iesr_inp, diis_ethr, diis_wthr, diis_delt, &
etot_conv_thr, diis_g0chmix, diis_nchmix, diis_g1chmix, xc_type, &
narray_inp, rmxxc_inp, empty_states_maxstep, empty_states_delt, &
empty_states_maxstep, empty_states_delt, &
empty_states_emass, empty_states_ethr, vhnr_inp, vhiunit_inp, &
vhrmin_inp, tvhmean_inp, tscra_inp, scradir, max_seconds, wk, &
outdir, prefix, k3, nk3, k1, k2, woptical, noptical, boptical, &
@ -1485,7 +1485,7 @@ MODULE input_fpmd
IF( empty_states_ethr > 0.d0 ) ethr_emp_inp = empty_states_ethr
CALL empty_setup( empty_states_maxstep, delt_emp_inp, ethr_emp_inp )
CALL exch_corr_setup(xc_type, narray_inp, rmxxc_inp)
CALL exch_corr_setup( xc_type )
CALL check_stop_init( max_seconds )
!
scradir_ = TRIM( scradir )

475
CPV/lda.f90
View File

@ -1,475 +0,0 @@
!
! Copyright (C) 2002 FPMD group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
MODULE local_density_approximation
USE kinds
IMPLICIT NONE
SAVE
INTEGER :: narray
REAL(dbl) :: ddro
REAL(dbl) :: rmxxc
REAL(dbl), ALLOCATABLE :: eexc(:)
REAL(dbl), ALLOCATABLE :: vvxc(:)
PUBLIC :: lda_setup
CONTAINS
SUBROUTINE lda_setup(narray_inp, rmxxc_inp)
INTEGER, INTENT(IN) :: narray_inp
REAL(dbl), INTENT(IN) :: rmxxc_inp
narray = narray_inp
rmxxc = rmxxc_inp
CALL genxc()
RETURN
END SUBROUTINE lda_setup
SUBROUTINE lda_print_info(iunit)
INTEGER, INTENT(IN) :: iunit
WRITE(iunit,10) RMXXC,NARRAY
10 FORMAT( 3X,'Using Local Density Approximation',/ &
,3X,'functional type : SLATER, PERDEW & ZUNGER',/ &
,3X,'Using an interpolating table with : ',/ &
,3X,' Density maximum value = ',E10.4,/ &
,3X,' No. of table entries = ',I6)
RETURN
END SUBROUTINE lda_print_info
SUBROUTINE allocate_lda
ALLOCATE(EEXC(0:NARRAY))
ALLOCATE(VVXC(0:NARRAY))
RETURN
END SUBROUTINE allocate_lda
SUBROUTINE deallocate_lda
IF( ALLOCATED( EEXC ) ) DEALLOCATE(EEXC)
IF( ALLOCATED( VVXC ) ) DEALLOCATE(VVXC)
RETURN
END SUBROUTINE deallocate_lda
SUBROUTINE vofxc_lda(rhoe, vpot, sxc, vxc)
USE io_global, ONLY: stdout
REAL(dbl) :: RHOE(:,:,:)
REAL(dbl) :: vpot(:,:,:)
REAL(dbl) :: sxc, vxc
REAL(dbl) :: ratio, dd, ddi, ee, vxc1, exc1, OBYDRO, roe
INTEGER :: i, j, k, i1
LOGICAL :: tmax, tmin
sxc = 0.0d0
vxc = 0.0d0
tmax = .FALSE.
tmin = .FALSE.
OBYDRO = 1.D0/DDRO
DO k = 1, SIZE(rhoe,3)
DO j = 1, SIZE(rhoe,2)
DO i = 1, SIZE(rhoe,1)
ROE = RHOE(i,j,k)
RATIO = OBYDRO * ROE
I1 = INT(RATIO)
DDI = REAL(I1)
DD = RATIO - DDI
EE = 1.D0 - DD
IF(I1 .GE. NARRAY) THEN
I1 = NARRAY - 1
tmax = .TRUE.
END IF
IF(I1 .LT. 0) then
I1 = 0
tmin = .TRUE.
END IF
VXC1 = VVXC(I1)*EE + VVXC(I1+1)*DD
EXC1 = EEXC(I1)*EE + EEXC(I1+1)*DD
SXC = SXC + EXC1*ROE
VXC = VXC + VXC1*ROE
VPOT(i,j,k) = VXC1
END DO
END DO
END DO
IF( tmax ) THEN
WRITE( stdout, fmt= &
'(" ** vofxc_lda: WARNING charge values excede the table max value")')
END IF
IF( tmin ) THEN
WRITE( stdout, fmt= &
'(" ** vofxc_lda: WARNING negative charge values")')
END IF
RETURN
END SUBROUTINE vofxc_lda
!
!---------------------------------------------------------------
REAL(dbl) FUNCTION ECCA(RS,IFLG)
!---------------------------------------------------------------
!
REAL(dbl), INTENT(IN) :: RS
INTEGER :: IFLG
REAL(dbl) :: A(2),B(2),C(2),D(2),G(2),B1(2),B2(2), RSL, RSQ
DATA A/0.0622D0,0.0311D0/, B/-0.096D0,-0.0538D0/, C/0.0040D0,0.0014D0/,&
D/-0.0232D0,-0.0096D0/, B1/1.0529D0,1.3981D0/, B2/0.3334D0,0.2611D0/,&
G/-0.2846D0,-0.1686D0/
IF(RS.LE.1.D0) THEN
RSL=LOG(RS)
ECCA=A(IFLG)*RSL+B(IFLG)+C(IFLG)*RS*RSL+D(IFLG)*RS
ELSE
RSQ=SQRT(RS)
ECCA=G(IFLG)/(1.D0+B1(IFLG)*RSQ+B2(IFLG)*RS)
END IF
RETURN
END FUNCTION ECCA
!
!---------------------------------------------------------------
REAL(dbl) FUNCTION ECGL(RS,IFLG)
!---------------------------------------------------------------
!
REAL(dbl), INTENT(IN) :: RS
INTEGER :: IFLG
REAL(dbl) :: C(2), R(2), X
REAL(dbl), PARAMETER :: THIRD=1.D0/3.D0
DATA C/0.0666D0,0.0406D0/, R/11.4D0,15.9D0/
X=RS/R(IFLG)
ECGL=-C(IFLG) * ( (1.D0+X**3)*LOG(1.D0+1.D0/X) - THIRD +X*(0.5D0-X) )
RETURN
END FUNCTION ECGL
!
!---------------------------------------------------------------
REAL(dbl) FUNCTION ECHL(RS,IFLG)
!---------------------------------------------------------------
!
REAL(dbl), INTENT(IN) :: RS
INTEGER :: IFLG
REAL(dbl) :: C(2), R(2), X
REAL(dbl), PARAMETER :: THIRD=1.D0/3.D0
DATA C/0.045D0,0.0225D0/, R/21.D0,52.917D0/
X=RS/R(IFLG)
ECHL=-C(IFLG) * ( (1.D0+X**3)*LOG(1.D0+1.D0/X) - THIRD +X*(0.5D0-X) )
RETURN
END FUNCTION ECHL
!
!--------------------------------------------------------------
REAL(dbl) FUNCTION EXC(RHO,ZETA,IXC)
!---------------------------------------------------------------
!
! INTERPOLATION FORMULA FOR THE EXCHANGE-CORRELATION ENERGY
! DENSITY:
! IXC = 0 EXCHANGE ONLY (KOHN & SHAM)
! +/- 1 GUNNARSON & LUNDQVIST
! +/- 2 HEDIN & LUNDQVIST
! +/- 3 CEPERLEY & ALDER
! > 0 EXCHANGE + CORRELATION
! < 0 CORRELATION ONLY
!
IMPLICIT REAL(dbl) (A-H,O-Z)
REAL(dbl) :: rho, zeta
INTEGER :: ixc
INTEGER :: ic, ifl
REAL(dbl), PARAMETER :: PI34 = 0.75D0 / 3.141592653589793D+00
REAL(dbl), PARAMETER :: ALPHA=0.5210617612D+00
REAL(dbl), PARAMETER :: THIRD=1.D0/3.D0
REAL(dbl), PARAMETER :: TW43=2.519842099D0
REAL(dbl), PARAMETER :: SMALL=1.0D-10
IC=IABS(IXC)
IF(IC.GT.3) STOP ' IXC NOT ALLOWED IN EXC'
EXC=0.D0
IF(RHO.LE.SMALL) RETURN
RS=(PI34/RHO)**THIRD
IFL=IC+1
IF(IXC.GE.0) EXC=-2.D0*PI34/RS/ALPHA
GO TO (40,10,20,30) IFL
10 EXC=EXC+ECGL(RS,1)
GO TO 40
20 EXC=EXC+ECHL(RS,1)
GO TO 40
30 EXC=EXC+ECCA(RS,1)
40 IF(ZETA.EQ.0.D0) RETURN
EXCF=0.D0
IF(IXC.GE.0) EXCF=-TW43*PI34/RS/ALPHA
GO TO (41,11,21,31) IFL
11 EXCF=EXCF+ECGL(RS,2)
GO TO 41
21 EXCF=EXCF+ECHL(RS,2)
GO TO 41
31 EXCF=EXCF+ECCA(RS,2)
41 EXC=EXC+FZ(ZETA,RS)*(EXCF-EXC)
RETURN
END FUNCTION EXC
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION FZ(ZETA,RS)
!---------------------------------------------------------------
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
INTEGER IFLAG
SAVE IFLAG
IFLAG = 0
IF(IFLAG.EQ.0) THEN
FZ=FZ0(ZETA)
ELSE
FZ=FZX(ZETA,RS)
END IF
RETURN
END FUNCTION FZ
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION FZ0(ZETA)
!---------------------------------------------------------------
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
REAL(dbl), PARAMETER :: F43= 4.D0/3.D0
REAL(dbl), PARAMETER :: DEN= 5.198420997897450D-01
ZS=ABS(ZETA)
IF( ZS-1.D0 .GT. 1.0D-10 ) STOP ' ZETA OUT OF RANGE IN FZ0'
ZS=MIN(ZS,1.D0)
IF(ZS.EQ.0.D0) THEN
FZ0=0.D0
ELSE IF(ZS.EQ.1.D0) THEN
FZ0=1.D0
ELSE
FZ0=( (1.D0+ZS)**F43 + (1.D0-ZS)**F43 -2.D0 ) / DEN
END IF
RETURN
END FUNCTION FZ0
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION FZX(ZETA,RS)
!---------------------------------------------------------------
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
REAL(dbl) :: RSX(9),A0(9),A1(9),A2(9),A3(9),B0(9),B1(9),B2(9) &
,B3(9),D0(9),D1(9),D2(9),D3(9)
REAL(dbl), PARAMETER :: FZ2=1.709920933D0
REAL(dbl), PARAMETER :: RSXM=15.D0
DATA RSX / 5.0000D-01, 1.0000D+00, 2.0000D+00, 3.0000D+00, &
4.0000D+00, 5.0000D+00, 6.0000D+00, 7.5000D+00, 1.0000D+01/ &
,A0 / 9.7980D+01, 7.7100D+01, 5.7860D+01, 4.7620D+01, &
4.0920D+01, 3.6090D+01, 3.2390D+01, 2.8180D+01, 2.3290D+01/ &
,A1 /-5.3220D+01,-3.1618D+01,-1.2133D+01,-8.2905D+00, &
-5.5253D+00,-4.1984D+00,-3.2712D+00,-2.4165D+00,-1.5615D+00/ &
,A2 / 2.5555D+01, 1.7649D+01, 1.8363D+00, 2.0062D+00, &
7.5897D-01, 5.6795D-01, 3.5924D-01, 2.1057D-01, 1.3144D-01/ &
,A3 /-5.2708D+00,-5.2708D+00, 5.6645D-02,-4.1574D-01, &
-6.3674D-02,-6.9568D-02,-3.3039D-02,-1.0550D-02,-1.0550D-02/ &
,B0 / 2.3870D-01, 1.9430D-01, 1.4660D-01, 1.1920D-01, &
1.0080D-01, 8.7500D-02, 7.7400D-02, 6.6100D-02, 5.3300D-02/ &
,B1 /-1.0887D-01,-7.0850D-02,-3.3049D-02,-2.2253D-02, &
-1.5338D-02,-1.1496D-02,-8.8781D-03,-6.4154D-03,-4.0142D-03/ &
,B2 / 4.4399D-02, 3.1650D-02, 6.1515D-03, 4.6443D-03, &
2.2714D-03, 1.5702D-03, 1.0478D-03, 5.9401D-04, 3.6650D-04/ &
,B3 /-8.4994D-03,-8.4994D-03,-5.0241D-04,-7.9097D-04, &
-2.3372D-04,-1.7416D-04,-1.0083D-04,-3.0335D-05,-3.0335D-05/
DATA D0 / 7.0977D+01, 5.3849D+01, 3.8797D+01, 3.1170D+01, &
2.6345D+01, 2.2952D+01, 2.0408D+01, 1.7570D+01, 1.4348D+01/ &
,D1 /-4.4080D+01,-2.5573D+01,-9.0941D+00,-6.0874D+00, &
-3.9121D+00,-2.9180D+00,-2.2268D+00,-1.6129D+00,-1.0129D+00/ &
,D2 / 2.1929D+01, 1.5084D+01, 1.3947D+00, 1.6120D+00, &
5.6329D-01, 4.3084D-01, 2.6036D-01, 1.4891D-01, 9.1108D-02/ &
,D3 /-4.5632D+00,-4.5632D+00, 7.2443D-02,-3.4957D-01, &
-4.4150D-02,-5.6827D-02,-2.4766D-02,-7.7072D-03,-7.7072D-03/
IF(RS.LT.RSX(1).OR.RS.GT.RSXM) THEN
FZX=FZ0(ZETA)
ELSE
I=1
10 I=I+1
IF(RSX(I).LT.RS) GO TO 10
I=I-1
DRS=RS-RSX(I)
A=A0(I)+DRS*(A1(I)+DRS*(A2(I)+DRS*A3(I)))
B=B0(I)+DRS*(B1(I)+DRS*(B2(I)+DRS*B3(I)))
D=D0(I)+DRS*(D1(I)+DRS*(D2(I)+DRS*D3(I)))
DRPA=A*FZ0(ZETA)*(1.D0+B*ZETA**4)/FZ2
FZX=DRPA/D
END IF
RETURN
END FUNCTION FZX
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION VCCA(RS,IFLG)
!---------------------------------------------------------------
!
! CEPERLEY & ALDER'S CORRELATION POTENTIAL.
! AFTER J.P. PERDEW & A. ZUNGER PRB 23, 5048 (1981)
! IFLG = 1: PARAMAGNETIC RESULTS
! = 2: FERROMAGNETIC RESULTS
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
REAL(dbl) :: A(2),B(2),C(2),D(2),BT1(2),BT2(2)
REAL(dbl) :: X76, X43, AP, BP, CP, DP, AF, BF, CF, DF
REAL(dbl) :: BP1, CP1, DP1, BF1, CF1, DF1
PARAMETER(X76=7.D0/6.D0, X43=4.D0/3.D0, &
AP=0.03110D0*2.D0, BP=-0.0480D0*2.D0, &
CP=0.0020D0*2.D0, DP=-0.0116D0*2.D0, &
AF=0.01555D0*2.D0, BF=-0.0269D0*2.D0, &
CF=0.0007D0*2.D0, DF=-0.0048D0*2.D0, &
BP1=BP-AP/3.D0, CP1=2.D0*CP/3.D0, DP1=(2.D0*DP-CP)/3.D0, &
BF1=BF-AF/3.D0, CF1=2.D0*CF/3.D0, DF1=(2.D0*DF-CF)/3.D0 )
DATA A/AP ,AF /, B/BP1,BF1/, C/CP1,CF1/, D/DP1,DF1/, &
BT1/1.0529D0,1.3981D0/, BT2/0.3334D0,0.2611D0/
!
IF(RS.LE.1.D0) THEN
RSL=LOG(RS)
VCCA=A(IFLG)*RSL+B(IFLG)+C(IFLG)*RS*RSL+D(IFLG)*RS
ELSE
RSQ=SQRT(RS)
VCCA=ECCA(RS,IFLG)*(1.D0+X76*BT1(IFLG)*RSQ+X43*BT2(IFLG)*RS)/ &
(1.D0+ BT1(IFLG)*RSQ+ BT2(IFLG)*RS)
END IF
RETURN
END FUNCTION VCCA
!
!***************************************************************
!***** *****
!***** THIS SECTION CONTAINS SEVERAL PARAMETRIZATIONS *****
!***** FOR THE EXCHANGE-CORRELATION ENERGY AND POTENTIAL *****
!***** *****
!***************************************************************
!
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION VCGL(RS,IFLG)
!---------------------------------------------------------------
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
REAL(dbl), PARAMETER :: CP=0.0666D0
REAL(dbl), PARAMETER :: RP=11.4D0
VCGL=-CP*LOG(1.D0+RP/RS)
RETURN
END FUNCTION VCGL
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION VCHL(RS,IFLG)
!---------------------------------------------------------------
!
! INTERPOLATION FORMULA FOR THE CORRELATION POTENTIAL AFTER
! L. HEDIN & S. LUNDQVIST J.PHYS. C4, 2064 (1971).
! IFLG=1 : POTENTIEL DE CORRELATION PARAMAGNETIQUE;
! IFLG=2 : POTENTIEL DE CORRELATION FERROMAGNETIQUE (MISAWA).
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
REAL(dbl) :: R(2),C(2)
DATA R/21.0D0,52.917D0/, C/0.045D0,0.0225D0/
X=RS/R(IFLG)
VCHL=-C(IFLG)*LOG(1.0D0+1.D0/X)
RETURN
END FUNCTION VCHL
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION VX(RS,IFLG)
!---------------------------------------------------------------
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
REAL(dbl) :: CX(2)
DATA CX/1.221774115D0, 1.539338926D0/
VX=-CX(IFLG)/RS
RETURN
END FUNCTION VX
!
!---------------------------------------------------------------
DOUBLE PRECISION FUNCTION VXC0(RHO,IXC)
!---------------------------------------------------------------
!
! INTERPOLATION FORMULA FOR THE EXCHANGE-CORRELATION POTENTIAL
! IN THE PARAMAGNETIC CASE:
! IXC = 0 EXCHANGE ONLY (KOHN & SHAM)
! +/- 1 GUNNARSON & LUNDQVIST
! +/- 2 HEDIN & LUNDQVIST
! +/- 3 CEPERLEY & ALDER
! > 0 EXCHANGE + CORRELATION
! < 0 CORRELATION ONLY
!
IMPLICIT REAL(dbl) (A-H,O-Z) , INTEGER (I-N)
REAL(dbl) :: PI34, ALPHA, THIRD, CX, SMALL
PARAMETER(PI34= 0.75D0 / 3.141592653589793D+00 &
,ALPHA=0.5210617612D+00, THIRD=1.D0/3.D0 &
,CX=8.D0*PI34*THIRD/ALPHA, SMALL=1.0D-10)
!
IC=IABS(IXC)
IF(IC.GT.3) STOP ' IXC NOT ALLOWED IN VXC0'
VXC0=0.D0
IF(RHO.LE.SMALL) RETURN
RS=(PI34/RHO)**THIRD
IFL=IC+1
IF(IXC.GE.0) VXC0=-CX/RS
GO TO (40,10,20,30) IFL
10 VXC0=VXC0+VCGL(RS,1)
GO TO 40
20 VXC0=VXC0+VCHL(RS,1)
GO TO 40
30 VXC0=VXC0+VCCA(RS,1)
40 RETURN
END FUNCTION VXC0
!
!===============================================================
!
SUBROUTINE XC(IRO,RO1,RO2,EX,EC,VX1,VX2,VC1,VC2,T58)
!
IMPLICIT REAL(dbl) (A-H,O-Z), INTEGER (I-N)
IF(IRO.NE.1) STOP ' WRONG IRO IN XC'
EX=EXC(RO1,0.D0,0)/2.D0
EC=EXC(RO1,0.D0,-3)/2.D0
VX1=VXC0(RO1,0)/2.D0
VX2=0.D0
VC1=VXC0(RO1,-3)/2.D0
VC2=0.D0
RETURN
END SUBROUTINE XC
!
!===============================================================
!
SUBROUTINE GENXC
!
!........THIS PROGRAM GENERATES THE FILE TABXC. THE FILE TABXC
!........CONTAINS A TABULATION OF XC ENERGIES AND POTENTIALS ON
!........A GRID CONSISTING OF NARRAY ELEMENTS UNIFORMLY SPACED
!........BETWEEN 0 AND RMAX.
!........RECOMMENDED VALUES FOR SILICON: RMAX=0.2, NARRAY=10000
!........RECOMMENDED VALUES FOR GALLIUM: RMAX=0.2, NARRAY=10000
!........RECOMMENDED VALUES FOR CARBON : RMAX=2.0, NARRAY=40000
!........IBM VERSION (2.AUG.85).
!
IMPLICIT NONE
REAL(dbl) RR,EX,EC,VX,VC,DUM
INTEGER I
call allocate_lda
! WRITE( stdout,10) RMXXC,NARRAY
DDRO=RMXXC/REAL(NARRAY-1)
! WRITE( stdout,40)
DO I=0,NARRAY-1
RR=REAL(I)*DDRO
CALL XC(1,RR,DUM,EX,EC,VX,DUM,VC,DUM,DUM)
EEXC(I)=EX+EC
VVXC(I)=VX+VC
! IF(MOD(I,100).EQ.0) THEN
! WRITE( stdout,30) RR,EX,EC,EEXC(I),VX,VC,VVXC(I)
! END IF
END DO
!
10 FORMAT(//' GENXC: RMAX, NARRAY = ',E10.4,3X,I5//)
40 FORMAT(3X,'RHO',8X,'EX',8X,'EC',8X,'EEXC',8X,'VX',8X,'VC', 8X,'VVXC')
30 FORMAT(1X,F9.6,6(1X,F9.6))
!
RETURN
END SUBROUTINE GENXC
END MODULE local_density_approximation

1124
CPV/lsd.f90
View File

File diff suppressed because it is too large Load Diff

2
CPV/main.f90
View File

@ -94,7 +94,6 @@
USE polarization, ONLY: deallocate_polarization, ddipole
USE energies, ONLY: dft_energy_type, debug_energies
USE recvecs_indexes, ONLY: deallocate_recvecs_indexes
USE exchange_correlation, ONLY: deallocate_lda
USE turbo, ONLY: tturbo, deallocate_turbo
USE nose_ions, ONLY: movenosep, nosep_velocity, update_nose_ions
USE nose_electrons, ONLY: movenosee, nosee_velocity, update_nose_electrons
@ -830,7 +829,6 @@
CALL deallocate_recvecs()
CALL deallocate_recvecs_indexes()
CALL deallocate_rvecs(gv)
CALL deallocate_lda( )
CALL deallocate_phfac(eigr)
CALL deallocate_electrons
IF(tdipole) THEN

52
CPV/potentials.f90
View File

@ -257,8 +257,7 @@
USE energies, ONLY: total_energy, dft_energy_type
USE stress, ONLY: pstress
USE gvecw, ONLY: tecfix
USE exchange_correlation, ONLY: tgc, vofxc_lda, vofxc_lsd, v2gc, &
exch_corr_energy
USE exchange_correlation, ONLY: tgc, exch_corr_energy
USE charge_density, ONLY: gradrho
USE non_local_core_correction, ONLY: add_core_charge, core_charge_forces
USE chi2, ONLY: tchi2, rhochi, allocate_chi2, deallocate_chi2
@ -315,15 +314,13 @@
REAL(dbl), ALLOCATABLE :: rhoetr(:,:,:,:)
REAL(dbl), ALLOCATABLE :: fion_vdw(:,:)
REAL(dbl), ALLOCATABLE :: grho(:,:,:,:,:)
REAL(dbl), ALLOCATABLE :: v2xc(:,:,:,:)
REAL(dbl), ALLOCATABLE :: v2xc2(:,:,:)
REAL(dbl), ALLOCATABLE :: v2xc(:,:,:,:,:)
REAL(dbl), ALLOCATABLE :: fion(:,:)
REAL(dbl), ALLOCATABLE :: self_rho(:,:,:,:)
REAL(dbl), ALLOCATABLE :: self_vpot(:,:,:,:)
REAL(dbl), ALLOCATABLE :: self_grho(:,:,:,:,:)
REAL(dbl), ALLOCATABLE :: self_v2xc(:,:,:,:)
REAL(dbl), ALLOCATABLE :: self_v2xc2(:,:,:)
REAL(dbl), ALLOCATABLE :: self_v2xc(:,:,:,:,:)
REAL(dbl) :: pail(3,3)
@ -395,12 +392,10 @@
IF(tgc) THEN
ALLOCATE( grho( nr1x, nr2x, nr3x, 3, nspin ) )
ALLOCATE( v2xc( nr1x, nr2x, nr3x, 2*nspin-1) )
ALLOCATE( v2xc2(nr1x, nr2x, nr3x) )
ALLOCATE( v2xc( nr1x, nr2x, nr3x, nspin, nspin) )
ELSE
ALLOCATE( grho( 1, 1, 1, 1, 1 ) )
ALLOCATE( v2xc( 1, 1, 1, 1 ) )
ALLOCATE( v2xc2(1, 1, 1) )
ALLOCATE( v2xc( 1, 1, 1, 1, 1 ) )
END IF
ALLOCATE( rhoeg(gv%ng_l, nspin) )
@ -418,13 +413,10 @@
ALLOCATE(self_grho( nr1x, nr2x, nr3x, 3, nspin ), STAT = ierr)
IF( ierr /= 0 ) CALL errore(' vofrhos ', ' allocating self_grho ', ierr)
ALLOCATE(self_v2xc( nr1x, nr2x, nr3x, 2*nspin-1), STAT = ierr)
ALLOCATE(self_v2xc( nr1x, nr2x, nr3x, nspin, nspin ), STAT = ierr)
IF( ierr /= 0 ) CALL errore(' vofrhos ', ' allocating self_v2xc ', ierr)
ALLOCATE(self_v2xc2( nr1x, nr2x, nr3x), STAT = ierr)
IF( ierr /= 0 ) CALL errore(' vofrhos ', ' allocating self_v2xc2 ', ierr)
self_v2xc = 0.D0
self_v2xc2 = 0.D0
END IF !on tgc
@ -509,8 +501,7 @@
! ... Self-interaction correction --- Excor Part
!In any case calculate the Excor part with rhoetr
CALL exch_corr_energy(rhoetr, rhoetg, grho, vpot, &
sxcp, vxc, v2xc, v2xc2, gv)
CALL exch_corr_energy(rhoetr, rhoetg, grho, vpot, sxcp, vxc, v2xc, gv)
IF( ttsic ) THEN
IF( ionode ) THEN
@ -542,12 +533,12 @@
self_grho(:,:,:,:,2) = 0.D0
ENDIF
CALL exch_corr_energy(self_rho, rhoetg, self_grho, self_vpot, &
self_sxcp, self_vxc, self_v2xc, self_v2xc2, gv)
self_sxcp, self_vxc, self_v2xc, gv)
vpot(:,:,:,1) = vpot(:,:,:,1) - self_vpot(:,:,:,1)
vpot(:,:,:,2) = vpot(:,:,:,2) + self_vpot(:,:,:,1)
IF (tgc) THEN
v2xc(:,:,:,1) = v2xc(:,:,:,1) - self_v2xc(:,:,:,1)
v2xc(:,:,:,2) = v2xc(:,:,:,2) + self_v2xc(:,:,:,1)
v2xc(:,:,:,1,1) = v2xc(:,:,:,1,1) - self_v2xc(:,:,:,1,1)
v2xc(:,:,:,2,2) = v2xc(:,:,:,2,2) + self_v2xc(:,:,:,1,1)
ENDIF
! write(stdout,*) 'da exc la parte da rhodwn rhodwn',self_sxcp
@ -575,15 +566,15 @@
! write(stdout,*)'DA XC SELF_RHO2',self_rho(:,:,:,2)
CALL exch_corr_energy(self_rho, rhoetg, self_grho, self_vpot, &
self_sxcp, self_vxc, self_v2xc, self_v2xc2, gv)
self_sxcp, self_vxc, self_v2xc, gv)
vpot (:,:,:,2) = self_vpot(:,:,:,2) + self_vpot(:,:,:,1)
vpot (:,:,:,1) = 0.d0
write(stdout,*) 'da XC with tgc: E_XC==',self_sxcp
IF (tgc) THEN
v2xc(:,:,:,1) = 0.d0
v2xc(:,:,:,2) = self_v2xc(:,:,:,2) + self_v2xc(:,:,:,1)
v2xc(:,:,:,1,1) = 0.d0
v2xc(:,:,:,2,2) = self_v2xc(:,:,:,2,2) + self_v2xc(:,:,:,1,1)
ENDIF
write(stdout,*) 'da exc la parte da rhodwn rhodwn',self_sxcp
edft%self_sxc= sxcp - self_sxcp
@ -607,19 +598,7 @@
END IF
IF ( ttstress ) THEN
strvxc = 0.0d0
DO ispin = 1, nspin
DO k = 1, nr3_l
DO j = 1, nr2_l
strvxc = strvxc + &
DDOT ( nr1_l, vpot(1,j,k,ispin), 1, rhoetr(1,j,k,ispin), 1 )
END DO
END DO
END DO
vxc = strvxc ! ... SUM ( u(r) * rho(r) )
strvxc = (edft%sxc - strvxc) * omega / REAL( nr1_g * nr2_g * nr3_g )
! strvxc = SUM(v2xc2) * omega / REAL(nr1_g*nr2_g*nr3_g)
strvxc = ( edft%sxc - vxc ) * omega / REAL( nr1_g * nr2_g * nr3_g )
END IF
IF( ANY( ps%tnlcc ) ) THEN
@ -698,7 +677,6 @@
IF( ALLOCATED( self_grho ) ) DEALLOCATE( self_grho )
IF( ALLOCATED( self_v2xc ) ) DEALLOCATE( self_v2xc )
IF( ALLOCATED( self_v2xc2 ) ) DEALLOCATE( self_v2xc2 )
IF( ALLOCATED( self_vpot ) ) DEALLOCATE( self_vpot )
IF( ALLOCATED( self_rho ) ) DEALLOCATE( self_rho )
@ -781,7 +759,7 @@
END IF
END IF
DEALLOCATE( rhoeg, rhoetg, rhoetr, grho, v2xc, v2xc2, fion )
DEALLOCATE( rhoeg, rhoetg, rhoetr, grho, v2xc, fion )
IF(tchi2) THEN
CALL deallocate_chi2

1
CPV/problem_size.f90
View File

@ -25,7 +25,6 @@
use brillouin, only: get_kpoints_number
use reciprocal_vectors, only: ngwx, ngmx, ngmt
use pseudopotential, only: lmax, ngh
use exchange_correlation, only: narray
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE fft_base, ONLY: dfftp, dffts

174
CPV/restart_sub.f90
View File

@ -16,24 +16,24 @@ MODULE from_restart_module
PUBLIC :: from_restart
INTERFACE from_restart
MODULE PROCEDURE from_restart_cp, from_restart_fpmd
MODULE PROCEDURE from_restart_cp, from_restart_fpmd, from_restart_sm
END INTERFACE
CONTAINS
!=----------------------------------------------------------------------------=!
SUBROUTINE from_restart_cp( sfac, eigr, ei1, ei2, ei3, bec, becdr, tfirst, eself, fion, &
taub, irb, eigrb, b1, b2, b3, nfi, rhog, rhor, rhos, rhoc, enl, ekin, stress, &
detot, enthal, etot, lambda, lambdam, lambdap, ema0bg, dbec, delt, &
bephi, becp, velh, dt2bye, iforce, fionm, nbeg, xnhe0, xnhem, vnhe, ekincm )
SUBROUTINE from_restart_cp( sfac, eigr, ei1, ei2, ei3, bec, becdr, tfirst, fion, &
taub, irb, eigrb, b1, b2, b3, nfi, rhog, rhor, rhos, rhoc, stress, &
detot, enthal, lambda, lambdam, lambdap, ema0bg, dbec, &
bephi, becp, velh, dt2bye, fionm, ekincm )
USE kinds, ONLY: dbl
USE control_flags, ONLY: tranp, trane, trhor, iprsta, tpre, tzeroc
USE control_flags, ONLY: tzerop, tzeroe, tfor, thdyn, lwf, tprnfor, tortho
USE control_flags, ONLY: amprp, taurdr, ortho_eps, ortho_max
USE control_flags, ONLY: amprp, taurdr, ortho_eps, ortho_max, nbeg
USE ions_positions, ONLY: taus, tau0, tausm, vels, velsm, ions_hmove
USE ions_base, ONLY: na, nsp, randpos, zv, ions_vel, pmass
USE ions_base, ONLY: na, nsp, randpos, zv, ions_vel, pmass, iforce
USE cell_base, ONLY: ainv, h, s_to_r, ibrav, omega, press, hold, r_to_s, deth
USE cell_base, ONLY: wmass, iforceh, cell_force, cell_hmove
USE efcalc, ONLY: ef_force
@ -49,33 +49,34 @@ CONTAINS
USE reciprocal_vectors, ONLY: gstart
USE wave_base, ONLY: wave_verlet
USE cvan, ONLY: nvb
USE ions_nose, ONLY: xnhp0, xnhpm, vnhp
USE cell_nose, ONLY: xnhh0, xnhhm, vnhh
USE ions_nose, ONLY: xnhp0, xnhpm, vnhp
USE cp_electronic_mass, ONLY: emass, emaec => emass_cutoff
USE efield_module, ONLY: efield_berry_setup, tefield
USE runcp_module, ONLY: runcp_uspp
USE wave_constrains, ONLY: interpolate_lambda
USE energies, ONLY: eself, enl, etot, ekin
USE time_step, ONLY: delt
use electrons_nose, only: xnhe0, xnhem, vnhe
USE cell_nose, ONLY: xnhh0, xnhhm, vnhh, cell_nosezero
COMPLEX(kind=8) :: eigr(:,:,:), ei1(:,:,:), ei2(:,:,:), ei3(:,:,:)
COMPLEX(kind=8) :: eigrb(:,:,:)
INTEGER :: irb(:,:,:)
REAL(kind=8) :: bec(:,:), fion(:,:), becdr(:,:,:), fionm(:,:)
REAL(kind=8) :: eself
REAL(kind=8) :: taub(:,:)
REAL(kind=8) :: b1(:), b2(:), b3(:)
INTEGER :: irb(:,:,:)
INTEGER :: nfi, iforce(:,:), nbeg
INTEGER :: nfi
LOGICAL :: tfirst
COMPLEX(kind=8) :: sfac(:,:)
COMPLEX(kind=8) :: rhog(:,:)
REAL(kind=8) :: rhor(:,:), rhos(:,:), rhoc(:), enl, ekin
REAL(kind=8) :: stress(:,:), detot(:,:), enthal, etot
REAL(kind=8) :: rhor(:,:), rhos(:,:), rhoc(:)
REAL(kind=8) :: stress(:,:), detot(:,:), enthal, ekincm
REAL(kind=8) :: lambda(:,:), lambdam(:,:), lambdap(:,:)
REAL(kind=8) :: ema0bg(:)
REAL(kind=8) :: dbec(:,:,:,:)
REAL(kind=8) :: delt
REAL(kind=8) :: bephi(:,:), becp(:,:)
REAL(kind=8) :: velh(:,:)
REAL(kind=8) :: dt2bye, xnhe0, xnhem, vnhe, ekincm
REAL(kind=8) :: velh(3,3)
REAL(kind=8) :: dt2bye
REAL(kind=8), ALLOCATABLE :: emadt2(:), emaver(:)
@ -103,7 +104,7 @@ CONTAINS
END IF
call s_to_r( taus, tau0, na, nsp, h )
!
if( trane .and. trhor ) then
call prefor( eigr, betae )
call graham( betae, bec, c0 )
@ -125,7 +126,7 @@ CONTAINS
lambdam = lambda
WRITE( stdout, * ) 'Electronic velocities set to zero'
end if
!
call phfac( tau0, ei1, ei2, ei3, eigr )
call strucf( ei1, ei2, ei3, sfac )
call formf( tfirst, eself )
@ -136,21 +137,23 @@ CONTAINS
call efield_berry_setup( eigr, tau0 )
end if
if ( tzerop .or. tzeroe .or. taurdr ) then
call initbox ( tau0, taub, irb )
call phbox( taub, eigrb )
!
if( lwf ) then
call get_wannier_center( tfirst, c0, bec, becdr, eigr, eigrb, taub, irb, ibrav, b1, b2, b3 )
end if
!
call rhoofr ( nfi, c0, irb, eigrb, bec, rhovan, rhor, rhog, rhos, enl, ekin )
!
! put core charge (if present) in rhoc(r)
!
!
! put core charge (if present) in rhoc(r)
!
if ( nlcc_any ) call set_cc( irb, eigrb, rhoc )
if( lwf ) then
call write_charge_and_exit( rhog )
call ef_tune( rhog, tau0 )
@ -159,61 +162,68 @@ CONTAINS
call vofrho( nfi, rhor, rhog, rhos, rhoc, tfirst, tlast, &
& ei1, ei2, ei3, irb, eigrb, sfac, tau0, fion )
call compute_stress( stress, detot, h, omega )
if( lwf ) then
call wf_options( tfirst, nfi, c0, rhovan, bec, becdr, eigr, eigrb, taub, irb, &
ibrav, b1, b2, b3, rhor, rhog, rhos, enl, ekin )
end if
call newd( rhor, irb, eigrb, rhovan, fion )
call prefor( eigr, betae )
!
if( tzeroe ) then
CALL runcp_uspp( nfi, fccc, ccc, ema0bg, dt2bye, rhos, bec, c0, cm, restart = .TRUE. )
end if
!
! nlfq needs deeq bec
!
!
! nlfq needs deeq bec
!
if( tfor .or. tprnfor ) call nlfq( c0, eigr, bec, becdr, fion )
!
if( tfor .or. thdyn ) then
CALL interpolate_lambda( lambdap, lambda, lambdam )
endif
!
! calphi calculates phi
! the electron mass rises with g**2
!
!
! calphi calculates phi
! the electron mass rises with g**2
!
call calphi( c0, ema0bg, bec, betae, phi )
!
! begin try and error loop (only one step!)
!
! nlfl and nlfh need: lambda (guessed) becdr
!
!
! begin try and error loop (only one step!)
!
! nlfl and nlfh need: lambda (guessed) becdr
!
if ( tfor .or. tprnfor ) call nlfl( bec, becdr, lambda, fion )
if( tpre ) then
call nlfh( bec, dbec, lambda )
endif
if( tortho ) then
call ortho( eigr, cm, phi, lambda, bigr, iter, dt2bye, ortho_eps, ortho_max, delt, bephi, becp )
else
call graham( betae, bec, cm )
endif
!
if ( tortho ) call updatc( dt2bye, lambda, phi, bephi, becp, bec, cm )
call calbec ( nvb+1, nsp, eigr, cm, bec )
if ( tpre ) call caldbec( 1, nsp, eigr, cm )
!
if( thdyn ) then
call cell_force( fcell, ainv, stress, omega, press, wmass )
call cell_hmove( h, hold, delt, iforceh, fcell )
call invmat( 3, h, ainv, deth )
endif
if( tfor ) then
if( lwf ) then
call ef_force( fion, na, nsp, zv )
@ -225,20 +235,22 @@ CONTAINS
if (tpre) call caldbec( 1, nsp, eigr, c0 )
endif
xnhp0=0.
xnhpm=0.
vnhp =0.
fionm=0.
xnhp0=0.0d0
xnhpm=0.0d0
vnhp =0.0d0
fionm=0.0d0
CALL ions_vel( vels, taus, tausm, na, nsp, delt )
xnhh0(:,:)=0.
xnhhm(:,:)=0.
vnhh (:,:) =0.
velh (:,:)=(h(:,:)-hold(:,:))/delt
!
! ======================================================
! kinetic energy of the electrons
! ======================================================
CALL cell_nosezero( vnhh, xnhh0, xnhhm )
velh = ( h - hold ) / delt
! ======================================================
! kinetic energy of the electrons
! ======================================================
call elec_fakekine( ekincm, ema0bg, emass, c0, cm, ngw, n, delt )
@ -250,12 +262,70 @@ CONTAINS
lambdam(:,:)=lambda(:,:)
end if
return
end subroutine
!=----------------------------------------------------------------------------=!
SUBROUTINE from_restart_sm &
( tfirst, taus, tau0, h, eigr, bec, c0, cm, ei1, ei2, ei3, sfac, eself )
USE kinds, ONLY: dbl
USE control_flags, ONLY: trane, trhor, iprsta, tpre
use uspp, only: betae => vkb
use ions_base, only: na, nsp
use electrons_base, only: n => nbnd
use io_global, only: stdout
use cell_base, only: s_to_r
use cpr_subroutines, only: print_atomic_var
IMPLICIT NONE
LOGICAL :: tfirst
REAL(dbl) :: taus( :, : ), tau0( :, : )
REAL(dbl) :: h(3,3)
COMPLEX(dbl) :: eigr( :, :, : )
REAL(dbl) :: bec( :, : )
COMPLEX(dbl) :: c0( :, : )
COMPLEX(dbl) :: cm( :, : )
COMPLEX(dbl) :: ei1( :, :, : )
COMPLEX(dbl) :: ei2( :, :, : )
COMPLEX(dbl) :: ei3( :, :, : )
COMPLEX(kind=8) :: sfac(:,:)
REAL(dbl) :: eself
INTEGER :: j
CALL s_to_r( taus, tau0, na, nsp, h )
!
if(trane.and.trhor) then
call prefor(eigr,betae)
call graham(betae, bec, c0 )
cm(:, 1:n) = c0(:, 1:n)
endif
!
if(iprsta.gt.2) then
call print_atomic_var( taus, na, nsp, ' read: taus ' )
WRITE( stdout,*) ' read: cell parameters h '
WRITE( stdout,*) (h(1,j),j=1,3)
WRITE( stdout,*) (h(2,j),j=1,3)
WRITE( stdout,*) (h(3,j),j=1,3)
endif
!
call phfac(tau0,ei1,ei2,ei3,eigr)
call strucf(ei1,ei2,ei3,sfac)
call formf(tfirst,eself)
call calbec (1,nsp,eigr,c0,bec)
if (tpre) call caldbec(1,nsp,eigr,c0)
RETURN
END SUBROUTINE
!=----------------------------------------------------------------------------=!
SUBROUTINE from_restart_fpmd( nfi, acc, gv, kp, ps, rhoe, desc, cm, c0, cdesc, &

65
CPV/smcp.f90
View File

@ -121,7 +121,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
use ions_positions, only: ions_hmove, ions_move
use ions_nose, only: gkbt, qnp, ions_nosevel, ions_noseupd, tempw
USE cell_base, ONLY: cell_kinene, cell_move, cell_gamma, cell_hmove
USE cell_nose, ONLY: cell_nosevel, cell_noseupd, qnh, temph
USE cell_nose, ONLY: cell_nosevel, cell_noseupd, qnh, temph, cell_nosezero
USE gvecw, ONLY: ecutw
USE gvecp, ONLY: ecutp
@ -151,6 +151,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
USE smd_rep
USE smd_ene
USE from_restart_module, ONLY: from_restart
!
!
implicit none
@ -583,6 +584,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
#endif
!
666 continue
!
!
@ -650,6 +652,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
!
!
!
INI_REP_LOOP : DO sm_k=1,smpm ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !
!
@ -905,29 +908,31 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
& tps, mat_z, f )
call from_restart( tfirst, rep(sm_k)%taus, rep(sm_k)%tau0, h, eigr, &
rep_el(sm_k)%bec, rep_el(sm_k)%c0, rep_el(sm_k)%cm, ei1, ei2, ei3, sfac, eself )
!
CALL s_to_r( rep(sm_k)%taus, rep(sm_k)%tau0, na, nsp, h )
!
if(trane.and.trhor) then
call prefor(eigr,betae)
call graham(betae,rep_el(sm_k)%bec,rep_el(sm_k)%c0)
rep_el(sm_k)%cm(:, 1:n)=rep_el(sm_k)%c0(:, 1:n)
endif
!
if(iprsta.gt.2) then
call print_atomic_var( rep(sm_k)%taus, na, nsp, ' read: taus ' )
WRITE( stdout,*) ' read: cell parameters h '
WRITE( stdout,*) (h(1,j),j=1,3)
WRITE( stdout,*) (h(2,j),j=1,3)
WRITE( stdout,*) (h(3,j),j=1,3)
endif
!
call phfac(rep(sm_k)%tau0,ei1,ei2,ei3,eigr)
call strucf(ei1,ei2,ei3,sfac)
call formf(tfirst,eself)
call calbec (1,nsp,eigr,rep_el(sm_k)%c0,rep_el(sm_k)%bec)
if (tpre) call caldbec(1,nsp,eigr,rep_el(sm_k)%c0)
! CALL s_to_r( rep(sm_k)%taus, rep(sm_k)%tau0, na, nsp, h )
!
! !
! if(trane.and.trhor) then
! call prefor(eigr,betae)
! call graham(betae,rep_el(sm_k)%bec,rep_el(sm_k)%c0)
! rep_el(sm_k)%cm(:, 1:n)=rep_el(sm_k)%c0(:, 1:n)
! endif
! !
! if(iprsta.gt.2) then
! call print_atomic_var( rep(sm_k)%taus, na, nsp, ' read: taus ' )
! WRITE( stdout,*) ' read: cell parameters h '
! WRITE( stdout,*) (h(1,j),j=1,3)
! WRITE( stdout,*) (h(2,j),j=1,3)
! WRITE( stdout,*) (h(3,j),j=1,3)
! endif
! !
! call phfac(rep(sm_k)%tau0,ei1,ei2,ei3,eigr)
! call strucf(ei1,ei2,ei3,sfac)
! call formf(tfirst,eself)
! call calbec (1,nsp,eigr,rep_el(sm_k)%c0,rep_el(sm_k)%bec)
! if (tpre) call caldbec(1,nsp,eigr,rep_el(sm_k)%c0)
!
end if ! <<<<<<<<<<<<<<<<<<<<<<<< !
@ -964,6 +969,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
ENDDO INI_REP_LOOP ! <<<<<<<<<<<<<<<<<<<<<<< !
!
!
!
IF(smlm .and. nbeg < 0) THEN
! .... temp assignment ...
@ -1013,10 +1019,9 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
vnhp(sm_k) =0.
rep(sm_k)%fionm=0.d0
CALL ions_vel( rep(sm_k)%vels, rep(sm_k)%taus, rep(sm_k)%tausm, na, nsp, delt )
xnhh0(:,:)=0.
xnhhm(:,:)=0.
vnhh (:,:) =0.
velh (:,:)=(h(:,:)-hold(:,:))/delt
CALL cell_nosezero( vnhh, xnhh0, xnhhm )
velh = ( h - hold ) / delt
!
! ======================================================
! kinetic energy of the electrons
@ -1298,6 +1303,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
!
ENDDO EL_REP_LOOP ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< !
!_________________________________________________________________________!
!
!
@ -1363,6 +1369,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
& 1x,i3,1x,f8.5, &
& 1x,f8.5, &
& 1x,E12.5,1x,E12.5)
!
!
!________________________________________________________________________!
@ -1535,6 +1542,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
!
ekinp(sm_k) = 0.d0
ekinpr(sm_k) = 0.d0
!
! ionic kinetic energy
!
@ -1803,11 +1811,13 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
END IF
END IF
tstop = tstop .OR. tconv
ENDDO POST_REP_LOOP ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< !
if( (nfi >= nomore) .OR. tstop ) EXIT MAIN_LOOP
END DO MAIN_LOOP
@ -1816,6 +1826,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
!============================= END of main LOOP ======================
!
!
! Here copy relevant physical quantities into the output arrays/variables
! for 0 and sm_p replicas.

25
CPV/stress.f90
View File

@ -90,7 +90,7 @@
! ... declare subroutine arguments
LOGICAL, INTENT(IN) :: prn, tgc
REAL(dbl) :: pail(:,:), desr(:), strvxc
REAL(dbl) :: grho(:,:,:,:,:), v2xc(:,:,:,:)
REAL(dbl) :: grho(:,:,:,:,:), v2xc(:,:,:,:,:)
COMPLEX(dbl) :: rhoeg(:,:), vxc(:,:)
COMPLEX(dbl), INTENT(IN) :: sfac(:,:)
REAL(dbl), INTENT(IN) :: occ(:,:,:)
@ -870,17 +870,14 @@
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nnr
REAL(dbl) :: v2xc(:,:,:,:)
REAL(dbl) :: v2xc(:,:,:,:,:)
REAL(dbl) :: grho(:,:,:,:,:)
REAL(dbl) :: gcpail(6)
REAL(dbl) :: omega
!
REAL(dbl) :: stre, grhoi, grhoj
INTEGER :: i, j, k, ipol, jpol, ic, nxl, nyl, nzl, is, js, nspin
INTEGER, DIMENSION(2,2), PARAMETER :: kk = reshape &
( (/ 1, 3, 3, 2 /), (/ 2, 2 /) )
INTEGER, DIMENSION(2,2), PARAMETER :: nn = reshape &
( (/ 1, 1, 2, 2 /), (/ 2, 2 /) )
INTEGER :: kk(2,2), nn(2,2)
! ...
nxl = SIZE(grho,1)
@ -888,17 +885,25 @@
nzl = SIZE(grho,3)
nspin = SIZE(grho,5)
kk(1,1)=1
kk(1,2)=3
kk(2,1)=3
kk(2,2)=2
nn(1,1)=1
nn(1,2)=2
nn(2,1)=1
nn(2,2)=2
DO ic = 1, 6
ipol = alpha(ic)
jpol = beta(ic)
stre = 0.0d0
DO is = 1, nspin
stre = 0.d0
DO js = 1, nspin
DO k = 1, nzl
DO j = 1, nyl
DO i = 1, nxl
stre = stre + v2xc(i,j,k,kk(is,js)) * &
grho(i,j,k,ipol,nn(is,js)) * grho(i,j,k,jpol,is)
stre = stre + v2xc(i,j,k,is,js) * grho(i,j,k,ipol,js) * grho(i,j,k,jpol,is)
END DO
END DO
END DO
@ -935,7 +940,7 @@
COMPLEX(dbl), INTENT(IN) :: sfac(:,:)
REAL(dbl) :: dexc(:), strvxc
REAL(dbl) :: grho(:,:,:,:,:)
REAL(dbl) :: v2xc(:,:,:,:)
REAL(dbl) :: v2xc(:,:,:,:,:)
REAL(dbl) :: GAgx_L(:,:)
REAL(dbl) :: rhocp(:,:)

2
D3/Makefile
View File

@ -256,7 +256,6 @@ PWOBJS = \
../PW/estimate.o \
../PW/ewald.o \
../PW/fftw.o \
../PW/functionals.o \
../PW/g_psi.o \
../PW/g_psi_mod.o \
../PW/gen_at_dj.o \
@ -285,7 +284,6 @@ PWOBJS = \
../PW/kpoint_grid.o \
../PW/lchk_tauxk.o \
../PW/linmin.o \
../PW/lsda_functionals.o \
../PW/make_pointlists.o \
../PW/mix_pot.o \
../PW/mix_rho.o \

2
Gamma/Makefile
View File

@ -79,7 +79,6 @@ PWOBJS = \
../PW/force_lc.o \
../PW/forces.o \
../PW/force_us.o \
../PW/functionals.o \
../PW/gen_at_dj.o \
../PW/gen_at_dy.o \
../PW/gen_us_dj.o \
@ -109,7 +108,6 @@ PWOBJS = \
../PW/kpoint_grid.o \
../PW/lchk_tauxk.o \
../PW/linmin.o \
../PW/lsda_functionals.o \
../PW/make_pointlists.o \
../PW/mix_pot.o \
../PW/mix_rho.o \

8
Modules/cell_base.f90
View File

@ -690,6 +690,14 @@ CONTAINS
return
end subroutine
subroutine cell_nosezero( vnhh, xnhh0, xnhhm )
real(dbl), intent(out) :: vnhh(3,3), xnhh0(3,3), xnhhm(3,3)
xnhh0=0.0d0
xnhhm=0.0d0
vnhh =0.0d0
return
end subroutine
subroutine cell_nosevel( vnhh, xnhh0, xnhhm, delt, velh, h, hold )
implicit none
real(kind=8), intent(inout) :: vnhh(3,3), velh(3,3)

2
Nmr/Makefile
View File

@ -214,7 +214,6 @@ PWOBJS = \
../PW/error_handler.o \
../PW/estimate.o \
../PW/fftw.o \
../PW/functionals.o \
../PW/g_psi.o \
../PW/g_psi_mod.o \
../PW/gen_at_dj.o \
@ -243,7 +242,6 @@ PWOBJS = \
../PW/kpoint_grid.o \
../PW/lchk_tauxk.o \
../PW/linmin.o \
../PW/lsda_functionals.o \
../PW/make_pointlists.o \
../PW/mix_pot.o \
../PW/mix_rho.o \

2
PH/Makefile
View File

@ -236,7 +236,6 @@ PWOBJS = \
../PW/force_lc.o \
../PW/force_us.o \
../PW/forces.o \
../PW/functionals.o \
../PW/g_psi.o \
../PW/g_psi_mod.o \
../PW/gen_at_dj.o \
@ -269,7 +268,6 @@ PWOBJS = \
../PW/kpoint_grid.o \
../PW/lchk_tauxk.o \
../PW/linmin.o \
../PW/lsda_functionals.o \
../PW/make_pointlists.o \
../PW/mix_pot.o \
../PW/mix_rho.o \

6
PP/Makefile
View File

@ -128,7 +128,6 @@ PWOBJS = \
../PW/ewald.o \
../PW/ewald_dipole.o \
../PW/fftw.o \
../PW/functionals.o \
../PW/g_psi.o \
../PW/g_psi_mod.o \
../PW/gen_at_dj.o \
@ -159,7 +158,6 @@ PWOBJS = \
../PW/kpoint_grid.o \
../PW/lchk_tauxk.o \
../PW/linmin.o \
../PW/lsda_functionals.o \
../PW/make_pointlists.o \
../PW/mix_pot.o \
../PW/mix_rho.o \
@ -291,6 +289,10 @@ pp.x : postproc.o $(PPOBJS)
$(LD) -o $@ postproc.o $(PPOBJS) $(PWOBJS) $(MODULES) $(LDFLAGS)
- ( cd ../bin ; ln -fs ../PP/$@ . )
xctest.x : xctest.o $(PPOBJS)
$(LD) -o $@ xctest.o $(PPOBJS) $(PWOBJS) $(MODULES) $(LDFLAGS)
- ( cd ../bin ; ln -fs ../PP/$@ . )
projwfc.x : projwfc.o $(PPOBJS)
$(LD) -o $@ projwfc.o $(PPOBJS) $(PWOBJS) $(MODULES) $(LDFLAGS)
- ( cd ../bin ; ln -fs ../PP/$@ . )

81
PP/xctest.f90 Normal file
View File

@ -0,0 +1,81 @@
program xctest
USE mp, ONLY: mp_start, mp_end
use kinds, only: dbl
use funct
implicit none
integer, parameter :: nnr = 1000
integer, parameter :: nspin = 2
CALL mp_start()
iexch=1
icorr=3
igcx=1
igcc=3
CALL test_xc( nnr, nspin )
CALL mp_end()
end program
subroutine test_xc( nnr,nspin )
use kinds, only: dbl
use funct
implicit none
integer, intent(in) :: nnr, nspin
real(dbl) :: rhor( nnr, nspin )
real(dbl) :: grhor( nnr, 3, nspin )
real(dbl) :: rhon( nnr, nspin )
real(dbl) :: grhon( nnr, 3, nspin )
real(dbl) :: exc, excn
integer :: ir, is, ipol
do is = 1, nspin
do ir = 1, nnr
rhor( ir, is ) = DBLE( is * ir ) / DBLE( nspin * nnr )
grhor( ir, 1, is ) = 0.13 * rhor( ir, is )
grhor( ir, 2, is ) = 0.43 * rhor( ir, is )
grhor( ir, 3, is ) = 0.73 * rhor( ir, is )
end do
end do
rhon = rhor
grhon = grhor
!
! original CP xc selection
!
if (iexch==1.and.icorr==1.and.igcx==0.and.igcc==0) then
! LDA (Perdew-Zunger)
call expxc(nnr,nspin,rhor,exc)
else if (iexch==1.and.icorr==4.and.igcx==2.and.igcc==2) then
! PW91
call ggapwold(nnr,nspin,grhor,rhor,exc)
else if (iexch==1.and.icorr==3.and.igcx==1.and.igcc==3) then
! BLYP
call ggablyp4(nnr,nspin,grhor,rhor,exc)
else if (iexch==1.and.icorr==4.and.igcx==3.and.igcc==4) then
! PBE
call ggapbe(nnr,nspin,grhor,rhor,exc)
else
call errore('exc-cor','no such exch-corr',1)
end if
!
! Wrapper to PW xc selection
!
call exch_corr_cp(nnr,nspin,grhon,rhon,excn)
!
write(6,*) 'EXC = ', exc, excn
write(6,*) 'VXC '
do is = 1, nspin
do ir = 1, nnr
WRITE(6,100) ir,is,rhor( ir, is ),rhon( ir, is )
end do
end do
do ipol = 1, 3
write(6,*) 'IPOL ', ipol
do is = 1, nspin
do ir = 1, nnr
WRITE(6,100) ir,is,grhor( ir, ipol, is ),grhon( ir, ipol, is )
end do
end do
end do
100 FORMAT( I5, I2, 1X, E15.8, 1X, E15.8 )
end subroutine

2
PW/Makefile
View File

@ -87,7 +87,6 @@ force_hub.o \
force_lc.o \
force_us.o \
forces.o \
functionals.o \
g_psi.o \
g_psi_mod.o \
gen_at_dj.o \
@ -121,7 +120,6 @@ iweights.o \
kpoint_grid.o \
lchk_tauxk.o \
linmin.o \
lsda_functionals.o \
make_pointlists.o \
mix_pot.o \
mix_rho.o \

2
PWCOND/Makefile
View File

@ -136,7 +136,6 @@ PWOBJS = \
../PW/estimate.o \
../PW/ewald.o \
../PW/fftw.o \
../PW/functionals.o \
../PW/g_psi.o \
../PW/g_psi_mod.o \
../PW/gen_at_dj.o \
@ -165,7 +164,6 @@ PWOBJS = \
../PW/kpoint_grid.o \
../PW/lchk_tauxk.o \
../PW/linmin.o \
../PW/lsda_functionals.o \
../PW/make_pointlists.o \
../PW/mix_pot.o \
../PW/mix_rho.o \

4
atomic/Makefile
View File

@ -86,9 +86,7 @@ MODULES = \
../Modules/io_global.o \
../Modules/kind.o \
../Modules/io_files.o \
../PW/bachel.o \
../PW/functionals.o \
../PW/lsda_functionals.o
../PW/bachel.o
all: ld1.x

3
flib/Makefile
View File

@ -9,6 +9,9 @@ atomic_number.o \
capital.o \
dost.o \
erf.o \
functionals.o \
lsda_functionals.o \
more_functionals.o \
iceil.o \
iglocal.o \
invmat.o \

0
PW/functionals.f90 → flib/functionals.f90
View File

2
PW/lsda_functionals.f90 → flib/lsda_functionals.f90
View File

@ -58,6 +58,8 @@ subroutine xc_spin (rho, zeta, ex, ec, vxup, vxdw, vcup, vcdw)
vcdw = 0.0d0
elseif (icorr == 1) then
call pz_spin (rs, zeta, ec, vcup, vcdw)
elseif (icorr == 3) then
call lsd_lyp (rho, zeta, ec, vcup, vcdw) ! from CP/FPMD (more_functionals)
elseif (icorr == 4) then
call pw_spin (rs, zeta, ec, vcup, vcdw)
else

1964
flib/more_functionals.f90 Normal file
View File

File diff suppressed because it is too large Load Diff