Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_xcwgn
8 : !-----------------------------------------------------------------------
9 : ! Called in case of icorr=1 : non-spinpolarized exchange-correlation
10 : ! from wigners interpolation formula.
11 :
12 : ! krla=1: Relativistic correction of exchange energy and potential
13 : ! related to Dirac kinetic energy, according to:
14 : ! A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
15 :
16 : ! be careful: calculation in rydberg!
17 :
18 : ! vxc calculates the XC-potential and
19 : ! exc calculates the XC-energy
20 :
21 : ! based on a subroutine by S. Bluegel; r.pentcheva 22.01.96
22 : !-----------------------------------------------------------------------
23 :
24 : USE m_constants, ONLY : pi_const
25 : USE m_relcor
26 : IMPLICIT NONE
27 :
28 : REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0
29 : REAL, PARAMETER, PRIVATE :: thrd = one/three , d_15 = 1.e-15
30 : REAL, PARAMETER, PRIVATE :: cex = 0.91633058742 ! 3/2 * ( 3/(2*pi) )^(2/3)
31 :
32 : CONTAINS
33 : !************************************************************************
34 0 : SUBROUTINE vxcwgn( &
35 : krla,jspins, &
36 0 : mgrid,ngrid,rh, &
37 0 : vx,vxc)
38 : !************************************************************************
39 :
40 : ! .. Scalar Arguments ..
41 : INTEGER, INTENT (IN) :: jspins
42 : INTEGER, INTENT (IN) :: krla ! run mode parameters
43 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
44 :
45 : ! .. Array Arguments ..
46 : REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
47 : REAL, INTENT (OUT) :: vx(mgrid,jspins) ! x potential
48 : REAL, INTENT (OUT) :: vxc(mgrid,jspins) ! xc potential
49 :
50 : ! .. Local Scalars ..
51 : REAL :: fothrd,vxp,vcp
52 : REAL :: rs, rho, thfpi, exp, ecp
53 : INTEGER :: ir
54 :
55 : ! .. Local Arrays ..
56 0 : REAL, ALLOCATABLE :: psi(:) ! relativistic exchange potential corr.
57 :
58 : !-----s Intrinsic Functions
59 : INTRINSIC alog,max
60 :
61 0 : thfpi = three / ( four * pi_const )
62 0 : fothrd = four/three
63 :
64 : !-----> evaluate relativistic corrections for exchange
65 :
66 0 : ALLOCATE ( psi(ngrid) )
67 : CALL relcor( &
68 : mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
69 0 : psi)
70 :
71 0 : DO ir = 1,ngrid
72 0 : IF (jspins == 1) THEN
73 0 : rho = max(d_15,rh(ir,1))
74 : ELSE
75 0 : rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
76 : ENDIF
77 0 : rs = (thfpi/rho)**thrd
78 0 : exp = - cex / rs ! exchange energy = -0.9163306/rs
79 0 : ecp = - 0.88 / (rs + 7.8) ! correlation energy = -0.88/(rs + 7.8)
80 :
81 0 : vxp = fothrd * exp ! exchange potential = 4/3 exp
82 0 : vcp = fothrd * ecp + 2.288 / ( rs + 7.8 ) ** 2
83 :
84 0 : vxc(ir,1) = vcp + vxp * psi(ir)
85 :
86 0 : vx(ir,1) = vxp * psi(ir)
87 : ENDDO
88 0 : IF ( jspins == 2 ) THEN ! spinpolarized calculation
89 0 : DO ir = 1,ngrid
90 0 : vxc(ir,jspins) = vxc(ir,1)
91 :
92 0 : vx(ir,jspins) = vx(ir,1)
93 : ENDDO
94 : ENDIF
95 :
96 0 : DEALLOCATE (psi)
97 0 : RETURN
98 :
99 : END SUBROUTINE vxcwgn
100 : !***********************************************************************
101 0 : SUBROUTINE excwgn( &
102 : iofile,krla,jspins, &
103 0 : mgrid,ngrid,rh, &
104 0 : exc)
105 : !***********************************************************************
106 :
107 : ! .. Scalar Arguments ..
108 : INTEGER, INTENT (IN) :: jspins
109 : INTEGER, INTENT (IN) :: krla ! run mode parameters
110 : INTEGER, INTENT (IN) :: iofile ! file number for read and write
111 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
112 :
113 : ! .. Array Arguments ..
114 : REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
115 : REAL, INTENT (OUT) :: exc(mgrid) ! xc energy
116 :
117 : ! .. Local Arrays ..
118 0 : REAL, ALLOCATABLE :: phi(:) ! relativistic exchange energy correct.
119 :
120 : REAL :: rs, rho, thfpi, exp, ecp
121 : INTEGER :: ir
122 :
123 0 : thfpi = three / ( four * pi_const )
124 :
125 0 : ALLOCATE ( phi(ngrid) )
126 : CALL relcor( &
127 : mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
128 0 : phi)
129 :
130 0 : DO ir = 1,ngrid
131 0 : IF (jspins == 1) THEN
132 0 : rho = max(d_15,rh(ir,1))
133 : ELSE
134 0 : rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
135 : ENDIF
136 0 : rs = (thfpi/rho)**thrd
137 0 : exp = - cex / rs ! exchange energy = -0.9163306/rs
138 0 : ecp = - 0.88 / (rs + 7.8) ! correlation energy = -0.88/(rs + 7.8)
139 :
140 0 : exc(ir) = ecp + exp * phi(ir)
141 : ENDDO
142 0 : IF (jspins == 2) THEN
143 0 : WRITE (iofile,'('' WARNING: Wigner correlation ! Only applicable for non-spinpolarized calculations'')')
144 : ENDIF
145 :
146 0 : DEALLOCATE (phi)
147 0 : RETURN
148 :
149 : END SUBROUTINE excwgn
150 : END MODULE m_xcwgn
|
