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