TABLE OF CONTENTS


ABINIT/amalgam [ Functions ]

[ Top ] [ 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

PARENTS

      evdw_wannier

CHILDREN

SOURCE

1470  subroutine amalgam(amagr,ngr,nsppol,nw,mwan,ord,nwan,vdw_nfrag,wanncent,wannspr)
1471 
1472 
1473 !This section has been created automatically by the script Abilint (TD).
1474 !Do not modify the following lines by hand.
1475 #undef ABI_FUNC
1476 #define ABI_FUNC 'amalgam'
1477 !End of the abilint section
1478 
1479  implicit none
1480  !Arguments
1481  integer,intent(in) :: nsppol,mwan,vdw_nfrag
1482  integer,intent(in) :: ord(mwan,nsppol),nwan(nsppol)
1483  real(dp),intent(in):: wanncent(3,mwan,nsppol),wannspr(mwan,nsppol)
1484  integer,intent(out):: ngr
1485  integer,intent(out):: nw(nsppol,mwan/2),amagr(mwan,nsppol,mwan/2)
1486  !local variables
1487  integer :: dimen,ii,igr,isppol,iw,iwan,jj,jsppol,jwan,ll
1488  real(dp):: dis, dnrm2
1489  real(dp),allocatable :: tmp(:)
1490 ! *************************************************************************
1491 
1492  ABI_ALLOCATE(tmp,(3))
1493 
1494 !Selecting pairs of MLWFs satisfying amalgamation criteria
1495  write(std_out,*) 'Searching for MLWFs close enough to amalgamate...',ch10
1496 
1497  dimen = iabs(vdw_nfrag)
1498 
1499 !Grouping MLWFs and amalgamation
1500 
1501  ngr = 0
1502  amagr(:,:,:) = 0
1503  nw(:,:) = 0
1504 
1505  do ll = 1 , dimen
1506    do isppol = 1 , nsppol
1507      jsppol = isppol
1508      do iwan = 2 , nwan(isppol)
1509        do jwan = 1 , iwan-1
1510 
1511          if (ord(iwan,isppol)==ll .and. ord(jwan,jsppol)==ll ) then
1512 
1513            tmp(:) = wanncent(:,iwan,isppol) - wanncent(:,jwan,jsppol)
1514            dis = dnrm2(3,tmp,1)
1515 
1516 !           dis=sqrt( dot_product(wanncent(:,iwan,isppol),wanncent(:,iwan,isppol)) &
1517 !&           + dot_product(wanncent(:,jwan,jsppol),wanncent(:,jwan,jsppol))&
1518 !&           - 2*(dot_product(wanncent(:,iwan,isppol),wanncent(:,jwan,jsppol))) )
1519 
1520            if ( dis <= (wannspr(iwan,isppol) + wannspr(jwan,jsppol)) / three ) then
1521 
1522              if ( all(amagr(:,isppol,:) /= iwan) .and. &
1523 &             all(amagr(:,jsppol,:) /= jwan) ) then
1524 
1525                ngr = ngr + 1
1526                amagr(1,isppol,ngr) = jwan
1527                amagr(2,jsppol,ngr) = iwan
1528                nw(isppol,ngr) = 2
1529                cycle
1530 
1531              end if
1532 
1533              if  ( any(amagr(:,isppol,:) == iwan) .and. &
1534 &             any(amagr(:,jsppol,:) == jwan) ) cycle
1535 
1536              do igr = 1 , mwan/2
1537                do iw = 1 , mwan
1538 
1539                  if ( amagr(iw,isppol,igr) ==  jwan .and. &
1540 &                 all(amagr(:,isppol,igr) /= iwan) ) then
1541                    nw(isppol,igr) = nw(isppol,igr) + 1
1542                    amagr(nw(isppol,igr),isppol,igr) = iwan
1543                    cycle
1544                  end if
1545 
1546                  if ( amagr(iw,isppol,igr) ==  iwan .and. &
1547 &                 all(amagr(:,isppol,igr) /= jwan) ) then
1548                    nw(isppol,igr) = nw(isppol,igr) + 1
1549                    amagr(nw(isppol,igr),isppol,igr) = jwan
1550                    cycle
1551                  end if
1552 
1553                end do
1554              end do
1555 
1556            end if  !if dis < (wannspr(iwan,isppol) + wannspr(jwan,jsppol))/three
1557          end if  !if (ord(iwan,isppol)==ll .and. ord(jwan,jsppol)==ll )
1558        end do  !jwan
1559      end do  !iwan
1560    end do  !isppol
1561 
1562 
1563    if (nsppol == 2) then
1564      isppol = 1
1565      jsppol = 2
1566      do iwan = 1 , nwan(isppol)
1567        do jwan = 1 , nwan(jsppol)
1568 
1569          if (ord(iwan,isppol)==ll .and. ord(jwan,jsppol)==ll ) then
1570 
1571            tmp(:) =  wanncent(:,iwan,isppol) - wanncent(:,jwan,jsppol)
1572            dis = dnrm2(3,tmp,1)
1573 
1574 !           dis=sqrt( dot_product(wanncent(:,iwan,isppol),wanncent(:,iwan,isppol)) &
1575 !&           + dot_product(wanncent(:,jwan,jsppol),wanncent(:,jwan,jsppol))&
1576 !&           - 2*(dot_product(wanncent(:,iwan,isppol),wanncent(:,jwan,jsppol))) )
1577 
1578            if ( dis <= (wannspr(iwan,isppol) + wannspr(jwan,jsppol)) / three ) then
1579 
1580              if ( all(amagr(:,isppol,:) /= iwan) .and. &
1581 &             all(amagr(:,jsppol,:) /= jwan) ) then
1582 
1583                ngr = ngr + 1
1584                amagr(1,isppol,ngr) = iwan
1585                amagr(1,jsppol,ngr) = jwan
1586                nw(isppol,ngr) = nw(isppol,ngr) + 1
1587                nw(jsppol,ngr) = nw(jsppol,ngr) + 1
1588                cycle
1589 
1590              end if
1591 
1592              if  ( any(amagr(:,isppol,:) == iwan) .and. &
1593 &             any(amagr(:,jsppol,:) == jwan) ) cycle
1594 
1595              do igr = 1 , mwan/2
1596                do iw = 1 , mwan
1597 
1598                  if ( amagr(iw,jsppol,igr) ==  jwan .and. &
1599 &                 all(amagr(:,isppol,igr) /= iwan) ) then
1600                    nw(isppol,igr) = nw(isppol,igr) + 1
1601                    amagr(nw(isppol,igr),isppol,igr) = iwan
1602                    cycle
1603                  end if
1604 
1605                  if ( amagr(iw,isppol,igr) ==  iwan .and. &
1606 &                 all(amagr(:,jsppol,igr) /= jwan) ) then
1607                    nw(jsppol,igr) = nw(jsppol,igr) + 1
1608                    amagr(nw(jsppol,igr),jsppol,igr) = jwan
1609                    cycle
1610                  end if
1611 
1612                end do
1613              end do
1614 
1615            end if
1616 
1617          end if
1618 
1619        end do
1620      end do
1621    end if !if (nsppol == 2)
1622 
1623  end do !ll
1624 
1625  write(std_out,*) 'Number of amalgamation groups:',ngr,ch10
1626  if(ngr/=0)then
1627    do ii = 1 , ngr
1628      do isppol = 1 ,nsppol
1629        write(std_out,*) 'Number of MLWFs in group',ii,':',nw(isppol,ii),ch10
1630        write(std_out,*) 'MLWFs in group',ii,': WFindex,spin,group ',ch10
1631        do jj = 1, nw(isppol,ii)
1632          write(std_out,*) amagr(jj,isppol,ii),isppol,ii,ch10
1633        end do
1634      end do
1635    end do
1636  end if
1637 
1638 !DEBUG
1639 !write(std_out,*)' amalgam : exit '
1640 !write(std_out,*)' nw =',nw
1641 !call flush
1642 !ENDDEBUG
1643 
1644  ABI_DEALLOCATE(tmp)
1645  end subroutine amalgam

ABINIT/evdw_wannier [ Functions ]

[ Top ] [ 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.

PARENTS

      mlwfovlp

CHILDREN

SOURCE

  85  subroutine evdw_wannier(csix,corrvdw,origmwan,natom,nsppol,orignwan,tdocc_wan,vdw_nfrag,&
  86 & vdw_supercell,vdw_typfrag,vdw_xc,rprimd,wann_centres,wann_spreads,xcart)
  87 
  88 
  89 !This section has been created automatically by the script Abilint (TD).
  90 !Do not modify the following lines by hand.
  91 #undef ABI_FUNC
  92 #define ABI_FUNC 'evdw_wannier'
  93 !End of the abilint section
  94 
  95  implicit none
  96 
  97 !Arguments ------------------------------------
  98  integer , intent(in)  :: origmwan,nsppol,natom,orignwan(nsppol)
  99  integer , intent(in)  :: vdw_nfrag,vdw_supercell(3),vdw_typfrag(natom),vdw_xc
 100  real(dp), intent(in)  :: rprimd(3,3),wann_centres(3,origmwan,nsppol),wann_spreads(origmwan,nsppol)
 101  real(dp), intent(in)  :: xcart(3,natom)
 102  real(dp), intent(out) :: corrvdw
 103  real(dp), intent(out) :: csix(origmwan,origmwan,nsppol,nsppol)
 104  real(dpc), intent(in) :: tdocc_wan(origmwan,nsppol)
 105 
 106 !Local variables-------------------------------
 107  integer  :: ier,igr,icx,icy,icz,ii,inx,iny,inz,isppol,iwan
 108  integer  :: jwan,jj,ll,mm,mwan,nc,ngr,tmp_mwan,mwan_half
 109  integer, allocatable:: amagr(:,:,:),inwan(:,:),nw(:,:),nwan(:),npwf(:),ord(:,:)
 110  integer, allocatable:: tmp_nwan(:)
 111  real(dp) :: dnrm2,fij,rij,rij_c(3),fu,shift,erfValue
 112  real(dp), parameter :: a = 20.d0 !Parameter related to the damping function.
 113  real(dp), parameter :: gama = 4.5d0/(sqrt3**3) !alpha=gama*S**3.
 114  real(dp), parameter :: gama1 = 0.88d0 !alpha=gama*S**3.
 115  real(dp), parameter :: zeta = 1.30d0 !polar=zeta*(Z/omega**2).
 116  real(dp), parameter :: beta = 1.39d0 !    .
 117  real(dp), allocatable:: amawf(:,:),amaspr(:),amaocc(:),dcenters(:,:,:),rc(:,:)
 118  real(dp), allocatable:: tmp_cent(:,:,:),tmp_spr(:,:),tmp_occ(:,:)
 119  real(dp), allocatable:: rv(:,:),wanncent(:,:,:),wannspr(:,:),wc_rec(:,:,:),xi(:,:)
 120  real(dp), allocatable:: c_QHO(:,:),Tij_dip(:,:),polar(:),omega(:),eigv(:),zhpev2(:)
 121  real(dpc), allocatable :: newocc_wan(:,:)
 122  complex(dpc), allocatable :: eigvec(:,:),matrx(:),zhpev1(:)
 123  character(len=500) :: message                   ! to be uncommented, if needed
 124 ! *************************************************************************
 125 
 126 !Determining presence p-like MLWFs see J.Chem.Phys.135:154105 (2011) [[cite:Andrinopoulos2011]]
 127  ABI_ALLOCATE(npwf,(nsppol))
 128  ABI_ALLOCATE(inwan,(origmwan,nsppol))
 129  ABI_ALLOCATE(nwan,(nsppol))
 130 
 131  ll = 0
 132  npwf(:) = 0
 133  inwan(:,:) = 0
 134  do jj=1,nsppol
 135    do iwan=1,orignwan(jj)
 136      if(tdocc_wan(iwan,jj)*nsppol<=1.50d0) then
 137        npwf(jj) = npwf(jj) + 1
 138        ll = ll+1
 139        inwan(ll,jj) = iwan
 140      end if
 141    end do
 142  end do
 143 
 144  write(std_out,*) ch10,'Number of p-like MLWFs per spin pol:',ch10
 145  write(std_out,*) (npwf(ii),ii=1,nsppol), ch10
 146 
 147  mwan=origmwan+(sum(npwf(:))) !two new MLWFs per p-like MLWF
 148  nwan(:)=orignwan(:)+npwf(:)
 149 
 150 
 151  ABI_ALLOCATE(wanncent,(3,mwan,nsppol))
 152  ABI_ALLOCATE(wannspr,(mwan,nsppol))
 153  ABI_ALLOCATE(wc_rec,(3,mwan,nsppol))
 154  ABI_ALLOCATE(ord,(mwan,nsppol))
 155  ABI_ALLOCATE(newocc_wan,(mwan,nsppol))
 156 
 157  wanncent(:,:,:) = zero
 158  wannspr(:,:) = zero
 159  wc_rec(:,:,:) = zero
 160  newocc_wan(:,:) = zero
 161  ord(:,:) = zero
 162 
 163 !The vdW correction is calculated in atomic units:
 164  do ii=1,nsppol
 165    do iwan=1,orignwan(ii)
 166 !    converting to bohr**2 and then squared
 167      wanncent(:,iwan,ii)=wann_centres(:,iwan,ii)/Bohr_Ang
 168 !    write(std_out,*) "spread of WF",i, "=", wann_spreads(i)
 169      wannspr(iwan,ii)=sqrt(wann_spreads(iwan,ii)/Bohr_Ang**2)
 170      newocc_wan(iwan,ii)=tdocc_wan(iwan,ii)
 171    end do
 172  end do
 173 
 174 !write(std_out,*) 'Number of MLWFs:',ch10
 175 !do ii=1,nsppol
 176 !write(std_out,*) 'nsppol=',ii, 'nwan(nsppol)=',nwan(nsppol),ch10
 177 !end do
 178 
 179  write(std_out,*) 'Original Wannier centres and spreads:',ch10
 180  do ii=1,nsppol
 181    write(std_out,*) 'nsppol=',ii,ch10
 182    do iwan=1,orignwan(ii)
 183      write(std_out,*) (wanncent(jj,iwan,ii),jj=1,3), wannspr(iwan,ii),ch10
 184    end do
 185  end do
 186 
 187 !Translate MLWFs to the original unit cell if vdw_nfrag > 0 :
 188 
 189  if(vdw_nfrag>0)then
 190    do jj=1,nsppol
 191      call xcart2xred(orignwan(jj),rprimd,wanncent(:,1:orignwan(jj),jj), &
 192 &     wc_rec(:,1:orignwan(jj),jj))
 193 !    got centers in reduced coor
 194      do iwan=1,orignwan(jj)
 195        do ii=1,3
 196          if(wc_rec(ii,iwan,jj)<zero) then
 197            shift=REAL(CEILING(ABS(wc_rec(ii,iwan,jj))),dp)
 198            wc_rec(ii,iwan,jj) = wc_rec(ii,iwan,jj)+shift
 199          end if
 200          if(wc_rec(ii,iwan,jj)>one) then
 201            shift=-REAL(INT(wc_rec(ii,iwan,jj)),dp)
 202            wc_rec(ii,iwan,jj) = wc_rec(ii,iwan,jj)+shift
 203          end if
 204        end do
 205      end do
 206      call xred2xcart(orignwan(jj),rprimd,wanncent(:,1:orignwan(jj),jj), &
 207 &     wc_rec(:,1:orignwan(jj),jj))
 208    end do
 209 
 210 !  ====================================================================
 211 
 212    write(std_out,*) ch10,'Wannier centres translated to unit cell and spr:',ch10
 213    do jj=1,nsppol
 214      write(std_out,*) 'nsppol=',jj,ch10
 215      do iwan=1,orignwan(jj)
 216        write(std_out,*) (wanncent(ii,iwan,jj),ii=1,3), wannspr(iwan,jj)
 217      end do
 218    end do
 219  end if !vdw_nfrag>0
 220 
 221 !Spliting of p-like into 2 s-like MLWFs
 222 !Eqs. (22) and (23) of J.Chem.Phys.135:154105 (2011) [[cite:Andrinopoulos2011]]
 223 
 224  if ( any (npwf(:)/=0) ) then
 225 
 226    write(std_out,*) 'Indexes of p-like MLWFs and its spin:'
 227 
 228    do isppol=1,nsppol
 229      do jj=1,npwf(isppol)
 230 
 231        write(std_out,*) inwan(jj,isppol),isppol
 232 
 233        wanncent(1:2,orignwan(isppol)+jj,isppol) = wanncent(1:2,inwan(jj,isppol),isppol)
 234 
 235        wanncent(3,orignwan(isppol)+jj,isppol) = wanncent(3,inwan(jj,isppol),isppol)  &
 236 &       + 15.d0*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt(30.d0))
 237 
 238        wanncent(3,inwan(jj,isppol),isppol) = wanncent(3,inwan(jj,isppol),isppol)  &
 239 &       - 15.d0*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt(30.d0))
 240 
 241        wannspr(orignwan(isppol)+jj,isppol) = seven*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt2)
 242 
 243        wannspr(inwan(jj,isppol),isppol) = seven*wannspr(inwan(jj,isppol),isppol) / (eight*sqrt2)
 244 
 245        newocc_wan(orignwan(isppol)+jj,isppol) = tdocc_wan(inwan(jj,isppol),isppol) / two
 246 
 247        newocc_wan(inwan(jj,isppol),isppol) = tdocc_wan(inwan(jj,isppol),isppol) / two
 248 
 249      end do
 250    end do
 251 
 252    write(std_out,*) ch10,'Wannier centres and spreads after splitting of p-like MLWFs:',ch10
 253    do isppol=1,nsppol
 254      write(std_out,*) 'nsppol=',isppol,ch10
 255      do iwan=1,nwan(isppol)
 256        write(std_out,*) (wanncent(jj,iwan,isppol),jj=1,3), wannspr(iwan,isppol)
 257      end do
 258    end do
 259 
 260  end if ! any(npwf(:)/=0)
 261 
 262 !Asign each MLWFs to one fragment, the same as their nearest atom:
 263 
 264  call order_wannier(mwan,natom,nwan,nsppol,ord,vdw_typfrag,wanncent,xcart)
 265 
 266  write(std_out,*) ch10,'Wannier centres and fragments',ch10
 267  do ll=1,abs(vdw_nfrag)
 268    write(std_out,*) 'MLWF centers in fragment',ll,ch10
 269    do jj=1,nsppol
 270      do iwan=1,nwan(jj)
 271        if (ord(iwan,jj)==ll) then
 272          write(std_out,*) 'X', (Bohr_Ang*wanncent(ii,iwan,jj),ii=1,3),iwan,jj
 273        end if
 274      end do
 275    end do
 276  end do
 277 
 278  write(std_out,*) ch10,'Occupation Matrix diagonal terms:',ch10
 279  do ll=1,abs(vdw_nfrag)
 280    write(std_out,*) 'For MLWF centers in fragment',ll,ch10
 281    do jj=1,nsppol
 282      do iwan=1,nwan(jj)
 283        if (ord(iwan,jj)==ll) then
 284          write(std_out,*) newocc_wan(iwan,jj),ch10
 285        end if
 286      end do
 287    end do
 288  end do
 289 
 290 !Amalgamation of close MLWFs, see J.Chem.Phys.135:154105 (2011) [[cite:Andrinopoulos2011]]
 291 
 292  if (all(npwf(:)==0).and.vdw_xc/=14) then !amalgamation is done only if no p-like
 293 
 294    mwan_half=mwan/2
 295    ABI_ALLOCATE(amagr,(mwan,nsppol,mwan_half))
 296    ABI_ALLOCATE(nw,(nsppol,mwan_half))
 297    nw=0
 298 
 299    call amalgam(amagr,ngr,nsppol,nw,mwan,ord,nwan,vdw_nfrag,wanncent,wannspr)
 300 
 301 !  Calculating amalgamated centres, spreads and occupancies if any:
 302 
 303    if( any(nw(:,:) /= 0) ) then
 304 
 305      ABI_ALLOCATE(amawf,(3,ngr))
 306      ABI_ALLOCATE(amaspr,(ngr))
 307      ABI_ALLOCATE(amaocc,(ngr))
 308 
 309      amawf(:,:) = 0
 310      amaspr(:) = 0
 311      amaocc(:) = 0
 312 
 313      do igr = 1 , ngr
 314        do isppol =  1 , nsppol
 315          do ii = 1 , nw(isppol,igr)
 316 
 317            amawf(:,igr) =  amawf(:,igr) + wanncent(:,amagr(ii,isppol,igr),isppol)
 318            amaspr(igr)  =  amaspr(igr) + wannspr(amagr(ii,isppol,igr),isppol)
 319            amaocc(igr)  =  amaocc(igr) + newocc_wan(amagr(ii,isppol,igr),isppol)
 320 
 321          end do
 322        end do
 323 
 324        amawf(:,igr) = amawf(:,igr) / real(sum(nw(1:nsppol,igr)),dp )
 325        amaspr(igr)  = amaspr(igr) / real(sum(nw(1:nsppol,igr)),dp )
 326 
 327      end do
 328 
 329      write(std_out,*) ch10,'Amalgamated MLWFs Centres, Spreads and Occupancies:',ch10
 330      do igr = 1 , ngr
 331        write(std_out,*) (amawf(ii,igr),ii=1,3),amaspr(igr),amaocc(igr)
 332      end do
 333 
 334 !    Redefining centres, spreads and occps arrays:
 335      ABI_ALLOCATE(tmp_nwan,(nsppol))
 336 
 337      tmp_nwan(:) = nwan(:) - sum(nw(:,1:ngr))
 338      tmp_mwan = maxval(tmp_nwan(:))
 339 
 340      ABI_ALLOCATE(tmp_cent,(3,tmp_mwan,nsppol))
 341      ABI_ALLOCATE(tmp_spr,(tmp_mwan,nsppol))
 342      ABI_ALLOCATE(tmp_occ,(tmp_mwan,nsppol))
 343 
 344      tmp_cent(:,:,:) = zero
 345      tmp_spr(:,:) = zero
 346      tmp_occ(:,:) = zero
 347 
 348      do isppol = 1 , nsppol
 349        ii = 0
 350        do iwan = 1 , nwan(isppol)
 351 
 352          if ( any(amagr(:,isppol,:) == iwan) ) cycle
 353 
 354          ii = ii + 1
 355          tmp_cent(:,ii,isppol) = wanncent(:,iwan,isppol)
 356          tmp_spr(ii,isppol) = wannspr(iwan,isppol)
 357          tmp_occ(ii,isppol) = newocc_wan(iwan,isppol)
 358 
 359        end do
 360      end do
 361 
 362 !    Redefining wanncent, wannspr, newocc_wan:
 363 !    Even if amalgamation occurs with MLWFs of different spins
 364 !    the new WF are gathered with isppol=1 functions...
 365 
 366      nwan(1) = nwan(1) - sum(nw(1,1:ngr)) + ngr
 367 
 368      if (nsppol == 2) then
 369        nwan(2) = nwan(2) - sum(nw(2,1:ngr))
 370      end if
 371 
 372      mwan = maxval(nwan(:))
 373 
 374      do isppol = 1 , nsppol
 375        do iwan = 1 , tmp_nwan(isppol)
 376 
 377          wanncent(:,iwan,isppol) = tmp_cent(:,iwan,isppol)
 378          wannspr(iwan,isppol) = tmp_spr(iwan,isppol)
 379          newocc_wan(iwan,isppol) = tmp_occ(iwan,isppol)
 380 
 381        end do
 382      end do
 383 
 384      do igr = 1 , ngr
 385 
 386        wanncent(:,tmp_nwan(1)+igr,1) = amawf(:,igr)
 387        wannspr(tmp_nwan(1)+igr,1) = amaspr(igr)
 388        newocc_wan(tmp_nwan(1)+igr,1) = amaocc(igr)
 389 
 390      end do
 391 
 392 !    Ordering again:
 393 !    Asign each MLWFs to one fragment, the same as their nearest atom:
 394 
 395      call order_wannier(mwan,natom,nwan,nsppol,ord,vdw_typfrag,wanncent,xcart)
 396 
 397 
 398      write(std_out,*) ch10,'Full set of Wannier functions and spreads'
 399      write(std_out,*) 'after both splitting of p-like WFs and amalgamation',ch10
 400 
 401      do ll=1,abs(vdw_nfrag)
 402        write(std_out,*) 'MLWF centers and spreads in fragment',ll,ch10
 403        do jj=1,nsppol
 404          do iwan=1,nwan(jj)
 405            if (ord(iwan,jj)==ll) then
 406              write(std_out,*) 'X', (Bohr_Ang*wanncent(ii,iwan,jj),ii=1,3),Bohr_Ang*wannspr(iwan,jj)
 407            end if
 408          end do
 409        end do
 410      end do
 411 
 412    end if ! any(nw(:,:) /= 0)
 413  end if ! all((npwf(:)==0).and.vdw_xc/=14)
 414 
 415 !vdW-WF VERSION 1
 416 
 417  if(vdw_xc==10) then
 418 
 419    ABI_ALLOCATE(dcenters,(3,mwan,nsppol))
 420    ABI_ALLOCATE(rc,(mwan,nsppol))
 421    ABI_ALLOCATE(rv,(mwan,nsppol))
 422 !  Calculate intermediate quantities
 423    do jj=1,nsppol
 424      do iwan=1, nwan(jj)
 425        rc(iwan,jj)= three*(0.769d0+half*dlog(wannspr(iwan,jj)))
 426 !      rv(iwan,jj)= (1.475d0-half_sqrt3*dlog(wannspr(iwan,jj)))*wannspr(iwan,jj)
 427 !      r_v suggested in JPhysChemA 113:5224 [[cite:Silvestrelli2009]]
 428        rv(iwan,jj)= (rc(iwan,jj)*wannspr(iwan,jj))/sqrt3 
 429      end do
 430    end do
 431    corrvdw=0.0d0  !Initializing the vdW correction energy.
 432 
 433    do ii=1,nsppol
 434      do jj=1,nsppol
 435        do iwan=1,nwan(ii)
 436          do jwan=1,nwan(jj)
 437 
 438            call getFu(wannspr(iwan,ii),wannspr(jwan,jj),rc(iwan,ii),rc(jwan,jj),&
 439 &           newocc_wan(iwan,ii),newocc_wan(jwan,jj),fu)
 440 
 441            csix(iwan,jwan,ii,jj)=( ( ((wannspr(iwan,ii))**1.5d0)*&
 442 &           (wannspr(jwan,jj)**three))/(two*(three**1.25d0) ) )*fu
 443 
 444          end do
 445        end do
 446      end do
 447    end do
 448 
 449 !  if (nsppol == 1) then
 450 !  csix(:,:,:,:)=sqrt2*csix(:,:,:,:)  !For non spin polarized systems
 451 !  end if
 452 
 453 
 454 !  DEBUG
 455    write(std_out,*) ch10,'C6ij coefficients matrix',ch10
 456    do ii=1,nsppol
 457      do jj=1,nsppol
 458        do iwan=1,nwan(ii)
 459          write(std_out,*) (csix(iwan,jwan,ii,jj),jwan=1,nwan(jj)),ch10
 460        end do
 461      end do
 462    end do
 463 !  END DEBUG
 464 
 465 !  test   k=0
 466    do ii=1,nsppol
 467      do iwan=1,nwan(ii)
 468        do inx=-abs(vdw_supercell(1)),abs(vdw_supercell(1))
 469          do iny=-abs(vdw_supercell(2)),abs(vdw_supercell(2))
 470            do inz=-abs(vdw_supercell(3)),abs(vdw_supercell(3))
 471              do jj=1,nsppol
 472                do jwan=1,nwan(jj)
 473 
 474                  if(inx==0.and.iny==0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 475 !                This avoids intrafragment vdW interactions.
 476                  if(vdw_supercell(1)<=0.and.inx==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 477                  if(vdw_supercell(2)<=0.and.iny==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 478                  if(vdw_supercell(3)<=0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 479 !                Last three conditions allow proper treatment of layered systems.
 480 
 481                  dcenters(:,jwan,jj) = (real(inx,dp))*rprimd(:,1)+(real(iny,dp))*rprimd(:,2)+&
 482 &                 (real(inz,dp))*rprimd(:,3)+wanncent(:,jwan,jj)
 483                  rij=sqrt((dcenters(1,jwan,jj)-wanncent(1,iwan,ii))**2+&
 484 &                 (dcenters(2,jwan,jj)-wanncent(2,iwan,ii))**2+&
 485 &                 (dcenters(3,jwan,jj)-wanncent(3,iwan,ii))**2)
 486 
 487                  fij=one/(one+exp(-a*(rij/(rv(iwan,ii)+rv(jwan,jj))-one))) !Damping function.
 488 
 489                  corrvdw = corrvdw - csix(iwan,jwan,ii,jj)*fij/(two*(rij**6)) !making the sum of eq(4) of
 490 !                JPhysChemA 113:5224-5234 [[cite:Silvestrelli2009]]. Each term is divided by two because
 491 !                we are counting twice within the unit cell, also the
 492 !                interactions with neighbor cells are properly acounted for in
 493 !                this way.
 494 
 495 !                write(std_out,*) 'i=',iwan, 'j=',jwan, 'C6ij=', csix(iwan,jwan)
 496 !                write(std_out,*) 'inx=',inx, "iny=",iny, "inz=",inz, "Evdw=",&
 497 !                & -(csix(iwan,jwan)*fij/(two*rij**6))*Ha_ev*ten**3
 498 !                write(std_out,*) 'rnl=',rnl
 499                end do
 500              end do
 501            end do
 502          end do
 503        end do
 504      end do
 505    end do
 506 
 507    ABI_DEALLOCATE(dcenters)
 508    ABI_DEALLOCATE(rc)
 509    ABI_DEALLOCATE(rv)
 510 
 511    write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
 512 &   ' vdw_xc : ',10,ch10,&
 513 &   ' van der Waals correction(Ha):',   corrvdw,ch10,&
 514 &   ' van der Waals correction(eV):',   corrvdw*Ha_ev,ch10
 515    call wrtout(std_out,message,'COLL')
 516    call wrtout(ab_out,message,'COLL')
 517 
 518  end if
 519 
 520 !vdW-WF VERSION 2: Phys. Rev. B. 85:073101 (2012) [[cite:Ambrosetti2012]]
 521 
 522  if (vdw_xc==11) then
 523 
 524    ABI_ALLOCATE(dcenters,(3,mwan,nsppol))
 525    ABI_ALLOCATE(rv,(mwan,nsppol))
 526    ABI_ALLOCATE(xi,(mwan,nsppol))
 527 
 528 !  Calculate intermediate quantities
 529    do jj=1,nsppol
 530      do iwan=1, nwan(jj)
 531        rv(iwan,jj)= ( (1.20d0/Bohr_Ang)*wannspr(iwan,jj) )/sqrt3
 532        write(std_out,*) 'rv(iwan,jj)=',rv(iwan,jj),ch10
 533      end do
 534    end do
 535 
 536 !  C6 coefficients between WF
 537    csix(:,:,:,:) = 0.0d0
 538    corrvdw = 0.0d0
 539 
 540    call ovlp_wann(mwan,nwan,nsppol,ord,wanncent,wannspr,xi)
 541 
 542 !  DEBUG
 543    write(std_out,*)ch10,'xi(iwan,isspol)=',ch10
 544    do jj=1,nsppol
 545      write(std_out,*) (xi(iwan,jj),iwan=1,nwan(jj))
 546    end do
 547 !  END DEBUG
 548 
 549    do ii=1,nsppol
 550      do jj=1,nsppol
 551        do iwan=1,nwan(ii)
 552          do jwan=1,nwan(jj)
 553 
 554            csix(iwan,jwan,ii,jj)=onehalf*( (wannspr(iwan,ii)*wannspr(jwan,jj))**three )*&
 555 &           ((xi(iwan,ii)*xi(jwan,jj))*gama**onehalf)/( sqrt(xi(iwan,ii))*&
 556 &           wannspr(iwan,ii)**onehalf + sqrt(xi(jwan,jj))*wannspr(jwan,jj)**onehalf )
 557 
 558          end do
 559        end do
 560      end do
 561    end do
 562 
 563 !  if (nsppol == 1) then
 564 !  csix(:,:,:,:)=sqrt2*csix(:,:,:,:)  !For non spin polarized systems
 565 !  end if
 566 
 567 !  DEBUG
 568    write(std_out,*) ch10,'C6ij coefficients:',ch10
 569    do ii=1,nsppol
 570      do jj=1,nsppol
 571        do iwan=1,nwan(ii)
 572          write(std_out,*) (csix(iwan,jwan,ii,jj),jwan=1,nwan(jj))
 573        end do
 574      end do
 575    end do
 576 !  END DEBUG
 577    do ii=1,nsppol
 578      do iwan=1,nwan(ii)
 579        do inx=-abs(vdw_supercell(1)),abs(vdw_supercell(1))
 580          do iny=-abs(vdw_supercell(2)),abs(vdw_supercell(2))
 581            do inz=-abs(vdw_supercell(3)),abs(vdw_supercell(3))
 582              do jj=1,nsppol
 583                do jwan=1,nwan(jj)
 584 
 585                  if(inx==0.and.iny==0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 586 !                This avoids intrafragment vdW interactions.
 587                  if(vdw_supercell(1)<=0.and.inx==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 588                  if(vdw_supercell(2)<=0.and.iny==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 589                  if(vdw_supercell(3)<=0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 590 !                Last three conditions allow proper treatment of layered systems.
 591 
 592                  dcenters(:,jwan,jj) = (real(inx,dp))*rprimd(:,1)+(real(iny,dp))*rprimd(:,2)+&
 593 &                 (real(inz,dp))*rprimd(:,3)+wanncent(:,jwan,jj)
 594                  rij=sqrt((dcenters(1,jwan,jj)-wanncent(1,iwan,ii))**2+&
 595 &                 (dcenters(2,jwan,jj)-wanncent(2,iwan,ii))**2+&
 596 &                 (dcenters(3,jwan,jj)-wanncent(3,iwan,ii))**2)
 597 
 598                  fij=one/(one+exp(-a*(rij/(rv(iwan,ii)+rv(jwan,jj))-one))) !Damping function.
 599 !                DEBUG
 600 !                write(std_out,*) 'f_i,j=',fij,ch10
 601 !                END DEBUG
 602                  corrvdw = corrvdw - csix(iwan,jwan,ii,jj)*fij/(two*(rij**6)) !making the sum of eq(4) of
 603 !                JPhysChemA 113:5224-5234 [[cite:Silvestrelli2009]]
 604                end do
 605              end do
 606            end do
 607          end do
 608        end do
 609      end do
 610    end do
 611 
 612    write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
 613 &   ' vdw_xc : ',11,ch10,&
 614 &   ' van der Waals correction(Ha):',   corrvdw,ch10,&
 615 &   ' van der Waals correction(eV):',   corrvdw*Ha_ev,ch10
 616    call wrtout(std_out,message,'COLL')
 617    call wrtout(ab_out,message,'COLL')
 618 
 619    ABI_DEALLOCATE(dcenters)
 620    ABI_DEALLOCATE(rv)
 621    ABI_DEALLOCATE(xi)
 622  end if
 623 
 624 !vdW-WF VERSION 3 (Using the long range limit of VV10 functional)
 625 
 626  if(vdw_xc==12) then
 627 
 628    ABI_ALLOCATE(dcenters,(3,mwan,nsppol))
 629    ABI_ALLOCATE(rc,(mwan,nsppol))
 630    ABI_ALLOCATE(rv,(mwan,nsppol))
 631 !  Calculate intermediate quantities
 632    do jj=1,nsppol
 633      do iwan=1, nwan(jj)
 634 !      rc(iwan,jj)= three*(0.769d0+half*dlog(wannspr(iwan,jj))) !from Silvestrelli see above.
 635        rc(iwan,jj)=three*wannspr(iwan,jj) !integral cutoff
 636 !      rv(iwan,jj)= (1.475d0-half_sqrt3*dlog(wannspr(iwan,jj)))*wannspr(iwan,jj)
 637        rv(iwan,jj)= wannspr(iwan,jj)*sqrt3*(0.769d0+half*dlog(wannspr(iwan,jj)))
 638 !      r_v suggested in JPhysChemA 113:5224 [[cite:Silvestrelli2009]]
 639      end do
 640    end do
 641    corrvdw=0.0d0  !Initializing the vdW correction energy.
 642 
 643    do ii=1,nsppol
 644      do jj=1,nsppol
 645        do iwan=1,nwan(ii)
 646          do jwan=1,nwan(jj)
 647 
 648            call vv10limit(wannspr(iwan,ii),wannspr(jwan,jj),rc(iwan,ii),rc(jwan,jj),fu)
 649 
 650            csix(iwan,jwan,ii,jj)=(1296.0d0/( (wannspr(iwan,ii)*wannspr(jwan,jj) )**3))*fu
 651 !          vv10limit needs revision as an error regarding occupations has been included
 652 !          currently we are calculating it with 1 electron per MLWF and there is an error four-->two
 653          end do
 654        end do
 655      end do
 656    end do
 657 
 658 !  if (nsppol == 1) then
 659 !  csix(:,:,:,:)=sqrt2*csix(:,:,:,:)  !For non spin polarized systems
 660 !  end if
 661 
 662 
 663 !  DEBUG
 664 
 665    write(std_out,*) ch10,'C6ij coefficients matrix',ch10
 666    do ii=1,nsppol
 667      do jj=1,nsppol
 668        do iwan=1,nwan(ii)
 669          write(std_out,*) (csix(iwan,jwan,ii,jj),jwan=1,nwan(jj)),ch10
 670        end do
 671      end do
 672    end do
 673 !  END DEBUG
 674 
 675 !  test   k=0
 676    do ii=1,nsppol
 677      do iwan=1,nwan(ii)
 678        do inx=-abs(vdw_supercell(1)),abs(vdw_supercell(1))
 679          do iny=-abs(vdw_supercell(2)),abs(vdw_supercell(2))
 680            do inz=-abs(vdw_supercell(3)),abs(vdw_supercell(3))
 681              do jj=1,nsppol
 682                do jwan=1,nwan(jj)
 683 
 684                  if(inx==0.and.iny==0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 685 !                This avoids intrafragment vdW interactions.
 686                  if(vdw_supercell(1)<=0.and.inx==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 687                  if(vdw_supercell(2)<=0.and.iny==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 688                  if(vdw_supercell(3)<=0.and.inz==0.and.ord(jwan,jj)==ord(iwan,ii)) cycle
 689 !                Last three conditions allow proper treatment of layered systems.
 690 
 691                  dcenters(:,jwan,jj) = (real(inx,dp))*rprimd(:,1)+(real(iny,dp))*rprimd(:,2)+&
 692 &                 (real(inz,dp))*rprimd(:,3)+wanncent(:,jwan,jj)
 693                  rij=sqrt((dcenters(1,jwan,jj)-wanncent(1,iwan,ii))**2+&
 694 &                 (dcenters(2,jwan,jj)-wanncent(2,iwan,ii))**2+&
 695 &                 (dcenters(3,jwan,jj)-wanncent(3,iwan,ii))**2)
 696 
 697                  fij=one/(one+exp(-a*(rij/(rv(iwan,ii)+rv(jwan,jj))-one))) !Damping function.
 698 
 699                  corrvdw = corrvdw - csix(iwan,jwan,ii,jj)*fij/(two*(rij**6)) !making the sum of eq(4) of
 700 !                JPhysChemA 113:5224-5234 [[cite:Silvestrelli2009]]. Each term is divided by two because
 701 !                we are counting twice within the unit cell, also the
 702 !                interactions with neighbor cells are properly acounted for in
 703 !                this way.
 704 
 705 !                write(std_out,*) 'i=',iwan, 'j=',jwan, 'C6ij=', csix(iwan,jwan)
 706 !                write(std_out,*) 'inx=',inx, "iny=",iny, "inz=",inz, "Evdw=",&
 707 !                & -(csix(iwan,jwan)*fij/(two*rij**6))*Ha_ev*ten**3
 708 !                write(std_out,*) 'rnl=',rnl
 709                end do
 710              end do
 711            end do
 712          end do
 713        end do
 714      end do
 715    end do
 716 
 717    ABI_DEALLOCATE(dcenters)
 718    ABI_DEALLOCATE(rc)
 719    ABI_DEALLOCATE(rv)
 720 
 721    write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
 722 &   ' vdw_xc : ',12,ch10,&
 723 &   ' van der Waals correction(Ha):',   corrvdw,ch10,&
 724 &   ' van der Waals correction(eV):',   corrvdw*Ha_ev,ch10
 725    call wrtout(std_out,message,'COLL')
 726    call wrtout(ab_out,message,'COLL')
 727 
 728  end if
 729 
 730 
 731 !vdW-QHO-WF method.
 732 
 733  if(vdw_xc==14) then
 734 
 735 ! There is no need of building the full set of MLWFs corresponding to the vdw_supercell
 736 ! since the matrix elements can be computed on the fly by translating the MLWFs centers.
 737 ! The polarizability and QHO frequencies are obteined for the MLWFs in the unit cell:
 738 
 739    ABI_ALLOCATE(polar,(mwan))
 740    ABI_ALLOCATE(omega,(mwan))
 741 
 742    corrvdw=zero
 743    fu=zero
 744 
 745    do isppol=1,nsppol
 746      polar=zero
 747      omega=zero
 748      do iwan=1,nwan(isppol)
 749        polar(iwan)=gama1*(wannspr(iwan,isppol)**3)
 750 ! assuming Z( not zeta) is the charge of a single Wannier function, 2 for non polarized and 1 for polarized)
 751        omega(iwan)=sqrt(zeta*(3-nsppol)/polar(iwan))
 752        fu=fu+omega(iwan)
 753      end do
 754 !DEBUG
 755      write(std_out,*) 'Unit cell non interacting QHO energy:',ch10
 756      write(std_out,*) (1.5d0)*fu,ch10
 757 !ENDDEBUG
 758 
 759 !  Total number of unit cells considered:
 760      nc=(2*abs(vdw_supercell(1))+1)*(2*abs(vdw_supercell(2))+1)*(2*abs(vdw_supercell(3))+1)
 761 !DEBUG
 762      write(std_out,*) 'Evaluation of vdW energy from ',nc,' unit cells.',ch10
 763 
 764      write(std_out,*) 'VdW supercell non interacting QHO energy:',ch10
 765      fu=nc*fu
 766      write(std_out,*) (1.5d0)*fu,ch10
 767 !ENDDEBUG
 768      ABI_ALLOCATE(c_QHO,(3*mwan*nc,3*mwan*nc))
 769      ABI_ALLOCATE(Tij_dip,(3*mwan*nc,3*mwan*nc))
 770      ABI_ALLOCATE(dcenters,(3,mwan,nsppol))
 771 
 772      c_QHO(:,:)=zero
 773      Tij_dip(:,:)=zero
 774      inx=0
 775      iny=0
 776 
 777    !writing matrix diagonal terms
 778      do ll=1,nc
 779        do iwan=1,nwan(isppol)
 780          do ii=1,3
 781            inx=inx+1
 782            c_QHO(inx,inx)=omega(iwan)*omega(iwan)
 783          end do
 784        end do
 785      end do
 786 
 787    !writing terms for interactions from each cell to the central unit cell
 788    ! icx, icy, icz labels the cells in the supercell defined by vdw_supercell
 789    ! while inx and iny label QHO matrix elements.
 790    ! iny should start from (nc/2)*3*mwan + 1 --> r=r_central-r_cell, and r_cell=displaced positions
 791    ! inx starts from 1 up to (nc/2)*3*mwan +1, this includes central cell intra-interactions.
 792 
 793      inx=0
 794      do icz=-abs(vdw_supercell(3)),0
 795        if (icz==0) then
 796          mm=0
 797        else
 798          mm=abs(vdw_supercell(2))
 799        end if
 800        do icy=-abs(vdw_supercell(2)),mm
 801          if (icy==0) then
 802            ll=0
 803          else
 804            ll=abs(vdw_supercell(1))
 805          end if
 806          do icx=-abs(vdw_supercell(1)),ll
 807            do iwan=1,nwan(isppol)
 808              do ii=1,3
 809                inx=inx+1
 810        ! loop over the MLWFs in the 'central' unit cell:
 811                iny=(nc/2)*3*mwan
 812                do jwan=1,nwan(isppol)
 813                  do jj=1,3
 814                    iny=iny+1
 815 
 816                    if (inx==iny) cycle !in order to avoid digonal terms which  were already computed
 817 
 818                    dcenters(:,iwan,isppol) = (real(icx,dp))*rprimd(:,1)+(real(icy,dp))*rprimd(:,2)+&
 819 &                   (real(icz,dp))*rprimd(:,3)+wanncent(:,iwan,isppol)
 820 
 821                    rij_c = -dcenters(:,iwan,isppol)+wanncent(:,jwan,isppol)
 822                    rij = dnrm2(3,rij_c,1)
 823 !                  rij=sqrt(dot_product(rij_c,rij_c))
 824                    if (rij==zero) cycle
 825 !DEBUG
 826 !    write(std_out,*) 'rij=',rij,' inx=',inx,' iny=',iny, ch10
 827 !ENDDEBUG
 828 ! This corresponds to beta*sigma_ij in the original paper:
 829                    fij=beta*sqrt(wannspr(iwan,isppol)*wannspr(iwan,isppol)+wannspr(jwan,isppol)*wannspr(jwan,isppol))
 830                    erfValue = abi_derf(rij/fij)
 831 
 832                    if (ii==jj) then
 833                      ll=1
 834                    else
 835                      ll=0
 836                    end if ! ii==jj
 837 
 838                    Tij_dip(inx,iny)=-((3*rij_c(ii)*rij_c(jj)-rij*rij*ll)/rij**5)* &
 839 &                   (erfValue-two*rij*exp(-((rij/fij)**2))/(sqrt(pi)*fij)) + &
 840 &                   2.*two*rij_c(ii)*rij_c(jj)*exp(-((rij/fij)**2))/(sqrt(pi)*fij*fij*fij*rij*rij)
 841 
 842                    c_QHO(inx,iny)=omega(iwan)*omega(jwan)*sqrt(polar(iwan)*polar(jwan))*Tij_dip(inx,iny)
 843 
 844                  end do !jj=1,3
 845                end do  !jwan=1,nwan
 846              end do   !ii=1,3
 847            end do    !iwan=1,nwan
 848          end do     !icx=-abs(vdw_supercell(1)),ll
 849        end do      !icy=-abs(vdw_supercell(2)),mm
 850      end do       !icz=-abs(vdw_supercell(3)),0
 851 
 852 
 853    !writing terms for interactions from the central unit cell to each cell
 854    ! icx, icy, icz labels the cells in the supercell defined by vdw_supercell
 855    ! while inx and iny label QHO matrix elements.
 856    ! inx should start from (nc/2)*3*mwan + 1 --> r=-r_central+r_cell, and r_cell=displaced positions
 857    ! iny starts from  (nc/2)*3*mwan+3*man+1   to avoid central cell intra interactions which were built
 858    ! before
 859 
 860      iny=(nc/2)*3*mwan+3*mwan
 861 
 862      do icz=0,abs(vdw_supercell(3))
 863        if (icz==0) then
 864          mm=1
 865        else
 866          mm=-abs(vdw_supercell(2))
 867        end if
 868        do icy=mm,abs(vdw_supercell(2))
 869          if (icy==0) then
 870            ll=1
 871          else
 872            ll=-abs(vdw_supercell(1))
 873          end if
 874          do icx=ll,abs(vdw_supercell(1))
 875            do iwan=1,nwan(isppol)
 876              do ii=1,3
 877                iny=iny+1
 878        ! loop over the MLWFs in the 'central' unit cell:
 879                inx=(nc/2)*3*mwan
 880                do jwan=1,nwan(isppol)
 881                  do jj=1,3
 882                    inx=inx+1
 883 
 884                    if (inx==iny) cycle !in order to avoid digonal terms which  were already computed
 885 
 886 
 887                    dcenters(:,iwan,isppol) = (real(icx,dp))*rprimd(:,1)+(real(icy,dp))*rprimd(:,2)+&
 888 &                   (real(icz,dp))*rprimd(:,3)+wanncent(:,iwan,isppol)
 889 
 890                    rij_c = dcenters(:,iwan,isppol)-wanncent(:,jwan,isppol)
 891                    rij = dnrm2(3,rij_c,1)
 892 !                  rij=sqrt(dot_product(rij_c,rij_c))
 893                    if(rij==zero) cycle
 894 !DEBUG
 895 !    write(std_out,*) 'rij=',rij,' inx=',inx,' iny=',iny, ch10
 896 !ENDDEBUG
 897 ! This corresponds to beta*sigma_ij in the original paper:
 898                    fij=beta*sqrt(wannspr(iwan,isppol)*wannspr(iwan,isppol)+wannspr(jwan,isppol)*wannspr(jwan,isppol))
 899                    erfValue = abi_derf(rij/fij)
 900 
 901                    if (ii==jj) then
 902                      ll=1
 903                    else
 904                      ll=0
 905                    end if ! ii==jj
 906 
 907                    Tij_dip(inx,iny)=-((3*rij_c(ii)*rij_c(jj)-rij*rij*ll)/rij**5)* &
 908 &                   (erfValue-two*rij*exp(-((rij/fij)**2))/(sqrt(pi)*fij)) + &
 909 &                   2.*two*rij_c(ii)*rij_c(jj)*exp(-((rij/fij)**2))/(sqrt(pi)*fij*fij*fij*rij*rij)
 910 
 911                    c_QHO(inx,iny)=omega(iwan)*omega(jwan)*sqrt(polar(iwan)*polar(jwan))*Tij_dip(inx,iny)
 912 
 913                  end do !jj=1,3
 914                end do  !jwan=1,nwan
 915              end do   !ii=1,3
 916            end do    !iwan=1,nwan
 917          end do     !icx=-abs(vdw_supercell(1)),ll
 918        end do      !icy=-abs(vdw_supercell(2)),mm
 919      end do       !icz=-abs(vdw_supercell(3)),0
 920 
 921 
 922 ! Here we diagonalize the matrix c_QHO and the eigenvalues come back in vector eigv
 923      ABI_ALLOCATE(matrx,((3*mwan*nc*(3*mwan*nc+1))/2))
 924      ABI_ALLOCATE(eigv,(3*mwan*nc))
 925      ABI_ALLOCATE(eigvec,(3*mwan*nc,3*mwan*nc))
 926      ABI_ALLOCATE(zhpev1,(3*2*mwan*nc-1))
 927      ABI_ALLOCATE(zhpev2,(3*3*mwan*nc-2))
 928      matrx(:)=cmplx(zero,zero)
 929      do jj=1,3*mwan*nc
 930        do ii=1,jj
 931          matrx(ii+(jj-1)*jj/2)=cmplx(c_QHO(ii,jj),0.0d0)
 932        end do
 933      end do
 934 
 935 !DEBUG
 936 !   write(std_out,*) 'Printing the real part of elements in array matrx:',ch10
 937 !   do jj=1,3*mwan*nc*(3*mwan*nc+1)/2
 938 !    write(std_out,*) real(matrx(jj))
 939 !   enddo
 940 !ENDDEBUG
 941      call ZHPEV ('N','U',3*mwan*nc,matrx,eigv,eigvec,3*mwan*nc,zhpev1,zhpev2,ier)
 942 !DEBUG
 943      write(std_out,*) 'Last argument of ZHPEV: ier=',ch10
 944      write(std_out,*) ier,ch10
 945      write(std_out,*) 'List of c_QHO eigenvaules:',ch10
 946      do ll=1,3*mwan*nc
 947        write(std_out,*) eigv(ll)
 948      end do
 949 !ENDDEBUG
 950      if(ier/=0) then !vz_d
 951        MSG_ERROR('zhpev fails!') !vz_d
 952      end if !vz_d
 953 
 954      ABI_DEALLOCATE(matrx)
 955      ABI_DEALLOCATE(eigvec)
 956      ABI_DEALLOCATE(zhpev1)
 957      ABI_DEALLOCATE(zhpev2)
 958 
 959      do ii=1,3*mwan*nc  !3*nwan(isppol)
 960        corrvdw=corrvdw+sqrt(eigv(ii))
 961      end do
 962 
 963    end do  ! end isppol
 964 
 965 
 966    corrvdw=0.5*corrvdw
 967 
 968 !DEBUG
 969    write(std_out,*) 'Half the sum of interacting matrix eigenvalues square roots:',ch10
 970    write(std_out,*) corrvdw,ch10
 971 !ENDDEBUG
 972 
 973 
 974 
 975    corrvdw=corrvdw-1.5d0*fu
 976 
 977    ABI_DEALLOCATE(c_QHO)
 978    ABI_DEALLOCATE(Tij_dip)
 979    ABI_DEALLOCATE(dcenters)
 980    ABI_DEALLOCATE(eigv)
 981    ABI_DEALLOCATE(polar)
 982    ABI_DEALLOCATE(omega)
 983 
 984    write(message, '(2a,i2,2a,f12.6,2a,f12.6,a)' )ch10,&
 985 &   ' vdw_xc : ',14,ch10,&
 986 &   ' van der Waals correction(Ha):',   corrvdw,ch10,&
 987 &   ' van der Waals correction(eV):',   corrvdw*Ha_ev,ch10
 988    call wrtout(std_out,message,'COLL')
 989    call wrtout(ab_out,message,'COLL')
 990 
 991  end if ! vdw-QHO
 992 
 993  if(allocated(ord))then
 994    ABI_DEALLOCATE(ord)
 995  end if
 996  if(allocated(wanncent))then
 997    ABI_DEALLOCATE(wanncent)
 998  end if
 999  if(allocated(wannspr))then
1000    ABI_DEALLOCATE(wannspr)
1001  end if
1002  if(allocated(newocc_wan))then
1003    ABI_DEALLOCATE(newocc_wan)
1004  end if
1005  if(allocated(npwf))then
1006    ABI_DEALLOCATE(npwf)
1007  end if
1008  if (allocated(inwan))then
1009    ABI_DEALLOCATE(inwan)
1010  end if
1011  if(allocated(nw))then
1012    ABI_DEALLOCATE(nw)
1013  end if
1014  if(allocated(nwan))then
1015    ABI_DEALLOCATE(nwan)
1016  end if
1017  ABI_DEALLOCATE(wc_rec)
1018  if(allocated(amagr))then
1019    ABI_DEALLOCATE(amagr)
1020  end if
1021  if(allocated(amawf))then
1022    ABI_DEALLOCATE(amawf)
1023  end if
1024  if(allocated(Tij_dip))then
1025    ABI_DEALLOCATE(Tij_dip)
1026  end if
1027  if(allocated(c_QHO))then
1028    ABI_DEALLOCATE(c_QHO)
1029  end if
1030  if(allocated(amaspr))then
1031    ABI_DEALLOCATE(amaspr)
1032  end if
1033  if(allocated(amaocc))then
1034    ABI_DEALLOCATE(amaocc)
1035  end if
1036  if(allocated(tmp_cent))then
1037    ABI_DEALLOCATE(tmp_cent)
1038  end if
1039  if(allocated(tmp_nwan))then
1040    ABI_DEALLOCATE(tmp_nwan)
1041  end if
1042  if(allocated(tmp_spr))then
1043    ABI_DEALLOCATE(tmp_spr)
1044  end if
1045  if(allocated(tmp_occ))then
1046    ABI_DEALLOCATE(tmp_occ)
1047  end if
1048 
1049 
1050 end subroutine evdw_wannier

ABINIT/getFu [ Functions ]

[ Top ] [ Functions ]

NAME

 getFu

FUNCTION

  Performs double integral needed to evaluate C6
  coefficients. Eq. (9) in J.Phys.Chem. 113:5224 [[cite:Silvestrelli2009]]

INPUTS

OUTPUT

PARENTS

      evdw_wannier

CHILDREN

SOURCE

1072  subroutine getFu(sn,sl,rn,rl,occn,occl,fu) ! sn-->spread(n), sl-->spread(l), rn --> rc(n), rl --> rc(l)
1073 
1074 
1075 !This section has been created automatically by the script Abilint (TD).
1076 !Do not modify the following lines by hand.
1077 #undef ABI_FUNC
1078 #define ABI_FUNC 'getFu'
1079 !End of the abilint section
1080 
1081  implicit none
1082  real(dp),intent(in)::sn,sl,rn,rl,occn,occl
1083  real(dp),intent(out)::fu
1084  !local variables
1085  integer::nx,ny,ix,iy
1086  real(dp)::deltax,deltay
1087  real(dp)::beta,xc,yc,y,x
1088  real(dp),allocatable::arg1(:),res1(:),arg2(:),res2(:)
1089 
1090 ! *************************************************************************
1091 
1092  ny=100
1093  nx=100
1094  beta=(sn/sl)**(1.5d0)
1095  xc=rn
1096  yc=rl
1097  deltax=xc/(real(nx,dp)-1.d0)
1098  deltay=yc/(real(ny,dp)-1.d0)
1099 
1100  ABI_ALLOCATE(arg1,(ny))
1101  ABI_ALLOCATE(res1,(ny))
1102  ABI_ALLOCATE(arg2,(nx))
1103  ABI_ALLOCATE(res2,(nx))
1104 
1105  do ix=1,nx
1106 
1107    x=deltax*(real(ix,dp)-1.d0)
1108 
1109    do iy=1,ny
1110      y=deltay*(real(iy,dp)-1.d0)
1111      arg1(iy)=( (y**2.d0)*exp(-y) )/( (exp(-x)/(beta*sqrt(occn))) + exp(-y)/(sqrt(occl)) )
1112    end do
1113 
1114    call simpson_int(ny,deltay,arg1,res1)
1115    arg2(ix)=(x**2.d0)*exp(-x)*res1(ny)
1116 
1117  end do
1118 
1119  call simpson_int(nx,deltax,arg2,res2)
1120 
1121  Fu = res2(nx)
1122 
1123  ABI_DEALLOCATE(arg1)
1124  ABI_DEALLOCATE(res1)
1125  ABI_DEALLOCATE(arg2)
1126  ABI_DEALLOCATE(res2)
1127 end subroutine getFu

ABINIT/m_evdw_wannier [ Modules ]

[ Top ] [ Modules ]

NAME

 m_evdw_wannier

FUNCTION

COPYRIGHT

  Copyright (C) 2010-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_evdw_wannier
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32 
33  use m_special_funcs,   only : abi_derf
34  use m_numeric_tools,   only : simpson_int
35  use m_geometry,        only : xcart2xred, xred2xcart
36 
37  implicit none
38 
39  private

ABINIT/order_wannier [ Functions ]

[ Top ] [ 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

PARENTS

      evdw_wannier

CHILDREN

SOURCE

1150  subroutine order_wannier(mwan,natom,nwan,nsppol,ord,vdw_typfrag,wanncent,xcart)
1151 
1152 
1153 !This section has been created automatically by the script Abilint (TD).
1154 !Do not modify the following lines by hand.
1155 #undef ABI_FUNC
1156 #define ABI_FUNC 'order_wannier'
1157 !End of the abilint section
1158 
1159    implicit none
1160 !Arguments
1161    integer, intent(in)    :: mwan,natom,nsppol,nwan(nsppol),vdw_typfrag(natom) !vz_d
1162    integer, intent(inout) :: ord(mwan,nsppol)
1163    real(dp),intent(in)    :: wanncent(3,mwan,nsppol),xcart(3,natom)
1164 !Local variables
1165    integer :: ii,jj,ll
1166    real(dp):: dis,dnrm2,mindi
1167    real(dp), allocatable :: tmp(:)
1168 ! *************************************************************************
1169 
1170  ABI_ALLOCATE(tmp,(3))
1171 
1172  do ll=1,nsppol
1173    do ii=1,nwan(ll)
1174      tmp(:) = wanncent(:,ii,ll) - xcart(:,1)
1175      mindi = dnrm2(3,tmp,1)
1176 !     mindi=sqrt( dot_product(wanncent(:,ii,ll),wanncent(:,ii,ll))+dot_product(xcart(:,1),xcart(:,1))&
1177 !&     -2*(dot_product(wanncent(:,ii,ll),xcart(:,1))) )
1178      ord(ii,ll)=vdw_typfrag(1)
1179      do jj=2,natom
1180        tmp(:) = wanncent(:,ii,ll) - xcart(:,jj)
1181        dis = dnrm2(3,tmp,1)
1182 !       dis=sqrt( dot_product(wanncent(:,ii,ll),wanncent(:,ii,ll))+dot_product(xcart(:,jj),xcart(:,jj))&
1183 !&       -2*(dot_product(wanncent(:,ii,ll),xcart(:,jj))) )
1184        if(dis<=mindi) then
1185          mindi=dis
1186          ord(ii,ll)=vdw_typfrag(jj)
1187        end if
1188      end do
1189    end do
1190  end do
1191 
1192  ABI_DEALLOCATE(tmp)
1193 
1194  end subroutine order_wannier

ABINIT/ovlp_wann [ Functions ]

[ Top ] [ Functions ]

NAME

 ovlp_wann

FUNCTION

  Evaluate volumen reduction of MLWFs
  due to intrafragment overlapping

INPUTS

OUTPUT

PARENTS

      evdw_wannier

CHILDREN

SOURCE

1215  subroutine ovlp_wann(mwan,nwan,nsppol,ord,wanncent,wannspr,xi)
1216 
1217 
1218 !This section has been created automatically by the script Abilint (TD).
1219 !Do not modify the following lines by hand.
1220 #undef ABI_FUNC
1221 #define ABI_FUNC 'ovlp_wann'
1222 !End of the abilint section
1223 
1224    implicit none
1225 !Arguments
1226    integer, intent(in)  :: mwan,nsppol,nwan(nsppol),ord(mwan,nsppol) !vz_d
1227    real(dp),intent(in)  :: wanncent(3,mwan,nsppol),wannspr(mwan,nsppol)
1228    real(dp), intent(out)  :: xi(mwan,nsppol)
1229 !Local variables
1230    integer :: ii,iwan,ix,iy,iz,jj,jwan,neigh,steps
1231    real(dp):: dis,disi,discent,veff,vfree,dnrm2
1232    integer, allocatable :: intsec(:,:,:,:)
1233    real(dp), allocatable :: rpoint(:), tmp(:)
1234    real(dp), parameter :: delt = 0.05d0 !Bohr, spatial mesh (1D) step
1235 ! *************************************************************************
1236 
1237  ABI_ALLOCATE(intsec,(mwan,nsppol,mwan,nsppol))
1238  ABI_ALLOCATE(rpoint,(3))
1239  ABI_ALLOCATE(tmp,(3))
1240  intsec(:,:,:,:) = 0
1241  xi(:,:) = 0.0d0
1242 !detecting WF intersecting neighbors
1243  do ii=1,nsppol
1244    do iwan=1,nwan(ii)
1245      do jj=1,nsppol
1246        do jwan=1,nwan(jj)
1247          dis = 0.0d0
1248          if (ord(jwan,jj)==ord(iwan,ii)) then
1249 
1250 
1251            tmp(:) = wanncent(:,iwan,ii) - wanncent(:,jwan,jj)
1252            dis =  dnrm2(3,tmp,1)
1253 !           dis=sqrt(  dot_product(wanncent(:,iwan,ii),wanncent(:,iwan,ii))+&
1254 !&           dot_product(wanncent(:,jwan,jj),wanncent(:,jwan,jj))&
1255 !&           -2*( dot_product(wanncent(:,iwan,ii),wanncent(:,jwan,jj)) )  )
1256 
1257            if ( ii == jj ) then
1258              if ( dis<=(wannspr(iwan,ii)+wannspr(jwan,jj)).and.iwan/=jwan ) then
1259                intsec(iwan,ii,jwan,jj) = 1
1260              end if
1261            end if
1262            if ( ii /= jj) then
1263              if ( dis<=(wannspr(iwan,ii)+wannspr(jwan,jj)) ) then
1264                intsec(iwan,ii,jwan,jj) = 1
1265              end if
1266            end if
1267 
1268          end if
1269        end do
1270      end do
1271    end do
1272  end do
1273 
1274 !DEBUG
1275  write(std_out,*) 'intsec(iwan,ii,jwan,jj)=',ch10
1276  do ii=1,nsppol
1277    do iwan=1,nwan(ii)
1278      do jj=1,nsppol
1279        write(std_out,*) (intsec(iwan,ii,jwan,jj),jwan=1,nwan(jj)),ch10
1280      end do
1281    end do
1282  end do
1283 !END DEBUG
1284 !Determining both free and effective volumes.
1285 !Eqs (6) and (7) in PRB 85:073101 [[cite:Ambrosetti2012]].
1286 !Creation of grids around each WF centre.
1287 !Calculation of intersection volumes.
1288  do ii = 1,nsppol
1289    do iwan = 1,nwan(ii)
1290 !    Spatial meshes and volume parameters
1291      steps=NINT(wannspr(iwan,ii)/delt)
1292      vfree = 0
1293      veff = 0
1294      rpoint(:) = 0.0d0
1295      do iz=-steps,steps
1296        do iy=-steps,steps
1297          do ix=-steps,steps
1298            neigh = 0
1299            rpoint(1) = wanncent(1,iwan,ii) + ix*delt
1300            rpoint(2) = wanncent(2,iwan,ii) + iy*delt
1301            rpoint(3) = wanncent(3,iwan,ii) + iz*delt
1302 
1303            tmp(:) = wanncent(:,iwan,ii) - rpoint(:)
1304            discent = dnrm2(3,tmp,1)
1305 
1306 !           discent = sqrt( dot_product(wanncent(:,iwan,ii),wanncent(:,iwan,ii))&
1307 !&           +dot_product( rpoint(:),rpoint(:) )&
1308 !&           -2*( dot_product( wanncent(:,iwan,ii),rpoint(:) ) ) )
1309 
1310            if (discent > wannspr(iwan,ii)) cycle
1311            if (discent <= wannspr(iwan,ii)) then
1312 
1313              neigh = 1
1314              do jj = 1,nsppol
1315                do jwan = 1,nwan(jj)
1316                  if ( intsec(iwan,ii,jwan,jj) == 0 ) cycle
1317                  if ( intsec(iwan,ii,jwan,jj) == 1 ) then
1318 
1319                    tmp(:) = rpoint(:) - wanncent(:,jwan,jj)
1320                    disi = dnrm2(3,tmp,1)
1321 !                   disi = sqrt( dot_product(rpoint(:),rpoint(:))&
1322 !&                   +dot_product( wanncent(:,jwan,jj),wanncent(:,jwan,jj) )&
1323 !&                   -2*( dot_product(rpoint(:),wanncent(:,jwan,jj)) ) )
1324                    if (disi <= wannspr(jwan,jj)) then
1325                      neigh = neigh + 1
1326                    end if
1327                  end if
1328                end do
1329              end do
1330              if (nsppol==1) then
1331                veff = veff + 1/(real(neigh,dp)**2)
1332              end if
1333              if (nsppol==2) then
1334                veff = veff + 1/real(neigh,dp)
1335              end if
1336              vfree = vfree + 1/real(neigh,dp)
1337            end if
1338          end do
1339        end do
1340      end do
1341 !    write(std_out,*) 'iwan=',iwan,'ii=',ii,ch10
1342 !    write(std_out,*) 'vfree=',vfree,'neigh=',neigh,'veff=',veff,ch10
1343      xi(iwan,ii) = veff/vfree
1344 !    write(std_out,*) 'xi(iwan,ii)=',xi(iwan,ii),ch10
1345    end do
1346  end do
1347 
1348  ABI_DEALLOCATE(intsec)
1349  ABI_DEALLOCATE(rpoint)
1350  ABI_DEALLOCATE(tmp)
1351 
1352  end subroutine ovlp_wann

ABINIT/vv10limit [ Functions ]

[ Top ] [ 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

PARENTS

      evdw_wannier

CHILDREN

SOURCE

1376  subroutine vv10limit(sn,sl,rn,rl,fu) ! sn-->spread(n), sl-->spread(l), rn --> rc(n), rl --> rc(l)
1377 
1378 
1379 !This section has been created automatically by the script Abilint (TD).
1380 !Do not modify the following lines by hand.
1381 #undef ABI_FUNC
1382 #define ABI_FUNC 'vv10limit'
1383 !End of the abilint section
1384 
1385  implicit none
1386  real(dp),intent(in)::sn,sl,rn,rl
1387  real(dp),intent(out)::fu
1388  !local variables
1389  integer::nx,ny,ix,iy
1390  real(dp)::deltax,deltay,pown,powl
1391  real(dp)::xc,yc,y,x,wgn,wgl,wox,woy
1392  real(dp),parameter :: cons = 0.0093d0 !related to local band gap model, VV10
1393  real(dp),allocatable::arg1(:),res1(:),arg2(:),res2(:)
1394 ! *************************************************************************
1395 
1396  ny=1000
1397  nx=1000
1398 
1399  xc=rn
1400  yc=rl
1401  deltax=xc/(real(nx,dp)-1.d0)
1402  deltay=yc/(real(ny,dp)-1.d0)
1403 
1404  ABI_ALLOCATE(arg1,(ny))
1405  ABI_ALLOCATE(res1,(ny))
1406  ABI_ALLOCATE(arg2,(nx))
1407  ABI_ALLOCATE(res2,(nx))
1408 
1409  wgn = cons*( (18.0d0/(sn*sqrt3**three))**4 )
1410  pown = two*sqrt3/sn
1411  wgl = cons*( (18.0d0/(sl*sqrt3**three))**4 )
1412  powl = two*sqrt3/sl
1413 
1414  do ix=1,nx
1415 
1416    x = deltax*(real(ix,dp)-1.d0)
1417    wox = sqrt(wgn + (four*pown/sn**two)*exp(-pown*x))
1418 
1419    do iy=1,ny
1420 
1421      y = deltay*(real(iy,dp)-1.d0)
1422      woy = sqrt(wgl + (four*powl/sl**two)*exp(-powl*y))
1423 
1424      arg1(iy)=( (y**two)*exp(-powl*y) )/( woy*(wox+woy) )
1425 
1426    end do
1427 
1428    call simpson_int(ny,deltay,arg1,res1)
1429    arg2(ix)=(x**two)*exp(-pown*x)*res1(ny)/wox
1430 
1431  end do
1432 
1433  call simpson_int(nx,deltax,arg2,res2)
1434 
1435  fu = res2(nx)
1436 
1437 !DEBUG
1438  write(std_out,*) ch10,'Int argument',ch10
1439  do ix=1,nx
1440    write(std_out,*) deltax*(real(ix,dp)-1.d0), arg2(ix)
1441  end do
1442 !END DEBUG
1443 
1444  ABI_DEALLOCATE(arg1)
1445  ABI_DEALLOCATE(res1)
1446  ABI_DEALLOCATE(arg2)
1447  ABI_DEALLOCATE(res2)
1448 end subroutine vv10limit