## ABINIT/bracketing [ Functions ]

[ Top ] [ Functions ]

NAME

 bracketing


FUNCTION

 bracket a minimun of a function f


INPUTS

 dp_dum_vdp: the function of which the mimimum should be bracketted


OUTPUT

 b= last member of the bracketing triplet a < x < b
fa,fx,fb= value of the function at dp_dum_vdp(v(:)+y*grad(:))


SIDE EFFECTS

 v: the initial vector for the function (return unchanged)
grad: the direction on which the bracketting is to be performed (return unchanged)
a,x: two members of the bracketing triplet (see b)


PARENTS

      linmin


CHILDREN

SOURCE

3186 subroutine bracketing (nv1,nv2,dp_dum_v2dp,v,grad,a,x,b,fa,fx,fb)
3187
3188
3189 !This section has been created automatically by the script Abilint (TD).
3190 !Do not modify the following lines by hand.
3191 #undef ABI_FUNC
3192 #define ABI_FUNC 'bracketing'
3193 !End of the abilint section
3194
3195  implicit none
3196
3197 !Arguments ------------------------------------
3198 include "dummy_functions.inc"
3199 !scalars
3200  integer,intent(in) :: nv1,nv2
3201  real(dp),intent(inout) :: a,x
3202  real(dp),intent(out) :: b,fa,fb,fx
3203 !arrays
3205
3206 !Local variables-------------------------------
3207 !scalars
3208  real(dp),parameter :: maglimit=10000.0_dp
3209  real(dp) :: c,fu,q,r,u,ulim
3210
3211 ! *************************************************************************
3212
3215  if(fx > fa) then
3216   c=a
3217   a=x
3218   x=c
3219   c=fa
3220   fa=fx
3221   fx=c
3222  end if
3223  b=x+gold*(x-a)
3225  do
3226   if (fx <= fb) return
3227   r=(x-a)*(fx-fb)
3228   q=(x-b)*(fx-fa)
3229   u=x-((x-b)*q-(x-a)*r)/(two*sign(max(abs(q-r),smallest_real),q-r))
3230   ulim=x+maglimit*(b-x)
3231   if((x-u)*(u-b) > zero) then
3233    if(fu < fb) then
3234     a=x
3235     fa=fx
3236     x=u
3237     fx=fu
3238     return
3239    else if (fx < fu) then
3240     b=u
3241     fb=fu
3242     return
3243    end if
3244    u=b+gold*(b-x)
3246   else if((b-u)*(u-ulim) > zero) then
3248    if(fu<fb) then
3249     x=b
3250     b=u
3251     u=b+gold*(b-x)
3252     fx=fb
3253     fb=fu
3255    end if
3256   else if((u-ulim)*(ulim-b) >= zero) then
3257    u=ulim
3259   else
3260    u=b+gold*(b-x)
3262   end if
3263   a=x
3264   x=b
3265   b=u
3266   fa=fx
3267   fx=fb
3268   fb=fu
3269  end do
3270
3271 end subroutine bracketing


## ABINIT/brent [ Functions ]

[ Top ] [ Functions ]

NAME

 brent


FUNCTION

 minimizes a function along a line


INPUTS

 dp_dum_vdp: function  to be minimized (return a dp from a vector of dp)
vdp_dum_vdp: derivative of the function (return a vector of dp from a vector of dp)
itmax: number of iterations allowed
tol: tolerance on error. It depend on the precision of the numbers
(usualy chosen as sqrt(max precision available with your floating point reresentation))
ax,xx,bx: a bracketing triplet around the minimum to be find


OUTPUT

 xmin: value such that dp_dum_vdp(v(:)+xmin*grad(:)) is minimum


SIDE EFFECTS

 grad(:): direction along which the minimization is performed
v(:): starting and ending point of the minimization


PARENTS

 linmin


CHILDREN

 dotproduct


SOURCE

3304 function brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,itmax,v,grad,ax,xx,bx,tol,xmin)
3305
3306
3307 !This section has been created automatically by the script Abilint (TD).
3308 !Do not modify the following lines by hand.
3309 #undef ABI_FUNC
3310 #define ABI_FUNC 'brent'
3311 !End of the abilint section
3312
3313  implicit none
3314
3315 !Arguments ------------------------------------
3316 include "dummy_functions.inc"
3317 !scalars
3318  integer,intent(in) :: itmax,nv1,nv2
3319  real(dp) :: brent
3320  real(dp),intent(in) :: ax,bx,tol,xx
3321  real(dp),intent(out) :: xmin
3322 !arrays
3324
3325 !Local variables-------------------------------
3326 !scalars
3327  integer :: iter
3328  real(dp) :: a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,vv,w
3329  real(dp) :: x,xm,zeps
3330  logical :: ok1,ok2,ok3,ok4
3331
3332 !************************************************************************
3333  zeps=epsilon(ax*real(1e-2,dp))
3334  a=min(ax,bx)
3335  b=max(ax,bx)
3336  vv=xx
3337  w=xx
3338  x=xx
3339  e=zero
3341  fv=fx
3342  fw=fx
3343 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of
3345 !but for instance renormilizing the density if brent is used on a density...
3346 !vp(:,:) = v(:,:)
3350  dv=dx
3351  dw=dx
3352  do iter=1,itmax
3353   xm=half*(a+b)
3354   tol1=tol*abs(x)+zeps
3355   tol2=two*tol1
3356   if(abs(x-xm) <= (tol2-half*(b-a))) then
3357    exit
3358   end if
3359   if(abs(e) > tol1) then
3360    d1=two*(b-a)
3361    d2=d1
3362    if(dw /= dx) d1=(w-x)*dx/(dx-dw)
3363    if(dv /= dx) d2=(vv-x)*dx/(dx-dv)
3364    u1=x+d1
3365    u2=x+d2
3366    ok1=((a-u1)*(u1-b)>zero).and.(dx*d1<=zero)
3367    ok2=((a-u2)*(u2-b)>zero).and.(dx*d2<=zero)
3368    olde=e
3369    e=d
3370    if(ok1.or.ok2) then
3371     if(ok1.and.ok2) then
3372      d=merge(d1,d2,abs(d1)<abs(d2))
3373     else
3374      d=merge(d1,d2,ok1)
3375     end if
3376     if(abs(d)<=abs(half*olde)) then
3377      u=x+d
3378      if(((u-a)<tol2).or.((b-u)<tol2)) d=sign(tol1,xm-x)
3379     else
3380      e=merge(a,b,dx>=zero)-x
3381      d=half*e
3382     end if
3383    else
3384     e=merge(a,b,dx>=zero)-x
3385     d=half*e
3386    end if
3387   else
3388    e=merge(a,b,dx>=zero)-x
3389    d=half*e
3390   end if
3391
3392   if(abs(d) >=tol1)then
3393    u=x+d
3395   else
3396    u=x+sign(tol1,d)
3398    if(fu>fx) then
3399     exit
3400    end if
3401   end if
3403   if(fu<=fx)then
3404    if(u>=x)then
3405     a=x
3406    else
3407     b=x
3408    end if
3409    vv=w
3410    fv=fw
3411    dv=dw
3412    w=x
3413    fw=fx
3414    dw=dx
3415    x=u
3416    dx=du
3417    fx=fu
3418   else
3419    if(u<x) then
3420     a=u
3421    else
3422     b=u
3423    end if
3424    ok3=(w==x).or.(fu.le.fw)
3425    ok4=(vv==w).or.(vv==x).or.(fu.lt.fv)
3426    if(ok3) then
3427     vv=w
3428     fv=fw
3429     dv=dw
3430     w=u
3431     fw=fu
3432     dw=du
3433    else if( ok4 ) then
3434     vv=u
3435     fv=fu
3436     dv=du
3437    end if
3438   end if
3439  end do
3440  xmin=x
3441 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of
3443 !but for instance renormilizing the density if brent is used on a density...
3445  brent=fx
3446
3447 end function brent


## ABINIT/cgpr [ Functions ]

[ Top ] [ Functions ]

NAME

 cgpr


FUNCTION

 perform Polak-Ribiere conjugate gradient on a function f
implementation based on the cg recipe of "numerical recipe"


INPUTS

 dp_dum_vdp: function  to be minimized (return a dp from a vector of dp)
vdp_dum_vdp: derivative of f
dtol: precision precision required for the minimization
itmax: number of iterations allowed (each linmin will be done with at max 10 times
this number


OUTPUT

 fmin: value of f at the minimum
lastdelta: absolute value of the last delta between steps


SIDE EFFECTS

 v: vector on which minimization is to be performed, starting point
and resulting min


PARENTS

      prcrskerker1,prcrskerker2


CHILDREN

      linmin


SOURCE

3021 subroutine cgpr(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,dtol,itmax,v,fmin,delta)
3022
3023
3024 !This section has been created automatically by the script Abilint (TD).
3025 !Do not modify the following lines by hand.
3026 #undef ABI_FUNC
3027 #define ABI_FUNC 'cgpr'
3028 !End of the abilint section
3029
3030  implicit none
3031
3032 !Arguments ------------------------------------
3033 include "dummy_functions.inc"
3034 !scalars
3035  integer,intent(in) :: itmax,nv1,nv2
3036  real(dp),intent(in) :: dtol
3037  real(dp),intent(out) :: delta,fmin
3038 !arrays
3039  real(dp),intent(inout) :: v(nv1,nv2)
3040
3041 !Local variables-------------------------------
3042 !scalars
3043  integer :: iiter
3044  real(dp) :: fv,gam,gscal,gscal2,sto
3045 !arrays
3047 !no_abirules
3048
3049 !************************************************************************
3050  fv = dp_dum_v2dp(nv1,nv2,v(:,:))
3054  do iiter=1,itmax
3056 ! return if the min is reached
3057   sto=dtol*(abs(fmin)+abs(fv)+tol14)
3058   delta=abs(fv-fmin)
3059   delta=abs(delta)
3060   if((delta.lt.sto).or.(iiter==itmax)) then
3061 !  DEBUG
3062 !  write(std_out,*) 'cgpr (01cg) : stop cond for cgpr:',sto,'delta:',delta,'fv:',fv,'fmin:',fmin
3063 !  ENDDEBUG
3064    return
3065   end if
3066 ! a new step
3067   fv=fmin
3072   gam=gscal2/gscal
3076 ! DEBUG
3077 ! write(std_out,*) 'cgpr (01cg) :================================================================================='
3078 ! write(std_out,*) 'cgpr (01cg) : step',iiter,'delta:',delta ,'fv',fv,'fmin',fmin
3079 ! write(std_out,*) 'cgpr (01cg) :================================================================================='
3080 ! ENDDEBUG
3081  end do
3082
3083 end subroutine cgpr


## ABINIT/dielmt [ Functions ]

[ Top ] [ Functions ]

NAME

 dielmt


FUNCTION

 Compute dielectric matrix from susceptibility matrix
Diagonalize it, then invert it.


INPUTS

  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
npwdiel=size of the dielinv and susmat arrays.
nspden=number of spin-density components
occopt=option for occupancies
prtvol=control print volume and debugging output
susmat(2,npwdiel,nspden,npwdiel,nspden)=
the susceptibility (or density-density response) matrix in reciprocal space


OUTPUT

  dielinv(2,npwdiel,(nspden+4)/3,npwdiel,(nspden+4)/3)=inverse of the (non-hermitian)
TC dielectric matrix in reciprocal space.


NOTES

 Warning : will not work in the spin-polarized, metallic case.
Output (not cleaned)
!!! Spin behaviour is not obvious !!!


TODO

 Write equation below (hermitian matrix)


PARENTS

      prcref,prcref_PMA


CHILDREN

      timab,wrtout,zhpev


SOURCE

1626 subroutine dielmt(dielinv,gmet,kg_diel,npwdiel,nspden,occopt,prtvol,susmat)
1627
1628
1629 !This section has been created automatically by the script Abilint (TD).
1630 !Do not modify the following lines by hand.
1631 #undef ABI_FUNC
1632 #define ABI_FUNC 'dielmt'
1633 !End of the abilint section
1634
1635  implicit none
1636
1637 !Arguments ------------------------------------
1638 !scalars
1639  integer,intent(in) :: npwdiel,nspden,occopt,prtvol
1640 !arrays
1641  integer,intent(in) :: kg_diel(3,npwdiel)
1642  real(dp),intent(in) :: gmet(3,3)
1643  real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
1644  real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden)
1645
1646 !Local variables-------------------------------
1647 !scalars
1648  integer :: ieig,ier,ii,index,ipw,ipw1,ipw2,isp,jj,npwsp
1649  real(dp) :: ai1,ai2,ar1,ar2,eiginv,gfact,gfactinv,gred1,gred2,gred3,gsquar
1650  real(dp) :: tpisq
1651  character(len=500) :: message
1652 !arrays
1653  real(dp) :: tsec(2)
1654  real(dp),allocatable :: dielh(:),dielmat(:,:,:,:,:),dielvec(:,:,:)
1655  real(dp),allocatable :: eig_diel(:),zhpev1(:,:),zhpev2(:)
1656 !no_abirules
1657 !integer :: ipw3
1658 !real(dp) :: elementi,elementr
1659
1660 ! *************************************************************************
1661
1662 !DEBUG
1663 !write(std_out,*)' dielmt : enter '
1664 !if(.true.)stop
1665 !ENDDEBUG
1666
1667 !tpisq is (2 Pi) **2:
1668  tpisq=(two_pi)**2
1669
1670  call timab(90,1,tsec)
1671
1672 !-Compute now the hermitian dielectric matrix------------------------------
1673 !Following remarks are only valid within RPA approximation (Kxc=0):
1674
1675 !for the spin-unpolarized case, 1 - 4pi (1/G) chi0(G,Gp) (1/Gp)
1676
1677 !for the spin-polarized case,
1678 !( 1  0 ) - 4pi ( 1/G  1/G )   ( chi0 upup  chi0 updn )   ( 1/Gp 1/Gp )
1679 !( 0  1 )       ( 1/G  1/G )   ( chi0 dnup  chi0 dndn )   ( 1/Gp 1/Gp )
1680 !which is equal to
1681 !( 1  0 ) - 4pi (1/G  0 ) (chi0 upup+dndn+updn+dnup  chi0 upup+dndn+updn+dnup) (1/Gp 0  )
1682 !( 0  1 )       ( 0  1/G) (chi0 upup+dndn+updn+dnup  chi0 upup+dndn+updn+dnup) ( 0  1/Gp)
1683 !So, if spin-polarized, sum all spin contributions
1684 !Note: chi0 updn = chi0 dnup = zero for non-metallic systems
1685
1686 !In the case of non-collinear magnetism, within RPA, this is the same because:
1687 !chi0_(s1,s2),(s3,s4) = delta_s1,s2 * delta_s3,s4 * chi0_(s1,s1),(s3,s3)
1688 !Only chi_upup,upup, chi_dndn,dndn, chi_upup,dndn and chi_dndn,upup
1689 !have to be taken into account (stored, susmat(:,ipw1,1:2,ipw2,1:2)
1690
1691  ABI_ALLOCATE(dielmat,(2,npwdiel,min(nspden,2),npwdiel,min(nspden,2)))
1692
1693  if(nspden/=1)then
1694    if (occopt<3) then
1695      do ipw2=1,npwdiel
1696        do ipw1=1,npwdiel
1697          dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2)
1698          dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2)
1699        end do
1700      end do
1701    else
1702      do ipw2=1,npwdiel
1703        do ipw1=1,npwdiel
1704          dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2)+susmat(1,ipw1,1,ipw2,2)+susmat(1,ipw1,2,ipw2,1)
1705          dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2)+susmat(2,ipw1,1,ipw2,2)+susmat(2,ipw1,2,ipw2,1)
1706        end do
1707      end do
1708    end if
1709  else
1710    do ipw2=1,npwdiel
1711      do ipw1=1,npwdiel
1712        dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)
1713        dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)
1714      end do
1715    end do
1716  end if
1717 !Compute 1/G factors and include them in the dielectric matrix
1718  do ipw1=1,npwdiel
1719    gred1=dble(kg_diel(1,ipw1))
1720    gred2=dble(kg_diel(2,ipw1))
1721    gred3=dble(kg_diel(3,ipw1))
1722    gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 &
1723 &   +two*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 +      &
1724 &   gmet(2,3)*gred2*gred3)                        )
1725 !  Distinguish G=0 from other elements
1726    if(gsquar>tol12)then
1727 !    !$gfact=\sqrt (4.0_dp \pi/gsquar/dble(nspden))$
1728      gfact=sqrt(four_pi/gsquar)
1729      do ipw2=1,npwdiel
1730 !      Must multiply both rows and columns, and also changes the sign
1731        dielmat(1,ipw2,1,ipw1,1)=-dielmat(1,ipw2,1,ipw1,1)*gfact
1732        dielmat(2,ipw2,1,ipw1,1)=-dielmat(2,ipw2,1,ipw1,1)*gfact
1733        dielmat(1,ipw1,1,ipw2,1)= dielmat(1,ipw1,1,ipw2,1)*gfact
1734        dielmat(2,ipw1,1,ipw2,1)= dielmat(2,ipw1,1,ipw2,1)*gfact
1735      end do
1736    else
1737 !    Zero the G=0 elements, head and wings
1738      do ipw2=1,npwdiel
1739        dielmat(1,ipw2,1,ipw1,1)=zero
1740        dielmat(2,ipw2,1,ipw1,1)=zero
1741        dielmat(1,ipw1,1,ipw2,1)=zero
1742        dielmat(2,ipw1,1,ipw2,1)=zero
1743      end do
1744    end if
1745  end do
1746
1747 !Complete the matrix in the spin-polarized case
1748 !should this be nspden==2??
1749  if(nspden/=1)then
1750    do ipw1=1,npwdiel
1751      do ipw2=1,npwdiel
1752        dielmat(1,ipw1,1,ipw2,2)=dielmat(1,ipw1,1,ipw2,1)
1753        dielmat(2,ipw1,1,ipw2,2)=dielmat(2,ipw1,1,ipw2,1)
1754        dielmat(1,ipw1,2,ipw2,1)=dielmat(1,ipw1,1,ipw2,1)
1755        dielmat(2,ipw1,2,ipw2,1)=dielmat(2,ipw1,1,ipw2,1)
1756        dielmat(1,ipw1,2,ipw2,2)=dielmat(1,ipw1,1,ipw2,1)
1757        dielmat(2,ipw1,2,ipw2,2)=dielmat(2,ipw1,1,ipw2,1)
1758      end do
1759    end do
1760  end if
1761
1762 !DEBUG
1763 !write(std_out,*)' dielmt : make dielmat equal to identity matrix '
1764 !do ipw1=1,npwdiel
1765 !do ipw2=1,npwdiel
1766 !dielmat(1,ipw1,1,ipw2,1)=0.0_dp
1767 !dielmat(2,ipw1,1,ipw2,1)=0.0_dp
1768 !end do
1769 !end do
1770 !ENDDEBUG
1771
1773  do isp=1,min(nspden,2)
1774    do ipw=1,npwdiel
1775      dielmat(1,ipw,isp,ipw,isp)=one+dielmat(1,ipw,isp,ipw,isp)
1776    end do
1777  end do
1778
1779 !-The hermitian dielectric matrix is computed ------------------------------
1780 !-Now, diagonalize it ------------------------------------------------------
1781
1782 !In RPA, everything is projected on the spin-symmetrized
1783 !space. This was coded here (for the time being).
1784
1785 !Diagonalize the hermitian dielectric matrix
1786
1787 !npwsp=npwdiel*nspden
1788  npwsp=npwdiel
1789
1790  ABI_ALLOCATE(dielh,(npwsp*(npwsp+1)))
1791  ABI_ALLOCATE(dielvec,(2,npwsp,npwsp))
1792  ABI_ALLOCATE(eig_diel,(npwsp))
1793  ABI_ALLOCATE(zhpev1,(2,2*npwsp-1))
1794  ABI_ALLOCATE(zhpev2,(3*npwsp-2))
1795  ier=0
1796 !Store the dielectric matrix in proper mode before calling zhpev
1797  index=1
1798  do ii=1,npwdiel
1799    do jj=1,ii
1800      dielh(index  )=dielmat(1,jj,1,ii,1)
1801      dielh(index+1)=dielmat(2,jj,1,ii,1)
1802      index=index+2
1803    end do
1804  end do
1805 !If spin-polarized and non RPA, need to store other parts of the matrix
1806 !if(nspden/=1)then
1807 !do ii=1,npwdiel
1808 !Here, spin-flip contribution
1809 !do jj=1,npwdiel
1810 !dielh(index  )=dielmat(1,jj,1,ii,2)
1811 !dielh(index+1)=dielmat(2,jj,1,ii,2)
1812 !index=index+2
1813 !end do
1814 !Here spin down-spin down upper matrix
1815 !do jj=1,ii
1816 !dielh(index  )=dielmat(1,jj,2,ii,2)
1817 !dielh(index+1)=dielmat(2,jj,2,ii,2)
1818 !index=index+2
1819 !end do
1820 !end do
1821 !end if
1822
1823  call ZHPEV ('V','U',npwsp,dielh,eig_diel,dielvec,npwdiel,zhpev1,&
1824 & zhpev2,ier)
1825  ABI_DEALLOCATE(zhpev1)
1826  ABI_DEALLOCATE(zhpev2)
1827
1828  if(prtvol>=10)then
1829    write(message, '(a,a,a,5es12.4)' )ch10,&
1830 &   ' Five largest eigenvalues of the hermitian RPA dielectric matrix:',&
1831 &   ch10,eig_diel(npwdiel:npwdiel-4:-1)
1832    call wrtout(ab_out,message,'COLL')
1833  end if
1834
1835  write(message, '(a,a)' )ch10,&
1836 & ' dielmt : 15 largest eigenvalues of the hermitian RPA dielectric matrix'
1837  call wrtout(std_out,message,'COLL')
1838  write(message, '(a,5es12.5)' )'  1-5  :',eig_diel(npwdiel:npwdiel-4:-1)
1839  call wrtout(std_out,message,'COLL')
1840  write(message, '(a,5es12.5)' )'  6-10 :',eig_diel(npwdiel-5:npwdiel-9:-1)
1841  call wrtout(std_out,message,'COLL')
1842  write(message, '(a,5es12.5)' )'  11-15:',eig_diel(npwdiel-10:npwdiel-14:-1)
1843  call wrtout(std_out,message,'COLL')
1844  write(message, '(a,a)' )ch10,&
1845 & ' dielmt : 5 smallest eigenvalues of the hermitian RPA dielectric matrix'
1846  call wrtout(std_out,message,'COLL')
1847  write(message, '(a,5es12.5)' )'  1-5  :',eig_diel(1:5)
1848  call wrtout(std_out,message,'COLL')
1849
1850 !Invert the hermitian dielectric matrix,
1851 !Should use a BLAS call !
1852  do ipw2=1,npwdiel
1853    do ipw1=ipw2,npwdiel
1854      dielinv(1,ipw1,1,ipw2,1)=zero
1855      dielinv(2,ipw1,1,ipw2,1)=zero
1856    end do
1857  end do
1858  do ieig=1,npwdiel
1859    eiginv=one/eig_diel(ieig)
1860    do ipw2=1,npwdiel
1861      do ipw1=ipw2,npwdiel
1862        ar1=dielvec(1,ipw1,ieig)
1863        ai1=dielvec(2,ipw1,ieig)
1864        ar2=dielvec(1,ipw2,ieig)
1865        ai2=dielvec(2,ipw2,ieig)
1866        dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)+&
1867 &       (ar1*ar2+ai1*ai2)*eiginv
1868        dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)+&
1869 &       (ai1*ar2-ar1*ai2)*eiginv
1870      end do
1871    end do
1872  end do
1873  do ipw2=1,npwdiel-1
1874    do ipw1=ipw2+1,npwdiel
1875      dielinv(1,ipw2,1,ipw1,1)= dielinv(1,ipw1,1,ipw2,1)
1876      dielinv(2,ipw2,1,ipw1,1)=-dielinv(2,ipw1,1,ipw2,1)
1877    end do
1878  end do
1879
1880  ABI_DEALLOCATE(dielh)
1881  ABI_DEALLOCATE(dielvec)
1882  ABI_DEALLOCATE(eig_diel)
1883
1884 !DEBUG
1885 !Checks whether the inverse of the hermitian dielectric matrix
1886 !has been correctly generated
1887 !do ipw1=1,npwdiel
1888 !do ipw2=1,npwdiel
1889 !elementr=0.0_dp
1890 !elementi=0.0_dp
1891 !do ipw3=1,npwdiel
1892 !elementr=elementr+dielinv(1,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1)&
1893 !&                    -dielinv(2,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1)
1894 !elementi=elementi+dielinv(1,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1)&
1895 !&                    +dielinv(2,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1)
1896 !end do
1897 !if(elementr**2+elementi**2 > 1.0d-12)then
1898 !if( ipw1 /= ipw2 .or. &
1899 !&        ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then
1900 !write(std_out,*)' dielmt : the inversion procedure is not correct '
1901 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2
1902 !write(std_out,*)' elementr,elementi=',elementr,elementi
1903 !stop
1904 !end if
1905 !end if
1906 !end do
1907 !end do
1908 !write(std_out,*)'dielmt : matrix has been inverted successfully '
1909 !stop
1910 !ENDDEBUG
1911
1912 !Then get the inverse of the asymmetric
1913 !dielectric matrix, as required for the preconditioning.
1914
1915 !Inverse of the dielectric matrix : ( 1 - 4pi (1/G^2) chi0(G,Gp) )^(-1)
1916 !In dielinv there is now (1 - 4pi (1/G) chi0(G,Gp) (1/Gp) )^(-1)
1917 !So, evaluate dielinv_after(G,Gp) =
1918 !(4pi/G^2)^(1/2) dielinv_before(G,Gp) (4pi/Gp^2)^(-1/2)
1919 !In RPA, can focus on the spin-averaged quantities
1920  do ipw1=1,npwdiel
1921    gred1=dble(kg_diel(1,ipw1))
1922    gred2=dble(kg_diel(2,ipw1))
1923    gred3=dble(kg_diel(3,ipw1))
1924    gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 &
1925 &   +two*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 +      &
1926 &   gmet(2,3)*gred2*gred3)                        )
1927 !  Distinguish G=0 from other elements
1928    if(gsquar>tol12)then
1929      gfact=sqrt(four_pi/gsquar)
1930      gfactinv=one/gfact
1931      do ipw2=1,npwdiel
1932 !      Must multiply both rows and columns
1933        dielinv(1,ipw2,1,ipw1,1)=dielinv(1,ipw2,1,ipw1,1)*gfactinv
1934        dielinv(2,ipw2,1,ipw1,1)=dielinv(2,ipw2,1,ipw1,1)*gfactinv
1935        dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)*gfact
1936        dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)*gfact
1937      end do
1938    else
1939 !    Zero the G=0 elements, head
1940      do ipw2=1,npwdiel
1941        if (ipw2/=ipw1) dielinv(1:2,ipw1,1,ipw2,1)=zero
1942      end do
1943    end if
1944  end do
1945
1946  ABI_DEALLOCATE(dielmat)
1947
1948  call timab(90,2,tsec)
1949
1950 end subroutine dielmt


## ABINIT/dieltcel [ Functions ]

[ Top ] [ Functions ]

NAME

 dieltcel


FUNCTION

 Compute either test charge or electronic dielectric matrices
from susceptibility matrix
Diagonalize it, then invert it.


INPUTS

  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
kxc(nfft,nkxc)=exchange-correlation kernel,
needed if the electronic dielectric matrix is computed
nfft=(effective) number of FFT grid points (for this processor)
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
nkxc=second dimension of the array kxc, see rhotoxc.f for a description
npwdiel=size of the dielinv and susmat arrays.
nspden=number of spin-density components
occopt=option for occupancies
option=1 for Test Charge dielectric matrix, 2 for electronic dielectric matrix
prtvol=control print volume and debugging output
susmat(2,npwdiel,nspden,npwdiel,nspden)=
the susceptibility (or density-density response) matrix in reciprocal space


OUTPUT

  dielinv(2,npwdiel,nspden,npwdiel,nspden)=inverse of the (non-hermitian)
TC dielectric matrix in reciprocal space.


NOTES

 Output (not cleaned)
!!! Spin behaviour is not obvious !!!
Will not work in the spin-polarized, metallic case.


PARENTS

      prcref,prcref_PMA


CHILDREN

      destroy_mpi_enreg,fourdp,init_distribfft_seq,initmpi_seq,timab,wrtout
zhpev


SOURCE

1997 subroutine dieltcel(dielinv,gmet,kg_diel,kxc,&
1998 &  nfft,ngfft,nkxc,npwdiel,nspden,occopt,option,paral_kgb,prtvol,susmat)
1999
2000
2001 !This section has been created automatically by the script Abilint (TD).
2002 !Do not modify the following lines by hand.
2003 #undef ABI_FUNC
2004 #define ABI_FUNC 'dieltcel'
2005 !End of the abilint section
2006
2007  implicit none
2008
2009 !Arguments ------------------------------------
2010 !scalars
2011  integer,intent(in) :: nfft,nkxc,npwdiel,nspden,occopt,option,paral_kgb
2012  integer,intent(in) :: prtvol
2013 !arrays
2014  integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18)
2015  real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc)
2016  real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
2017  real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden)
2018
2019 !Local variables-------------------------------
2020 !scalars
2021  integer :: i1,i2,i3,ieig,ier,ifft,ii,index,ipw0,ipw1,ipw2,ispden,j1
2022  integer :: j2,j3,jj,k1,k2,k3,n1,n2,n3
2023  real(dp) :: ai,ai2,ar,ar2,eiginv,gred1,gred2,gred3,gsquar,si
2024  real(dp) :: sr,tpisq
2025  character(len=500) :: message
2026  type(MPI_type) :: mpi_enreg_seq
2027 !arrays
2028  real(dp) :: tsec(2)
2029  real(dp),allocatable :: eig_msusinvsqr(:),eig_msussqr(:)
2030  real(dp),allocatable :: eig_sus(:),eig_sym(:),invsqrsus(:,:,:,:,:)
2031  real(dp),allocatable :: khxc(:,:,:,:,:),kxcg(:,:),sqrsus(:,:,:,:,:),sush(:)
2032  real(dp),allocatable :: susvec(:,:,:),symdielmat(:,:,:,:,:),symh(:)
2033  real(dp),allocatable :: symvec(:,:,:,:,:),wkxc(:),work(:,:,:,:,:)
2034  real(dp),allocatable :: work2(:,:,:,:,:),zhpev1(:,:),zhpev2(:)
2035 !no_abirules
2036 !integer :: ipw3
2037 !real(dp) :: elementi,elementr
2038 !DEBUG
2039 !Used to moderate divergence effect near rho=0 of the Kxc
2040 !this limit value is truly empirical (exprmt on small Sr cell).
2041 !real(dp) :: kxc_min=-200.0
2042 !ENDDEBUG
2043
2044 ! *************************************************************************
2045
2046  call timab(96,1,tsec)
2047
2048 !tpisq is (2 Pi) **2:
2049  tpisq=(two_pi)**2
2050
2051  if(nspden/=1 .and. (occopt>=3 .and. occopt<=8) )then
2052    write(message, '(a,a,a)' )&
2053 &   'In the present version of the code, one cannot produce',ch10,&
2054 &   'the dielectric matrix in the metallic, spin-polarized case.'
2055    MSG_BUG(message)
2056  end if
2057
2058  if(nspden==4)then
2059    write(message,'(a,a,a)')&
2060 &   'In the present version of the code, one cannot produce',ch10,&
2061 &   'the dielectric matrix in the non-collinear spin-polarized case.'
2062    MSG_ERROR(message)
2063  end if
2064
2065
2066 !-Diagonalize the susceptibility matrix
2067
2068  ABI_ALLOCATE(sush,(npwdiel*(npwdiel+1)))
2069  ABI_ALLOCATE(susvec,(2,npwdiel,npwdiel))
2070  ABI_ALLOCATE(eig_msusinvsqr,(npwdiel))
2071  ABI_ALLOCATE(eig_msussqr,(npwdiel))
2072  ABI_ALLOCATE(eig_sus,(npwdiel))
2073  ABI_ALLOCATE(zhpev1,(2,2*npwdiel-1))
2074  ABI_ALLOCATE(zhpev2,(3*npwdiel-2))
2075  ABI_ALLOCATE(work,(2,npwdiel,nspden,npwdiel,nspden))
2076  ABI_ALLOCATE(work2,(2,npwdiel,nspden,npwdiel,nspden))
2077  ABI_ALLOCATE(sqrsus,(2,npwdiel,nspden,npwdiel,nspden))
2078  ABI_ALLOCATE(invsqrsus,(2,npwdiel,nspden,npwdiel,nspden))
2079
2080 !At some time, should take care of different spin channels
2081  do ispden=1,nspden
2082
2083    if(nspden/=1)then
2084      MSG_ERROR('dieltcel : stop, nspden/=1')
2085    end if
2086
2087 !  Store the susceptibility matrix in proper mode before calling zhpev
2088    index=1
2089    do ii=1,npwdiel
2090      do jj=1,ii
2091        sush(index  )=susmat(1,jj,1,ii,1)
2092        sush(index+1)=susmat(2,jj,1,ii,1)
2093        index=index+2
2094      end do
2095    end do
2096
2097    ier=0
2098    call ZHPEV ('V','U',npwdiel,sush,eig_sus,susvec,npwdiel,zhpev1,zhpev2,ier)
2099
2100 !  DEBUG
2101 !  write(std_out,*)' dieltcel : print eigenvalues of the susceptibility matrix'
2102 !  do ii=1,npwdiel
2103 !  write(std_out,'(i5,es16.6)' )ii,eig_sus(ii)
2104 !  end do
2105 !  ENDDEBUG
2106
2107    do ii=1,npwdiel
2108      if(-eig_sus(ii)>1.d-12)then
2109        eig_msussqr(ii)=sqrt(-eig_sus(ii))
2110        eig_msusinvsqr(ii)=1._dp/eig_msussqr(ii)
2111      else if(-eig_sus(ii)< -1.d-12)then
2112        message = "Found positive eigenvalue of susceptibility matrix."
2113        MSG_BUG(message)
2114      else
2115 !      Set the eigenvalue corresponding to a constant potential change to 1,
2116 !      while it will be set to zero in Khx.
2117        eig_msussqr(ii)=1._dp
2118        eig_msusinvsqr(ii)=1._dp
2119      end if
2120    end do
2121
2122 !  Compute square root of minus susceptibility matrix
2123 !  and inverse square root of minus susceptibility matrix
2124    do ii=1,npwdiel
2125      work(:,:,1,ii,1)=susvec(:,:,ii)*eig_msussqr(ii)
2126      work2(:,:,1,ii,1)=susvec(:,:,ii)*eig_msusinvsqr(ii)
2127    end do
2128    do ipw2=1,npwdiel
2129      do ipw1=ipw2,npwdiel
2130        ar=0._dp ; ai=0._dp ; ar2=0._dp ; ai2=0._dp
2131        do ii=1,npwdiel
2132          sr=susvec(1,ipw2,ii) ; si=susvec(2,ipw2,ii)
2133          ar =ar  +work(1,ipw1,1,ii,1)*sr  +work(2,ipw1,1,ii,1)*si
2134          ai =ai  +work(2,ipw1,1,ii,1)*sr  -work(1,ipw1,1,ii,1)*si
2135          ar2=ar2 +work2(1,ipw1,1,ii,1)*sr +work2(2,ipw1,1,ii,1)*si
2136          ai2=ai2 +work2(2,ipw1,1,ii,1)*sr -work2(1,ipw1,1,ii,1)*si
2137        end do
2138        sqrsus(1,ipw1,1,ipw2,1)=ar
2139        sqrsus(2,ipw1,1,ipw2,1)=ai
2140        invsqrsus(1,ipw1,1,ipw2,1)=ar2
2141        invsqrsus(2,ipw1,1,ipw2,1)=ai2
2142        if(ipw1/=ipw2)then
2143          sqrsus(1,ipw2,1,ipw1,1)=ar
2144          sqrsus(2,ipw2,1,ipw1,1)=-ai
2145          invsqrsus(1,ipw2,1,ipw1,1)=ar2
2146          invsqrsus(2,ipw2,1,ipw1,1)=-ai2
2147        end if
2148      end do
2149    end do
2150
2151 !  DEBUG
2152 !  Checks whether sqrsus and invsqrsus are inverse of each other.
2153 !  do ipw1=1,npwdiel
2154 !  do ipw2=1,npwdiel
2155 !  elementr=0.0_dp
2156 !  elementi=0.0_dp
2157 !  do ipw3=1,npwdiel
2158 !  elementr=elementr+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1)&
2159 !  &                    -sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1)
2160 !  elementi=elementi+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1)&
2161 !  &                    +sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1)
2162 !  end do
2163 !  if(elementr**2+elementi**2 > 1.0d-12)then
2164 !  if( ipw1 /= ipw2 .or. &
2165 !  &        ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then
2166 !  write(std_out,*)' dieltcel : sqrsus and invsqrsus are not (pseudo)',&
2167 !  &        'inverse of each other'
2168 !  write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2
2169 !  write(std_out,*)' elementr,elementi=',elementr,elementi
2170 !  stop
2171 !  end if
2172 !  end if
2173 !  end do
2174 !  end do
2175 !  ENDDEBUG
2176
2177 !  End loop over spins
2178  end do
2179
2180  ABI_DEALLOCATE(eig_msusinvsqr)
2181  ABI_DEALLOCATE(eig_msussqr)
2182  ABI_DEALLOCATE(eig_sus)
2183  ABI_DEALLOCATE(sush)
2184  ABI_DEALLOCATE(susvec)
2185
2186 !-Compute the Hxc kernel
2187
2188  ABI_ALLOCATE(khxc,(2,npwdiel,nspden,npwdiel,nspden))
2189  ABI_ALLOCATE(symdielmat,(2,npwdiel,nspden,npwdiel,nspden))
2190
2191  khxc(:,:,:,:,:)=0.0_dp
2192
2193 !Compute Hartree kernel
2194  do ipw1=1,npwdiel
2195    gred1=dble(kg_diel(1,ipw1))
2196    gred2=dble(kg_diel(2,ipw1))
2197    gred3=dble(kg_diel(3,ipw1))
2198    gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 &
2199 &   +2.0_dp*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 +      &
2200 &   gmet(2,3)*gred2*gred3)                        )
2201 !  Distinguish G=0 from other elements
2202    if(gsquar>1.0d-12)then
2203      khxc(1,ipw1,1,ipw1,1)= 4.0_dp*pi/gsquar
2204    else
2205 !    G=0
2206      ipw0=ipw1
2207    end if
2208  end do
2209
2210 !Eventually add the xc part
2211  if(option>=2)then
2212
2213    ABI_ALLOCATE(wkxc,(nfft))
2214    ABI_ALLOCATE(kxcg,(2,nfft))
2215    wkxc(:)=kxc(:,1)
2216 !  DEBUG
2217 !  Used to moderate divergenc effect near rho=0 of the Kxc (see above).
2218 !  wkxc(:)=merge(kxc(:,1), kxc_min, kxc(:,1) > kxc_min)
2219 !  ENDDEBUG
2220    call initmpi_seq(mpi_enreg_seq)
2221    call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all')
2222    call fourdp(1,kxcg,wkxc,-1,mpi_enreg_seq,nfft,ngfft,paral_kgb,0) ! trsfrm R to G
2223    call destroy_mpi_enreg(mpi_enreg_seq)
2224
2225 !  Compute difference in G vectors
2226    n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
2227    do ipw2=1,npwdiel
2228      if(ipw2/=ipw0)then
2229
2230        j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2)
2231 !      Fills diagonal
2232        khxc(1,ipw2,1,ipw2,1)=khxc(1,ipw2,1,ipw2,1)+kxcg(1,1)
2233        khxc(2,ipw2,1,ipw2,1)=khxc(2,ipw2,1,ipw2,1)+kxcg(2,1)
2234
2235        if(ipw2/=npwdiel)then
2236 !        Fills off-diagonal part of the matrix, except G=0
2237          do ipw1=ipw2+1,npwdiel
2238            if(ipw1/=ipw0)then
2239              i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1)
2240 !            Use of two mod calls handles both i1-j1>=ndiel1 AND i1-j1<0
2241              k1=mod(n1+mod(i1-j1,n1),n1)
2242              k2=mod(n2+mod(i2-j2,n2),n2)
2243              k3=mod(n3+mod(i3-j3,n3),n3)
2244              ifft=k1+1+n1*(k2+n2*k3)
2245 !            The signs of imaginary contributions have been checked
2246              khxc(1,ipw1,1,ipw2,1)=kxcg(1,ifft)
2247              khxc(2,ipw1,1,ipw2,1)=kxcg(2,ifft)
2248              khxc(1,ipw2,1,ipw1,1)=kxcg(1,ifft)
2249              khxc(2,ipw2,1,ipw1,1)=-kxcg(2,ifft)
2250            end if
2251          end do
2252        end if
2253
2254      end if
2255    end do
2256
2257    ABI_DEALLOCATE(wkxc)
2258    ABI_DEALLOCATE(kxcg)
2259
2260 !  Endif option 2
2261  end if
2262
2263 !Now, get the symmetric dielectric matrix
2264 !Premultiplication by square root of minus susceptibility matrix
2265  do ipw2=1,npwdiel
2266    do ipw1=1,npwdiel
2267      ar=0._dp ; ai=0._dp
2268      do ii=1,npwdiel
2269        ar=ar+sqrsus(1,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) &
2270 &       -sqrsus(2,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1)
2271        ai=ai+sqrsus(2,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) &
2272 &       +sqrsus(1,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1)
2273      end do
2274      work(1,ipw1,1,ipw2,1)=ar
2275      work(2,ipw1,1,ipw2,1)=ai
2276    end do
2277  end do
2278 !Postmultiplication by square root of minus susceptibility matrix
2279  do ipw2=1,npwdiel
2280 !  do ipw1=ipw2,npwdiel
2281    do ipw1=1,npwdiel
2282      ar=0._dp ; ai=0._dp
2283      do ii=1,npwdiel
2284        ar=ar+work(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2285 &       -work(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2286        ai=ai+work(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2287 &       +work(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2288      end do
2289      symdielmat(1,ipw1,1,ipw2,1)=ar
2290      symdielmat(2,ipw1,1,ipw2,1)=ai
2291 !    if(ipw1/=ipw2)then
2292 !    symdielmat(1,ipw2,1,ipw1,1)=ar
2293 !    symdielmat(2,ipw2,1,ipw1,1)=-ai
2294 !    end if
2295    end do
2296 !  Add the unity matrix
2297    symdielmat(1,ipw2,1,ipw2,1)=1._dp+symdielmat(1,ipw2,1,ipw2,1)
2298  end do
2299
2300  ABI_DEALLOCATE(khxc)
2301
2302  ABI_ALLOCATE(symh,(npwdiel*(npwdiel+1)))
2303  ABI_ALLOCATE(symvec,(2,npwdiel,nspden,npwdiel,nspden))
2304  ABI_ALLOCATE(eig_sym,(npwdiel))
2305
2306 !Store the symmetrized dielectric matrix in proper mode before calling zhpev
2307  index=1
2308  do ii=1,npwdiel
2309    do jj=1,ii
2310      symh(index  )=symdielmat(1,jj,1,ii,1)
2311      symh(index+1)=symdielmat(2,jj,1,ii,1)
2312      index=index+2
2313    end do
2314  end do
2315
2316  ier=0
2317  call ZHPEV ('V','U',npwdiel,symh,eig_sym,symvec,npwdiel,zhpev1,&
2318 & zhpev2,ier)
2319
2320  if(prtvol>=10)then
2321    write(message, '(a,a,a,5es12.4)' )ch10,&
2322 &   ' Five largest eigenvalues of the symmetrized dielectric matrix:',&
2323 &   ch10,eig_sym(npwdiel:npwdiel-4:-1)
2324    call wrtout(ab_out,message,'COLL')
2325  end if
2326
2327  write(message, '(a,a)' )ch10,&
2328 & ' dieltcel : 15 largest eigenvalues of the symmetrized dielectric matrix'
2329  call wrtout(std_out,message,'COLL')
2330  write(message, '(a,5es12.5)' )'  1-5  :',eig_sym(npwdiel:npwdiel-4:-1)
2331  call wrtout(std_out,message,'COLL')
2332  write(message, '(a,5es12.5)' )'  6-10 :',eig_sym(npwdiel-5:npwdiel-9:-1)
2333  call wrtout(std_out,message,'COLL')
2334  write(message, '(a,5es12.5)' )'  11-15:',eig_sym(npwdiel-10:npwdiel-14:-1)
2335  call wrtout(std_out,message,'COLL')
2336  write(message, '(a,a)' )ch10,&
2337 & ' dieltcel : 5 smallest eigenvalues of the symmetrized dielectric matrix'
2338  call wrtout(std_out,message,'COLL')
2339  write(message, '(a,5es12.5)' )'  1-5  :',eig_sym(1:5)
2340  call wrtout(std_out,message,'COLL')
2341
2342 !Invert the hermitian dielectric matrix,
2343  work(:,:,:,:,:)=0.0_dp
2344  do ieig=1,npwdiel
2345    eiginv=1.0_dp/eig_sym(ieig)
2346    do ipw2=1,npwdiel
2347 !    do ipw1=ipw2,npwdiel
2348      do ipw1=1,npwdiel
2349        work(1,ipw1,1,ipw2,1)=work(1,ipw1,1,ipw2,1)+&
2350 &       (symvec(1,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)+ &
2351 &       symvec(2,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv
2352        work(2,ipw1,1,ipw2,1)=work(2,ipw1,1,ipw2,1)+&
2353 &       (symvec(2,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)- &
2354 &       symvec(1,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv
2355      end do
2356    end do
2357  end do
2358 !if(npwdiel>1)then
2359 !do ipw2=2,npwdiel
2360 !do ipw1=1,ipw2-1
2361 !work(1,ipw1,1,ipw2,1)= work(1,ipw2,1,ipw1,1)
2362 !work(2,ipw1,1,ipw2,1)=-work(2,ipw2,1,ipw1,1)
2363 !end do
2364 !end do
2365 !end if
2366
2367  ABI_DEALLOCATE(eig_sym)
2368  ABI_DEALLOCATE(symh)
2369  ABI_DEALLOCATE(symvec)
2370
2371 !DEBUG
2372 !Checks whether the inverse of the symmetric dielectric matrix
2373 !has been correctly generated
2374 !do ipw1=1,npwdiel
2375 !do ipw2=1,npwdiel
2376 !elementr=0.0_dp
2377 !elementi=0.0_dp
2378 !do ipw3=1,npwdiel
2379 !elementr=elementr+work(1,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1)&
2380 !&                    -work(2,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1)
2381 !elementi=elementi+work(1,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1)&
2382 !&                    +work(2,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1)
2383 !end do
2384 !if(elementr**2+elementi**2 > 1.0d-12)then
2385 !if( ipw1 /= ipw2 .or. &
2386 !&        ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then
2387 !write(std_out,*)' dieltcel : the inversion procedure is not correct '
2388 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2
2389 !write(std_out,*)' elementr,elementi=',elementr,elementi
2390 !stop
2391 !end if
2392 !end if
2393 !end do
2394 !end do
2395 !write(std_out,*)'dieltcel : matrix has been inverted successfully '
2396 !ENDDEBUG
2397
2398  ABI_DEALLOCATE(symdielmat)
2399
2400 !Then get the inverse of the asymmetric
2401 !dielectric matrix, as required for the preconditioning.
2402 !Premultiplication by square root of minus susceptibility matrix
2403  do ipw2=1,npwdiel
2404    do ipw1=1,npwdiel
2405      ar=0._dp ; ai=0._dp
2406      do ii=1,npwdiel
2407        ar=ar+invsqrsus(1,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) &
2408 &       -invsqrsus(2,ipw1,1,ii,1)*work(2,ii,1,ipw2,1)
2409        ai=ai+invsqrsus(2,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) &
2410 &       +invsqrsus(1,ipw1,1,ii,1)*work(2,ii,1,ipw2,1)
2411      end do
2412      work2(1,ipw1,1,ipw2,1)=ar
2413      work2(2,ipw1,1,ipw2,1)=ai
2414    end do
2415  end do
2416 !Postmultiplication by square root of minus susceptibility matrix
2417  do ipw2=1,npwdiel
2418    do ipw1=1,npwdiel
2419      ar=0._dp ; ai=0._dp
2420      do ii=1,npwdiel
2421        ar=ar+work2(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2422 &       -work2(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2423        ai=ai+work2(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2424 &       +work2(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2425      end do
2426      dielinv(1,ipw1,1,ipw2,1)=ar
2427      dielinv(2,ipw1,1,ipw2,1)=ai
2428    end do
2429  end do
2430
2431  ABI_DEALLOCATE(invsqrsus)
2432  ABI_DEALLOCATE(sqrsus)
2433  ABI_DEALLOCATE(work)
2434  ABI_DEALLOCATE(work2)
2435  ABI_DEALLOCATE(zhpev1)
2436  ABI_DEALLOCATE(zhpev2)
2437
2438  call timab(96,2,tsec)
2439
2440 end subroutine dieltcel


## ABINIT/linmin [ Functions ]

[ Top ] [ Functions ]

NAME

 linmin


FUNCTION

 minimizes a function along a gradient line:
first bracket the minimum then perform the minimization


 Copyright (C) 1998-2018 ABINIT group (DCA, XG, MT)
GNU General Public License, see ~ABINIT/COPYING
or http://www.gnu.org/copyleft/gpl.txt .
For the initials of contributors, see ~ABINIT/Infos/contributors .


INPUTS

 dp_dum_vdp: function  to be minimized (return a dp from a vector of dp)
vdp_dum_vdp: derivative of f


OUTPUT

 fmin: minimun value reached for dp_dum_vdp


SIDE EFFECTS

 grad: the gradient line along which the minimization is performed (not changed)
v: the starting and then ending point of the minimization


PARENTS

      cgpr


CHILDREN

      bracketing


SOURCE

3120 subroutine linmin(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,v,grad,fmin)
3121
3122
3123 !This section has been created automatically by the script Abilint (TD).
3124 !Do not modify the following lines by hand.
3125 #undef ABI_FUNC
3126 #define ABI_FUNC 'linmin'
3127 !End of the abilint section
3128
3129  implicit none
3130
3131 !Arguments ------------------------------------
3132 include "dummy_functions.inc"
3133 !scalars
3134  integer,intent(in) :: nv1,nv2
3135  real(dp),intent(out) :: fmin
3136 !arrays
3138
3139 !Local variables-------------------------------
3140 !scalars
3141  real(dp),parameter :: maglimit=10000.0_dp,tol=tol8*tol8*tol8
3142  real(dp) :: a,b,fa,fb,fx,x,xmin
3143 !no_abirules
3144
3145 !************************************************************************
3146  a=zero
3147  x=ninth*real(1e-4,dp)
3149 !DEBUG
3150 !write(std_out,*) 'linmin (01cg) : linmin after bracketing'
3151 !write(std_out,*) 'linmin (01cg) : point',a,x,b,'value',fa,fx,fb
3152 !ENDDEBUG
3154
3155 end subroutine linmin


## ABINIT/m_prcref [ Modules ]

[ Top ] [ Modules ]

NAME

  m_prcref


FUNCTION

  Routines to precondition residual potential (or density) and forces.


  Copyright (C) 1998-2018 ABINIT group (DCA, XG, MT, PMA)
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_prcref
28
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_errors
34  use m_abicore
35  use m_xmpi
36  use m_xcdata
37
38  use m_frskerker1
39  use m_frskerker2
40  use mod_prc_memory
41
42  use m_time,     only : timab
43  use m_numeric_tools, only : dotproduct
44  use m_geometry, only : xcart2xred, metric
45  use m_cgtools,  only : dotprod_vn, mean_fftr
46  use m_mpinfo,   only : ptabs_fourdp, destroy_mpi_enreg, initmpi_seq
47  use m_pawtab,   only : pawtab_type
48  use m_pawrhoij, only : pawrhoij_type
49  use m_fftcore,  only : kgindex
50  use m_fft,      only : zerosym, indirect_parallel_fourier, fourdp
51  use m_kg,       only : getph
52  use m_dtset,    only : testsusmat
53  use m_spacepar, only : hartre, laplacian
54  use m_distribfft, only : init_distribfft_seq
55  use m_forces,     only : fresid
56  use m_atm2fft,    only : atm2fft
57  use m_rhotoxc,    only : rhotoxc
58  use m_mklocl,     only : mklocl
59  use m_mkcore,     only : mkcore
60
61
62  implicit none
63
64  private


## ABINIT/moddiel [ Functions ]

[ Top ] [ Functions ]

NAME

 moddiel


FUNCTION

 Precondition the residual, using a model dielectric function.
When cplex=1, assume q=(0 0 0), and vresid and vrespc will be REAL
When cplex=2, q must be taken into account, and vresid and vrespc will be COMPLEX


INPUTS

  cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX
dielar(7)=input parameters for dielectric matrix:
diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
nfft=(effective) number of FFT grid points (for this processor)
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
nspden=number of spin-density components
optreal=1 if residual potential is in REAL space, 2 if it is in RECIPROCAL SPACE
optres=0: the array vresid contains a potential residual
1: the array vresid contains a density residual
qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
rprimd(3,3)=dimensional primitive translations in real space (bohr)
vresid(cplex*nfft,nspden)=residual density/potential in REAL space       (if optreal==1)
residual density/potential in RECIPROCAL space (if optreal==2)


OUTPUT

  vrespc(cplex*nfft,nspden)=preconditioned residual of the density/potential
in REAL space       if optreal==1
in RECIPROCAL space if optreal==2


SIDE EFFECTS

NOTES

 optreal==2 is not compatible with cplex==1


PARENTS

      dfpt_newvtr,prcref,prcref_PMA


CHILDREN

      fourdp,metric,ptabs_fourdp


SOURCE

1379 subroutine moddiel(cplex,dielar,mpi_enreg,nfft,ngfft,nspden,optreal,optres,paral_kgb,qphon,rprimd,vresid,vrespc)
1380
1381
1382 !This section has been created automatically by the script Abilint (TD).
1383 !Do not modify the following lines by hand.
1384 #undef ABI_FUNC
1385 #define ABI_FUNC 'moddiel'
1386 !End of the abilint section
1387
1388  implicit none
1389
1390 !Arguments-------------------------------
1391 !scalars
1392  integer,intent(in) :: cplex,nfft,nspden,optreal,optres,paral_kgb
1393  type(MPI_type),intent(in) :: mpi_enreg
1394 !arrays
1395  integer,intent(in) :: ngfft(18)
1396  real(dp),intent(in) :: dielar(7),qphon(3),rprimd(3,3)
1397  real(dp),intent(in) :: vresid(cplex*nfft,nspden)
1398  real(dp),intent(out) :: vrespc(cplex*nfft,nspden)
1399
1400 !Local variables-------------------------------
1401 !scalars
1402  integer,parameter :: im=2,re=1
1403  integer :: i1,i2,i23,i3,ifft,ig,ii,ii1,ing,ispden,me_fft,mg,n1,n2,n3,nproc_fft
1404  integer :: nspden_eff,qeq0
1405  logical :: magn_precon
1406  real(dp) :: dielng,diemac,diemac_inv,diemix,diemixmag,diemix_eff,factor,gqg2p3,gqgm12,gqgm13
1407  real(dp) :: gqgm23,gs,gs2,gs3,l2g2,length2,ucvol
1408  character(len=500) :: message
1409 !arrays
1410  integer :: id(3)
1411  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1412  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1413  real(dp) :: gmet(3,3),gprimd(3,3),potg0(4),rmet(3,3)
1414  real(dp),allocatable :: gq(:,:),work1(:,:),work2(:)
1415
1416 ! *************************************************************************
1417
1418 !Check that cplex has an allowed value
1419  if(cplex/=1 .and. cplex/=2)then
1420    write(message,'(a,i0,a,a)')&
1421 &   '  From the calling routine, cplex=',cplex,ch10,&
1422 &   '  but the only value allowed are 1 and 2.'
1423    MSG_BUG(message)
1424  end if
1425
1426  if(cplex==1.and.optreal==2)then
1427    MSG_BUG('When optreal=2, cplex must be 2.')
1428  end if
1429
1430 !This is to allow q=0
1431  qeq0=0
1432  if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
1433
1434 !If cplex=1 then qphon should be 0 0 0
1435  if (cplex==1.and. qeq0/=1) then
1436    write(message,'(a,3e12.4,a)' )' cplex=1 but qphon=',qphon,' qphon should be 0 0 0.'
1437    MSG_BUG(message)
1438  end if
1439
1440  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1441  me_fft=ngfft(11)
1442  nproc_fft=ngfft(10)
1443
1444 !Get the distrib associated with this fft_grid
1445  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1446
1447 !Compute different geometric tensor, as well as ucvol, from rprimd
1448  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1449
1450  dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ; diemixmag=dielar(7)
1451
1452  magn_precon=(diemixmag>=zero) ! Set to true if magnetization has to be preconditionned
1453  diemixmag=abs(diemixmag)
1454
1455 !DEBUG
1456 !write(std_out,*)' moddiel : diemac, diemix, diemixmag =',diemac,diemix,diemixmag
1457 !ENDDEBUG
1458
1459  if(abs(diemac-1.0_dp)<1.0d-6)then
1460
1461 !  Here, simple mixing is required, through macroscopic
1462 !  dielectric constant set to 1.0_dp .
1463    vrespc(:,1)=diemix*vresid(:,1)
1464    if (nspden/=1) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden)
1465  else
1466
1467 !  Magnetization is not preconditionned
1468    if (optres==1.and.nspden>1.and.(.not.magn_precon)) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden)
1469
1470 !  Here, model dielectric function (G-diagonal operator)
1471
1472    length2=(two_pi*dielng)**2
1473    diemac_inv=1.0_dp/diemac
1474    ABI_ALLOCATE(work1,(2,nfft))
1475    if (optreal==1) then
1476      ABI_ALLOCATE(work2,(cplex*nfft))
1477    end if
1478
1479 !  In order to speed the routine, precompute the components of g
1480    mg=maxval(ngfft)
1481    ABI_ALLOCATE(gq,(3,mg))
1482    do ii=1,3
1483      id(ii)=ngfft(ii)/2+2
1484      do ing=1,ngfft(ii)
1485        ig=ing-(ing/id(ii))*ngfft(ii)-1
1486        gq(ii,ing)=ig+qphon(ii)
1487      end do
1488    end do
1489
1490 !  Do-loop on spins
1491 !  Note XG 010922 : I doubt the preconditioner is OK for the magnetization
1492    nspden_eff=nspden;if (optres==1.and.(.not.magn_precon)) nspden_eff=1
1493    do ispden=1,nspden_eff
1494
1495      diemix_eff=diemix;if (ispden>1) diemix_eff=diemixmag
1496
1497 !    Do fft from real space (work2) to G space (work1)
1498      if (optreal==1) then
1499        work2(:)=vresid(:,ispden)
1500        call fourdp(cplex,work1,work2,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1501      else
1502 !      work1(:,:)=reshape(vresid(:,ispden),(/2,nfft/))
1503 !      Reshape function does not work with big arrays for some compilers
1504        do ifft=1,nfft
1505          work1(1,ifft)=vresid(2*ifft-1,ispden)
1506          work1(2,ifft)=vresid(2*ifft  ,ispden)
1507        end do
1508      end if
1509
1510 !    Store G=0 value
1511      potg0(ispden)=work1(re,1)
1512
1513 !    Triple loop, for the three dimensions
1514      do i3=1,n3
1515 !      Precompute some products that do not depend on i2 and i1
1516        gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
1517        gqgm23=gq(3,i3)*gmet(2,3)*2
1518        gqgm13=gq(3,i3)*gmet(1,3)*2
1519        do i2=1,n2
1520          if (fftn2_distrib(i2)==me_fft) then
1521            gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
1522            gqgm12=gq(2,i2)*gmet(1,2)*2
1523            gqg2p3=gqgm13+gqgm12
1524            i23=n1*( ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
1525
1526 !          Do the test that eliminates the Gamma point outside
1527 !          of the inner loop
1528            ii1=1
1529            if(i2 == 1 .and. i3 == 1 .and. qeq0==1)then
1530 !            if(i23==0 .and. qeq0==1)then: this changes with the number of fft procs...
1531 !            and seems to be wrong.Pls check
1532              ii1=2
1533            end if
1534
1535 !          Here, unlike in hartre.f, the G=0 term is not eliminated, albeit
1536 !          not changed.
1537            do i1=ii1,n1
1538
1539 !            One obtains the square of the norm of q+G (defined by i1,i2,i3)
1540              gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
1541              ifft=i1+i23
1542
1543              l2g2=length2*gs
1544 !            The model dielectric function is now computed
1545              factor = (l2g2+diemac_inv)/(l2g2+1.0_dp) * diemix_eff
1546              work1(re,ifft)=work1(re,ifft)*factor
1547              work1(im,ifft)=work1(im,ifft)*factor
1548
1549            end do
1550          end if
1551        end do
1552      end do
1553
1554 !    Might get rid of the G=0 term
1555 !    if(qeq0==1)then
1556 !    work1(re,1)=0.0_dp
1557 !    work1(im,1)=0.0_dp
1558 !    end if
1559
1560 !    Fourier transform
1561      if (optreal==1) then
1562        call fourdp(cplex,work1,work2,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1563        vrespc(:,ispden)=work2(:)
1564      else
1565 !      vrespc(:,ispden)=reshape(work1(:,:),(/nfft*2/))
1566 !      Reshape function does not work with big arrays for some compilers
1567        do ifft=1,nfft
1568          vrespc(2*ifft-1,ispden)=work1(1,ifft)
1569          vrespc(2*ifft  ,ispden)=work1(2,ifft)
1570        end do
1571      end if
1572
1573 !    End of loop on spin polarizations
1574    end do
1575
1576    ABI_DEALLOCATE(gq)
1577    ABI_DEALLOCATE(work1)
1578    if (optreal==1) then
1579      ABI_DEALLOCATE(work2)
1580    end if
1581
1582 !  End condition diemac/=1.0
1583  end if
1584
1585 end subroutine moddiel


## ABINIT/prcref [ Functions ]

[ Top ] [ Functions ]

NAME

 prcref


FUNCTION

 Compute preconditioned residual potential (or density) and forces.
iprcel, densfor_pred and iprcfc govern the choice of the preconditioner.
1) Preconditioning of the forces (residual has already been included)
using the approximate force constant matrix. Get proposed
change of atomic positions.
2) Precondition the residual, get first part of proposed trial
potential change.
3) PAW only: precondition the rhoij residuals (simple preconditionning)
4) Take into account the proposed change of atomic positions to
modify the proposed trial potential change.

NOTE
This routine is almost similar to prcref_PMA.F90 which is employed in
case of potential mixing. Yet it has undergone strong changes simultaneously
from two different sources at the same time which resulted in a splitting.


INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
dielar(7)=input parameters for dielectric matrix:
diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
dielstrt=number of the step at which the dielectric preconditioning begins.
dtset <type(dataset_type)>=all input variables in this dataset
| densfor_pred= not yet used here
| iprcel= governs the preconditioning of the potential residual
|    0 => simple model dielectric matrix, described by the
|              parameters dielng, diemac, diemix and diemixmag contained in dielar.
|    between 21 and 39 => until istep=dielstart, same as iprcel=0, then uses
|              the RPA dielectric matrix (routine dielmt)
|    between 41 and 49 => uses the RPA dielectric matrix (routine dielmt).
|    between 51 and 59 => uses the RPA dielectric matrix (routine dieltcel).
|    between 61 and 69 => uses the electronic dielectric matr (routine dieltcel).
|    between 71 and 79 => uses the real-space preconditioner based on Kerker prc (prcrskerkerN)
|    between 81 and 99 => reserved for futur version of the real-space preconditioner
|    between 141 and 169 -> same as between 41 and 69 but with a different periodicity: modulo(iprcel modulo (10))
| iprcfc= governs the preconditioning of the forces
|         0 => hessian is the identity matrix
|         1 => hessian is 0.5 times the identity matrix
|         2 => hessian is 0.25 times the identity matrix
| ixc=exchange-correlation choice parameter.
| natom=number of atoms
| nspden=number of spin-density components
| occopt=option for occupancies
| prtvol=control print volume and debugging
| typat(natom)=integer type for each atom in cell
etotal=total ennergy
fcart(3,natom)=cartesian forces (hartree/bohr)
ffttomix(nfft*(1-nfftprc/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
gmet(3,3)=metrix tensor in G space in Bohr**-2.
gsqcut=cutoff on (k+G)^2 (bohr^-2)
istep= number of the step in the SCF cycle
kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
mgfft=maximum size of 1D FFTs
moved_atm_inside= if 1, then the preconditioned forces
as well as the preconditioned potential residual must be computed;
otherwise, compute only the preconditioned potential residual.
my_natom=number of atoms treated by current processor
nattyp(ntypat)=number of atoms of each type in cell.
nfft=number of fft grid points
nfftprc=size of FFT grid on which the potential residual will be preconditionned
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
ngfftprc(18)=contain all needed information about 3D FFT for the grid corresponding to nfftprc
nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description
npawmix=-PAW only- number of spherical part elements to be mixed
npwdiel=number of planewaves for dielectric matrix
ntypat=number of types of atoms in cell.
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
optreal=1 if residual potential is is REAL space, 2 if it is in RECIPROCAL SPACE
optres=0: the array vresid contains a potential residual
1: the array vresid contains a density residual
pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
Use here rhoij residuals (and gradients)
pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rhog(2,nfft)=array for electron density in reciprocal space
rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
rprimd(3,3)=dimensional primitive translations in real space (bohr)
susmat(2,npwdiel,nspden,npwdiel,nspden)=
the susceptibility (or density-density response) matrix in reciprocal space
vresid(optreal*nfftprc,nspden)=residual potential
vxc(nfft,nspden)=exchange-correlation potential (hartree)
vhartr(nfft)=array for holding Hartree potential
vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
vpsp(nfft)=array for holding local psp
xred(3,natom)=reduced dimensionless atomic coordinates


OUTPUT

  dtn_pc(3,natom)=preconditioned change of atomic position,
in reduced coordinates
vrespc(optreal*nfftprc,nspden)=preconditioned residual of the potential
==== if psps%usepaw==1
rhoijrespc(npawmix)= preconditionned rhoij residuals at output

SIDE EFFECT
dielinv(2,npwdiel,nspden,npwdiel,nspden)=
inverse of the dielectric matrix in rec. space
kxc(nfft,nkxc)=exchange-correlation kernel,
needed if the electronic dielectric matrix is computed
===== if densfor_pred==3 .and. moved_atm_inside==1 =====
ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases


PARENTS

      newrho


CHILDREN

      atm2fft,dielmt,dieltcel,fourdp,fresid,getph,hartre
indirect_parallel_fourier,kgindex,mean_fftr,metric,mkcore,mklocl
moddiel,prcrskerker1,prcrskerker2,rhotoxc,testsusmat,xcart2xred
xcdata_init,xmpi_sum,zerosym


SOURCE

195 subroutine prcref(atindx,dielar,dielinv,&
196 &  dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,&
197 &  istep,kg_diel,kxc,&
198 &  mgfft,moved_atm_inside,mpi_enreg,my_natom,&
199 &  nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
200 &  optreal,optres,pawrhoij,pawtab,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,&
201 &  susmat,vhartr,vpsp,vresid,vrespc,vxc,wvl,wvl_den,xred)
202
203
204 !This section has been created automatically by the script Abilint (TD).
205 !Do not modify the following lines by hand.
206 #undef ABI_FUNC
207 #define ABI_FUNC 'prcref'
208 !End of the abilint section
209
210  implicit none
211
212 !Arguments-------------------------------
213 !scalars
214  integer,intent(in) :: dielstrt,istep,my_natom,mgfft,moved_atm_inside,n1xccc
215  integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat,optreal,optres
216  real(dp),intent(in) :: etotal,gsqcut
217  type(MPI_type),intent(in) :: mpi_enreg
218  type(dataset_type),intent(in) :: dtset
219  type(pseudopotential_type),intent(in) :: psps
220  type(wvl_internal_type), intent(in) :: wvl
221  type(wvl_denspot_type), intent(inout) :: wvl_den
222 !arrays
223  integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft))
224  integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18)
225  real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),rhog(2,nfft)
226  real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3)
227  real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
228  real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden)
229  real(dp),intent(in) :: vxc(nfft,dtset%nspden)
230  real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
231  real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc)
232  real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft)
233  real(dp),intent(inout) :: xred(3,dtset%natom)
234  real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix)
235  real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden)
236  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
237  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
238
239 !Local variables-------------------------------
240 !scalars
241  integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1
242  integer :: ipw2,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm
243  integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method
244  real(dp) :: ai,ar,diemix,diemixmag,eei,enxc
245  real(dp) :: mixfac
246  real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg
247  logical :: computediel
248  logical :: non_magnetic_xc
249  character(len=500) :: message
250  type(xcdata_type) :: xcdata
251 !arrays
252  integer :: qprtrb(3)
253  integer,allocatable :: indpw_prc(:)
254  real(dp) :: dummy6(6),dummy7(6),gprimd(3,3),qphon(3),rmet(3,3),strsxc(6)
255  real(dp) :: vmean(dtset%nspden),vprtrb(2)
256  real(dp),allocatable :: dummy(:),dummy1(:),dummy2(:),dummy3(:),dummy4(:),dummy5(:),dummy8(:),dummy9(:)
257  real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:)
258  real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:)
259  real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:),rhor_new(:,:)
260  real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:)
261  real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:)
262  real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:)
264
265 ! *************************************************************************
266
267 !Compute different geometric tensor, as well as ucvol, from rprimd
268  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
269
270 ! Initialise non_magnetic_xc for rhohxc
271  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
272
273 !1) Eventually take care of the forces
274
275  if(moved_atm_inside==1)then
276    ABI_ALLOCATE(fcart_pc,(3,dtset%natom))
277
278    if(dtset%iprcfc==0)then
279      fcart_pc(:,:)=fcart(:,:)
280    else
281      fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:)
282    end if
283
284 !  Compute preconditioned delta xred from preconditioned fcart and rprimd
285    call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc)
286
287    ABI_DEALLOCATE(fcart_pc)
288  end if
289
290 !#######################################################################
291
292 !2) Take care of the potential residual
293
294 !Compute the residuals corresponding to the solution
295 !of an approximate realspace dielectric function according
296 !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]]
297  if(dtset%iprcel>=71.and.dtset%iprcel<=79) then
298    if (nfft==nfftprc) then
299      if (dtset%iprcel<=78) then
300        call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, &
301 &       gprimd,vresid,vrespc,rhor(:,1))
302      else
303        call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, &
304 &       vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol)
305      end if
306    else
307 !    If preconditionning has to be done on a coarse grid,
308 !    has to transfer several arrays
309      ABI_ALLOCATE(work1,(nfftprc,dtset%nspden))
310      ABI_ALLOCATE(work3,(nfftprc,dtset%nspden))
311      ABI_ALLOCATE(work,(2*nfftprc))
312      do ispden=1,dtset%nspden
313        work(:)=vresid(:,ispden)
314        call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
315      end do
316      ABI_DEALLOCATE(work)
317      if (dtset%iprcel<=78) then
318        ABI_ALLOCATE(rhog_wk,(2,nfftprc))
319        rhog_wk(:,:)=zero
320        if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then
321          nfftot=PRODUCT(ngfft(1:3))
322          call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,&
323 &         ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot)
324        else
325          do ii=1,nfft
326            if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii)
327          end do
328        end if
329        call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),&
330 &       comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
331        ABI_ALLOCATE(work,(nfftprc))
332        call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
333        call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, &
334 &       gprimd,work1,work3,work)
335        ABI_DEALLOCATE(work)
336      else
337        call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, &
338 &       work1,work3,dtset%natom,xred,mpi_enreg,ucvol)
339      end if
340      do ispden=1,dtset%nspden
341        call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
342      end do
343      ABI_DEALLOCATE(work1)
344      ABI_DEALLOCATE(work3)
345    end if
346
347  else
348
349    if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then
350      cplex=optreal
351      qphon(:)=zero
352 !    Simple scalar multiplication, or model dielectric function
353      call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,dtset%paral_kgb,qphon,rprimd,vresid,vrespc)
354
355 !    Use the inverse dielectric matrix in a small G sphere
356    else if( (istep>=dielstrt .and. dtset%iprcel>=21) .or. modulo(dtset%iprcel,100)>=41 )then
357
358 !    With dielop=1, the matrices will be computed when istep=dielstrt
359 !    With dielop=2, the matrices will be computed when istep=dielstrt and 1
360      dielop=1
361      if(modulo(dtset%iprcel,100)>=41)dielop=2
362      call testsusmat(computediel,dielop,dielstrt,dtset,istep) !test if the matrix is to be computed
363      if(computediel) then
364 !      Compute the inverse dielectric matrix from the susceptibility matrix
365 !      There are two routines for the RPA matrix, while for the electronic
366 !      dielectric matrix, only dieltcel will do the work
367        if(modulo(dtset%iprcel,100)<=49)then
368          call dielmt(dielinv,gmet,kg_diel,&
369 &         npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat)
370        else
371          option=1
372          if(modulo(dtset%iprcel,100)>=61)option=2
373          call dieltcel(dielinv,gmet,kg_diel,kxc,&
374 &         nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%paral_kgb,dtset%prtvol,susmat)
375        end if
376      end if
377
378      ABI_ALLOCATE(work1,(2,nfftprc))
379      ABI_ALLOCATE(work2,(optreal*nfftprc))
380
381 !    Presently, one uses the inverse of the RPA dielectric matrix,
382 !    for which spin must be averaged.
383
384 !    Do fft from real space (work2) to G space (work1)
385      if (optreal==1) then
386        work2(:)=vresid(:,1)
387 !      Must average over spins in the case of a potential residual
388        if(dtset%nspden/=1.and.optres==0)work2(:)=(work2(:)+vresid(:,2))*half
389        call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
390      else
391        work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/))
392        if(dtset%nspden/=1.and.optres==0)work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half
393      end if
394
395 !    Multiply by restricted inverse of dielectric matrix.
396 !    Must first copy relevant elements of work1 to a npwdiel-dimensioned array,
397 !    then zero work1, operate with the dielinv matrix, and store in work1.
398
399      ABI_ALLOCATE(vres_diel,(2,npwdiel))
400      ABI_ALLOCATE(indpw_prc,(npwdiel))
404
405      do ipw1=1,npwdiel
407          vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1))
408          vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1))
409        end if
410      end do
411      work1(:,:)=zero
412      do ipw1=1,npwdiel
413        ar=zero ; ai=zero
414
415 !      Use inverse of dielectric matrix (potential mixing)
416        if (optres==0) then
417          do ipw2=1,npwdiel
419              ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
420 &             -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
421              ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
422 &             +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
423            end if
424          end do
425        else
426 !        Use symetric of inverse of dielectric matrix (density mixing)
427          do ipw2=1,npwdiel
429              ar=ar+dielinv(1,ipw2,1,ipw1,1)*vres_diel(1,ipw2) &
430 &             +dielinv(2,ipw2,1,ipw1,1)*vres_diel(2,ipw2)
431              ai=ai-dielinv(2,ipw2,1,ipw1,1)*vres_diel(1,ipw2) &
432 &             +dielinv(1,ipw2,1,ipw1,1)*vres_diel(2,ipw2)
433            end if
434          end do
435        end if
436 !      Must be careful not to count the diagonal 1 twice : it is added later,
437 !      so must be subtracted now.
438        call xmpi_sum(ar,mpi_enreg%comm_fft,ier)
439        call xmpi_sum(ai,mpi_enreg%comm_fft,ier)
441          work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1)
442          work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1)
444      end do ! ipw1
445      ABI_DEALLOCATE(vres_diel)
446      ABI_DEALLOCATE(indpw_prc)
448 !    Fourier transform
449      if (optreal==1) then
450        call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
451      else
452        work2(:)=reshape(work1(:,:),(/nfftprc*2/))
453      end if
454
455 !    Add to get the preconditioned vresid, must be careful about spins.
456      if(dtset%iprcel>=30)then
457        diemix=dielar(4);diemixmag=abs(dielar(7))
458        vrespc(:,1)=diemix*(vresid(:,1)+work2(:))
459        if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:))
460        if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=diemixmag*vresid(:,3:4)
461        if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=diemixmag*vresid(:,2:dtset%nspden)
462      else
463        vrespc(:,1)=vresid(:,1)+work2(:)
464        if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=vresid(:,2)+work2(:)
465        if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=vresid(:,3:4)
466        if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=vresid(:,2:dtset%nspden)
467      end if
468
469      ABI_DEALLOCATE(work1)
470      ABI_DEALLOCATE(work2)
471
472 !    Other choice ?
473    else
474      write(message, '(a,i3,a,a,a,a)' )&
475 &     'From the calling routine, iprcel=',dtset%iprcel,ch10,&
476 &     'The only allowed values are 0 or larger than 20.',ch10,&
477 &     'Action: correct your input file.'
478      MSG_ERROR(message)
479    end if
480  end if
481 !#######################################################################
482
483 !3) PAW only : precondition the rhoij quantities (augmentation
484 !occupancies) residuals. Use a simple preconditionning
485 !with the same mixing factor as the model dielectric function.
486
487  if (psps%usepaw==1.and.my_natom>0) then
488    if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then
489      mixfac=one;mixfacmag=one
490    else
491      mixfac=dielar(4);mixfacmag=abs(dielar(7))
492    end if
493    if (pawrhoij(1)%cplex==1) then
494      index=0
495      do iatom=1,my_natom
496        do ispden=1,pawrhoij(iatom)%nspden
497          mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
498          do kmix=1,pawrhoij(iatom)%lmnmix_sz
499            index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix)
500            rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden)
501          end do
502        end do
503      end do
504    else
505      index=-1
506      do iatom=1,my_natom
507        do ispden=1,pawrhoij(iatom)%nspden
508          mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
509          do kmix=1,pawrhoij(iatom)%lmnmix_sz
510            index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
511            rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
512          end do
513        end do
514      end do
515    end if
516  end if
517
518 !#######################################################################
519
520 !4) Take care of the change of atomic positions
521 !Note : this part is very demanding on memory...
522 !however, since this algorithm is still in development,
523 !it was NOT included in the estimation provided by memory.f
524  if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then
525
526 !  Not yet compatible with resid given in reciprocal space
527    if (optreal/=1) then
528      write(message, '(5a)' )&
529 &     'From the calling routine, densfor_pred=3',ch10,&
530 &     'You cannot use residuals in reciprocal space.',ch10,&
531 &     'Action: correct your input file.'
532      MSG_ERROR(message)
533    end if
534 !  Not compatible with non-collinear magnetism
535    if(dtset%nspden==4)then
536      MSG_ERROR('densfor_pred=3 does not work for nspden=4 !')
537    end if
538
539    n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
540    nfftot=PRODUCT(ngfft(1:3))
541
542    if (optres==0) then  ! Array vresid contains a potential residual
543 !    -----------------------------------------------------------------
544
545 !    First subtract the current local, hartree and exchange correlation potentials
546      do ispden=1,min(dtset%nspden,2)
547        vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden)
548      end do
549      if (dtset%nspden==4) then
550        do ispden=3,4
551          vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden)
552        end do
553      end if
554
555 !    Compute the modified density, in rhor_wk
556      option=2
557      ABI_ALLOCATE(gresid,(3,dtset%natom))
558      ABI_ALLOCATE(grxc,(3,dtset%natom))
559      ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden))
560      ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden))
561      ABI_ALLOCATE(xred_wk,(3,dtset%natom))
562      xred_wk(:,:)=xred(:,:)+dtn_pc(:,:)
563      call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
564 &     ntypat,option,pawtab,rhor,rprimd,&
565 &     ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp)
566
567 !    Compute up+down rhog_wk(G) by fft
568      ABI_ALLOCATE(work,(nfft))
569      ABI_ALLOCATE(rhog_wk,(2,nfft))
570      work(:)=rhor_wk(:,1)
571      call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
572      ABI_DEALLOCATE(work)
573
574 !    Compute structure factor phases for new atomic pos:
575      call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk)
576
577 !    Compute local ionic pseudopotential vpsp:
578 !    and core electron density xccc3d, if needed.
579      n3xccc=0;if (n1xccc/=0) n3xccc=nfft
580      ABI_ALLOCATE(xccc3d,(n3xccc))
581      ABI_ALLOCATE(vpsp_wk,(nfft))
582      vprtrb(1:2)=zero
583
584 !    Determine by which method the local ionic potential and/or
585 !    the pseudo core charge density contributions have to be computed
586 !    Local ionic potential:
587 !     Method 1: PAW
588 !     Method 2: Norm-conserving PP, icoulomb>0, wavelets
589      vloc_method=1;if (psps%usepaw==0) vloc_method=2
590      if (dtset%icoulomb>0) vloc_method=2
591      if (psps%usewvl==1) vloc_method=2
592 !    Pseudo core charge density:
593 !     Method 1: PAW, nc_xccc_gspace
594 !     Method 2: Norm-conserving PP, wavelets
595      coredens_method=1;if (psps%usepaw==0) coredens_method=2
596      if (psps%nc_xccc_gspace==1) coredens_method=1
597      if (psps%nc_xccc_gspace==0) coredens_method=2
598      if (psps%usewvl==1) coredens_method=2
599
600 !    Local ionic potential and/or pseudo core charge by method 1
601      if (vloc_method==1.or.coredens_method==1) then
602        optv=0;if (vloc_method==1) optv=1
603        optn=0;if (coredens_method==1) optn=n3xccc/nfft
604        optatm=1;optdyfr=0;optgr=0;optstr=0;optn2=1;opteltfr=0
605 !      Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused...
606        call atm2fft(atindx,xccc3d,vpsp,dummy,dummy2,dummy9,dummy1,gmet,gprimd,dummy3,dummy4,gsqcut,&
607 &       mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,&
608 &       optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dummy5,dummy6,dummy7,&
609 &       ucvol,psps%usepaw,dummy8,dummy8,dummy8,vprtrb,psps%vlspl,&
610 &       comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
611 &       paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
612      end if
613
614 !    Local ionic potential by method 2
615      if (vloc_method==2) then
616        option=1
617        ABI_ALLOCATE(dyfrlo_indx,(3,3,dtset%natom))
618        ABI_ALLOCATE(grtn_indx,(3,dtset%natom))
619        call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,&
620 &       mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,&
621 &       ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,&
622 &       ucvol,vprtrb,vpsp_wk,wvl,wvl_den,xred)
623        ABI_DEALLOCATE(dyfrlo_indx)
624        ABI_DEALLOCATE(grtn_indx)
625      end if
626
627 !    Pseudo core electron density by method 2
628      if (coredens_method==2.and.n1xccc/=0) then
629        option=1
630        ABI_ALLOCATE(dyfrx2,(3,3,dtset%natom))
631        ABI_ALLOCATE(grxc_indx,(3,dtset%natom))
632        call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
633 &       n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,&
634 &       psps%xccc1d,xccc3d,xred_wk)
635        ABI_DEALLOCATE(dyfrx2)
636        ABI_DEALLOCATE(grxc_indx)
637      end if
638
639 !    Compute Hartree+xc potentials
640      ABI_ALLOCATE(vxc_wk,(nfft,dtset%nspden))
641      ABI_ALLOCATE(vhartr_wk,(nfft))
642      option=1
643
644      call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog_wk,rprimd,vhartr_wk)
645
646 !    Prepare the call to rhotoxc
647      call xcdata_init(xcdata,dtset=dtset)
648      nk3xc=1
649      call rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft,&
650 &     work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor_wk,rprimd,strsxc,1,&
651 &     vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk)
652      ABI_DEALLOCATE(xccc3d)
653
654 !    Sum all contributions
655      do ispden=1,min(dtset%nspden,2)
656        do ifft=1,nfft
657          vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden)
658        end do
659      end do
660      if (dtset%nspden==4) then
661        do ispden=3,4
662          do ifft=1,nfft
663            vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden)
664          end do
665        end do
666      end if
667      call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
668      if(dtset%nspden==2) then
669        vmean(1)=half*(vmean(1)+vmean(2))
670        vmean(2)=vmean(1)
671      end if
672      do ispden=1,dtset%nspden
673        vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden)
674      end do
675      ABI_DEALLOCATE(gresid)
676      ABI_DEALLOCATE(grxc)
677      ABI_DEALLOCATE(rhog_wk)
678      ABI_DEALLOCATE(rhor_wk)
679      ABI_DEALLOCATE(rhor_wk0)
680      ABI_DEALLOCATE(xred_wk)
681      ABI_DEALLOCATE(vhartr_wk)
682      ABI_DEALLOCATE(vpsp_wk)
683      ABI_DEALLOCATE(vxc_wk)
684
685    else                 ! Array vresid contains a density residual
686 !    -----------------------------------------------------------------
687
688 !    Only have to compute the modified preconditionned density residual
689      option=2
690      ABI_ALLOCATE(gresid,(3,dtset%natom))
691      ABI_ALLOCATE(grxc,(3,dtset%natom))
692      ABI_ALLOCATE(rhor_new,(nfft,dtset%nspden))
693      ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden))
694      ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden))
695      ABI_ALLOCATE(xred_wk,(3,dtset%natom))
696      xred_wk(:,:)=xred(:,:)+dtn_pc(:,:)
697      rhor_new(:,1)=rhor(:,1)+vrespc(:,1)
698      if (dtset%nspden==2) then
699        rhor_new(:,1)=rhor_new(:,1)+vrespc(:,2)
700        rhor_new(:,2)=rhor(:,2)+vrespc(:,1)
701      end if
702      call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
703 &     ntypat,option,pawtab,rhor,rprimd,&
704 &     ucvol,rhor_wk0,xred_wk,xred,psps%znuclpsp)
705      call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
706 &     ntypat,option,pawtab,rhor_new,rprimd,&
707 &     ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp)
708      vrespc(:,1)=rhor_wk(:,dtset%nspden)-rhor_wk0(:,dtset%nspden)
709      if (dtset%nspden==2) vrespc(:,2)=rhor_wk(:,1)-rhor_wk0(:,1)-vrespc(:,1)
710      ABI_DEALLOCATE(gresid)
711      ABI_DEALLOCATE(grxc)
712      ABI_DEALLOCATE(rhor_new)
713      ABI_DEALLOCATE(rhor_wk)
714      ABI_DEALLOCATE(rhor_wk0)
715      ABI_DEALLOCATE(xred_wk)
716    end if
717
718  end if
719
720 end subroutine prcref


## ABINIT/prcref_PMA [ Functions ]

[ Top ] [ Functions ]

NAME

 prcref_PMA


FUNCTION

 Compute preconditioned residual potential (or density) and forces.
iprcel, densfor_pred and iprcfc govern the choice of the preconditioner.
1) Preconditioning of the forces (residual has already been included)
using the approximate force constant matrix. Get proposed
change of atomic positions.
2) Precondition the residual, get first part of proposed trial
potential change.
3) PAW only: precondition the rhoij residuals (simple preconditionning)
4) Take into account the proposed change of atomic positions to
modify the proposed trial potential change.

NOTE
This routine is almost similar to prcref.F90 which is employed in
case of density mixing. Yet it has undergone strong changes simultaneously
from two different sources at the same time which resulted in a splitting.


INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
dielar(7)=input parameters for dielectric matrix:
diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
dielstrt=number of the step at which the dielectric preconditioning begins.
dtset <type(dataset_type)>=all input variables in this dataset
| densfor_pred= not yet used here
| iprcel= governs the preconditioning of the potential residual
|    0 => simple model dielectric matrix, described by the
|              parameters dielng, diemac, diemix and diemixmag contained in dielar.
|    between 21 and 39 => until istep=dielstart, same as iprcel=0, then uses
|              the RPA dielectric matrix (routine dielmt)
|    between 41 and 49 => uses the RPA dielectric matrix (routine dielmt).
|    between 51 and 59 => uses the RPA dielectric matrix (routine dieltcel).
|    between 61 and 69 => uses the electronic dielectric matr (routine dieltcel).
|    between 71 and 79 => uses the real-space preconditioner based on Kerker prc (prcrskerkerN)
|    between 81 and 99 => reserved for futur version of the real-space preconditioner
|    between 141 and 169 -> same as between 41 and 69 but with a different periodicity: modulo(iprcel modulo (10))
| iprcfc= governs the preconditioning of the forces
|         0 => hessian is the identity matrix
|         1 => hessian is 0.5 times the identity matrix
|         2 => hessian is 0.25 times the identity matrix
| ixc=exchange-correlation choice parameter.
| natom=number of atoms
| nspden=number of spin-density components
| occopt=option for occupancies
| prtvol=control print volume and debugging
| typat(natom)=integer type for each atom in cell
fcart(3,natom)=cartesian forces (hartree/bohr)
ffttomix(nfft*(1-nfftprc/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
gmet(3,3)=metrix tensor in G space in Bohr**-2.
gsqcut=cutoff on (k+G)^2 (bohr^-2)
istep= number of the step in the SCF cycle
kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
mgfft=maximum size of 1D FFTs
moved_atm_inside= if 1, then the preconditioned forces
as well as the preconditioned potential residual must be computed;
otherwise, compute only the preconditioned potential residual.
my_natom=number of atoms treated by current processor
nattyp(ntypat)=number of atoms of each type in cell.
nfft=number of fft grid points
nfftprc=size of FFT grid on which the potential residual will be preconditionned
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
ngfftprc(18)=contain all needed information about 3D FFT for the grid corresponding to nfftprc
nkxc=second dimension of the array kxc, see rhotoxc.f for a description
npawmix=-PAW only- number of spherical part elements to be mixed
npwdiel=number of planewaves for dielectric matrix
ntypat=number of types of atoms in cell.
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
optreal=1 if residual potential is is REAL space, 2 if it is in RECIPROCAL SPACE
optres=0: the array vresid contains a potential residual
1: the array vresid contains a density residual
pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
Use here rhoij residuals (and gradients)
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rhog(2,nfft)=array for electron density in reciprocal space
rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
rprimd(3,3)=dimensional primitive translations in real space (bohr)
susmat(2,npwdiel,nspden,npwdiel,nspden)=
the susceptibility (or density-density response) matrix in reciprocal space
vresid(optreal*nfftprc,nspden)=residual potential
vxc(nfft,nspden)=exchange-correlation potential (hartree)
vhartr(nfft)=array for holding Hartree potential
vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
vpsp(nfft)=array for holding local psp
xred(3,natom)=reduced dimensionless atomic coordinates

etotal
pawtab


OUTPUT

  dtn_pc(3,natom)=preconditioned change of atomic position,
in reduced coordinates
vrespc(optreal*nfftprc,nspden)=preconditioned residual of the potential
==== if psps%usepaw==1
rhoijrespc(npawmix)= preconditionned rhoij residuals at output

SIDE EFFECT
dielinv(2,npwdiel,nspden,npwdiel,nspden)=
inverse of the dielectric matrix in rec. space
kxc(nfft,nkxc)=exchange-correlation kernel,
needed if the electronic dielectric matrix is computed
===== if densfor_pred==3 .and. moved_atm_inside==1 =====
ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases


PARENTS

      newvtr


CHILDREN

      atm2fft,dielmt,dieltcel,fourdp,fresid,getph,hartre
indirect_parallel_fourier,kgindex,mean_fftr,metric,mkcore,mklocl
moddiel,prcrskerker1,prcrskerker2,rhotoxc,testsusmat,xcart2xred
xcdata_init,xmpi_sum,zerosym


SOURCE

 844   subroutine prcref_PMA(atindx,dielar,dielinv,&
845 &  dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,&
846 &  istep,kg_diel,kxc,&
847 &  mgfft,moved_atm_inside,mpi_enreg,my_natom,&
848 &  nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
849 &  optreal,optres,pawrhoij,ph1d,psps,rhog, rhoijrespc,rhor,rprimd,&
850 &  susmat,vhartr,vpsp,vresid,vrespc,vxc,xred,&
851 &  etotal,pawtab,wvl)
852
853
854 !This section has been created automatically by the script Abilint (TD).
855 !Do not modify the following lines by hand.
856 #undef ABI_FUNC
857 #define ABI_FUNC 'prcref_PMA'
858 !End of the abilint section
859
860  implicit none
861
862 !Arguments-------------------------------
863 !variables used for tfvw
864 !scalars
865  integer,intent(in) :: dielstrt,istep,mgfft,moved_atm_inside,my_natom,n1xccc
866  integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat
867  integer,intent(in) :: optreal,optres
868  real(dp),intent(in) :: etotal,gsqcut
869  type(MPI_type),intent(in) :: mpi_enreg
870  type(dataset_type),intent(in) :: dtset
871  type(pseudopotential_type),intent(in) :: psps
872  type(wvl_data), intent(inout) :: wvl
873 !arrays
874  integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft))
875  integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18)
876  real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom)
877  real(dp),intent(in) :: rhog(2,nfft)
878  real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3)
879  real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
880  real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden)
881  real(dp),intent(in) :: vxc(nfft,dtset%nspden)
882  real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
883  real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc)
884  real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft)
885  real(dp),intent(inout) :: xred(3,dtset%natom)
886  real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix)
887  real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden)
888  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
889  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
890
891 !Local variables-------------------------------
892 !scalars
893  integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1
894  integer :: ipw2,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm
895  integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method
896  real(dp) :: ai,ar,diemix,diemixmag,eei,enxc
897  real(dp) :: mixfac
898  real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg
899  logical :: computediel
900  logical :: non_magnetic_xc
901  character(len=500) :: message
902  type(xcdata_type) :: xcdata
903 !arrays
904  integer :: qprtrb(3)
905  integer,allocatable :: indpw_prc(:)
906  real(dp) :: dummy6(6),gprimd(3,3),qphon(3),rmet(3,3),strsxc(6)
907  real(dp) :: vmean(dtset%nspden),vprtrb(2)
908  real(dp),allocatable :: dummy_in(:)
909  real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0)
910  real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:)
911  real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:)
912  real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:)
913  real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:)
914  real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:)
915  real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:)
917
918 ! *************************************************************************
919
920  if(optres==1)then
921    MSG_ERROR('density mixing (optres=1) not admitted!')
922  end if
923
924 !Compute different geometric tensor, as well as ucvol, from rprimd
925  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
926
927 ! Initialise non_magnetic_xc for rhohxc
928  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
929
930 !1) Eventually take care of the forces
931
932  if(moved_atm_inside==1)then
933    ABI_ALLOCATE(fcart_pc,(3,dtset%natom))
934
935    if(dtset%iprcfc==0)then
936      fcart_pc(:,:)=fcart(:,:)
937    else
938      fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:)
939    end if
940
941 !  Compute preconditioned delta xred from preconditioned fcart and rprimd
942    call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc)
943
944    ABI_DEALLOCATE(fcart_pc)
945  end if
946
947 !#######################################################################
948
949 !2) Take care of the potential residual
950
951 !Compute the residuals corresponding to the solution
952 !of an approximate realspace dielectric function according
953 !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]]
954  if(dtset%iprcel>=71.and.dtset%iprcel<=79) then
955    if (nfft==nfftprc) then
956      if (dtset%iprcel<=78) then
957        call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, &
958 &       gprimd,vresid,vrespc,rhor(:,1))
959      else
960        call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, &
961 &       vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol)
962      end if
963    else
964 !    If preconditionning has to be done on a coarse grid,
965 !    has to transfer several arrays
966      ABI_ALLOCATE(work1,(nfftprc,dtset%nspden))
967      ABI_ALLOCATE(work3,(nfftprc,dtset%nspden))
968      ABI_ALLOCATE(work,(2*nfftprc))
969      do ispden=1,dtset%nspden
970        work(:)=vresid(:,ispden)
971        call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
972      end do
973      ABI_DEALLOCATE(work)
974      if (dtset%iprcel<=78) then
975        ABI_ALLOCATE(rhog_wk,(2,nfftprc))
976        rhog_wk(:,:)=zero
977        if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then
978          nfftot=PRODUCT(ngfft(1:3))
979          call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,&
980 &         ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot)
981        else
982          do ii=1,nfft
983            if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii)
984          end do
985        end if
986        call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),&
987 &       comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
988        ABI_ALLOCATE(work,(nfftprc))
989        call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
990        call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, &
991 &       gprimd,work1,work3,work)
992        ABI_DEALLOCATE(work)
993      else
994        call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, &
995 &       work1,work3,dtset%natom,xred,mpi_enreg,ucvol)
996      end if
997      do ispden=1,dtset%nspden
998        call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
999      end do
1000      ABI_DEALLOCATE(work1)
1001      ABI_DEALLOCATE(work3)
1002    end if
1003
1004  else
1005
1006    if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then
1007      cplex=optreal
1008      qphon(:)=zero
1009 !    Simple scalar multiplication, or model dielectric function
1010      call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,dtset%paral_kgb,qphon,rprimd,vresid,vrespc)
1011
1012 !    Use the inverse dielectric matrix in a small G sphere
1013    else if( (istep>=dielstrt .and. dtset%iprcel>=21) .or. modulo(dtset%iprcel,100)>=41 )then
1014
1015 !    Wnith dielop=1, the matrices will be computed when istep=dielstrt
1016 !    With dielop=2, the matrices will be computed when istep=dielstrt and 1
1017      dielop=1
1018      if(modulo(dtset%iprcel,100)>=41)dielop=2
1019      call testsusmat(computediel,dielop,dielstrt,dtset,istep) !test if the matrix is to be computed
1020      if(computediel) then
1021 !      Compute the inverse dielectric matrix from the susceptibility matrix
1022 !      There are two routines for the RPA matrix, while for the electronic
1023 !      dielectric matrix, only dieltcel will do the work
1024        if(modulo(dtset%iprcel,100)<=49)then
1025          call dielmt(dielinv,gmet,kg_diel,&
1026 &         npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat)
1027        else
1028          option=1
1029          if(modulo(dtset%iprcel,100)>=61)option=2
1030          call dieltcel(dielinv,gmet,kg_diel,kxc,&
1031 &         nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%paral_kgb,dtset%prtvol,susmat)
1032        end if
1033      end if
1034
1035      ABI_ALLOCATE(work1,(2,nfftprc))
1036      ABI_ALLOCATE(work2,(optreal*nfftprc))
1037
1038 !    Presently, one uses the inverse of the RPA dielectric matrix,
1039 !    for which spin must be averaged.
1040
1041 !    Do fft from real space (work2) to G space (work1)
1042      if (optreal==1) then
1043        work2(:)=vresid(:,1)
1044 !      Must average over spins if needed.
1045        if(dtset%nspden/=1)work2(:)=(work2(:)+vresid(:,2))*half
1046        call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
1047      else
1048        work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/))
1049        if (dtset%nspden/=1) work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half
1050      end if
1051
1052 !    Multiply by restricted inverse of dielectric matrix.
1053 !    Must first copy relevant elements of work1 to a npwdiel-dimensioned array,
1054 !    then zero work1, operate with the dielinv matrix, and store in work1.
1055
1056      ABI_ALLOCATE(vres_diel,(2,npwdiel))
1057      ABI_ALLOCATE(indpw_prc,(npwdiel))
1061      do ipw1=1,npwdiel
1063          vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1))
1064          vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1))
1065        end if
1066      end do
1067
1068      work1(:,:)=zero
1069      do ipw1=1,npwdiel
1070        ar=zero ; ai=zero
1071
1072 !      Use inverse of dielectric matrix (potential mixing)
1073        do ipw2=1,npwdiel
1075            ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
1076 &           -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
1077            ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
1078 &           +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
1079          end if
1080        end do
1081 !      Must be careful not to count the diagonal 1 twice : it is added later,
1082 !      so must be subtracted now.
1083        call xmpi_sum(ar,mpi_enreg%comm_fft,ier)
1084        call xmpi_sum(ai,mpi_enreg%comm_fft,ier)
1086          work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1)
1087          work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1)
1089      end do ! ipw1
1090      ABI_DEALLOCATE(vres_diel)
1091      ABI_DEALLOCATE(indpw_prc)
1093
1094 !    Fourier transform
1095      if (optreal==1) then
1096        call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0)
1097      else
1098        work2(:)=reshape(work1(:,:),(/nfftprc*2/))
1099      end if
1100
1101 !    Add to get the preconditioned vresid, must be careful about spins.
1102      if(dtset%iprcel>=30)then
1103        diemix=dielar(4);diemixmag=abs(dielar(7))
1104        vrespc(:,1)=diemix*(vresid(:,1)+work2(:))
1105        if(dtset%nspden/=1)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:))
1106        if(dtset%nspden==4)vrespc(:,3:4)=diemixmag*vresid(:,3:4)
1107      else
1108        vrespc(:,1)=vresid(:,1)+work2(:)
1109        if(dtset%nspden/=1)vrespc(:,2)=vresid(:,2)+work2(:)
1110        if(dtset%nspden==4)vrespc(:,3:4)=vresid(:,3:4)
1111      end if
1112
1113      ABI_DEALLOCATE(work1)
1114      ABI_DEALLOCATE(work2)
1115
1116 !    Other choice ?
1117    else
1118      write(message, '(a,i0,a,a,a,a)' )&
1119 &     'From the calling routine, iprcel= ',dtset%iprcel,ch10,&
1120 &     'The only allowed values are 0 or larger than 20.',ch10,&
1121 &     'Action: correct your input file.'
1122      MSG_ERROR(message)
1123    end if
1124  end if
1125 !#######################################################################
1126
1127 !3) PAW only : precondition the rhoij quantities (augmentation
1128 !occupancies) residuals. Use a simple preconditionning
1129 !with the same mixing factor as the model dielectric function.
1130
1131  if (psps%usepaw==1.and.my_natom>0) then
1132    if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then
1133      mixfac=one;mixfacmag=one
1134    else
1135      mixfac=dielar(4);mixfacmag=abs(dielar(7))
1136    end if
1137    if (pawrhoij(1)%cplex==1) then
1138      index=0
1139      do iatom=1,my_natom
1140        do ispden=1,pawrhoij(iatom)%nspden
1141          mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
1142          do kmix=1,pawrhoij(iatom)%lmnmix_sz
1143            index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix)
1144            rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden)
1145          end do
1146        end do
1147      end do
1148    else
1149      index=-1
1150      do iatom=1,my_natom
1151        do ispden=1,pawrhoij(iatom)%nspden
1152          mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
1153          do kmix=1,pawrhoij(iatom)%lmnmix_sz
1154            index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
1155            rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
1156          end do
1157        end do
1158      end do
1159    end if
1160  end if
1161 !#######################################################################
1162
1163 !4) Take care of the change of atomic positions
1164 !Note : this part is very demanding on memory...
1165 !however, since this algorithm is still in development,
1166 !it was NOT included in the estimation provided by memory.f
1167  if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then
1168
1169 !  Not yet compatible with resid given in reciprocal space
1170    if (optreal/=1) then
1171      write(message, '(5a)' )&
1172 &     'From the calling routine, densfor_pred=3',ch10,&
1173 &     'You cannot use residuals in reciprocal space.',ch10,&
1174 &     'Action: correct your input file.'
1175      MSG_ERROR(message)
1176    end if
1177
1178 !  Not compatible with non-collinear magnetism
1179    if(dtset%nspden==4)then
1180      MSG_ERROR('densfor_pred=3 does not work for nspden=4!')
1181    end if
1182
1183    n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1184    nfftot=PRODUCT(ngfft(1:3))
1185
1186 !  First subtract the current local, hartree and exchange correlation potentials
1187    do ispden=1,min(dtset%nspden,2)
1188      vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden)
1189    end do
1190    if (dtset%nspden==4) then
1191      do ispden=3,4
1192        vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden)
1193      end do
1194    end if
1195
1196 !  Compute the modified density, in rhor_wk
1197    option=2
1198    ABI_ALLOCATE(gresid,(3,dtset%natom))
1199    ABI_ALLOCATE(grxc,(3,dtset%natom))
1200    ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden))
1201    ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden))
1202    ABI_ALLOCATE(xred_wk,(3,dtset%natom))
1203    xred_wk(:,:)=xred(:,:)+dtn_pc(:,:)
1204
1205    call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
1206 &   ntypat,option,pawtab,rhor,rprimd,&
1207 &   ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp)
1208
1209 !  Compute up+down rhog_wk(G) by fft
1210    ABI_ALLOCATE(work,(nfft))
1211    ABI_ALLOCATE(rhog_wk,(2,nfft))
1212    work(:)=rhor_wk(:,1)
1213    call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
1214    ABI_DEALLOCATE(work)
1215
1216 !  Compute structure factor phases for new atomic pos:
1217    call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk)
1218
1219 !  Compute local ionic pseudopotential vpsp:
1220 !  and core electron density xccc3d, if needed.
1221    n3xccc=0;if (n1xccc/=0) n3xccc=nfft
1222    ABI_ALLOCATE(xccc3d,(n3xccc))
1223    ABI_ALLOCATE(vpsp_wk,(nfft))
1224    vprtrb(1:2)=zero
1225
1226 !  Determine by which method the local ionic potential and/or
1227 !  the pseudo core charge density have to be computed
1228 !  Local ionic potential:
1229 !   Method 1: PAW
1230 !   Method 2: Norm-conserving PP, icoulomb>0, wavelets
1231    vloc_method=1;if (psps%usepaw==0) vloc_method=2
1232    if (dtset%icoulomb>0) vloc_method=2
1233    if (psps%usewvl==1) vloc_method=2
1234 !  Pseudo core charge density:
1235 !   Method 1: PAW, nc_xccc_gspace
1236 !   Method 2: Norm-conserving PP, wavelets
1237    coredens_method=1;if (psps%usepaw==0) coredens_method=2
1238    if (psps%nc_xccc_gspace==1) coredens_method=1
1239    if (psps%nc_xccc_gspace==0) coredens_method=2
1240    if (psps%usewvl==1) coredens_method=2
1241
1242 !  Local ionic potential and/or pseudo core charge by method 1
1243    if (vloc_method==1.or.coredens_method==1) then
1244      optv=0;if (vloc_method==1) optv=1
1245      optn=0;if (coredens_method==1) optn=n3xccc/nfft
1246      optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1
1247 !    Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused...
1248      call atm2fft(atindx,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,gmet,&
1249 &     gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,&
1250 &     nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
1251 &     psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dummy_in,dummy_out6,dummy_out7,ucvol,&
1252 &     psps%usepaw,dummy_in,dummy_in,dummy_in,vprtrb,psps%vlspl,&
1253 &     comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
1254 &     paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
1255    end if
1256
1257 !  Local ionic potential by method 2
1258    if (vloc_method==2) then
1259      option=1
1260      ABI_ALLOCATE(dyfrlo_indx,(3,3,dtset%natom))
1261      ABI_ALLOCATE(grtn_indx,(3,dtset%natom))
1262      call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,&
1263 &     mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,&
1264 &     ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,&
1265 &     ucvol,vprtrb,vpsp_wk,wvl%descr,wvl%den,xred)
1266      ABI_DEALLOCATE(dyfrlo_indx)
1267      ABI_DEALLOCATE(grtn_indx)
1268    end if
1269
1270 !  Pseudo core electron density by method 2
1271    if (coredens_method==2.and.n1xccc/=0) then
1272      option=1
1273      ABI_ALLOCATE(dyfrx2,(3,3,dtset%natom))
1274      ABI_ALLOCATE(grxc_indx,(3,dtset%natom))
1275      call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
1276 &     n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,&
1277 &     psps%xccc1d,xccc3d,xred_wk)
1278      ABI_DEALLOCATE(dyfrx2)
1279      ABI_DEALLOCATE(grxc_indx)
1280    end if
1281
1282 !  Compute Hartree+xc potentials
1283    ABI_ALLOCATE(vxc_wk,(nfft,dtset%nspden))
1284    ABI_ALLOCATE(vhartr_wk,(nfft))
1285    option=1
1286
1287    call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog_wk,rprimd,vhartr_wk)
1288
1289 !  Prepare the call to rhotoxc
1290    call xcdata_init(xcdata,dtset=dtset)
1291    nk3xc=1
1292    ABI_ALLOCATE(work,(0))
1293    call rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft,&
1294 &   work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor_wk,rprimd,strsxc,1,&
1295 &   vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk)
1296    ABI_DEALLOCATE(work)
1297    ABI_DEALLOCATE(xccc3d)
1298
1299 !  Sum all contributions
1300    do ispden=1,min(dtset%nspden,2)
1301      do ifft=1,nfft
1302        vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden)
1303      end do
1304    end do
1305    if (dtset%nspden==4) then
1306      do ispden=3,4
1307        do ifft=1,nfft
1308          vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden)
1309        end do
1310      end do
1311    end if
1312    call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1313    if(dtset%nspden==2) then
1314      vmean(1)=half*(vmean(1)+vmean(2))
1315      vmean(2)=vmean(1)
1316    end if
1317    do ispden=1,dtset%nspden
1318      vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden)
1319    end do
1320    ABI_DEALLOCATE(gresid)
1321    ABI_DEALLOCATE(grxc)
1322    ABI_DEALLOCATE(rhog_wk)
1323    ABI_DEALLOCATE(rhor_wk)
1324    ABI_DEALLOCATE(rhor_wk0)
1325    ABI_DEALLOCATE(xred_wk)
1326    ABI_DEALLOCATE(vhartr_wk)
1327    ABI_DEALLOCATE(vpsp_wk)
1328    ABI_DEALLOCATE(vxc_wk)
1329
1330  end if
1331
1332 end subroutine prcref_PMA


## ABINIT/prcrskerker1 [ Functions ]

[ Top ] [ Functions ]

NAME

 prcrskerker1


FUNCTION

 preconditionning by a real-space conjugate gradient on residual
using a model dielectric function in real space


INPUTS

  nfft=number of fft grid points
nspden=number of spin-density components
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
dielar(7)=input parameters for dielectric matrix:
diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
gprimd(3,3)=dimensional primitive translations in fourier space (bohr**-1)
rprimd(3,3)=dimensional primitive translations in real space (bohr)
vresid(nfft,nspden)=residual potential
base(nfft) = real space function used as a basis to guess a fine dielectric funtion
see the calling routine to know the content


OUTPUT

  vrespc(nfft,nspden)=preconditioned residual of the potential


WARNINGS

 This is experimental code : input, ouptput, results and any other feature may vary greatly.


NOTES

  needs severe cleaning and this is abuse of modules as common blocks...


PARENTS

      prcref,prcref_PMA


CHILDREN

      cgpr,frskerker1__end,frskerker1__init,laplacian,prc_mem_init


SOURCE

2480 subroutine prcrskerker1(dtset,mpi_enreg,nfft,nspden,ngfft,dielar,etotal,gprimd,vresid,vrespc,base)
2481
2482
2483 !This section has been created automatically by the script Abilint (TD).
2484 !Do not modify the following lines by hand.
2485 #undef ABI_FUNC
2486 #define ABI_FUNC 'prcrskerker1'
2487 !End of the abilint section
2488
2489  implicit none
2490
2491 !Arguments ------------------------------------
2492 !scalars
2493  integer,intent(in) :: nfft,nspden
2494  real(dp) :: etotal
2495  type(MPI_type),intent(in) :: mpi_enreg
2496  type(dataset_type),intent(in) :: dtset
2497 !arrays
2498  integer,intent(in) :: ngfft(18)
2499  real(dp),intent(in) :: base(nfft),dielar(7),gprimd(3,3)
2500  real(dp),intent(in) :: vresid(nfft,nspden)
2501  real(dp),intent(out) :: vrespc(nfft,nspden)
2502
2503 !Local variables-------------------------------
2504 !scalars
2505  integer ::  ifft,ispden,n1,n2,n3
2506  real(dp) :: base_delta,base_max,base_min,dielng,diemac,diemix
2507  real(dp) :: diemixmag
2508  real(dp) :: rdummy1,rdummy2
2509  logical ::  new_prc_func
2510 !arrays
2511  real(dp) :: deltaW(nfft,nspden)
2512  real(dp) :: g2cart(nfft)
2513  real(dp) :: mat(nfft,nspden)
2514
2515 ! *************************************************************************
2516
2517 !DEBUG
2518 !write(std_out,*)' prckerker1 : enter '
2519 !ENDDEBUG
2520 !if(cycle==0) then
2521  call prc_mem_init(nfft)
2522
2523  if(cycle==0) then
2524    new_prc_func=.TRUE.
2525    energy_min=etotal
2526  else if(etotal < energy_min) then
2527    new_prc_func=.TRUE.
2528    energy_min=etotal
2529  else
2530    new_prc_func=.FALSE.
2531  end if
2532
2533
2534  dielng=dielar(2)
2535  diemac=dielar(3)
2536  diemix=dielar(4)
2537  diemixmag=dielar(7)
2538 !******************************************************************
2539 !compute the diemac(r)                                          **
2540 !******************************************************************
2541 !this task will be devoted to a general function later
2542  n1=ngfft(1)
2543  n2=ngfft(2)
2544  n3=ngfft(3)
2545 !base_cp=base
2546  if(new_prc_func) then
2547    base_min=base(1)
2548    base_max=base(1)
2549    do ifft=1,nfft
2550      base_min = min(base_min,base(ifft))
2551      base_max = max(base_max,base(ifft))
2552    end do
2553    base_delta = base_max - base_min
2554 !  if(cycle.lt.2) then
2555    rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one)
2556 !  else
2557 !  rdiemac(:) = rdiemac(:)*0.5_dp+0.5_dp*(((base(:)-base_min) / (base_delta) ) *(diemac-one) + one)
2558 !  end if
2559 !  if(cycle==0) rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one)
2560 !  rdiemac(:) = exp(((base(:)-base_min) / (base_delta) *log(diemac)))
2561  end if
2562  cycle=cycle+1
2563 !if(cycle==5) cycle=0
2564 !end if
2565 !******************************************************************
2566 !compute deltaW                                                 **
2567 !******************************************************************
2568  vrespc=vresid !starting point
2569  call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,dtset%paral_kgb,rdfuncr=vrespc,laplacerdfuncr=deltaW,g2cart_out=g2cart) ! put the laplacian of the residuals into deltaW
2570 !call laplacian(vrespc,buffer,ngfft,gprimd) ! put the laplacian of the residuals into deltaW
2571 !do ifft=1,nfft
2572 !if (buffer(ifft,1)/=deltaW(ifft,1)) then
2573 !stop
2574 !end if
2575 !end do
2576  deltaW(:,1)= diemix*(((one/rdiemac(:))*vresid(:,1))-(((dielng)**2)*deltaW(:,1)))
2577  if (nspden>1.and.(diemixmag>=zero)) then
2578    do ispden=2,nspden
2579      deltaW(:,ispden)= abs(diemixmag)*(((one/rdiemac(:))*vresid(:,ispden))-(((dielng)**2)*deltaW(:,ispden)))
2580    end do
2581  end if
2582 !call random_number(deltaW)
2583 !call random_number(vrespc)
2584 !******************************************************************
2585 !Finding the preconditionned residuals which minimizes          **
2586 !half*(vrespc*(1-dielng2/4pi2 nabla2) vrespc) - vrespc * deltaW **
2587 !***********************************************************************
2588  vrespc(:,1)=diemix*vrespc(:,1)
2589  if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*vrespc(:,2:nspden)
2590 !buffer=vrespc
2591
2592
2593 !==============================================================================
2594 !==============================================================================
2595 !! Original loop
2596 !==============================================================================
2597 !==============================================================================
2598
2599  call frskerker1__init(dtset,mpi_enreg,nfft,ngfft,nspden,dielng,deltaW,gprimd,mat,g2cart)
2600
2601 !call cgpr(pf_rscgres,dpf_rscgres,newvres,real(1e-40,dp),700,vrespc,rdummy1,rdummy2)
2602 !rdummy1 = pf_rscgres(nfft,nspden,vrespc)
2603  call cgpr(nfft,nspden,frskerker1__pf,frskerker1__dpf,frskerker1__newvres,&
2604 & real(1e-10,dp),700,vrespc,rdummy1,rdummy2)
2605  call frskerker1__end()
2606
2607 !==============================================================================
2608 !==============================================================================
2609 !! Original loop end
2610 !==============================================================================
2611 !==============================================================================
2612
2613
2614 !cplex=1
2615 !qphon(:)=zero
2616 !call moddiel(cplex,dielar,nfft,ngfft,nspden,1,0,qphon,rprimd,vresid,buffer)
2617 !c1=0
2618 !do ifft=1,nfft,1
2619 !if((abs(buffer(ifft,1)-vrespc(ifft,1))/(abs(buffer(ifft,1)+vrespc(ifft,1))*half)) > 5e-3) then
2620 !c1=c1+1
2621 !end if
2622 !end do
2623 !call laplacian(vrespc,buffer,ngfft,gprimd)
2624 !buffer=vrespc(:,:)-buffer(:,:)*dielng**2
2625 !c2=0
2626 !do ifft=1,nfft,1
2627 !if((abs(buffer(ifft,1)-deltaW(ifft,1))/(abs(buffer(ifft,1)+deltaW(ifft,1))*half)) > 5e-3) then
2628 !c2=c2+1
2629 !end if
2630 !end do
2631 !!!  !stop
2632 !call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,&
2633 !& g2cart_out=g2cart)
2634
2635 !vrespc=vresid
2636 !do ispden=1,nspden
2637 !call fourdp(1, gvrespc(:,:,ispden), vrespc(:,ispden),-1,mpi_enreg,nfft,ngfft,0)
2638 !end do
2639 !filtering
2640 !do ispden=1,nspden
2641 !do ifft=1,nfft
2642 !!    gvrespc(:,ifft,ispden)=(one-exp(-g2cart(ifft)*15.0_dp))*gvrespc(:,ifft,ispden)
2643 !!      gvrespc(:,ifft,ispden)=(exp(-g2cart(ifft)*10.0_dp))*gvrespc(:,ifft,ispden)
2644 !!      gvrespc(:,ifft,ispden)=(one-one/(exp(-0.002_dp/g2cart(ifft)**2)+one))*gvrespc(:,ifft,ispden)
2645 !gvrespc(:,ifft,ispden)=(two-2_dp/(exp(-0.008_dp/(g2cart(ifft)+0.0012_dp))+one))*gvrespc(:,ifft,ispden)
2646 !gvrespc(:,ifft,ispden)=min(one,(sqrt(g2cart(ifft)/0.006_dp))**(one))*gvrespc(:,ifft,ispden)
2647 !end do
2648 !end do
2649 !change resulting potential to real space
2650 !do ispden=1,nspden
2651 !call fourdp(1,gvrespc(:,:,ispden),vrespc(:,ispden),1,mpi_enreg,nfft,ngfft,0)
2652 !end do
2653 !vrespc=vrespc*diemix
2654 !maxg2=g2cart(1)
2655 !ming2=g2cart(5)
2656 !do ifft=1,nfft
2657 !maxg2=max(g2cart(ifft),maxg2)
2658 !if(g2cart(ifft) .gt. zero) ming2=min(g2cart(ifft),ming2)
2659 !end do
2660 !stop
2661
2662 !DEBUG
2663 !write(std_out,*)' prckerker1 : exit '
2664 !ENDDEBUG
2665
2666 end subroutine prcrskerker1


## ABINIT/prcrskerker2 [ Functions ]

[ Top ] [ Functions ]

NAME

 prcrskerker2


FUNCTION

 preconditionning by a real-space conjugate gradient on residual
using a model dielectric function in real space
differing from prcrskerker1 by the
use of a linear response approach


INPUTS

  nfft=number of fft grid points
nspden=number of spin-density components
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
dielar(7)=input parameters for dielectric matrix:
diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
gprimd(3,3)=dimensional primitive translations in fourier space (bohr**-1)
rprimd(3,3)=dimensional primitive translations in real space (bohr)
vresid(nfft,nspden)=residual potential


OUTPUT

  vrespc(nfft,nspden)=preconditioned residual of the potential


WARNINGS

 This is experimental code : input, ouptput, results and any other feature may vary greatly.


NOTES

PARENTS

      prcref,prcref_PMA


CHILDREN

      cgpr,dotprod_vn,frskerker2__end,frskerker2__init,laplacian,ptabs_fourdp


SOURCE

2705 subroutine prcrskerker2(dtset,nfft,nspden,ngfft,dielar,gprimd,rprimd,vresid,vrespc,natom,xred,mpi_enreg,ucvol)
2706
2707
2708 !This section has been created automatically by the script Abilint (TD).
2709 !Do not modify the following lines by hand.
2710 #undef ABI_FUNC
2711 #define ABI_FUNC 'prcrskerker2'
2712 !End of the abilint section
2713
2714  implicit none
2715
2716 !Arguments ------------------------------------
2717 !scalars
2718  integer,intent(in) :: natom,nfft,nspden
2719  real(dp),intent(in) :: ucvol
2720  type(MPI_type),intent(in) :: mpi_enreg
2721  type(dataset_type),intent(in) :: dtset
2722 !arrays
2723  integer,intent(in) :: ngfft(18)
2724  real(dp),intent(in) :: dielar(7),gprimd(3,3),rprimd(3,3),vresid(nfft,nspden)
2725  real(dp),intent(in) :: xred(3,natom)
2726  real(dp),intent(out) :: vrespc(nfft,nspden)
2727
2728 !Local variables-------------------------------
2729   !logical,save ::ok=.FALSE.
2730 !scalars
2731  integer :: cplex,i1,i2,i3,iatom,iatom27,ifft,ispden,n1,n2,n3,natom27,nfftotf
2732  integer :: option
2733  real(dp),save :: lastp1=one,lastp2=one
2734  real(dp) :: C1,C2,DE,core,dielng,diemac,diemix,diemixmag,doti,dr,l1,l2,l3,l4,r
2735  real(dp) :: rdummy1,rdummy2,rmin,xr,y,yr,zr
2736 !arrays
2737  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
2738  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
2739  real(dp) :: V1(nfft,nspden),V2(nfft,nspden),buffer(nfft,nspden)
2740  real(dp) :: deltaW(nfft,nspden)
2741  real(dp) :: mat(nfft,nspden)
2742  real(dp) :: rdielng(nfft),rdiemac(nfft),xcart(3,natom)
2743  real(dp) :: xcart27(3,natom*27)
2744
2745 ! *************************************************************************
2746
2747  dielng=dielar(2)
2748  diemac=dielar(3)
2749  diemix=dielar(4)
2750  diemixmag=dielar(7)
2751 !******************************************************************
2752 !compute the diemac(r)                                          **
2753 !******************************************************************
2754 !this task will be devoted to a general function later
2755  n1=ngfft(1)
2756  n2=ngfft(2)
2757  n3=ngfft(3)
2758  nfftotf=n1*n2*n3
2759 !if(.not.ok) then
2760  xcart(1,:)=xred(1,:)*rprimd(1,1)+xred(2,:)*rprimd(1,2)+xred(3,:)*rprimd(1,3)
2761  xcart(2,:)=xred(1,:)*rprimd(2,1)+xred(2,:)*rprimd(2,2)+xred(3,:)*rprimd(2,3)
2762  xcart(3,:)=xred(1,:)*rprimd(3,1)+xred(2,:)*rprimd(3,2)+xred(3,:)*rprimd(3,3)
2763
2764  iatom27=0
2765  do i1=-1,1
2766    do i2=-1,1
2767      do i3=-1,1
2768        do iatom=1,natom
2769          iatom27=iatom27+1
2770          xcart27(:,iatom27)=xcart(:,iatom)+rprimd(:,1)*i1+rprimd(:,2)*i2+rprimd(:,3)*i3
2771        end do
2772      end do
2773    end do
2774  end do
2775
2776 !stop
2777  natom27=27*natom
2778
2779  l1=0.34580850339844665
2780 !l2=0.5123510203906797 !0.41242551019533985
2781 !l3=0.8001489796093203 !0.90007448980466009
2782
2783  l2=0.41242551019533985
2784  l3=0.90007448980466009
2785  l4=0.9666914966015534
2786
2787
2788  l1=0.31387233559896449
2789  l2=0.35828367346355994
2790  l3=0.9333829932031068
2791  l4=0.9777943310677023
2792
2793  l1=3.5
2794  l2=11.5
2795  l3=2.5
2796  l4=6.5
2797 !l1=30. !cellules pleines
2798
2799  rdielng=zero
2800  core=1. !(value of Er at the core of atoms)
2801  dr=2.65 ! radius of atoms=2.65165
2802  y=1. ! value of Er in the empty region
2803
2804 !Get the distrib associated with this fft_grid
2805  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2806
2807  do i3=1,n3
2808    ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft)
2809    do i2=1,n2
2810      if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
2811        do i1=1,n1
2812          ifft=ifft+1
2813 !        !!!!!!!!!!!!!!!!!!!!!!!!!
2814 !        ! calculation of the simplest part void/metal
2815 !        !!              x=real(real(i3,dp)/real(n3,dp),dp)
2816 !        !!              !x=i3/n3
2817 !        !!              if(x < l1) then
2818 !        !!                 rdiemac(ifft)=diemac
2819 !        !!                 rdielng(ifft)=dielng
2820 !        !!              else if(x < l2) then
2821 !        !!                 xp=(l2-x)/(l2-l1)
2822 !        !!                 rdiemac(ifft)=y+(diemac-y)&
2823 !        !!                      & * (1.-(1.-xp)**4)**4
2824 !        !!                 rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4
2825 !        !!              else if(x < l3) then
2826 !        !!                 rdiemac(ifft)=y
2827 !        !!                 rdielng(ifft)=zero
2828 !        !!              else if(x < l4) then
2829 !        !!                 xp=(l3-x)/(l3-l4)
2830 !        !!                 rdiemac(ifft)=y+(diemac-y)&
2831 !        !!                      & * (1.-(1.-xp)**4)**4
2832 !        !!                 rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4
2833 !        !!              else
2834 !        !!                 rdiemac(ifft)=diemac
2835 !        !!                 rdielng(ifft)=dielng
2836 !        !!              end if
2837 !        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2838 !        !!!! calculation of atomic core dielectric
2839 !        !!              rmin=1e16
2840 !        !!              xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)&
2841 !        !!                   &+real((i3-1),dp)/real(n3,dp)*rprimd(1,3)
2842 !        !!              yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)&
2843 !        !!                   &+real((i3-1),dp)/real(n3,dp)*rprimd(2,3)
2844 !        !!              zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)&
2845 !        !!                   &+real((i3-1),dp)/real(n3,dp)*rprimd(3,3)
2846 !        !!              do iatom=1,natom27
2847 !        !!                 r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2
2848 !        !!                 if (r<rmin) then
2849 !        !!                    rmin=r
2850 !        !!                 end if
2851 !        !!              end do
2852 !        !!              if(rmin < dr**2) then
2853 !        !!                 rdiemac(ifft)=min(rdiemac(ifft),core+(diemac-core)*(1.-(1.-sqrt(rmin)/dr)**2)**2)
2854 !        !!                 rdielng(ifft)=dielng-dielng*(1.-(1.-sqrt(rmin)/dr)**4)**4
2855 !        !!              else
2856 !        !!                 rdiemac(ifft)=min(rdiemac(ifft),diemac)
2857 !        !!              end if
2858          rmin=1e16
2859          xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)&
2860 &         +real((i3-1),dp)/real(n3,dp)*rprimd(1,3)
2861          yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)&
2862 &         +real((i3-1),dp)/real(n3,dp)*rprimd(2,3)
2863          zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)&
2864 &         +real((i3-1),dp)/real(n3,dp)*rprimd(3,3)
2865
2866          rdiemac(ifft)=y
2867          rdielng(ifft)=zero
2868          do iatom=1,natom27
2869            r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2
2870            if (r<rmin) then
2871              rmin=r
2872            end if
2873            if(r < l1) then
2874              rdiemac(ifft)= rdiemac(ifft) +  0.7_dp * (diemac-y)
2875            else if(r < l2) then
2876              rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y)*(one-((sqrt(r)-l1)/(l2-l1))**2)**2
2877            else
2878              rdiemac(ifft)=rdiemac(ifft)
2879            end if
2880            if(r < l3) then
2881              rdielng(ifft)= rdielng(ifft) +  0.5_dp * (dielng)
2882            else if(r < l4) then
2883              rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng)  *(one-((sqrt(r)-l3)/(l4-l3))**2)**2
2884            end if
2885          end do
2886
2887          rdielng(ifft)=min(rdielng(ifft),dielng)
2888 !        rdielng(ifft)=dielng
2889          rdiemac(ifft)=min(rdiemac(ifft),diemac)
2890          rdiemac(ifft)=diemac
2891        end do
2892      end if
2893    end do
2894  end do
2895 !rdielng(:)=dielng
2896
2897 !****************************************************************************************
2898 !****************************************************************************************
2899 !****************************************************************************************
2900 !****************************************************************************************
2901 !******************************************************************
2902 !compute V1
2903 !******************************************************************
2904  V1=vresid
2905  call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,dtset%paral_kgb,rdfuncr=V1,laplacerdfuncr=deltaW)
2906  deltaW(:,1)= (((one/rdiemac(:))*V1(:,1))-(((rdielng(:))**2)*deltaW(:,1)))
2907 !deltaW(:,1)= -diemix*(((rdielng(:))**2)*deltaW(:,ispden))
2908  if (nspden>1) then
2909    do ispden=2,nspden
2910      deltaW(:,ispden)= (((one/rdiemac(:))*V1(:,ispden))-(((rdielng(:))**2)*deltaW(:,ispden)))
2911 !    deltaW(:,ispden)= -abs(diemixmag)*(((rdielng(:))**2)*deltaW(:,ispden))
2912    end do
2913  end if
2914  call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat)
2915  call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,&
2916 & frskerker2__newvres2,lastp1*real(1e-6 ,dp),700,V1,rdummy1,rdummy2)
2917  lastp1=min(abs(rdummy1),1e-6_dp)
2918  call frskerker2__end()
2919
2920 !******************************************************************
2921 !compute V2
2922 !******************************************************************
2923  V2=vresid
2924  do ispden=1,nspden
2925    deltaW(:,ispden)= (rdielng(:)**2)
2926  end do
2927  call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat)
2928  call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,&
2929 & frskerker2__newvres2,lastp2*real(1e-6,dp),700,V2,rdummy1,rdummy2)
2930  lastp2=min(abs(rdummy1),1e-6_dp)
2931  call frskerker2__end()
2932
2933
2934 !******************************************************************
2935 !compute C1, C2 & DE
2936 !******************************************************************
2937  cplex=1;
2938  option=1;
2939  call dotprod_vn(cplex,& !complex density/pot
2940 &rdielng,&          !the density
2941 &DE,&  !resulting dorproduct integrated over r  ! here DE is used has a buffer
2942 &doti,&          !imaginary part of the integral
2943 &size(rdielng,1),&          !number of localy(cpu) attributed grid point
2944 &nfftotf,&        !real total number of grid point
2945 &nspden,&        !nspden
2946 &option,&        !1=compute only the real part 2=compute also the imaginary part
2947 &rdielng,&          !the potential
2948 &ucvol,&          !cell volume
2949 &mpi_comm_sphgrid=mpi_enreg%comm_fft)
2950  do ispden=1,nspden
2951    buffer(:,ispden)=rdielng(:)*V1(:,ispden)
2952  end do
2953  call dotprod_vn(cplex,& !complex density/pot
2954 &rdielng,&          !the density
2955 &C1,&  !resulting dorproduct integrated over r  ! here DE is used has a buffer
2956 &doti,&          !imaginary part of the integral
2957 &size(rdielng,1),&          !number of localy(cpu) attributed grid point
2958 &nfftotf,&        !real total number of grid point
2959 &nspden,&        !nspden
2960 &option,&        !1=compute only the real part 2=compute also the imaginary part
2961 &buffer,&          !the potential
2962 &ucvol,&         !cell volume
2963 &mpi_comm_sphgrid=mpi_enreg%comm_fft)
2964  do ispden=1,nspden
2965    buffer(:,ispden)=rdielng(:)*V2(:,ispden)
2966  end do
2967  call dotprod_vn(cplex,& !complex density/pot
2968 &rdielng,&          !the density
2969 &C2,&  !resulting dorproduct integrated over r  ! here DE is used has a buffer
2970 &doti,&          !imaginary part of the integral
2971 &size(rdielng,1),&          !number of localy(cpu) attributed grid point
2972 &nfftotf,&        !real total number of grid point
2973 &nspden,&        !nspden
2974 &option,&        !1=compute only the real part 2=compute also the imaginary part
2975 &buffer,&          !the potential
2976 &ucvol,&         !cell volume
2977 &mpi_comm_sphgrid=mpi_enreg%comm_fft)
2978  C1=C1/DE
2979  C2=C2/DE
2980  DE=C1/(one-C2)
2981
2982 !******************************************************************
2983 !compute the new preconditionned residuals
2984 !******************************************************************
2985  vrespc(:,1)=diemix*(V1(:,1)+DE*V2(:,1))
2986  if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*(V1(:,2:nspden)+DE*V2(:,2:nspden))
2987
2988 end subroutine prcrskerker2