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