Branch data Line data Source code
1 : : !!****m* ABINIT/m_evdw_wannier
2 : : !! NAME
3 : : !! m_evdw_wannier
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2010-2026 ABINIT group (CE, TR, AR)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_evdw_wannier
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : :
28 : : use m_special_funcs, only : abi_derf
29 : : use m_numeric_tools, only : simpson_int
30 : : use m_geometry, only : xcart2xred, xred2xcart
31 : :
32 : : implicit none
33 : :
34 : : private
35 : : !!***
36 : :
37 : : public :: evdw_wannier
38 : : !!***
39 : :
40 : : contains
41 : : !!***
42 : :
43 : : !!****f* ABINIT/evdw_wannier
44 : : !! NAME
45 : : !! evdw_wannier
46 : : !!
47 : : !! FUNCTION
48 : : !! FIXME: Evaluates the van der Waals correlation energy using maximally
49 : : !! localized Wannier functions (MLWF) as proposed by:
50 : : !! P. L. Silvestrelli in PRL 100:053002 (2008) [[cite:Sivestrelli2008]] vdw_xc=10 and
51 : : !! A. Ambrosetti and P. L. Silvestrelli in PRB 85:073101 (2012) [[cite:Ambrosetti2012]] vdw_xc=11.
52 : : !! P. L. Silvestrelli in J.Chem.Phys. 139:054106 (2013) [[cite:Silvestrelli2013]] vdw_xc=14.
53 : : !!
54 : : !! INPUTS
55 : : !! nsppol = Spin polarization.
56 : : !! nwan(nsppol) = Total number of MLWF in the system per spin component.
57 : : !! origmwan = max[nwan(nsppol)] from mlwfovlp.F90.
58 : : !! tdocc_wan = MLWFs occupation matrix diagonal terms
59 : : !! vdw_nfrag = Number of vdW interating fragments in the unit cell.
60 : : !! vdw_supercell(3) = Distance along each rprimd components for
61 : : !! which vdW interactions between MLWF will be taken into account.
62 : : !! vdw_typfrag(natom) = Fragment to which each atom belongs to.
63 : : !! vdw_xc = vdW-WF version.
64 : : !! rprimd = Real space primitive translations.
65 : : !! wann_centres(3,origmwan,nsppol) = The centers of MLWFs in a.u.
66 : : !! wann_spreads(origmwan,nsppol) = Spread of the MLWFs, in Ang**2. (from wannier90).
67 : : !! xcart = Coordinates of unit cell atoms in atomic units.
68 : : !!
69 : : !! OUTPUT
70 : : !! csix(origmwan,origmwan,nsppol,nsppol) = dispersion coefficient between each pair of MLWF.
71 : : !! corrvdw = van der Waals correction to the energy.
72 : : !!
73 : : !! SOURCE
74 : :
75 : 3 : subroutine evdw_wannier(csix,corrvdw,origmwan,natom,nsppol,orignwan,tdocc_wan,vdw_nfrag,&
76 : 3 : & vdw_supercell,vdw_typfrag,vdw_xc,rprimd,wann_centres,wann_spreads,xcart)
77 : :
78 : : !Arguments ------------------------------------
79 : : integer , intent(in) :: origmwan,nsppol,natom,orignwan(nsppol)
80 : : integer , intent(in) :: vdw_nfrag,vdw_supercell(3),vdw_typfrag(natom),vdw_xc
81 : : real(dp), intent(in) :: rprimd(3,3),wann_centres(3,origmwan,nsppol),wann_spreads(origmwan,nsppol)
82 : : real(dp), intent(in) :: xcart(3,natom)
83 : : real(dp), intent(out) :: corrvdw
84 : : real(dp), intent(out) :: csix(origmwan,origmwan,nsppol,nsppol)
85 : : real(dp), intent(in) :: tdocc_wan(origmwan,nsppol)
86 : :
87 : : !Local variables-------------------------------
88 : : integer :: ier,igr,icx,icy,icz,ii,inx,iny,inz,isppol,iwan
89 : : integer :: jwan,jj,ll,mm,mwan,nc,ngr,tmp_mwan,mwan_half
90 : 3 : integer, allocatable:: amagr(:,:,:),inwan(:,:),nw(:,:),nwan(:),npwf(:),ord(:,:)
91 : 3 : integer, allocatable:: tmp_nwan(:)
92 : : real(dp) :: dnrm2,fij,rij,rij_c(3),fu,shift,erfValue
93 : : real(dp), parameter :: a = 20.d0 !Parameter related to the damping function.
94 : : real(dp), parameter :: gama = 4.5d0/(sqrt3**3) !alpha=gama*S**3.
95 : : real(dp), parameter :: gama1 = 0.88d0 !alpha=gama*S**3.
96 : : real(dp), parameter :: zeta = 1.30d0 !polar=zeta*(Z/omega**2).
97 : : real(dp), parameter :: beta = 1.39d0 ! .
98 : 3 : real(dp), allocatable:: amawf(:,:),amaspr(:),amaocc(:),dcenters(:,:,:),rc(:,:)
99 : 3 : real(dp), allocatable:: tmp_cent(:,:,:),tmp_spr(:,:),tmp_occ(:,:)
100 : 3 : real(dp), allocatable:: rv(:,:),wanncent(:,:,:),wannspr(:,:),wc_rec(:,:,:),xi(:,:)
101 : 3 : real(dp), allocatable:: c_QHO(:,:),Tij_dip(:,:),polar(:),omega(:),eigv(:),zhpev2(:)
102 : 3 : real(dp), allocatable :: newocc_wan(:,:)
103 : 3 : complex(dp), allocatable :: eigvec(:,:),matrx(:),zhpev1(:)
104 : : character(len=500) :: message ! to be uncommented, if needed
105 : : ! *************************************************************************
106 : :
107 : : !Determining presence p-like MLWFs see J.Chem.Phys.135:154105 (2011) [[cite:Andrinopoulos2011]]
108 [ - + + - : 9 : ABI_MALLOC(npwf,(nsppol))
+ - + - -
+ - + ]
109 [ - + + - : 12 : ABI_MALLOC(inwan,(origmwan,nsppol))
- + + - +
- + - - +
- + ]
110 [ - + + - : 6 : ABI_MALLOC(nwan,(nsppol))
- + - + ]
111 : :
112 : 3 : ll = 0
113 [ + + ]: 6 : npwf(:) = 0
114 [ + + + + ]: 18 : inwan(:,:) = 0
115 [ + + ]: 6 : do jj=1,nsppol
116 [ + + ]: 18 : do iwan=1,orignwan(jj)
117 [ - + ]: 15 : if(tdocc_wan(iwan,jj)*nsppol<=1.50d0) then
118 : 0 : npwf(jj) = npwf(jj) + 1
119 : 0 : ll = ll+1
120 : 0 : inwan(ll,jj) = iwan
121 : : end if
122 : : end do
123 : : end do
124 : :
125 : 3 : write(std_out,*) ch10,'Number of p-like MLWFs per spin pol:',ch10
126 : 3 : write(std_out,*) (npwf(ii),ii=1,nsppol), ch10
127 : :
128 [ + + ]: 6 : mwan=origmwan+(sum(npwf(:))) !two new MLWFs per p-like MLWF
129 [ + + ]: 6 : nwan(:)=orignwan(:)+npwf(:)
130 : :
131 : :
132 [ - + + - : 12 : ABI_MALLOC(wanncent,(3,mwan,nsppol))
- + + - +
- + - - +
- + ]
133 [ - + + - : 12 : ABI_MALLOC(wannspr,(mwan,nsppol))
- + + - +
- + - - +
- + ]
134 [ - + + - : 9 : ABI_MALLOC(wc_rec,(3,mwan,nsppol))
- + + - -
+ - + ]
135 [ - + + - : 12 : ABI_MALLOC(ord,(mwan,nsppol))
- + + - +
- + - - +
- + ]
136 [ - + + - : 9 : ABI_MALLOC(newocc_wan,(mwan,nsppol))
- + + - -
+ - + ]
137 : :
138 [ + + + + : 54 : wanncent(:,:,:) = zero
+ + ]
139 [ + + + + ]: 18 : wannspr(:,:) = zero
140 [ + + + + : 54 : wc_rec(:,:,:) = zero
+ + ]
141 [ + + + + ]: 18 : newocc_wan(:,:) = zero
142 [ + + + + ]: 18 : ord(:,:) = zero
143 : :
144 : : !The vdW correction is calculated in atomic units:
145 [ + + ]: 6 : do ii=1,nsppol
146 [ + + ]: 18 : do iwan=1,orignwan(ii)
147 : : ! converting to bohr**2 and then squared
148 [ + + ]: 48 : wanncent(:,iwan,ii)=wann_centres(:,iwan,ii)/Bohr_Ang
149 : : ! write(std_out,*) "spread of WF",i, "=", wann_spreads(i)
150 : 12 : wannspr(iwan,ii)=sqrt(wann_spreads(iwan,ii)/Bohr_Ang**2)
151 : 15 : newocc_wan(iwan,ii)=tdocc_wan(iwan,ii)
152 : : end do
153 : : end do
154 : :
155 : : !write(std_out,*) 'Number of MLWFs:',ch10
156 : : !do ii=1,nsppol
157 : : !write(std_out,*) 'nsppol=',ii, 'nwan(nsppol)=',nwan(nsppol),ch10
158 : : !end do
159 : :
160 : 3 : write(std_out,*) 'Original Wannier centres and spreads:',ch10
161 [ + + ]: 6 : do ii=1,nsppol
162 : 3 : write(std_out,*) 'nsppol=',ii,ch10
163 [ + + ]: 18 : do iwan=1,orignwan(ii)
164 : 15 : write(std_out,*) (wanncent(jj,iwan,ii),jj=1,3), wannspr(iwan,ii),ch10
165 : : end do
166 : : end do
167 : :
168 : : !Translate MLWFs to the original unit cell if vdw_nfrag > 0 :
169 : :
170 [ - + ]: 3 : if(vdw_nfrag>0)then
171 [ # # ]: 0 : do jj=1,nsppol
172 : : call xcart2xred(orignwan(jj),rprimd,wanncent(:,1:orignwan(jj),jj), &
173 : 0 : & wc_rec(:,1:orignwan(jj),jj))
174 : : ! got centers in reduced coor
175 [ # # ]: 0 : do iwan=1,orignwan(jj)
176 [ # # ]: 0 : do ii=1,3
177 [ # # ]: 0 : if(wc_rec(ii,iwan,jj)<zero) then
178 [ # # ]: 0 : shift=REAL(CEILING(ABS(wc_rec(ii,iwan,jj))),dp)
179 : 0 : wc_rec(ii,iwan,jj) = wc_rec(ii,iwan,jj)+shift
180 : : end if
181 [ # # ]: 0 : if(wc_rec(ii,iwan,jj)>one) then
182 : 0 : shift=-REAL(INT(wc_rec(ii,iwan,jj)),dp)
183 : 0 : wc_rec(ii,iwan,jj) = wc_rec(ii,iwan,jj)+shift
184 : : end if
185 : : end do
186 : : end do
187 : : call xred2xcart(orignwan(jj),rprimd,wanncent(:,1:orignwan(jj),jj), &
188 : 0 : & wc_rec(:,1:orignwan(jj),jj))
189 : : end do
190 : :
191 : : ! ====================================================================
192 : :
193 : 0 : write(std_out,*) ch10,'Wannier centres translated to unit cell and spr:',ch10
194 [ # # ]: 0 : do jj=1,nsppol
195 : 0 : write(std_out,*) 'nsppol=',jj,ch10
196 [ # # ]: 0 : do iwan=1,orignwan(jj)
197 : 0 : write(std_out,*) (wanncent(ii,iwan,jj),ii=1,3), wannspr(iwan,jj)
198 : : end do
199 : : end do
200 : : end if !vdw_nfrag>0
201 : :
202 : : !Spliting of p-like into 2 s-like MLWFs
203 : : !Eqs. (22) and (23) of J.Chem.Phys.135:154105 (2011) [[cite:Andrinopoulos2011]]
204 : :
205 [ + + + - : 6 : if ( any (npwf(:)/=0) ) then
- + ]
206 : :
207 : 0 : write(std_out,*) 'Indexes of p-like MLWFs and its spin:'
208 : :
209 [ # # ]: 0 : do isppol=1,nsppol
210 [ # # ]: 0 : do jj=1,npwf(isppol)
211 : :
212 : 0 : write(std_out,*) inwan(jj,isppol),isppol
213 : :
214 [ # # ]: 0 : wanncent(1:2,orignwan(isppol)+jj,isppol) = wanncent(1:2,inwan(jj,isppol),isppol)
215 : :
216 : : wanncent(3,orignwan(isppol)+jj,isppol) = wanncent(3,inwan(jj,isppol),isppol) &
217 : 0 : & + 15.d0*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt(30.d0))
218 : :
219 : : wanncent(3,inwan(jj,isppol),isppol) = wanncent(3,inwan(jj,isppol),isppol) &
220 : 0 : & - 15.d0*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt(30.d0))
221 : :
222 : 0 : wannspr(orignwan(isppol)+jj,isppol) = seven*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt2)
223 : :
224 : 0 : wannspr(inwan(jj,isppol),isppol) = seven*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt2)
225 : :
226 : 0 : newocc_wan(orignwan(isppol)+jj,isppol) = tdocc_wan(inwan(jj,isppol),isppol) / two
227 : :
228 : 0 : newocc_wan(inwan(jj,isppol),isppol) = tdocc_wan(inwan(jj,isppol),isppol) / two
229 : :
230 : : end do
231 : : end do
232 : :
233 : 0 : write(std_out,*) ch10,'Wannier centres and spreads after splitting of p-like MLWFs:',ch10
234 [ # # ]: 0 : do isppol=1,nsppol
235 : 0 : write(std_out,*) 'nsppol=',isppol,ch10
236 [ # # ]: 0 : do iwan=1,nwan(isppol)
237 : 0 : write(std_out,*) (wanncent(jj,iwan,isppol),jj=1,3), wannspr(iwan,isppol)
238 : : end do
239 : : end do
240 : :
241 : : end if ! any(npwf(:)/=0)
242 : :
243 : : !Asign each MLWFs to one fragment, the same as their nearest atom:
244 : :
245 : 3 : call order_wannier(mwan,natom,nwan,nsppol,ord,vdw_typfrag,wanncent,xcart)
246 : :
247 : 3 : write(std_out,*) ch10,'Wannier centres and fragments',ch10
248 [ + + ]: 6 : do ll=1,abs(vdw_nfrag)
249 : 3 : write(std_out,*) 'MLWF centers in fragment',ll,ch10
250 [ + + ]: 9 : do jj=1,nsppol
251 [ + + ]: 18 : do iwan=1,nwan(jj)
252 [ + - ]: 15 : if (ord(iwan,jj)==ll) then
253 [ + - + + ]: 48 : write(std_out,*) 'X', (Bohr_Ang*wanncent(ii,iwan,jj),ii=1,3),iwan,jj
254 : : end if
255 : : end do
256 : : end do
257 : : end do
258 : :
259 : 3 : write(std_out,*) ch10,'Occupation Matrix diagonal terms:',ch10
260 [ + + ]: 6 : do ll=1,abs(vdw_nfrag)
261 : 3 : write(std_out,*) 'For MLWF centers in fragment',ll,ch10
262 [ + + ]: 9 : do jj=1,nsppol
263 [ + + ]: 18 : do iwan=1,nwan(jj)
264 [ + - ]: 15 : if (ord(iwan,jj)==ll) then
265 : 12 : write(std_out,*) newocc_wan(iwan,jj),ch10
266 : : end if
267 : : end do
268 : : end do
269 : : end do
270 : :
271 : : !Amalgamation of close MLWFs, see J.Chem.Phys.135:154105 (2011) [[cite:Andrinopoulos2011]]
272 : :
273 [ + + + - : 6 : if (all(npwf(:)==0).and.vdw_xc/=14) then !amalgamation is done only if no p-like
+ - + + ]
274 : :
275 : 2 : mwan_half=mwan/2
276 [ - + + - : 10 : ABI_MALLOC(amagr,(mwan,nsppol,mwan_half))
- + + - -
+ + - + -
+ - - + -
+ ]
277 [ - + + - : 8 : ABI_MALLOC(nw,(nsppol,mwan_half))
- + + - +
- + - - +
- + ]
278 [ + + + + ]: 10 : nw=0
279 : :
280 : 2 : call amalgam(amagr,ngr,nsppol,nw,mwan,ord,nwan,vdw_nfrag,wanncent,wannspr)
281 : :
282 : : ! Calculating amalgamated centres, spreads and occupancies if any:
283 : :
284 [ + + + + : 10 : if( any(nw(:,:) /= 0) ) then
+ - - + ]
285 : :
286 [ # # # # : 0 : ABI_MALLOC(amawf,(3,ngr))
# # # # #
# # # ]
287 [ # # # # : 0 : ABI_MALLOC(amaspr,(ngr))
# # # # #
# # # ]
288 [ # # # # : 0 : ABI_MALLOC(amaocc,(ngr))
# # # # ]
289 : :
290 [ # # # # ]: 0 : amawf(:,:) = 0
291 [ # # ]: 0 : amaspr(:) = 0
292 [ # # ]: 0 : amaocc(:) = 0
293 : :
294 [ # # ]: 0 : do igr = 1 , ngr
295 [ # # ]: 0 : do isppol = 1 , nsppol
296 [ # # ]: 0 : do ii = 1 , nw(isppol,igr)
297 : :
298 [ # # ]: 0 : amawf(:,igr) = amawf(:,igr) + wanncent(:,amagr(ii,isppol,igr),isppol)
299 : 0 : amaspr(igr) = amaspr(igr) + wannspr(amagr(ii,isppol,igr),isppol)
300 : 0 : amaocc(igr) = amaocc(igr) + newocc_wan(amagr(ii,isppol,igr),isppol)
301 : :
302 : : end do
303 : : end do
304 : :
305 [ # # # # ]: 0 : amawf(:,igr) = amawf(:,igr) / real(sum(nw(1:nsppol,igr)),dp )
306 [ # # ]: 0 : amaspr(igr) = amaspr(igr) / real(sum(nw(1:nsppol,igr)),dp )
307 : :
308 : : end do
309 : :
310 : 0 : write(std_out,*) ch10,'Amalgamated MLWFs Centres, Spreads and Occupancies:',ch10
311 [ # # ]: 0 : do igr = 1 , ngr
312 : 0 : write(std_out,*) (amawf(ii,igr),ii=1,3),amaspr(igr),amaocc(igr)
313 : : end do
314 : :
315 : : ! Redefining centres, spreads and occps arrays:
316 [ # # # # : 0 : ABI_MALLOC(tmp_nwan,(nsppol))
# # # # ]
317 : :
318 [ # # # # : 0 : tmp_nwan(:) = nwan(:) - sum(nw(:,1:ngr))
# # ]
319 [ # # ]: 0 : tmp_mwan = maxval(tmp_nwan(:))
320 : :
321 [ # # # # : 0 : ABI_MALLOC(tmp_cent,(3,tmp_mwan,nsppol))
# # # # #
# # # # #
# # ]
322 [ # # # # : 0 : ABI_MALLOC(tmp_spr,(tmp_mwan,nsppol))
# # # # #
# # # # #
# # ]
323 [ # # # # : 0 : ABI_MALLOC(tmp_occ,(tmp_mwan,nsppol))
# # # # #
# # # ]
324 : :
325 [ # # # # : 0 : tmp_cent(:,:,:) = zero
# # ]
326 [ # # # # ]: 0 : tmp_spr(:,:) = zero
327 [ # # # # ]: 0 : tmp_occ(:,:) = zero
328 : :
329 [ # # ]: 0 : do isppol = 1 , nsppol
330 : 0 : ii = 0
331 [ # # ]: 0 : do iwan = 1 , nwan(isppol)
332 : :
333 [ # # # # : 0 : if ( any(amagr(:,isppol,:) == iwan) ) cycle
# # # # ]
334 : :
335 : 0 : ii = ii + 1
336 [ # # ]: 0 : tmp_cent(:,ii,isppol) = wanncent(:,iwan,isppol)
337 : 0 : tmp_spr(ii,isppol) = wannspr(iwan,isppol)
338 : 0 : tmp_occ(ii,isppol) = newocc_wan(iwan,isppol)
339 : :
340 : : end do
341 : : end do
342 : :
343 : : ! Redefining wanncent, wannspr, newocc_wan:
344 : : ! Even if amalgamation occurs with MLWFs of different spins
345 : : ! the new WF are gathered with isppol=1 functions...
346 : :
347 [ # # ]: 0 : nwan(1) = nwan(1) - sum(nw(1,1:ngr)) + ngr
348 : :
349 [ # # ]: 0 : if (nsppol == 2) then
350 [ # # ]: 0 : nwan(2) = nwan(2) - sum(nw(2,1:ngr))
351 : : end if
352 : :
353 [ # # ]: 0 : mwan = maxval(nwan(:))
354 : :
355 [ # # ]: 0 : do isppol = 1 , nsppol
356 [ # # ]: 0 : do iwan = 1 , tmp_nwan(isppol)
357 : :
358 [ # # ]: 0 : wanncent(:,iwan,isppol) = tmp_cent(:,iwan,isppol)
359 : 0 : wannspr(iwan,isppol) = tmp_spr(iwan,isppol)
360 : 0 : newocc_wan(iwan,isppol) = tmp_occ(iwan,isppol)
361 : :
362 : : end do
363 : : end do
364 : :
365 [ # # ]: 0 : do igr = 1 , ngr
366 : :
367 [ # # ]: 0 : wanncent(:,tmp_nwan(1)+igr,1) = amawf(:,igr)
368 : 0 : wannspr(tmp_nwan(1)+igr,1) = amaspr(igr)
369 : 0 : newocc_wan(tmp_nwan(1)+igr,1) = amaocc(igr)
370 : :
371 : : end do
372 : :
373 : : ! Ordering again:
374 : : ! Asign each MLWFs to one fragment, the same as their nearest atom:
375 : :
376 : 0 : call order_wannier(mwan,natom,nwan,nsppol,ord,vdw_typfrag,wanncent,xcart)
377 : :
378 : :
379 : 0 : write(std_out,*) ch10,'Full set of Wannier functions and spreads'
380 : 0 : write(std_out,*) 'after both splitting of p-like WFs and amalgamation',ch10
381 : :
382 [ # # ]: 0 : do ll=1,abs(vdw_nfrag)
383 : 0 : write(std_out,*) 'MLWF centers and spreads in fragment',ll,ch10
384 [ # # ]: 0 : do jj=1,nsppol
385 [ # # ]: 0 : do iwan=1,nwan(jj)
386 [ # # ]: 0 : if (ord(iwan,jj)==ll) then
387 [ # # # # ]: 0 : write(std_out,*) 'X', (Bohr_Ang*wanncent(ii,iwan,jj),ii=1,3),Bohr_Ang*wannspr(iwan,jj)
388 : : end if
389 : : end do
390 : : end do
391 : : end do
392 : :
393 : : end if ! any(nw(:,:) /= 0)
394 : : end if ! all((npwf(:)==0).and.vdw_xc/=14)
395 : :
396 : : !vdW-WF VERSION 1
397 : :
398 [ + + ]: 3 : if(vdw_xc==10) then
399 : :
400 [ - + + - : 4 : ABI_MALLOC(dcenters,(3,mwan,nsppol))
- + + - +
- + - - +
- + ]
401 [ - + + - : 4 : ABI_MALLOC(rc,(mwan,nsppol))
- + + - +
- + - - +
- + ]
402 [ - + + - : 3 : ABI_MALLOC(rv,(mwan,nsppol))
- + + - -
+ - + ]
403 : : ! Calculate intermediate quantities
404 [ + + ]: 2 : do jj=1,nsppol
405 [ + + ]: 6 : do iwan=1, nwan(jj)
406 : 4 : rc(iwan,jj)= three*(0.769d0+half*dlog(wannspr(iwan,jj)))
407 : : ! rv(iwan,jj)= (1.475d0-half_sqrt3*dlog(wannspr(iwan,jj)))*wannspr(iwan,jj)
408 : : ! r_v suggested in JPhysChemA 113:5224 [[cite:Silvestrelli2009]]
409 : 5 : rv(iwan,jj)= (rc(iwan,jj)*wannspr(iwan,jj))/sqrt3
410 : : end do
411 : : end do
412 : 1 : corrvdw=0.0d0 !Initializing the vdW correction energy.
413 : :
414 [ + + ]: 2 : do ii=1,nsppol
415 [ + + ]: 3 : do jj=1,nsppol
416 [ + + ]: 6 : do iwan=1,nwan(ii)
417 [ + + ]: 21 : do jwan=1,nwan(jj)
418 : :
419 : : call getFu(wannspr(iwan,ii),wannspr(jwan,jj),rc(iwan,ii),rc(jwan,jj),&
420 : 16 : & newocc_wan(iwan,ii),newocc_wan(jwan,jj),fu)
421 : :
422 : : csix(iwan,jwan,ii,jj)=( ( ((wannspr(iwan,ii))**1.5d0)*&
423 : 20 : & (wannspr(jwan,jj)**three))/(two*(three**1.25d0) ) )*fu
424 : :
425 : : end do
426 : : end do
427 : : end do
428 : : end do
429 : :
430 : : ! if (nsppol == 1) then
431 : : ! csix(:,:,:,:)=sqrt2*csix(:,:,:,:) !For non spin polarized systems
432 : : ! end if
433 : :
434 : :
435 : : ! DEBUG
436 : 1 : write(std_out,*) ch10,'C6ij coefficients matrix',ch10
437 [ + + ]: 2 : do ii=1,nsppol
438 [ + + ]: 3 : do jj=1,nsppol
439 [ + + ]: 6 : do iwan=1,nwan(ii)
440 : 5 : write(std_out,*) (csix(iwan,jwan,ii,jj),jwan=1,nwan(jj)),ch10
441 : : end do
442 : : end do
443 : : end do
444 : : ! END DEBUG
445 : :
446 : : ! test k=0
447 [ + + ]: 2 : do ii=1,nsppol
448 [ + + ]: 6 : do iwan=1,nwan(ii)
449 [ + + ]: 25 : do inx=-abs(vdw_supercell(1)),abs(vdw_supercell(1))
450 [ + + ]: 124 : do iny=-abs(vdw_supercell(2)),abs(vdw_supercell(2))
451 [ + + ]: 620 : do inz=-abs(vdw_supercell(3)),abs(vdw_supercell(3))
452 [ + + ]: 1100 : do jj=1,nsppol
453 [ + + ]: 3000 : do jwan=1,nwan(jj)
454 : :
455 [ + + + - ]: 2000 : if(inx==0.and.iny==0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
456 : : ! This avoids intrafragment vdW interactions.
457 [ - + - - : 1984 : if(vdw_supercell(1)<=0.and.inx==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
- - ]
458 [ - + - - : 1984 : if(vdw_supercell(2)<=0.and.iny==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
- - ]
459 [ - + - - : 1984 : if(vdw_supercell(3)<=0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
- - ]
460 : : ! Last three conditions allow proper treatment of layered systems.
461 : :
462 : : dcenters(:,jwan,jj) = (real(inx,dp))*rprimd(:,1)+(real(iny,dp))*rprimd(:,2)+&
463 [ + + ]: 7936 : & (real(inz,dp))*rprimd(:,3)+wanncent(:,jwan,jj)
464 : : rij=sqrt((dcenters(1,jwan,jj)-wanncent(1,iwan,ii))**2+&
465 : : & (dcenters(2,jwan,jj)-wanncent(2,iwan,ii))**2+&
466 : 1984 : & (dcenters(3,jwan,jj)-wanncent(3,iwan,ii))**2)
467 : :
468 : 1984 : fij=one/(one+exp(-a*(rij/(rv(iwan,ii)+rv(jwan,jj))-one))) !Damping function.
469 : :
470 : 2500 : corrvdw = corrvdw - csix(iwan,jwan,ii,jj)*fij/(two*(rij**6)) !making the sum of eq(4) of
471 : : ! JPhysChemA 113:5224-5234 [[cite:Silvestrelli2009]]. Each term is divided by two because
472 : : ! we are counting twice within the unit cell, also the
473 : : ! interactions with neighbor cells are properly acounted for in
474 : : ! this way.
475 : :
476 : : ! write(std_out,*) 'i=',iwan, 'j=',jwan, 'C6ij=', csix(iwan,jwan)
477 : : ! write(std_out,*) 'inx=',inx, "iny=",iny, "inz=",inz, "Evdw=",&
478 : : ! & -(csix(iwan,jwan)*fij/(two*rij**6))*Ha_ev*ten**3
479 : : ! write(std_out,*) 'rnl=',rnl
480 : : end do
481 : : end do
482 : : end do
483 : : end do
484 : : end do
485 : : end do
486 : : end do
487 : :
488 [ - + ]: 1 : ABI_FREE(dcenters)
489 [ - + ]: 1 : ABI_FREE(rc)
490 [ - + ]: 1 : ABI_FREE(rv)
491 : :
492 : 1 : write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
493 : 1 : & ' vdw_xc : ',10,ch10,&
494 : 1 : & ' van der Waals correction(Ha):', corrvdw,ch10,&
495 : 2 : & ' van der Waals correction(eV):', corrvdw*Ha_ev,ch10
496 : 1 : call wrtout(std_out,message,'COLL')
497 : 1 : call wrtout(ab_out,message,'COLL')
498 : :
499 : : end if
500 : :
501 : : !vdW-WF VERSION 2: Phys. Rev. B. 85:073101 (2012) [[cite:Ambrosetti2012]]
502 : :
503 [ + + ]: 3 : if (vdw_xc==11) then
504 : :
505 [ - + + - : 4 : ABI_MALLOC(dcenters,(3,mwan,nsppol))
- + + - +
- + - - +
- + ]
506 [ - + + - : 4 : ABI_MALLOC(rv,(mwan,nsppol))
- + + - +
- + - - +
- + ]
507 [ - + + - : 3 : ABI_MALLOC(xi,(mwan,nsppol))
- + + - -
+ - + ]
508 : :
509 : : ! Calculate intermediate quantities
510 [ + + ]: 2 : do jj=1,nsppol
511 [ + + ]: 6 : do iwan=1, nwan(jj)
512 : 4 : rv(iwan,jj)= ( (1.20d0/Bohr_Ang)*wannspr(iwan,jj) )/sqrt3
513 : 5 : write(std_out,*) 'rv(iwan,jj)=',rv(iwan,jj),ch10
514 : : end do
515 : : end do
516 : :
517 : : ! C6 coefficients between WF
518 [ + + + + : 23 : csix(:,:,:,:) = 0.0d0
+ + + + ]
519 : 1 : corrvdw = 0.0d0
520 : :
521 : 1 : call ovlp_wann(mwan,nwan,nsppol,ord,wanncent,wannspr,xi)
522 : :
523 : : ! DEBUG
524 : 1 : write(std_out,*)ch10,'xi(iwan,isspol)=',ch10
525 [ + + ]: 2 : do jj=1,nsppol
526 : 2 : write(std_out,*) (xi(iwan,jj),iwan=1,nwan(jj))
527 : : end do
528 : : ! END DEBUG
529 : :
530 [ + + ]: 2 : do ii=1,nsppol
531 [ + + ]: 3 : do jj=1,nsppol
532 [ + + ]: 6 : do iwan=1,nwan(ii)
533 [ + + ]: 21 : do jwan=1,nwan(jj)
534 : :
535 : : csix(iwan,jwan,ii,jj)=onehalf*( (wannspr(iwan,ii)*wannspr(jwan,jj))**three )*&
536 : : & ((xi(iwan,ii)*xi(jwan,jj))*gama**onehalf)/( sqrt(xi(iwan,ii))*&
537 : 20 : & wannspr(iwan,ii)**onehalf + sqrt(xi(jwan,jj))*wannspr(jwan,jj)**onehalf )
538 : :
539 : : end do
540 : : end do
541 : : end do
542 : : end do
543 : :
544 : : ! if (nsppol == 1) then
545 : : ! csix(:,:,:,:)=sqrt2*csix(:,:,:,:) !For non spin polarized systems
546 : : ! end if
547 : :
548 : : ! DEBUG
549 : 1 : write(std_out,*) ch10,'C6ij coefficients:',ch10
550 [ + + ]: 2 : do ii=1,nsppol
551 [ + + ]: 3 : do jj=1,nsppol
552 [ + + ]: 6 : do iwan=1,nwan(ii)
553 : 5 : write(std_out,*) (csix(iwan,jwan,ii,jj),jwan=1,nwan(jj))
554 : : end do
555 : : end do
556 : : end do
557 : : ! END DEBUG
558 [ + + ]: 2 : do ii=1,nsppol
559 [ + + ]: 6 : do iwan=1,nwan(ii)
560 [ + + ]: 33 : do inx=-abs(vdw_supercell(1)),abs(vdw_supercell(1))
561 [ + + ]: 228 : do iny=-abs(vdw_supercell(2)),abs(vdw_supercell(2))
562 [ + + ]: 1596 : do inz=-abs(vdw_supercell(3)),abs(vdw_supercell(3))
563 [ + + ]: 2940 : do jj=1,nsppol
564 [ + + ]: 8232 : do jwan=1,nwan(jj)
565 : :
566 [ + + + - ]: 5488 : if(inx==0.and.iny==0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
567 : : ! This avoids intrafragment vdW interactions.
568 [ - + - - : 5472 : if(vdw_supercell(1)<=0.and.inx==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
- - ]
569 [ - + - - : 5472 : if(vdw_supercell(2)<=0.and.iny==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
- - ]
570 [ - + - - : 5472 : if(vdw_supercell(3)<=0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
- - ]
571 : : ! Last three conditions allow proper treatment of layered systems.
572 : :
573 : : dcenters(:,jwan,jj) = (real(inx,dp))*rprimd(:,1)+(real(iny,dp))*rprimd(:,2)+&
574 [ + + ]: 21888 : & (real(inz,dp))*rprimd(:,3)+wanncent(:,jwan,jj)
575 : : rij=sqrt((dcenters(1,jwan,jj)-wanncent(1,iwan,ii))**2+&
576 : : & (dcenters(2,jwan,jj)-wanncent(2,iwan,ii))**2+&
577 : 5472 : & (dcenters(3,jwan,jj)-wanncent(3,iwan,ii))**2)
578 : :
579 : 5472 : fij=one/(one+exp(-a*(rij/(rv(iwan,ii)+rv(jwan,jj))-one))) !Damping function.
580 : : ! DEBUG
581 : : ! write(std_out,*) 'f_i,j=',fij,ch10
582 : : ! END DEBUG
583 : 6860 : corrvdw = corrvdw - csix(iwan,jwan,ii,jj)*fij/(two*(rij**6)) !making the sum of eq(4) of
584 : : ! JPhysChemA 113:5224-5234 [[cite:Silvestrelli2009]]
585 : : end do
586 : : end do
587 : : end do
588 : : end do
589 : : end do
590 : : end do
591 : : end do
592 : :
593 : 1 : write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
594 : 1 : & ' vdw_xc : ',11,ch10,&
595 : 1 : & ' van der Waals correction(Ha):', corrvdw,ch10,&
596 : 2 : & ' van der Waals correction(eV):', corrvdw*Ha_ev,ch10
597 : 1 : call wrtout(std_out,message,'COLL')
598 : 1 : call wrtout(ab_out,message,'COLL')
599 : :
600 [ - + ]: 1 : ABI_FREE(dcenters)
601 [ - + ]: 1 : ABI_FREE(rv)
602 [ - + ]: 1 : ABI_FREE(xi)
603 : : end if
604 : :
605 : : !vdW-WF VERSION 3 (Using the long range limit of VV10 functional)
606 : :
607 [ - + ]: 3 : if(vdw_xc==12) then
608 : :
609 [ # # # # : 0 : ABI_MALLOC(dcenters,(3,mwan,nsppol))
# # # # #
# # # # #
# # ]
610 [ # # # # : 0 : ABI_MALLOC(rc,(mwan,nsppol))
# # # # #
# # # # #
# # ]
611 [ # # # # : 0 : ABI_MALLOC(rv,(mwan,nsppol))
# # # # #
# # # ]
612 : : ! Calculate intermediate quantities
613 [ # # ]: 0 : do jj=1,nsppol
614 [ # # ]: 0 : do iwan=1, nwan(jj)
615 : : ! rc(iwan,jj)= three*(0.769d0+half*dlog(wannspr(iwan,jj))) !from Silvestrelli see above.
616 : 0 : rc(iwan,jj)=three*wannspr(iwan,jj) !integral cutoff
617 : : ! rv(iwan,jj)= (1.475d0-half_sqrt3*dlog(wannspr(iwan,jj)))*wannspr(iwan,jj)
618 : 0 : rv(iwan,jj)= wannspr(iwan,jj)*sqrt3*(0.769d0+half*dlog(wannspr(iwan,jj)))
619 : : ! r_v suggested in JPhysChemA 113:5224 [[cite:Silvestrelli2009]]
620 : : end do
621 : : end do
622 : 0 : corrvdw=0.0d0 !Initializing the vdW correction energy.
623 : :
624 [ # # ]: 0 : do ii=1,nsppol
625 [ # # ]: 0 : do jj=1,nsppol
626 [ # # ]: 0 : do iwan=1,nwan(ii)
627 [ # # ]: 0 : do jwan=1,nwan(jj)
628 : :
629 : 0 : call vv10limit(wannspr(iwan,ii),wannspr(jwan,jj),rc(iwan,ii),rc(jwan,jj),fu)
630 : :
631 : 0 : csix(iwan,jwan,ii,jj)=(1296.0d0/( (wannspr(iwan,ii)*wannspr(jwan,jj) )**3))*fu
632 : : ! vv10limit needs revision as an error regarding occupations has been included
633 : : ! currently we are calculating it with 1 electron per MLWF and there is an error four-->two
634 : : end do
635 : : end do
636 : : end do
637 : : end do
638 : :
639 : : ! if (nsppol == 1) then
640 : : ! csix(:,:,:,:)=sqrt2*csix(:,:,:,:) !For non spin polarized systems
641 : : ! end if
642 : :
643 : :
644 : : ! DEBUG
645 : :
646 : 0 : write(std_out,*) ch10,'C6ij coefficients matrix',ch10
647 [ # # ]: 0 : do ii=1,nsppol
648 [ # # ]: 0 : do jj=1,nsppol
649 [ # # ]: 0 : do iwan=1,nwan(ii)
650 : 0 : write(std_out,*) (csix(iwan,jwan,ii,jj),jwan=1,nwan(jj)),ch10
651 : : end do
652 : : end do
653 : : end do
654 : : ! END DEBUG
655 : :
656 : : ! test k=0
657 [ # # ]: 0 : do ii=1,nsppol
658 [ # # ]: 0 : do iwan=1,nwan(ii)
659 [ # # ]: 0 : do inx=-abs(vdw_supercell(1)),abs(vdw_supercell(1))
660 [ # # ]: 0 : do iny=-abs(vdw_supercell(2)),abs(vdw_supercell(2))
661 [ # # ]: 0 : do inz=-abs(vdw_supercell(3)),abs(vdw_supercell(3))
662 [ # # ]: 0 : do jj=1,nsppol
663 [ # # ]: 0 : do jwan=1,nwan(jj)
664 : :
665 [ # # # # ]: 0 : if(inx==0.and.iny==0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
666 : : ! This avoids intrafragment vdW interactions.
667 [ # # # # : 0 : if(vdw_supercell(1)<=0.and.inx==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
# # ]
668 [ # # # # : 0 : if(vdw_supercell(2)<=0.and.iny==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
# # ]
669 [ # # # # : 0 : if(vdw_supercell(3)<=0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
# # ]
670 : : ! Last three conditions allow proper treatment of layered systems.
671 : :
672 : : dcenters(:,jwan,jj) = (real(inx,dp))*rprimd(:,1)+(real(iny,dp))*rprimd(:,2)+&
673 [ # # ]: 0 : & (real(inz,dp))*rprimd(:,3)+wanncent(:,jwan,jj)
674 : : rij=sqrt((dcenters(1,jwan,jj)-wanncent(1,iwan,ii))**2+&
675 : : & (dcenters(2,jwan,jj)-wanncent(2,iwan,ii))**2+&
676 : 0 : & (dcenters(3,jwan,jj)-wanncent(3,iwan,ii))**2)
677 : :
678 : 0 : fij=one/(one+exp(-a*(rij/(rv(iwan,ii)+rv(jwan,jj))-one))) !Damping function.
679 : :
680 : 0 : corrvdw = corrvdw - csix(iwan,jwan,ii,jj)*fij/(two*(rij**6)) !making the sum of eq(4) of
681 : : ! JPhysChemA 113:5224-5234 [[cite:Silvestrelli2009]]. Each term is divided by two because
682 : : ! we are counting twice within the unit cell, also the
683 : : ! interactions with neighbor cells are properly acounted for in
684 : : ! this way.
685 : :
686 : : ! write(std_out,*) 'i=',iwan, 'j=',jwan, 'C6ij=', csix(iwan,jwan)
687 : : ! write(std_out,*) 'inx=',inx, "iny=",iny, "inz=",inz, "Evdw=",&
688 : : ! & -(csix(iwan,jwan)*fij/(two*rij**6))*Ha_ev*ten**3
689 : : ! write(std_out,*) 'rnl=',rnl
690 : : end do
691 : : end do
692 : : end do
693 : : end do
694 : : end do
695 : : end do
696 : : end do
697 : :
698 [ # # ]: 0 : ABI_FREE(dcenters)
699 [ # # ]: 0 : ABI_FREE(rc)
700 [ # # ]: 0 : ABI_FREE(rv)
701 : :
702 : 0 : write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
703 : 0 : & ' vdw_xc : ',12,ch10,&
704 : 0 : & ' van der Waals correction(Ha):', corrvdw,ch10,&
705 : 0 : & ' van der Waals correction(eV):', corrvdw*Ha_ev,ch10
706 : 0 : call wrtout(std_out,message,'COLL')
707 : 0 : call wrtout(ab_out,message,'COLL')
708 : :
709 : : end if
710 : :
711 : :
712 : : !vdW-QHO-WF method.
713 : :
714 [ + + ]: 3 : if(vdw_xc==14) then
715 : :
716 : : ! There is no need of building the full set of MLWFs corresponding to the vdw_supercell
717 : : ! since the matrix elements can be computed on the fly by translating the MLWFs centers.
718 : : ! The polarizability and QHO frequencies are obteined for the MLWFs in the unit cell:
719 : :
720 [ - + + - : 3 : ABI_MALLOC(polar,(mwan))
+ - + - -
+ - + ]
721 [ - + + - : 2 : ABI_MALLOC(omega,(mwan))
- + - + ]
722 : :
723 : 1 : corrvdw=zero
724 : 1 : fu=zero
725 : :
726 [ + + ]: 2 : do isppol=1,nsppol
727 [ + + ]: 5 : polar=zero
728 [ + + ]: 5 : omega=zero
729 [ + + ]: 5 : do iwan=1,nwan(isppol)
730 : 4 : polar(iwan)=gama1*(wannspr(iwan,isppol)**3)
731 : : ! assuming Z( not zeta) is the charge of a single Wannier function, 2 for non polarized and 1 for polarized)
732 : 4 : omega(iwan)=sqrt(zeta*(3-nsppol)/polar(iwan))
733 : 5 : fu=fu+omega(iwan)
734 : : end do
735 : : !DEBUG
736 : 1 : write(std_out,*) 'Unit cell non interacting QHO energy:',ch10
737 : 1 : write(std_out,*) (1.5d0)*fu,ch10
738 : : !ENDDEBUG
739 : :
740 : : ! Total number of unit cells considered:
741 : 1 : nc=(2*abs(vdw_supercell(1))+1)*(2*abs(vdw_supercell(2))+1)*(2*abs(vdw_supercell(3))+1)
742 : : !DEBUG
743 : 1 : write(std_out,*) 'Evaluation of vdW energy from ',nc,' unit cells.',ch10
744 : :
745 : 1 : write(std_out,*) 'VdW supercell non interacting QHO energy:',ch10
746 : 1 : fu=nc*fu
747 : 1 : write(std_out,*) (1.5d0)*fu,ch10
748 : : !ENDDEBUG
749 [ + - - + : 4 : ABI_MALLOC(c_QHO,(3*mwan*nc,3*mwan*nc))
+ - - + +
- + - + -
- + - + -
+ ]
750 [ - + + - : 3 : ABI_MALLOC(Tij_dip,(3*mwan*nc,3*mwan*nc))
- + + - -
+ - + -
+ ]
751 [ - + + - : 4 : ABI_MALLOC(dcenters,(3,mwan,nsppol))
- + + - +
- + - - +
- + - + ]
752 : :
753 [ + + + + ]: 105301 : c_QHO(:,:)=zero
754 [ + + + + ]: 105301 : Tij_dip(:,:)=zero
755 : 1 : inx=0
756 : 1 : iny=0
757 : :
758 : : !writing matrix diagonal terms
759 [ + + ]: 28 : do ll=1,nc
760 [ + + ]: 136 : do iwan=1,nwan(isppol)
761 [ + + ]: 459 : do ii=1,3
762 : 324 : inx=inx+1
763 : 432 : c_QHO(inx,inx)=omega(iwan)*omega(iwan)
764 : : end do
765 : : end do
766 : : end do
767 : :
768 : : !writing terms for interactions from each cell to the central unit cell
769 : : ! icx, icy, icz labels the cells in the supercell defined by vdw_supercell
770 : : ! while inx and iny label QHO matrix elements.
771 : : ! iny should start from (nc/2)*3*mwan + 1 --> r=r_central-r_cell, and r_cell=displaced positions
772 : : ! inx starts from 1 up to (nc/2)*3*mwan +1, this includes central cell intra-interactions.
773 : :
774 : 1 : inx=0
775 [ + + ]: 3 : do icz=-abs(vdw_supercell(3)),0
776 [ + + ]: 2 : if (icz==0) then
777 : : mm=0
778 : : else
779 : 1 : mm=abs(vdw_supercell(2))
780 : : end if
781 [ + + ]: 8 : do icy=-abs(vdw_supercell(2)),mm
782 [ + + ]: 5 : if (icy==0) then
783 : 2 : ll=0
784 : : else
785 : 3 : ll=abs(vdw_supercell(1))
786 : : end if
787 [ + + ]: 20 : do icx=-abs(vdw_supercell(1)),ll
788 [ + + ]: 70 : do iwan=1,nwan(isppol)
789 [ + + ]: 221 : do ii=1,3
790 : 156 : inx=inx+1
791 : : ! loop over the MLWFs in the 'central' unit cell:
792 : 156 : iny=(nc/2)*3*mwan
793 [ + + ]: 832 : do jwan=1,nwan(isppol)
794 [ + + ]: 2652 : do jj=1,3
795 : 1872 : iny=iny+1
796 : :
797 [ - + ]: 1872 : if (inx==iny) cycle !in order to avoid digonal terms which were already computed
798 : :
799 : : dcenters(:,iwan,isppol) = (real(icx,dp))*rprimd(:,1)+(real(icy,dp))*rprimd(:,2)+&
800 [ + + ]: 7488 : & (real(icz,dp))*rprimd(:,3)+wanncent(:,iwan,isppol)
801 : :
802 [ + + ]: 7488 : rij_c = -dcenters(:,iwan,isppol)+wanncent(:,jwan,isppol)
803 : 1872 : rij = dnrm2(3,rij_c,1)
804 : : ! rij=sqrt(dot_product(rij_c,rij_c))
805 [ + + ]: 1872 : if (rij==zero) cycle
806 : : !DEBUG
807 : : ! write(std_out,*) 'rij=',rij,' inx=',inx,' iny=',iny, ch10
808 : : !ENDDEBUG
809 : : ! This corresponds to beta*sigma_ij in the original paper:
810 : 1836 : fij=beta*sqrt(wannspr(iwan,isppol)*wannspr(iwan,isppol)+wannspr(jwan,isppol)*wannspr(jwan,isppol))
811 : 1836 : erfValue = abi_derf(rij/fij)
812 : :
813 [ + + ]: 1836 : if (ii==jj) then
814 : 612 : ll=1
815 : : else
816 : 1224 : ll=0
817 : : end if ! ii==jj
818 : :
819 : : Tij_dip(inx,iny)=-((3*rij_c(ii)*rij_c(jj)-rij*rij*ll)/rij**5)* &
820 : : & (erfValue-two*rij*exp(-((rij/fij)**2))/(sqrt(pi)*fij)) + &
821 : 1836 : & 2.*two*rij_c(ii)*rij_c(jj)*exp(-((rij/fij)**2))/(sqrt(pi)*fij*fij*fij*rij*rij)
822 : :
823 : 2496 : c_QHO(inx,iny)=omega(iwan)*omega(jwan)*sqrt(polar(iwan)*polar(jwan))*Tij_dip(inx,iny)
824 : :
825 : : end do !jj=1,3
826 : : end do !jwan=1,nwan
827 : : end do !ii=1,3
828 : : end do !iwan=1,nwan
829 : : end do !icx=-abs(vdw_supercell(1)),ll
830 : : end do !icy=-abs(vdw_supercell(2)),mm
831 : : end do !icz=-abs(vdw_supercell(3)),0
832 : :
833 : :
834 : : !writing terms for interactions from the central unit cell to each cell
835 : : ! icx, icy, icz labels the cells in the supercell defined by vdw_supercell
836 : : ! while inx and iny label QHO matrix elements.
837 : : ! inx should start from (nc/2)*3*mwan + 1 --> r=-r_central+r_cell, and r_cell=displaced positions
838 : : ! iny starts from (nc/2)*3*mwan+3*man+1 to avoid central cell intra interactions which were built
839 : : ! before
840 : :
841 : 1 : iny=(nc/2)*3*mwan+3*mwan
842 : :
843 [ + + ]: 3 : do icz=0,abs(vdw_supercell(3))
844 [ + + ]: 2 : if (icz==0) then
845 : : mm=1
846 : : else
847 : 1 : mm=-abs(vdw_supercell(2))
848 : : end if
849 [ + + ]: 7 : do icy=mm,abs(vdw_supercell(2))
850 [ + + ]: 4 : if (icy==0) then
851 : 1 : ll=1
852 : : else
853 : 3 : ll=-abs(vdw_supercell(1))
854 : : end if
855 [ + + ]: 16 : do icx=ll,abs(vdw_supercell(1))
856 [ + + ]: 54 : do iwan=1,nwan(isppol)
857 [ + + ]: 170 : do ii=1,3
858 : 120 : iny=iny+1
859 : : ! loop over the MLWFs in the 'central' unit cell:
860 : 120 : inx=(nc/2)*3*mwan
861 [ + + ]: 640 : do jwan=1,nwan(isppol)
862 [ + + ]: 2040 : do jj=1,3
863 : 1440 : inx=inx+1
864 : :
865 [ - + ]: 1440 : if (inx==iny) cycle !in order to avoid digonal terms which were already computed
866 : :
867 : :
868 : : dcenters(:,iwan,isppol) = (real(icx,dp))*rprimd(:,1)+(real(icy,dp))*rprimd(:,2)+&
869 [ + + ]: 5760 : & (real(icz,dp))*rprimd(:,3)+wanncent(:,iwan,isppol)
870 : :
871 [ + + ]: 5760 : rij_c = dcenters(:,iwan,isppol)-wanncent(:,jwan,isppol)
872 : 1440 : rij = dnrm2(3,rij_c,1)
873 : : ! rij=sqrt(dot_product(rij_c,rij_c))
874 [ - + ]: 1440 : if(rij==zero) cycle
875 : : !DEBUG
876 : : ! write(std_out,*) 'rij=',rij,' inx=',inx,' iny=',iny, ch10
877 : : !ENDDEBUG
878 : : ! This corresponds to beta*sigma_ij in the original paper:
879 : 1440 : fij=beta*sqrt(wannspr(iwan,isppol)*wannspr(iwan,isppol)+wannspr(jwan,isppol)*wannspr(jwan,isppol))
880 : 1440 : erfValue = abi_derf(rij/fij)
881 : :
882 [ + + ]: 1440 : if (ii==jj) then
883 : 480 : ll=1
884 : : else
885 : 960 : ll=0
886 : : end if ! ii==jj
887 : :
888 : : Tij_dip(inx,iny)=-((3*rij_c(ii)*rij_c(jj)-rij*rij*ll)/rij**5)* &
889 : : & (erfValue-two*rij*exp(-((rij/fij)**2))/(sqrt(pi)*fij)) + &
890 : 1440 : & 2.*two*rij_c(ii)*rij_c(jj)*exp(-((rij/fij)**2))/(sqrt(pi)*fij*fij*fij*rij*rij)
891 : :
892 : 1920 : c_QHO(inx,iny)=omega(iwan)*omega(jwan)*sqrt(polar(iwan)*polar(jwan))*Tij_dip(inx,iny)
893 : :
894 : : end do !jj=1,3
895 : : end do !jwan=1,nwan
896 : : end do !ii=1,3
897 : : end do !iwan=1,nwan
898 : : end do !icx=-abs(vdw_supercell(1)),ll
899 : : end do !icy=-abs(vdw_supercell(2)),mm
900 : : end do !icz=-abs(vdw_supercell(3)),0
901 : :
902 : :
903 : : ! Here we diagonalize the matrix c_QHO and the eigenvalues come back in vector eigv
904 [ + - - + : 3 : ABI_MALLOC(matrx,((3*mwan*nc*(3*mwan*nc+1))/2))
+ - + - +
- - + -
+ ]
905 [ - + + - : 3 : ABI_MALLOC(eigv,(3*mwan*nc))
+ - + - -
+ - + -
+ ]
906 [ - + + - : 4 : ABI_MALLOC(eigvec,(3*mwan*nc,3*mwan*nc))
- + + - +
- + - - +
- + ]
907 [ + - - + : 3 : ABI_MALLOC(zhpev1,(3*2*mwan*nc-1))
+ - + - +
- - + -
+ ]
908 [ + - - + : 3 : ABI_MALLOC(zhpev2,(3*3*mwan*nc-2))
+ - + - +
- - + -
+ ]
909 [ + + ]: 52651 : matrx(:)=cmplx(zero,zero)
910 [ + + ]: 325 : do jj=1,3*mwan*nc
911 [ + + ]: 52975 : do ii=1,jj
912 : 52974 : matrx(ii+(jj-1)*jj/2)=cmplx(c_QHO(ii,jj),0.0d0)
913 : : end do
914 : : end do
915 : :
916 : : !DEBUG
917 : : ! write(std_out,*) 'Printing the real part of elements in array matrx:',ch10
918 : : ! do jj=1,3*mwan*nc*(3*mwan*nc+1)/2
919 : : ! write(std_out,*) real(matrx(jj))
920 : : ! enddo
921 : : !ENDDEBUG
922 : 1 : call ZHPEV ('N','U',3*mwan*nc,matrx,eigv,eigvec,3*mwan*nc,zhpev1,zhpev2,ier)
923 : : !DEBUG
924 : 1 : write(std_out,*) 'Last argument of ZHPEV: ier=',ch10
925 : 1 : write(std_out,*) ier,ch10
926 : 1 : write(std_out,*) 'List of c_QHO eigenvaules:',ch10
927 [ + + ]: 325 : do ll=1,3*mwan*nc
928 : 325 : write(std_out,*) eigv(ll)
929 : : end do
930 : : !ENDDEBUG
931 [ - + ]: 1 : if(ier/=0) then !vz_d
932 : 0 : ABI_ERROR('zhpev fails!') !vz_d
933 : : end if !vz_d
934 : :
935 [ - + ]: 1 : ABI_FREE(matrx)
936 [ - + ]: 1 : ABI_FREE(eigvec)
937 [ - + ]: 1 : ABI_FREE(zhpev1)
938 [ - + ]: 1 : ABI_FREE(zhpev2)
939 : :
940 [ + + ]: 326 : do ii=1,3*mwan*nc !3*nwan(isppol)
941 : 325 : corrvdw=corrvdw+sqrt(eigv(ii))
942 : : end do
943 : :
944 : : end do ! end isppol
945 : :
946 : :
947 : 1 : corrvdw=0.5*corrvdw
948 : :
949 : : !DEBUG
950 : 1 : write(std_out,*) 'Half the sum of interacting matrix eigenvalues square roots:',ch10
951 : 1 : write(std_out,*) corrvdw,ch10
952 : : !ENDDEBUG
953 : :
954 : :
955 : :
956 : 1 : corrvdw=corrvdw-1.5d0*fu
957 : :
958 [ - + ]: 1 : ABI_FREE(c_QHO)
959 [ - + ]: 1 : ABI_FREE(Tij_dip)
960 [ - + ]: 1 : ABI_FREE(dcenters)
961 [ - + ]: 1 : ABI_FREE(eigv)
962 [ - + ]: 1 : ABI_FREE(polar)
963 [ - + ]: 1 : ABI_FREE(omega)
964 : :
965 : 1 : write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
966 : 1 : & ' vdw_xc : ',14,ch10,&
967 : 1 : & ' van der Waals correction(Ha):', corrvdw,ch10,&
968 : 2 : & ' van der Waals correction(eV):', corrvdw*Ha_ev,ch10
969 : 1 : call wrtout(std_out,message,'COLL')
970 : 1 : call wrtout(ab_out,message,'COLL')
971 : :
972 : : end if ! vdw-QHO
973 : :
974 [ + - ]: 3 : if(allocated(ord))then
975 : 3 : ABI_FREE(ord)
976 : : end if
977 [ + - ]: 3 : if(allocated(wanncent))then
978 : 3 : ABI_FREE(wanncent)
979 : : end if
980 [ + - ]: 3 : if(allocated(wannspr))then
981 : 3 : ABI_FREE(wannspr)
982 : : end if
983 [ + - ]: 3 : if(allocated(newocc_wan))then
984 : 3 : ABI_FREE(newocc_wan)
985 : : end if
986 [ + - ]: 3 : if(allocated(npwf))then
987 : 3 : ABI_FREE(npwf)
988 : : end if
989 [ + - ]: 3 : if (allocated(inwan))then
990 : 3 : ABI_FREE(inwan)
991 : : end if
992 [ + + ]: 3 : if(allocated(nw))then
993 : 2 : ABI_FREE(nw)
994 : : end if
995 [ + - ]: 3 : if(allocated(nwan))then
996 : 3 : ABI_FREE(nwan)
997 : : end if
998 [ - + ]: 3 : ABI_FREE(wc_rec)
999 [ + + ]: 3 : if(allocated(amagr))then
1000 : 2 : ABI_FREE(amagr)
1001 : : end if
1002 [ - + ]: 3 : if(allocated(amawf))then
1003 : 0 : ABI_FREE(amawf)
1004 : : end if
1005 : : if(allocated(Tij_dip))then
1006 : : ABI_FREE(Tij_dip)
1007 : : end if
1008 : : if(allocated(c_QHO))then
1009 : : ABI_FREE(c_QHO)
1010 : : end if
1011 [ - + ]: 3 : if(allocated(amaspr))then
1012 : 0 : ABI_FREE(amaspr)
1013 : : end if
1014 [ - + ]: 3 : if(allocated(amaocc))then
1015 : 0 : ABI_FREE(amaocc)
1016 : : end if
1017 [ - + ]: 3 : if(allocated(tmp_cent))then
1018 : 0 : ABI_FREE(tmp_cent)
1019 : : end if
1020 [ - + ]: 3 : if(allocated(tmp_nwan))then
1021 : 0 : ABI_FREE(tmp_nwan)
1022 : : end if
1023 [ - + ]: 3 : if(allocated(tmp_spr))then
1024 : 0 : ABI_FREE(tmp_spr)
1025 : : end if
1026 [ - + ]: 3 : if(allocated(tmp_occ))then
1027 : 0 : ABI_FREE(tmp_occ)
1028 : : end if
1029 : :
1030 : :
1031 : 3 : end subroutine evdw_wannier
1032 : : !!***
1033 : :
1034 : : !!****f* ABINIT/getFu
1035 : : !! NAME
1036 : : !! getFu
1037 : : !!
1038 : : !! FUNCTION
1039 : : !! Performs double integral needed to evaluate C6
1040 : : !! coefficients. Eq. (9) in J.Phys.Chem. 113:5224 [[cite:Silvestrelli2009]]
1041 : : !!
1042 : : !! INPUTS
1043 : : !!
1044 : : !! OUTPUT
1045 : : !!
1046 : : !! SOURCE
1047 : :
1048 : 16 : subroutine getFu(sn,sl,rn,rl,occn,occl,fu) ! sn-->spread(n), sl-->spread(l), rn --> rc(n), rl --> rc(l)
1049 : :
1050 : : real(dp),intent(in)::sn,sl,rn,rl,occn,occl
1051 : : real(dp),intent(out)::fu
1052 : : !local variables
1053 : : integer::nx,ny,ix,iy
1054 : : real(dp)::deltax,deltay
1055 : : real(dp)::beta,xc,yc,y,x
1056 : 16 : real(dp),allocatable::arg1(:),res1(:),arg2(:),res2(:)
1057 : :
1058 : : ! *************************************************************************
1059 : :
1060 : 16 : ny=100
1061 : 16 : nx=100
1062 : 16 : beta=(sn/sl)**(1.5d0)
1063 : 16 : xc=rn
1064 : 16 : yc=rl
1065 : 16 : deltax=xc/(real(nx,dp)-1.d0)
1066 : 16 : deltay=yc/(real(ny,dp)-1.d0)
1067 : :
1068 [ - + ]: 16 : ABI_MALLOC(arg1,(ny))
1069 [ - + ]: 16 : ABI_MALLOC(res1,(ny))
1070 [ - + ]: 16 : ABI_MALLOC(arg2,(nx))
1071 [ - + ]: 16 : ABI_MALLOC(res2,(nx))
1072 : :
1073 [ + + ]: 1616 : do ix=1,nx
1074 : :
1075 : 1600 : x=deltax*(real(ix,dp)-1.d0)
1076 : :
1077 [ + + ]: 161600 : do iy=1,ny
1078 : 160000 : y=deltay*(real(iy,dp)-1.d0)
1079 : 161600 : arg1(iy)=( (y**2.d0)*exp(-y) )/( (exp(-x)/(beta*sqrt(occn))) + exp(-y)/(sqrt(occl)) )
1080 : : end do
1081 : :
1082 : 1600 : call simpson_int(ny,deltay,arg1,res1)
1083 : 1616 : arg2(ix)=(x**2.d0)*exp(-x)*res1(ny)
1084 : :
1085 : : end do
1086 : :
1087 : 16 : call simpson_int(nx,deltax,arg2,res2)
1088 : :
1089 : 16 : Fu = res2(nx)
1090 : :
1091 [ - + ]: 16 : ABI_FREE(arg1)
1092 [ - + ]: 16 : ABI_FREE(res1)
1093 [ - + ]: 16 : ABI_FREE(arg2)
1094 : 16 : ABI_FREE(res2)
1095 : 16 : end subroutine getFu
1096 : : !!***
1097 : :
1098 : : !!****f* ABINIT/order_wannier
1099 : : !! NAME
1100 : : !! order_wannier
1101 : : !!
1102 : : !! FUNCTION
1103 : : !! Assign each MLWF with a corresponding fragment of atoms, according
1104 : : !! to vdw_typfrag array. Assignation is done by evaluating the distance
1105 : : !! from each MLWF center to the unit cell atoms. MLWFs belong to the
1106 : : !! same fragment as their nearest atom.
1107 : : !!
1108 : : !! INPUTS
1109 : : !!
1110 : : !! OUTPUT
1111 : : !!
1112 : : !! SOURCE
1113 : 3 : subroutine order_wannier(mwan,natom,nwan,nsppol,ord,vdw_typfrag,wanncent,xcart)
1114 : :
1115 : : !Arguments
1116 : : integer, intent(in) :: mwan,natom,nsppol,nwan(nsppol),vdw_typfrag(natom) !vz_d
1117 : : integer, intent(inout) :: ord(mwan,nsppol)
1118 : : real(dp),intent(in) :: wanncent(3,mwan,nsppol),xcart(3,natom)
1119 : : !Local variables
1120 : : integer :: ii,jj,ll
1121 : : real(dp):: dis,dnrm2,mindi
1122 : 3 : real(dp), allocatable :: tmp(:)
1123 : : ! *************************************************************************
1124 : :
1125 [ + - ]: 3 : ABI_MALLOC(tmp,(3))
1126 : :
1127 [ + + ]: 6 : do ll=1,nsppol
1128 [ + + ]: 18 : do ii=1,nwan(ll)
1129 [ + + ]: 48 : tmp(:) = wanncent(:,ii,ll) - xcart(:,1)
1130 : 12 : mindi = dnrm2(3,tmp,1)
1131 : : ! mindi=sqrt( dot_product(wanncent(:,ii,ll),wanncent(:,ii,ll))+dot_product(xcart(:,1),xcart(:,1))&
1132 : : !& -2*(dot_product(wanncent(:,ii,ll),xcart(:,1))) )
1133 : 12 : ord(ii,ll)=vdw_typfrag(1)
1134 [ - + ]: 15 : do jj=2,natom
1135 [ # # ]: 0 : tmp(:) = wanncent(:,ii,ll) - xcart(:,jj)
1136 : 0 : dis = dnrm2(3,tmp,1)
1137 : : ! dis=sqrt( dot_product(wanncent(:,ii,ll),wanncent(:,ii,ll))+dot_product(xcart(:,jj),xcart(:,jj))&
1138 : : !& -2*(dot_product(wanncent(:,ii,ll),xcart(:,jj))) )
1139 [ - - ]: 12 : if(dis<=mindi) then
1140 : 0 : mindi=dis
1141 : 0 : ord(ii,ll)=vdw_typfrag(jj)
1142 : : end if
1143 : : end do
1144 : : end do
1145 : : end do
1146 : :
1147 [ - + ]: 3 : ABI_FREE(tmp)
1148 : :
1149 : 3 : end subroutine order_wannier
1150 : : !!***
1151 : :
1152 : : !!****f* ABINIT/ovlp_wann
1153 : : !! NAME
1154 : : !! ovlp_wann
1155 : : !!
1156 : : !! FUNCTION
1157 : : !! Evaluate volumen reduction of MLWFs
1158 : : !! due to intrafragment overlapping
1159 : : !!
1160 : : !! INPUTS
1161 : : !!
1162 : : !! OUTPUT
1163 : : !!
1164 : : !! SOURCE
1165 : 1 : subroutine ovlp_wann(mwan,nwan,nsppol,ord,wanncent,wannspr,xi)
1166 : :
1167 : : !Arguments
1168 : : integer, intent(in) :: mwan,nsppol,nwan(nsppol),ord(mwan,nsppol) !vz_d
1169 : : real(dp),intent(in) :: wanncent(3,mwan,nsppol),wannspr(mwan,nsppol)
1170 : : real(dp), intent(out) :: xi(mwan,nsppol)
1171 : : !Local variables
1172 : : integer :: ii,iwan,ix,iy,iz,jj,jwan,neigh,steps
1173 : : real(dp):: dis,disi,discent,veff,vfree,dnrm2
1174 : 1 : integer, allocatable :: intsec(:,:,:,:)
1175 : 1 : real(dp), allocatable :: rpoint(:), tmp(:)
1176 : : real(dp), parameter :: delt = 0.05d0 !Bohr, spatial mesh (1D) step
1177 : : ! *************************************************************************
1178 : :
1179 [ - + + - : 6 : ABI_MALLOC(intsec,(mwan,nsppol,mwan,nsppol))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1180 [ - + ]: 1 : ABI_MALLOC(rpoint,(3))
1181 [ - + ]: 1 : ABI_MALLOC(tmp,(3))
1182 [ + + + + : 26 : intsec(:,:,:,:) = 0
+ + + + ]
1183 [ + + + + ]: 6 : xi(:,:) = 0.0d0
1184 : : !detecting WF intersecting neighbors
1185 [ + + ]: 2 : do ii=1,nsppol
1186 [ + + ]: 6 : do iwan=1,nwan(ii)
1187 [ + + ]: 9 : do jj=1,nsppol
1188 [ + + ]: 24 : do jwan=1,nwan(jj)
1189 : 16 : dis = 0.0d0
1190 [ + - ]: 20 : if (ord(jwan,jj)==ord(iwan,ii)) then
1191 : :
1192 : :
1193 [ + + ]: 64 : tmp(:) = wanncent(:,iwan,ii) - wanncent(:,jwan,jj)
1194 : 16 : dis = dnrm2(3,tmp,1)
1195 : : ! dis=sqrt( dot_product(wanncent(:,iwan,ii),wanncent(:,iwan,ii))+&
1196 : : !& dot_product(wanncent(:,jwan,jj),wanncent(:,jwan,jj))&
1197 : : !& -2*( dot_product(wanncent(:,iwan,ii),wanncent(:,jwan,jj)) ) )
1198 : :
1199 [ + - ]: 16 : if ( ii == jj ) then
1200 [ + - + + ]: 16 : if ( dis<=(wannspr(iwan,ii)+wannspr(jwan,jj)).and.iwan/=jwan ) then
1201 : 12 : intsec(iwan,ii,jwan,jj) = 1
1202 : : end if
1203 : : end if
1204 [ - + ]: 16 : if ( ii /= jj) then
1205 [ # # ]: 0 : if ( dis<=(wannspr(iwan,ii)+wannspr(jwan,jj)) ) then
1206 : 0 : intsec(iwan,ii,jwan,jj) = 1
1207 : : end if
1208 : : end if
1209 : :
1210 : : end if
1211 : : end do
1212 : : end do
1213 : : end do
1214 : : end do
1215 : :
1216 : : !DEBUG
1217 : 1 : write(std_out,*) 'intsec(iwan,ii,jwan,jj)=',ch10
1218 [ + + ]: 2 : do ii=1,nsppol
1219 [ + + ]: 6 : do iwan=1,nwan(ii)
1220 [ + + ]: 9 : do jj=1,nsppol
1221 : 8 : write(std_out,*) (intsec(iwan,ii,jwan,jj),jwan=1,nwan(jj)),ch10
1222 : : end do
1223 : : end do
1224 : : end do
1225 : : !END DEBUG
1226 : : !Determining both free and effective volumes.
1227 : : !Eqs (6) and (7) in PRB 85:073101 [[cite:Ambrosetti2012]].
1228 : : !Creation of grids around each WF centre.
1229 : : !Calculation of intersection volumes.
1230 [ + + ]: 2 : do ii = 1,nsppol
1231 [ + + ]: 6 : do iwan = 1,nwan(ii)
1232 : : ! Spatial meshes and volume parameters
1233 : 4 : steps=NINT(wannspr(iwan,ii)/delt)
1234 : 4 : vfree = 0
1235 : 4 : veff = 0
1236 [ + + ]: 16 : rpoint(:) = 0.0d0
1237 [ + + ]: 248 : do iz=-steps,steps
1238 [ + + ]: 15132 : do iy=-steps,steps
1239 [ + + ]: 923052 : do ix=-steps,steps
1240 : 907924 : neigh = 0
1241 : 907924 : rpoint(1) = wanncent(1,iwan,ii) + ix*delt
1242 : 907924 : rpoint(2) = wanncent(2,iwan,ii) + iy*delt
1243 : 907924 : rpoint(3) = wanncent(3,iwan,ii) + iz*delt
1244 : :
1245 [ + + ]: 3631696 : tmp(:) = wanncent(:,iwan,ii) - rpoint(:)
1246 : 907924 : discent = dnrm2(3,tmp,1)
1247 : :
1248 : : ! discent = sqrt( dot_product(wanncent(:,iwan,ii),wanncent(:,iwan,ii))&
1249 : : !& +dot_product( rpoint(:),rpoint(:) )&
1250 : : !& -2*( dot_product( wanncent(:,iwan,ii),rpoint(:) ) ) )
1251 : :
1252 [ + + ]: 907924 : if (discent > wannspr(iwan,ii)) cycle
1253 [ + - ]: 473040 : if (discent <= wannspr(iwan,ii)) then
1254 : :
1255 : : neigh = 1
1256 [ + + ]: 916312 : do jj = 1,nsppol
1257 [ + + ]: 2748936 : do jwan = 1,nwan(jj)
1258 [ + + ]: 1832624 : if ( intsec(iwan,ii,jwan,jj) == 0 ) cycle
1259 [ - + ]: 1832624 : if ( intsec(iwan,ii,jwan,jj) == 1 ) then
1260 : :
1261 [ + + ]: 5497872 : tmp(:) = rpoint(:) - wanncent(:,jwan,jj)
1262 : 1374468 : disi = dnrm2(3,tmp,1)
1263 : : ! disi = sqrt( dot_product(rpoint(:),rpoint(:))&
1264 : : !& +dot_product( wanncent(:,jwan,jj),wanncent(:,jwan,jj) )&
1265 : : !& -2*( dot_product(rpoint(:),wanncent(:,jwan,jj)) ) )
1266 [ + + ]: 1374468 : if (disi <= wannspr(jwan,jj)) then
1267 : 574372 : neigh = neigh + 1
1268 : : end if
1269 : : end if
1270 : : end do
1271 : : end do
1272 [ + - ]: 458156 : if (nsppol==1) then
1273 : 458156 : veff = veff + 1/(real(neigh,dp)**2)
1274 : : end if
1275 [ - + ]: 458156 : if (nsppol==2) then
1276 : 0 : veff = veff + 1/real(neigh,dp)
1277 : : end if
1278 : 458156 : vfree = vfree + 1/real(neigh,dp)
1279 : : end if
1280 : : end do
1281 : : end do
1282 : : end do
1283 : : ! write(std_out,*) 'iwan=',iwan,'ii=',ii,ch10
1284 : : ! write(std_out,*) 'vfree=',vfree,'neigh=',neigh,'veff=',veff,ch10
1285 : 5 : xi(iwan,ii) = veff/vfree
1286 : : ! write(std_out,*) 'xi(iwan,ii)=',xi(iwan,ii),ch10
1287 : : end do
1288 : : end do
1289 : :
1290 [ - + ]: 1 : ABI_FREE(intsec)
1291 [ - + ]: 1 : ABI_FREE(rpoint)
1292 [ - + ]: 1 : ABI_FREE(tmp)
1293 : :
1294 : 1 : end subroutine ovlp_wann
1295 : : !!***
1296 : :
1297 : : !!****f* ABINIT/vv10limit
1298 : : !! NAME
1299 : : !! vv10limit
1300 : : !!
1301 : : !! FUNCTION
1302 : : !! Performs double integral needed to evaluate C6
1303 : : !! coefficients from the long range limit of VV10
1304 : : !! functional (Phys. Rev. A. 81:062708 (2010)) [[cite:Vydrov2010]]
1305 : : !! as expressed in terms of MLWFs.
1306 : : !!
1307 : : !! INPUTS
1308 : : !!
1309 : : !! OUTPUT
1310 : : !!
1311 : : !! SOURCE
1312 : :
1313 : 0 : subroutine vv10limit(sn,sl,rn,rl,fu) ! sn-->spread(n), sl-->spread(l), rn --> rc(n), rl --> rc(l)
1314 : :
1315 : : real(dp),intent(in)::sn,sl,rn,rl
1316 : : real(dp),intent(out)::fu
1317 : : !local variables
1318 : : integer::nx,ny,ix,iy
1319 : : real(dp)::deltax,deltay,pown,powl
1320 : : real(dp)::xc,yc,y,x,wgn,wgl,wox,woy
1321 : : real(dp),parameter :: cons = 0.0093d0 !related to local band gap model, VV10
1322 : 0 : real(dp),allocatable::arg1(:),res1(:),arg2(:),res2(:)
1323 : : ! *************************************************************************
1324 : :
1325 : 0 : ny=1000
1326 : 0 : nx=1000
1327 : :
1328 : 0 : xc=rn
1329 : 0 : yc=rl
1330 : 0 : deltax=xc/(real(nx,dp)-1.d0)
1331 : 0 : deltay=yc/(real(ny,dp)-1.d0)
1332 : :
1333 [ # # ]: 0 : ABI_MALLOC(arg1,(ny))
1334 [ # # ]: 0 : ABI_MALLOC(res1,(ny))
1335 [ # # ]: 0 : ABI_MALLOC(arg2,(nx))
1336 [ # # ]: 0 : ABI_MALLOC(res2,(nx))
1337 : :
1338 : 0 : wgn = cons*( (18.0d0/(sn*sqrt3**three))**4 )
1339 : 0 : pown = two*sqrt3/sn
1340 : 0 : wgl = cons*( (18.0d0/(sl*sqrt3**three))**4 )
1341 : 0 : powl = two*sqrt3/sl
1342 : :
1343 [ # # ]: 0 : do ix=1,nx
1344 : :
1345 : 0 : x = deltax*(real(ix,dp)-1.d0)
1346 : 0 : wox = sqrt(wgn + (four*pown/sn**two)*exp(-pown*x))
1347 : :
1348 [ # # ]: 0 : do iy=1,ny
1349 : :
1350 : 0 : y = deltay*(real(iy,dp)-1.d0)
1351 : 0 : woy = sqrt(wgl + (four*powl/sl**two)*exp(-powl*y))
1352 : :
1353 : 0 : arg1(iy)=( (y**two)*exp(-powl*y) )/( woy*(wox+woy) )
1354 : :
1355 : : end do
1356 : :
1357 : 0 : call simpson_int(ny,deltay,arg1,res1)
1358 : 0 : arg2(ix)=(x**two)*exp(-pown*x)*res1(ny)/wox
1359 : :
1360 : : end do
1361 : :
1362 : 0 : call simpson_int(nx,deltax,arg2,res2)
1363 : :
1364 : 0 : fu = res2(nx)
1365 : :
1366 : : !DEBUG
1367 : 0 : write(std_out,*) ch10,'Int argument',ch10
1368 [ # # ]: 0 : do ix=1,nx
1369 : 0 : write(std_out,*) deltax*(real(ix,dp)-1.d0), arg2(ix)
1370 : : end do
1371 : : !END DEBUG
1372 : :
1373 [ # # ]: 0 : ABI_FREE(arg1)
1374 [ # # ]: 0 : ABI_FREE(res1)
1375 [ # # ]: 0 : ABI_FREE(arg2)
1376 : 0 : ABI_FREE(res2)
1377 : 0 : end subroutine vv10limit
1378 : : !!***
1379 : :
1380 : : !!****f* ABINIT/amalgam
1381 : : !! NAME
1382 : : !! amalgam
1383 : : !!
1384 : : !! FUNCTION
1385 : : !! Amalgamates MLWFs, which are close enough,
1386 : : !! into one MLWF as suggested in J.Chem.Phys.135:154105 (2011) [[cite:Andrinopoulos2011]]
1387 : : !!
1388 : : !! INPUTS
1389 : : !!
1390 : : !! OUTPUT
1391 : : !!
1392 : : !! SOURCE
1393 : :
1394 : 2 : subroutine amalgam(amagr,ngr,nsppol,nw,mwan,ord,nwan,vdw_nfrag,wanncent,wannspr)
1395 : :
1396 : : !Arguments
1397 : : integer,intent(in) :: nsppol,mwan,vdw_nfrag
1398 : : integer,intent(in) :: ord(mwan,nsppol),nwan(nsppol)
1399 : : real(dp),intent(in):: wanncent(3,mwan,nsppol),wannspr(mwan,nsppol)
1400 : : integer,intent(out):: ngr
1401 : : integer,intent(out):: nw(nsppol,mwan/2),amagr(mwan,nsppol,mwan/2)
1402 : : !local variables
1403 : : integer :: dimen,ii,igr,isppol,iw,iwan,jj,jsppol,jwan,ll
1404 : : real(dp):: dis, dnrm2
1405 : 2 : real(dp),allocatable :: tmp(:)
1406 : : ! *************************************************************************
1407 : :
1408 [ - + ]: 2 : ABI_MALLOC(tmp,(3))
1409 : :
1410 : : !Selecting pairs of MLWFs satisfying amalgamation criteria
1411 : 2 : write(std_out,*) 'Searching for MLWFs close enough to amalgamate...',ch10
1412 : :
1413 : 2 : dimen = iabs(vdw_nfrag)
1414 : :
1415 : : !Grouping MLWFs and amalgamation
1416 : :
1417 : 2 : ngr = 0
1418 [ + + + + : 26 : amagr(:,:,:) = 0
+ + ]
1419 [ + + + + ]: 10 : nw(:,:) = 0
1420 : :
1421 [ + + ]: 4 : do ll = 1 , dimen
1422 [ + + ]: 4 : do isppol = 1 , nsppol
1423 : 2 : jsppol = isppol
1424 [ + + ]: 10 : do iwan = 2 , nwan(isppol)
1425 [ + + ]: 20 : do jwan = 1 , iwan-1
1426 : :
1427 [ + - + - ]: 18 : if (ord(iwan,isppol)==ll .and. ord(jwan,jsppol)==ll ) then
1428 : :
1429 [ + + ]: 48 : tmp(:) = wanncent(:,iwan,isppol) - wanncent(:,jwan,jsppol)
1430 : 12 : dis = dnrm2(3,tmp,1)
1431 : :
1432 : : ! dis=sqrt( dot_product(wanncent(:,iwan,isppol),wanncent(:,iwan,isppol)) &
1433 : : !& + dot_product(wanncent(:,jwan,jsppol),wanncent(:,jwan,jsppol))&
1434 : : !& - 2*(dot_product(wanncent(:,iwan,isppol),wanncent(:,jwan,jsppol))) )
1435 : :
1436 [ - + ]: 12 : if ( dis <= (wannspr(iwan,isppol) + wannspr(jwan,jsppol)) / three ) then
1437 : :
1438 [ # # # # : 0 : if ( all(amagr(:,isppol,:) /= iwan) .and. &
# # # # #
# # # #
# ]
1439 : : & all(amagr(:,jsppol,:) /= jwan) ) then
1440 : :
1441 : 0 : ngr = ngr + 1
1442 : 0 : amagr(1,isppol,ngr) = jwan
1443 : 0 : amagr(2,jsppol,ngr) = iwan
1444 : 0 : nw(isppol,ngr) = 2
1445 : 0 : cycle
1446 : :
1447 : : end if
1448 : :
1449 [ # # # # : 0 : if ( any(amagr(:,isppol,:) == iwan) .and. &
# # # # #
# # # #
# ]
1450 : : & any(amagr(:,jsppol,:) == jwan) ) cycle
1451 : :
1452 [ # # ]: 0 : do igr = 1 , mwan/2
1453 [ # # ]: 0 : do iw = 1 , mwan
1454 : :
1455 [ # # # # : 0 : if ( amagr(iw,isppol,igr) == jwan .and. &
# # # # ]
1456 : : & all(amagr(:,isppol,igr) /= iwan) ) then
1457 : 0 : nw(isppol,igr) = nw(isppol,igr) + 1
1458 : 0 : amagr(nw(isppol,igr),isppol,igr) = iwan
1459 : 0 : cycle
1460 : : end if
1461 : :
1462 [ # # # # : 0 : if ( amagr(iw,isppol,igr) == iwan .and. &
# # # # ]
1463 : 0 : & all(amagr(:,isppol,igr) /= jwan) ) then
1464 : 0 : nw(isppol,igr) = nw(isppol,igr) + 1
1465 : 0 : amagr(nw(isppol,igr),isppol,igr) = jwan
1466 : 0 : cycle
1467 : : end if
1468 : :
1469 : : end do
1470 : : end do
1471 : :
1472 : : end if !if dis < (wannspr(iwan,isppol) + wannspr(jwan,jsppol))/three
1473 : : end if !if (ord(iwan,isppol)==ll .and. ord(jwan,jsppol)==ll )
1474 : : end do !jwan
1475 : : end do !iwan
1476 : : end do !isppol
1477 : :
1478 : :
1479 [ - + ]: 4 : if (nsppol == 2) then
1480 : 0 : isppol = 1
1481 : 0 : jsppol = 2
1482 [ # # ]: 0 : do iwan = 1 , nwan(isppol)
1483 [ # # ]: 0 : do jwan = 1 , nwan(jsppol)
1484 : :
1485 [ # # # # ]: 0 : if (ord(iwan,isppol)==ll .and. ord(jwan,jsppol)==ll ) then
1486 : :
1487 [ # # ]: 0 : tmp(:) = wanncent(:,iwan,isppol) - wanncent(:,jwan,jsppol)
1488 : 0 : dis = dnrm2(3,tmp,1)
1489 : :
1490 : : ! dis=sqrt( dot_product(wanncent(:,iwan,isppol),wanncent(:,iwan,isppol)) &
1491 : : !& + dot_product(wanncent(:,jwan,jsppol),wanncent(:,jwan,jsppol))&
1492 : : !& - 2*(dot_product(wanncent(:,iwan,isppol),wanncent(:,jwan,jsppol))) )
1493 : :
1494 [ # # ]: 0 : if ( dis <= (wannspr(iwan,isppol) + wannspr(jwan,jsppol)) / three ) then
1495 : :
1496 [ # # # # : 0 : if ( all(amagr(:,isppol,:) /= iwan) .and. &
# # # # #
# # # #
# ]
1497 : : & all(amagr(:,jsppol,:) /= jwan) ) then
1498 : :
1499 : 0 : ngr = ngr + 1
1500 : 0 : amagr(1,isppol,ngr) = iwan
1501 : 0 : amagr(1,jsppol,ngr) = jwan
1502 : 0 : nw(isppol,ngr) = nw(isppol,ngr) + 1
1503 : 0 : nw(jsppol,ngr) = nw(jsppol,ngr) + 1
1504 : 0 : cycle
1505 : :
1506 : : end if
1507 : :
1508 [ # # # # : 0 : if ( any(amagr(:,isppol,:) == iwan) .and. &
# # # # #
# # # #
# ]
1509 : : & any(amagr(:,jsppol,:) == jwan) ) cycle
1510 : :
1511 [ # # ]: 0 : do igr = 1 , mwan/2
1512 [ # # ]: 0 : do iw = 1 , mwan
1513 : :
1514 [ # # # # : 0 : if ( amagr(iw,jsppol,igr) == jwan .and. &
# # # # ]
1515 : : & all(amagr(:,isppol,igr) /= iwan) ) then
1516 : 0 : nw(isppol,igr) = nw(isppol,igr) + 1
1517 : 0 : amagr(nw(isppol,igr),isppol,igr) = iwan
1518 : 0 : cycle
1519 : : end if
1520 : :
1521 [ # # # # : 0 : if ( amagr(iw,isppol,igr) == iwan .and. &
# # # # ]
1522 : 0 : & all(amagr(:,jsppol,igr) /= jwan) ) then
1523 : 0 : nw(jsppol,igr) = nw(jsppol,igr) + 1
1524 : 0 : amagr(nw(jsppol,igr),jsppol,igr) = jwan
1525 : 0 : cycle
1526 : : end if
1527 : :
1528 : : end do
1529 : : end do
1530 : :
1531 : : end if
1532 : :
1533 : : end if
1534 : :
1535 : : end do
1536 : : end do
1537 : : end if !if (nsppol == 2)
1538 : :
1539 : : end do !ll
1540 : :
1541 : 2 : write(std_out,*) 'Number of amalgamation groups:',ngr,ch10
1542 [ - + ]: 2 : if(ngr/=0)then
1543 [ # # ]: 0 : do ii = 1 , ngr
1544 [ # # ]: 0 : do isppol = 1 ,nsppol
1545 : 0 : write(std_out,*) 'Number of MLWFs in group',ii,':',nw(isppol,ii),ch10
1546 : 0 : write(std_out,*) 'MLWFs in group',ii,': WFindex,spin,group ',ch10
1547 [ # # ]: 0 : do jj = 1, nw(isppol,ii)
1548 : 0 : write(std_out,*) amagr(jj,isppol,ii),isppol,ii,ch10
1549 : : end do
1550 : : end do
1551 : : end do
1552 : : end if
1553 : :
1554 : : !DEBUG
1555 : : !write(std_out,*)' amalgam : exit '
1556 : : !write(std_out,*)' nw =',nw
1557 : : !call flush
1558 : : !ENDDEBUG
1559 : :
1560 [ - + ]: 2 : ABI_FREE(tmp)
1561 : 2 : end subroutine amalgam
1562 : : !!***
1563 : :
1564 : : end module m_evdw_wannier
1565 : : !!***
|
