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 : MODULE m_vgen_xcpot
7 :
8 : USE m_juDFT
9 : #ifdef CPP_MPI
10 : use mpi
11 : #endif
12 :
13 : CONTAINS
14 :
15 674 : SUBROUTINE vgen_xcpot(hybdat, input, xcpot, atoms, sphhar, stars, vacuum, sym, &
16 : cell, fmpi, noco, den, denRot, EnergyDen, vTot, vx, vxc, exc, results, &
17 : den1Rot, den1Rotimag, dfptvTotimag, starsq)
18 : !! FLAPW potential generator
19 : !! Calculates the density-potential integrals needed for the total energy
20 : !! TE_VCOUL: charge density-coulomb potential integral
21 : !! TE_VEFF : charge density-effective potential integral
22 : !! TE_EXC : charge density-ex-corr.energy density integral
23 : !!
24 : !! DFPT case: Calculate Vxc1 instead of Vxc. For this we need fxc, the xc Kernel.
25 : !! The calculation changes dramatically, so we enter different subroutines for it.
26 : !! TODO: They only work for LDA right now. GGA requires some mathematical considerations
27 : !! and will have to use grad_rho1 as input as well.
28 :
29 : USE m_types
30 : USE m_constants
31 : USE m_intnv
32 : USE m_vmt_xc
33 : USE m_vvac_xc
34 : USE m_vis_xc
35 : USE m_checkdopall
36 : USE m_cdn_io
37 : USE m_convol
38 : USE m_intgr
39 : USE m_metagga
40 : USE m_dfpt_vmt_xc
41 : USE m_dfpt_vis_xc
42 : USE m_dfpt_vvac_xc
43 :
44 : IMPLICIT NONE
45 :
46 : CLASS(t_xcpot), INTENT(IN) :: xcpot
47 : TYPE(t_hybdat), INTENT(IN) :: hybdat
48 : TYPE(t_mpi), INTENT(IN) :: fmpi
49 :
50 : TYPE(t_input), INTENT(IN) :: input
51 : TYPE(t_vacuum), INTENT(IN) :: vacuum
52 : TYPE(t_noco), INTENT(IN) :: noco
53 : TYPE(t_sym), INTENT(IN) :: sym
54 : TYPE(t_stars), INTENT(IN) :: stars
55 : TYPE(t_cell), INTENT(IN) :: cell
56 : TYPE(t_sphhar), INTENT(IN) :: sphhar
57 : TYPE(t_atoms), INTENT(IN) :: atoms
58 : TYPE(t_potden), INTENT(IN) :: den, denRot, EnergyDen
59 : TYPE(t_potden), INTENT(INOUT) :: vTot, vx, vxc, exc
60 : TYPE(t_results), INTENT(INOUT), OPTIONAL :: results
61 : TYPE(t_potden), INTENT(IN), OPTIONAL :: den1Rot, den1Rotimag
62 : TYPE(t_potden), INTENT(INOUT), OPTIONAL :: dfptvTotimag
63 : TYPE(t_stars), INTENT(IN), OPTIONAL :: starsq
64 :
65 : ! Local type instances
66 674 : TYPE(t_potden) :: workDen, veff
67 674 : Type(t_kinED) :: kinED
68 674 : REAL, ALLOCATABLE :: rhoc(:,:,:),rhoc_vx(:)
69 674 : REAL, ALLOCATABLE :: tec(:,:), qintc(:,:)
70 : ! Local Scalars
71 : INTEGER :: ifftd2, ispin, i, iType
72 : REAL :: dpdot
73 : LOGICAL :: l_dfptvgen
74 : #ifdef CPP_MPI
75 : integer:: ierr
76 : #endif
77 :
78 674 : l_dfptvgen = PRESENT(starsq)
79 :
80 : call set_kinED(fmpi, sphhar, atoms, sym, xcpot, &
81 674 : input, noco, stars,vacuum , cell, Den, EnergyDen, vTot,kinED)
82 :
83 674 : IF (PRESENT(results)) THEN
84 674 : CALL veff%init(stars, atoms, sphhar, vacuum, noco, input%jspins, 1)
85 : #ifndef CPP_OLDINTEL
86 2022 : ALLOCATE (veff%pw_w, mold=veff%pw)
87 : #else
88 : ALLOCATE (veff%pw_w(size(veff%pw, 1), size(veff%pw, 2)))
89 : #endif
90 : ENDIF
91 :
92 : ! exchange correlation potential
93 :
94 : ! vacuum region
95 674 : IF (fmpi%irank == 0) THEN
96 337 : IF (input%film) THEN
97 37 : CALL timestart("Vxc in vacuum")
98 :
99 37 : ifftd2 = 9*stars%mx1*stars%mx2
100 :
101 37 : IF (.NOT. l_dfptvgen) THEN
102 37 : CALL vvac_xc(ifftd2, stars, vacuum, noco, cell, xcpot, input, Den, vTot, exc)
103 : ELSE
104 0 : CALL dfpt_vvac_xc(ifftd2, stars, starsq, vacuum, noco, cell,denRot, den1Rot, xcpot, input, vTot)
105 : END IF
106 37 : CALL timestop("Vxc in vacuum")
107 : END IF
108 :
109 : ! interstitial region
110 337 : CALL timestart("Vxc in interstitial")
111 337 : IF (.NOT.l_dfptvgen) THEN
112 337 : CALL vis_xc(stars, sym, cell, den, xcpot, input, noco, EnergyDen,kinED, vTot, vx, exc, vxc)
113 : ELSE
114 : ! TODO: This is different enough to warrant a separate subroutine, right?
115 0 : CALL dfpt_vis_xc(stars, starsq, sym, cell, denRot, den1Rot, xcpot, input, vTot)
116 : END IF
117 337 : CALL timestop("Vxc in interstitial")
118 : END IF !irank==0
119 :
120 : !
121 : ! ------------------------------------------
122 : ! ----> muffin tin spheres region
123 :
124 674 : IF (fmpi%irank == 0) THEN
125 337 : CALL timestart("Vxc in MT")
126 : END IF
127 :
128 674 : IF (.NOT.l_dfptvgen) THEN
129 : CALL vmt_xc(fmpi, sphhar, atoms, den, xcpot, input, sym, &
130 674 : EnergyDen,kinED, noco,vTot, vx, exc, vxc)
131 : ELSE
132 0 : CALL dfpt_vmt_xc(fmpi,sphhar,atoms,denRot,den1Rot,den1Rotimag,xcpot,input,sym,noco,vTot,dfptvTotimag)
133 : END IF
134 :
135 : ! add MT EXX potential to vr
136 674 : IF (fmpi%irank == 0) THEN
137 337 : CALL timestop("Vxc in MT")
138 :
139 : ! check continuity of total potential
140 337 : IF (input%vchk) CALL checkDOPAll(input, sphhar, stars, atoms, sym, vacuum, cell, vTot, 1)
141 :
142 : ! TOTAL
143 337 : IF (PRESENT(results)) THEN
144 : ! CALCULATE THE INTEGRAL OF n1*Veff1 + n2*Veff2
145 : ! Veff = Vcoulomb + Vxc
146 337 : IF (noco%l_noco) THEN
147 96 : workDen = denRot
148 : ELSE
149 241 : workden = den
150 : END IF
151 337 : veff = vTot
152 337 : IF (xcpot%is_hybrid() .AND. hybdat%l_subvxc) THEN
153 111 : DO ispin = 1, input%jspins
154 111 : CALL convol(stars, vx%pw_w(:, ispin), vx%pw(:, ispin))
155 : END DO
156 14414 : veff%pw = vTot%pw - xcpot%get_exchange_weight()*vx%pw
157 14414 : veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight()*vx%pw_w
158 403578 : veff%mt = vTot%mt - xcpot%get_exchange_weight()*vx%mt
159 : END IF
160 :
161 873 : DO ispin = 1, input%jspins
162 1193516 : DO i = 1, stars%ng3
163 1192643 : vx%pw(i, ispin) = vx%pw(i, ispin)/stars%nstr(i)
164 1193179 : vx%pw_w(i, ispin) = vx%pw_w(i, ispin)/stars%nstr(i)
165 : END DO
166 : END DO
167 :
168 337 : results%te_veff = 0.0
169 873 : DO ispin = 1, input%jspins
170 536 : WRITE (oUnit, FMT=8050) ispin
171 : 8050 FORMAT(/, 10x, 'density-effective potential integrals for spin ', i2,/)
172 873 : CALL int_nv(ispin, stars, vacuum, atoms, sphhar, cell, sym, input, veff, workden, results%te_veff)
173 : END DO
174 :
175 337 : IF (xcpot%is_hybrid() .AND. hybdat%l_subvxc) THEN
176 :
177 147 : ALLOCATE(rhoc(atoms%jmtd,atoms%ntype,input%jspins), rhoc_vx(atoms%jmtd))
178 98 : ALLOCATE(tec(atoms%ntype,input%jspins),qintc(atoms%ntype,input%jspins))
179 49 : CALL readCoreDensity(input,atoms,rhoc,tec,qintc)
180 49 : DEALLOCATE(tec,qintc)
181 :
182 111 : DO ispin = 1, input%jspins
183 209 : DO iType = 1,atoms%ntype
184 : rhoc_vx(:atoms%jri(iType)) = rhoc(:atoms%jri(iType),iType,ispin) * &
185 76848 : vx%mt(:atoms%jri(iType),0,iType,ispin) / sfp_const
186 98 : CALL intgr3(rhoc_vx,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),dpdot)
187 160 : results%te_veff = results%te_veff + xcpot%get_exchange_weight() * dpdot*atoms%neq(iType)
188 : END DO
189 : END DO
190 :
191 49 : DEALLOCATE(rhoc,rhoc_vx)
192 :
193 : END IF
194 :
195 337 : WRITE (oUnit, FMT=8060) results%te_veff
196 : 8060 FORMAT(/, 10x, 'total density-effective potential integral :', t40, ES20.10)
197 :
198 : ! CALCULATE THE INTEGRAL OF n*exc
199 :
200 : ! perform spin summation of charge densities for the calculation of Exc
201 337 : CALL workden%sum_both_spin()
202 :
203 337 : WRITE (oUnit, FMT=8070)
204 : 8070 FORMAT(/, 10x, 'charge density-energy density integrals',/)
205 :
206 337 : results%te_exc = 0.0
207 337 : CALL int_nv(1, stars, vacuum, atoms, sphhar, cell, sym, input, exc, workDen, results%te_exc)
208 337 : WRITE (oUnit, FMT=8080) results%te_exc
209 :
210 : 8080 FORMAT(/, 10x, 'total charge density-energy density integral :', t40, ES20.10)
211 : END IF
212 : END IF ! fmpi%irank == 0
213 :
214 674 : END SUBROUTINE vgen_xcpot
215 :
216 : END MODULE m_vgen_xcpot
|
