Line data Source code
1 : module m_slater
2 :
3 : use m_types
4 : use m_constants
5 : use m_juDFT
6 : use m_differ
7 : use m_intgr
8 : use m_xmlOutput
9 :
10 : implicit none
11 :
12 : contains
13 144 : subroutine slater(input, jspin, atoms, vr, l_write, slater_parameters)
14 :
15 : ! Calculate the slater integrals for the orbitals specified
16 : ! to be corrected with LDA+OP
17 : !
18 : ! Adapted from core/cored.F90
19 : ! Author: Henning Janssen 2021
20 :
21 : type(t_input), intent(in) :: input
22 : type(t_atoms), intent(in) :: atoms
23 : integer, intent(in) :: jspin
24 : real, intent(in) :: vr(:, :)
25 : logical, optional, intent(in) :: l_write
26 : real, allocatable, optional, intent(out) :: slater_parameters(:,:)
27 :
28 : real :: eig, fj, fl, fn, t2
29 : real :: d, dxx, rn, rnot, z, t1, rr, r, c
30 : integer :: i, n, ncmsh, l, ierr, kk
31 : integer :: i_opc, atomType, ipm
32 :
33 144 : real :: vrd(atoms%msh), f(0:3, 2, atoms%n_opc)
34 144 : real :: a(atoms%msh), b(atoms%msh), ain(atoms%msh), ahelp(atoms%msh)
35 : real :: lambda
36 : character(LEN=20) :: attributes(6)
37 :
38 : logical :: l_write_arg
39 :
40 144 : l_write_arg = .false.
41 144 : if (present(l_write)) l_write_arg = l_write
42 :
43 144 : if (l_write_arg) call openXMLElementNoAttributes('slaterIntegrals')
44 :
45 144 : c = c_light(1.0)
46 :
47 432 : do i_opc = 1, atoms%n_opc
48 :
49 288 : l = atoms%lda_opc(i_opc)%l
50 288 : n = atoms%lda_opc(i_opc)%n
51 288 : atomType = atoms%lda_opc(i_opc)%atomType
52 288 : lambda = 0.0 !Set for OPC since it does not matter there
53 :
54 288 : z = atoms%zatom(atomType)
55 288 : dxx = atoms%dx(atomType)
56 288 : rnot = atoms%rmsh(1, atomType)
57 288 : d = exp(atoms%dx(atomType))
58 288 : ncmsh = nint(log((atoms%rmt(atomType) + 10.0)/rnot)/dxx + 1)
59 288 : ncmsh = min(ncmsh, atoms%msh)
60 288 : rn = rnot*(d**(ncmsh - 1))
61 288 : if (l_write_arg) write (oUnit, fmt=8000) z, rnot, dxx, atoms%jri(atomType)
62 215424 : vrd(:atoms%jri(atomType)) = vr(:atoms%jri(atomType), atomType)
63 :
64 : !
65 288 : if (input%l_core_confpot) then
66 : !---> linear extension of the potential with slope t1 / a.u.
67 288 : t1 = 0.125
68 : t1 = max((vrd(atoms%jri(atomType)) - vrd(atoms%jri(atomType) - 1)*d)* &
69 288 : d/(atoms%rmt(atomType)**2*(d - 1)), t1)
70 288 : t2 = vrd(atoms%jri(atomType))/atoms%rmt(atomType) - atoms%rmt(atomType)*t1
71 288 : rr = atoms%rmt(atomType)
72 : else
73 0 : t2 = vrd(atoms%jri(atomType))/(atoms%jri(atomType) - ncmsh)
74 : end if
75 288 : if (atoms%jri(atomType) < ncmsh) then
76 30816 : do i = atoms%jri(atomType) + 1, ncmsh
77 30816 : if (input%l_core_confpot) then
78 30528 : rr = d*rr
79 30528 : vrd(i) = rr*(t2 + rr*t1)
80 : ! vrd(i) = 2*vrd(jri(jatom)) - rr*( t2 + rr*t1 )
81 : else
82 0 : vrd(i) = vrd(atoms%jri(atomType)) + t2*(i - atoms%jri(atomType))
83 : end if
84 : !
85 : end do
86 : end if
87 :
88 288 : fl = real(l)
89 288 : fn = real(n)
90 864 : do ipm = 1, 2
91 576 : fj = fl + (ipm - 1.5)
92 :
93 576 : eig = -2*(z/(fn + fl))**2
94 576 : call differ(fn, fl, fj, c, z, dxx, rnot, rn, d, ncmsh, vrd, eig, a, b, ierr)
95 576 : if (l_write_arg) write (oUnit, fmt=8010) fn, fl, fj, eig
96 :
97 576 : if (ierr /= 0) call juDFT_error("error in slater routine", calledby="slater")
98 :
99 3168 : do kk = 0, 2*l, 2 ! F0 for s, F0 + F2 for p etc.
100 : !lambda = 3.5449*sqrt((l+1.0)/100) ! screening (TF) sqrt(4pi N(ef))
101 : !IF (kk.GT.0) lambda = 2*lambda
102 1728 : r = rnot
103 1475712 : do i = 1, ncmsh
104 1473984 : ain(i) = a(i)**2*r**(-kk - 1) ! prepare inner integrand
105 : if (kk == 0) then
106 : ain(i) = ain(i)*exp(-r*lambda)
107 : end if
108 1475712 : r = r*d
109 : end do
110 1728 : call intgr1(ain, rnot, dxx, ncmsh, ahelp)
111 :
112 1728 : r = rnot
113 1473984 : do i = 1, ncmsh - 1
114 1472256 : ain(i) = a(i)**2*r**kk*(ahelp(ncmsh) - ahelp(i))
115 : if (kk == 0) then
116 : ain(i) = ain(i)*exp(r*lambda)
117 : end if
118 1473984 : r = r*d
119 : end do
120 2304 : call intgr0(ain, rnot, dxx, ncmsh - 1, f(kk/2, ipm, i_opc))
121 : end do
122 : end do
123 2592 : f(0:l, :, i_opc) = f(0:l, :, i_opc)*2
124 288 : if (l_write_arg) then
125 144 : write (oUnit, fmt=8020) n, l, jspin, lambda
126 576 : write (oUnit, '(12x,i3,f7.1,4f20.10)') l, fj, (f(kk, 1, i_opc)*hartree_to_ev_const, kk=0, l)
127 576 : write (oUnit, '(12x,i3,f7.1,4f20.10)') l, fj, (f(kk, 2, i_opc)*hartree_to_ev_const, kk=0, l)
128 : end if
129 :
130 144 : if (l_write_arg) then
131 1008 : attributes = ''
132 144 : write (attributes(1), '(i0)') atomType
133 144 : write (attributes(2), '(i0)') nint(z)
134 144 : write (attributes(3), '(i0)') jspin
135 144 : write (attributes(4), '(i0)') l
136 144 : write (attributes(5), '(i0)') n
137 144 : write (attributes(6), '(f10.7)') lambda
138 : call openXMLElementForm('slaterIntegral', (/'atomType ', 'atomicNumber ', 'spin ', 'l ', &
139 : 'n ', 'screening '/), &
140 1008 : attributes, reshape((/8, 12, 4, 1, 1, 9, 6, 3, 1, 1, 1, 18/), (/6, 2/)))
141 432 : do ipm = 1, 2
142 288 : fj = fl + (ipm - 1.5)
143 2016 : attributes = ''
144 288 : write (attributes(1), '(i0)') n
145 288 : write (attributes(2), '(i0)') nint(fl)
146 288 : write (attributes(3), '(f4.1)') fj
147 288 : write (attributes(4), '(a)') 'eV'
148 : call writeXMLElementPoly('state', (/'n ', 'l ', 'j ','units'/), &
149 2736 : attributes(1:4), contentList=f(0:l, ipm, i_opc)*hartree_to_ev_const)
150 : end do
151 144 : call closeXMLElement('slaterIntegral')
152 : end if
153 : end do
154 :
155 144 : if (l_write_arg) call closeXMLElement('slaterIntegrals')
156 :
157 144 : if (present(slater_parameters)) then
158 1248 : allocate(slater_parameters(0:lmaxU_const,atoms%n_opc), source=0.0)
159 288 : do i_opc = 1, atoms%n_opc
160 192 : l = atoms%lda_opc(i_opc)%l
161 : !Is simply averaging the right approach??
162 864 : slater_parameters(0:l, i_opc) = (f(0:l,1,i_opc) + f(0:l,2,i_opc))/2.0
163 : enddo
164 : endif
165 :
166 : 8000 format(/, /, 10x, 'z=', f4.0, 5x, 'r(1)=', e14.6, 5x, 'dx=', f9.6, 5x,&
167 : & 'm.t.index=', i4, /, 15x, 'n', 4x, 'l', 5x, 'j', 4x, 'energy')
168 : 8010 format(12x, 2f5.0, f6.1, f10.4, f12.4)
169 : 8020 format(/, /, 10x, 'Slater integrals: n=', i2, 5x, 'l=', i2, 5x, 'spin=', i2, 5x, 'screening=', f9.6, 5x,&
170 : & /, 14x, 'l', 4x, 'j', 9x, 'F(0,2,4,6)')
171 144 : end subroutine slater
172 : end module m_slater
|
