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 : !!***
|
