TABLE OF CONTENTS
- ABINIT/bracketing
- ABINIT/brent
- ABINIT/cgpr
- ABINIT/dielmt
- ABINIT/dieltcel
- ABINIT/linmin
- ABINIT/m_prcref
- ABINIT/moddiel
- ABINIT/prcref
- ABINIT/prcref_PMA
- ABINIT/prcrskerker1
- ABINIT/prcrskerker2
ABINIT/bracketing [ 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)
SOURCE
3029 subroutine bracketing (nv1,nv2,dp_dum_v2dp,v,grad,a,x,b,fa,fx,fb) 3030 3031 !Arguments ------------------------------------ 3032 include "dummy_functions.inc" 3033 !scalars 3034 integer,intent(in) :: nv1,nv2 3035 real(dp),intent(inout) :: a,x 3036 real(dp),intent(out) :: b,fa,fb,fx 3037 !arrays 3038 real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2) 3039 3040 !Local variables------------------------------- 3041 !scalars 3042 real(dp),parameter :: maglimit=10000.0_dp 3043 real(dp) :: c,fu,q,r,u,ulim 3044 3045 ! ************************************************************************* 3046 3047 fa=dp_dum_v2dp(nv1,nv2,v(:,:)+(a*grad(:,:))) 3048 fx=dp_dum_v2dp(nv1,nv2,(x*grad(:,:))+v(:,:)) 3049 if(fx > fa) then 3050 c=a 3051 a=x 3052 x=c 3053 c=fa 3054 fa=fx 3055 fx=c 3056 end if 3057 b=x+gold*(x-a) 3058 fb=dp_dum_v2dp(nv1,nv2,(b*grad(:,:))+v(:,:)) 3059 do 3060 if (fx <= fb) return 3061 r=(x-a)*(fx-fb) 3062 q=(x-b)*(fx-fa) 3063 u=x-((x-b)*q-(x-a)*r)/(two*sign(max(abs(q-r),smallest_real),q-r)) 3064 ulim=x+maglimit*(b-x) 3065 if((x-u)*(u-b) > zero) then 3066 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3067 if(fu < fb) then 3068 a=x 3069 fa=fx 3070 x=u 3071 fx=fu 3072 return 3073 else if (fx < fu) then 3074 b=u 3075 fb=fu 3076 return 3077 end if 3078 u=b+gold*(b-x) 3079 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3080 else if((b-u)*(u-ulim) > zero) then 3081 fu=dp_dum_v2dp(nv1,nv2,u*grad(:,:)+v(:,:)) 3082 if(fu<fb) then 3083 x=b 3084 b=u 3085 u=b+gold*(b-x) 3086 fx=fb 3087 fb=fu 3088 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3089 end if 3090 else if((u-ulim)*(ulim-b) >= zero) then 3091 u=ulim 3092 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3093 else 3094 u=b+gold*(b-x) 3095 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3096 end if 3097 a=x 3098 x=b 3099 b=u 3100 fa=fx 3101 fx=fb 3102 fb=fu 3103 end do 3104 3105 end subroutine bracketing
ABINIT/brent [ 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 brent: dp_dum_vdp(v(:)+xmin*grad(:))
SIDE EFFECTS
grad(:): direction along which the minimization is performed v(:): starting and ending point of the minimization
SOURCE
3132 function brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,itmax,v,grad,ax,xx,bx,tol,xmin) 3133 3134 !Arguments ------------------------------------ 3135 include "dummy_functions.inc" 3136 !scalars 3137 integer,intent(in) :: itmax,nv1,nv2 3138 real(dp) :: brent 3139 real(dp),intent(in) :: ax,bx,tol,xx 3140 real(dp),intent(out) :: xmin 3141 !arrays 3142 real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2) 3143 3144 !Local variables------------------------------- 3145 !scalars 3146 integer :: iter 3147 real(dp) :: a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,vv,w 3148 real(dp) :: x,xm,zeps 3149 logical :: ok1,ok2,ok3,ok4 3150 3151 !************************************************************************ 3152 zeps=epsilon(ax*real(1e-2,dp)) 3153 a=min(ax,bx) 3154 b=max(ax,bx) 3155 vv=xx 3156 w=xx 3157 x=xx 3158 e=zero 3159 fx=dp_dum_v2dp(nv1,nv2,x*grad(:,:)+v(:,:)) 3160 fv=fx 3161 fw=fx 3162 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of 3163 !v(:,:)=v(:,:)+(grad(:,:)*x) 3164 !but for instance renormilizing the density if brent is used on a density... 3165 !vp(:,:) = v(:,:) 3166 !sub_dum_dp_v2dp_v2dp(x,grad(:,:),vp(:,:) 3167 !dx=dotproduct(v2dp_dum_v2dp(vp(:,:)),grad(:,:)) 3168 dx=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,v(:,:)+x*grad(:,:)),grad(:,:)) 3169 dv=dx 3170 dw=dx 3171 do iter=1,itmax 3172 xm=half*(a+b) 3173 tol1=tol*abs(x)+zeps 3174 tol2=two*tol1 3175 if(abs(x-xm) <= (tol2-half*(b-a))) then 3176 exit 3177 end if 3178 if(abs(e) > tol1) then 3179 d1=two*(b-a) 3180 d2=d1 3181 if(dw /= dx) d1=(w-x)*dx/(dx-dw) 3182 if(dv /= dx) d2=(vv-x)*dx/(dx-dv) 3183 u1=x+d1 3184 u2=x+d2 3185 ok1=((a-u1)*(u1-b)>zero).and.(dx*d1<=zero) 3186 ok2=((a-u2)*(u2-b)>zero).and.(dx*d2<=zero) 3187 olde=e 3188 e=d 3189 if(ok1.or.ok2) then 3190 if(ok1.and.ok2) then 3191 d=merge(d1,d2,abs(d1)<abs(d2)) 3192 else 3193 d=merge(d1,d2,ok1) 3194 end if 3195 if(abs(d)<=abs(half*olde)) then 3196 u=x+d 3197 if(((u-a)<tol2).or.((b-u)<tol2)) d=sign(tol1,xm-x) 3198 else 3199 e=merge(a,b,dx>=zero)-x 3200 d=half*e 3201 end if 3202 else 3203 e=merge(a,b,dx>=zero)-x 3204 d=half*e 3205 end if 3206 else 3207 e=merge(a,b,dx>=zero)-x 3208 d=half*e 3209 end if 3210 3211 if(abs(d) >=tol1)then 3212 u=x+d 3213 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3214 else 3215 u=x+sign(tol1,d) 3216 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3217 if(fu>fx) then 3218 exit 3219 end if 3220 end if 3221 du=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)),grad(:,:)) 3222 if(fu<=fx)then 3223 if(u>=x)then 3224 a=x 3225 else 3226 b=x 3227 end if 3228 vv=w 3229 fv=fw 3230 dv=dw 3231 w=x 3232 fw=fx 3233 dw=dx 3234 x=u 3235 dx=du 3236 fx=fu 3237 else 3238 if(u<x) then 3239 a=u 3240 else 3241 b=u 3242 end if 3243 ok3=(w==x).or.(fu.le.fw) 3244 ok4=(vv==w).or.(vv==x).or.(fu.lt.fv) 3245 if(ok3) then 3246 vv=w 3247 fv=fw 3248 dv=dw 3249 w=u 3250 fw=fu 3251 dw=du 3252 else if( ok4 ) then 3253 vv=u 3254 fv=fu 3255 dv=du 3256 end if 3257 end if 3258 end do 3259 xmin=x 3260 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of 3261 !v(:,:)=v(:,:)+(grad(:,:)*x) 3262 !but for instance renormilizing the density if brent is used on a density... 3263 call sub_dum_dp_v2dp_v2dp(nv1,nv2,x,grad(:,:),v(:,:)) 3264 brent=fx 3265 3266 end function brent
ABINIT/cgpr [ 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
SOURCE
2893 subroutine cgpr(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,dtol,itmax,v,fmin,delta) 2894 2895 !Arguments ------------------------------------ 2896 include "dummy_functions.inc" 2897 !scalars 2898 integer,intent(in) :: itmax,nv1,nv2 2899 real(dp),intent(in) :: dtol 2900 real(dp),intent(out) :: delta,fmin 2901 !arrays 2902 real(dp),intent(inout) :: v(nv1,nv2) 2903 2904 !Local variables------------------------------- 2905 !scalars 2906 integer :: iiter 2907 real(dp) :: fv,gam,gscal,gscal2,sto 2908 !arrays 2909 real(dp) :: grad0(nv1,nv2),grad1(nv1,nv2),grad2(nv1,nv2),grad3(nv1,nv2) 2910 !no_abirules 2911 2912 !************************************************************************ 2913 fv = dp_dum_v2dp(nv1,nv2,v(:,:)) 2914 grad0(:,:) = -v2dp_dum_v2dp(nv1,nv2,v(:,:)) 2915 grad1(:,:) = grad0(:,:) 2916 grad2(:,:) = grad0(:,:) 2917 do iiter=1,itmax 2918 call linmin(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,v,grad0,fmin) 2919 ! return if the min is reached 2920 sto=dtol*(abs(fmin)+abs(fv)+tol14) 2921 delta=abs(fv-fmin) 2922 delta=abs(delta) 2923 if((delta.lt.sto).or.(iiter==itmax)) then 2924 ! DEBUG 2925 ! write(std_out,*) 'cgpr (01cg) : stop cond for cgpr:',sto,'delta:',delta,'fv:',fv,'fmin:',fmin 2926 ! ENDDEBUG 2927 return 2928 end if 2929 ! a new step 2930 fv=fmin 2931 grad0(:,:)=v2dp_dum_v2dp(nv1,nv2,v(:,:)) 2932 gscal=dotproduct(nv1,nv2,grad1(:,:),grad1(:,:)) 2933 grad3(:,:)=grad0(:,:)+grad1(:,:) 2934 gscal2=dotproduct(nv1,nv2,grad3(:,:),grad0(:,:)) 2935 gam=gscal2/gscal 2936 grad1(:,:)=-grad0(:,:) 2937 grad2(:,:)=grad1(:,:)+gam*grad2(:,:) 2938 grad0(:,:)=grad2(:,:) 2939 ! DEBUG 2940 ! write(std_out,*) 'cgpr (01cg) :=================================================================================' 2941 ! write(std_out,*) 'cgpr (01cg) : step',iiter,'delta:',delta ,'fv',fv,'fmin',fmin 2942 ! write(std_out,*) 'cgpr (01cg) :=================================================================================' 2943 ! ENDDEBUG 2944 end do 2945 2946 end subroutine cgpr
ABINIT/dielmt [ 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)
SOURCE
1560 subroutine dielmt(dielinv,gmet,kg_diel,npwdiel,nspden,occopt,prtvol,susmat) 1561 1562 !Arguments ------------------------------------ 1563 !scalars 1564 integer,intent(in) :: npwdiel,nspden,occopt,prtvol 1565 !arrays 1566 integer,intent(in) :: kg_diel(3,npwdiel) 1567 real(dp),intent(in) :: gmet(3,3) 1568 real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 1569 real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden) 1570 1571 !Local variables------------------------------- 1572 !scalars 1573 integer :: ieig,ier,ii,index,ipw,ipw1,ipw2,isp,jj,npwsp 1574 real(dp) :: ai1,ai2,ar1,ar2,eiginv,gfact,gfactinv,kg_red1,kg_red2,kg_red3,gsquar 1575 real(dp) :: tpisq 1576 character(len=500) :: message 1577 !arrays 1578 real(dp) :: tsec(2) 1579 real(dp),allocatable :: dielh(:),dielmat(:,:,:,:,:),dielvec(:,:,:) 1580 real(dp),allocatable :: eig_diel(:),zhpev1(:,:),zhpev2(:) 1581 !no_abirules 1582 !integer :: ipw3 1583 !real(dp) :: elementi,elementr 1584 1585 ! ************************************************************************* 1586 1587 !DEBUG 1588 !write(std_out,*)' dielmt : enter ' 1589 !if(.true.)stop 1590 !ENDDEBUG 1591 1592 !tpisq is (2 Pi) **2: 1593 tpisq=(two_pi)**2 1594 1595 call timab(90,1,tsec) 1596 1597 !-Compute now the hermitian dielectric matrix------------------------------ 1598 !Following remarks are only valid within RPA approximation (Kxc=0): 1599 1600 !for the spin-unpolarized case, 1 - 4pi (1/G) chi0(G,Gp) (1/Gp) 1601 1602 !for the spin-polarized case, 1603 !( 1 0 ) - 4pi ( 1/G 1/G ) ( chi0 upup chi0 updn ) ( 1/Gp 1/Gp ) 1604 !( 0 1 ) ( 1/G 1/G ) ( chi0 dnup chi0 dndn ) ( 1/Gp 1/Gp ) 1605 !which is equal to 1606 !( 1 0 ) - 4pi (1/G 0 ) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) (1/Gp 0 ) 1607 !( 0 1 ) ( 0 1/G) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) ( 0 1/Gp) 1608 !So, if spin-polarized, sum all spin contributions 1609 !Note: chi0 updn = chi0 dnup = zero for non-metallic systems 1610 1611 !In the case of non-collinear magnetism, within RPA, this is the same because: 1612 !chi0_(s1,s2),(s3,s4) = delta_s1,s2 * delta_s3,s4 * chi0_(s1,s1),(s3,s3) 1613 !Only chi_upup,upup, chi_dndn,dndn, chi_upup,dndn and chi_dndn,upup 1614 !have to be taken into account (stored, susmat(:,ipw1,1:2,ipw2,1:2) 1615 1616 ABI_MALLOC(dielmat,(2,npwdiel,min(nspden,2),npwdiel,min(nspden,2))) 1617 1618 if(nspden/=1)then 1619 if (occopt<3) then 1620 do ipw2=1,npwdiel 1621 do ipw1=1,npwdiel 1622 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2) 1623 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2) 1624 end do 1625 end do 1626 else 1627 do ipw2=1,npwdiel 1628 do ipw1=1,npwdiel 1629 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) 1630 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) 1631 end do 1632 end do 1633 end if 1634 else 1635 do ipw2=1,npwdiel 1636 do ipw1=1,npwdiel 1637 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1) 1638 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1) 1639 end do 1640 end do 1641 end if 1642 !Compute 1/G factors and include them in the dielectric matrix 1643 do ipw1=1,npwdiel 1644 kg_red1=dble(kg_diel(1,ipw1)) 1645 kg_red2=dble(kg_diel(2,ipw1)) 1646 kg_red3=dble(kg_diel(3,ipw1)) 1647 gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 & 1648 & +two*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + & 1649 & gmet(2,3)*kg_red2*kg_red3) ) 1650 ! Distinguish G=0 from other elements 1651 if(gsquar>tol12)then 1652 ! !$ gfact=\sqrt (4.0_dp \pi/gsquar/dble(nspden))$ 1653 gfact=sqrt(four_pi/gsquar) 1654 do ipw2=1,npwdiel 1655 ! Must multiply both rows and columns, and also changes the sign 1656 dielmat(1,ipw2,1,ipw1,1)=-dielmat(1,ipw2,1,ipw1,1)*gfact 1657 dielmat(2,ipw2,1,ipw1,1)=-dielmat(2,ipw2,1,ipw1,1)*gfact 1658 dielmat(1,ipw1,1,ipw2,1)= dielmat(1,ipw1,1,ipw2,1)*gfact 1659 dielmat(2,ipw1,1,ipw2,1)= dielmat(2,ipw1,1,ipw2,1)*gfact 1660 end do 1661 else 1662 ! Zero the G=0 elements, head and wings 1663 do ipw2=1,npwdiel 1664 dielmat(1,ipw2,1,ipw1,1)=zero 1665 dielmat(2,ipw2,1,ipw1,1)=zero 1666 dielmat(1,ipw1,1,ipw2,1)=zero 1667 dielmat(2,ipw1,1,ipw2,1)=zero 1668 end do 1669 end if 1670 end do 1671 1672 !Complete the matrix in the spin-polarized case 1673 !should this be nspden==2?? 1674 if(nspden/=1)then 1675 do ipw1=1,npwdiel 1676 do ipw2=1,npwdiel 1677 dielmat(1,ipw1,1,ipw2,2)=dielmat(1,ipw1,1,ipw2,1) 1678 dielmat(2,ipw1,1,ipw2,2)=dielmat(2,ipw1,1,ipw2,1) 1679 dielmat(1,ipw1,2,ipw2,1)=dielmat(1,ipw1,1,ipw2,1) 1680 dielmat(2,ipw1,2,ipw2,1)=dielmat(2,ipw1,1,ipw2,1) 1681 dielmat(1,ipw1,2,ipw2,2)=dielmat(1,ipw1,1,ipw2,1) 1682 dielmat(2,ipw1,2,ipw2,2)=dielmat(2,ipw1,1,ipw2,1) 1683 end do 1684 end do 1685 end if 1686 1687 !DEBUG 1688 !write(std_out,*)' dielmt : make dielmat equal to identity matrix ' 1689 !do ipw1=1,npwdiel 1690 !do ipw2=1,npwdiel 1691 !dielmat(1,ipw1,1,ipw2,1)=0.0_dp 1692 !dielmat(2,ipw1,1,ipw2,1)=0.0_dp 1693 !end do 1694 !end do 1695 !ENDDEBUG 1696 1697 !Add the diagonal part 1698 do isp=1,min(nspden,2) 1699 do ipw=1,npwdiel 1700 dielmat(1,ipw,isp,ipw,isp)=one+dielmat(1,ipw,isp,ipw,isp) 1701 end do 1702 end do 1703 1704 !-The hermitian dielectric matrix is computed ------------------------------ 1705 !-Now, diagonalize it ------------------------------------------------------ 1706 1707 !In RPA, everything is projected on the spin-symmetrized 1708 !space. This was coded here (for the time being). 1709 1710 !Diagonalize the hermitian dielectric matrix 1711 1712 !npwsp=npwdiel*nspden 1713 npwsp=npwdiel 1714 1715 ABI_MALLOC(dielh,(npwsp*(npwsp+1))) 1716 ABI_MALLOC(dielvec,(2,npwsp,npwsp)) 1717 ABI_MALLOC(eig_diel,(npwsp)) 1718 ABI_MALLOC(zhpev1,(2,2*npwsp-1)) 1719 ABI_MALLOC(zhpev2,(3*npwsp-2)) 1720 ier=0 1721 !Store the dielectric matrix in proper mode before calling zhpev 1722 index=1 1723 do ii=1,npwdiel 1724 do jj=1,ii 1725 dielh(index )=dielmat(1,jj,1,ii,1) 1726 dielh(index+1)=dielmat(2,jj,1,ii,1) 1727 index=index+2 1728 end do 1729 end do 1730 !If spin-polarized and non RPA, need to store other parts of the matrix 1731 !if(nspden/=1)then 1732 !do ii=1,npwdiel 1733 !Here, spin-flip contribution 1734 !do jj=1,npwdiel 1735 !dielh(index )=dielmat(1,jj,1,ii,2) 1736 !dielh(index+1)=dielmat(2,jj,1,ii,2) 1737 !index=index+2 1738 !end do 1739 !Here spin down-spin down upper matrix 1740 !do jj=1,ii 1741 !dielh(index )=dielmat(1,jj,2,ii,2) 1742 !dielh(index+1)=dielmat(2,jj,2,ii,2) 1743 !index=index+2 1744 !end do 1745 !end do 1746 !end if 1747 1748 call ZHPEV ('V','U',npwsp,dielh,eig_diel,dielvec,npwdiel,zhpev1,& 1749 & zhpev2,ier) 1750 ABI_FREE(zhpev1) 1751 ABI_FREE(zhpev2) 1752 1753 if(prtvol>=10)then 1754 write(message, '(a,a,a,5es12.4)' )ch10,& 1755 & ' Five largest eigenvalues of the hermitian RPA dielectric matrix:',& 1756 & ch10,eig_diel(npwdiel:npwdiel-4:-1) 1757 call wrtout(ab_out,message,'COLL') 1758 end if 1759 1760 write(message, '(a,a)' )ch10,& 1761 & ' dielmt : 15 largest eigenvalues of the hermitian RPA dielectric matrix' 1762 call wrtout(std_out,message,'COLL') 1763 write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(npwdiel:npwdiel-4:-1) 1764 call wrtout(std_out,message,'COLL') 1765 write(message, '(a,5es12.5)' )' 6-10 :',eig_diel(npwdiel-5:npwdiel-9:-1) 1766 call wrtout(std_out,message,'COLL') 1767 write(message, '(a,5es12.5)' )' 11-15:',eig_diel(npwdiel-10:npwdiel-14:-1) 1768 call wrtout(std_out,message,'COLL') 1769 write(message, '(a,a)' )ch10,& 1770 & ' dielmt : 5 smallest eigenvalues of the hermitian RPA dielectric matrix' 1771 call wrtout(std_out,message,'COLL') 1772 write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(1:5) 1773 call wrtout(std_out,message,'COLL') 1774 1775 !Invert the hermitian dielectric matrix, 1776 !Should use a BLAS call ! 1777 do ipw2=1,npwdiel 1778 do ipw1=ipw2,npwdiel 1779 dielinv(1,ipw1,1,ipw2,1)=zero 1780 dielinv(2,ipw1,1,ipw2,1)=zero 1781 end do 1782 end do 1783 do ieig=1,npwdiel 1784 eiginv=one/eig_diel(ieig) 1785 do ipw2=1,npwdiel 1786 do ipw1=ipw2,npwdiel 1787 ar1=dielvec(1,ipw1,ieig) 1788 ai1=dielvec(2,ipw1,ieig) 1789 ar2=dielvec(1,ipw2,ieig) 1790 ai2=dielvec(2,ipw2,ieig) 1791 dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)+& 1792 & (ar1*ar2+ai1*ai2)*eiginv 1793 dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)+& 1794 & (ai1*ar2-ar1*ai2)*eiginv 1795 end do 1796 end do 1797 end do 1798 do ipw2=1,npwdiel-1 1799 do ipw1=ipw2+1,npwdiel 1800 dielinv(1,ipw2,1,ipw1,1)= dielinv(1,ipw1,1,ipw2,1) 1801 dielinv(2,ipw2,1,ipw1,1)=-dielinv(2,ipw1,1,ipw2,1) 1802 end do 1803 end do 1804 1805 ABI_FREE(dielh) 1806 ABI_FREE(dielvec) 1807 ABI_FREE(eig_diel) 1808 1809 !DEBUG 1810 !Checks whether the inverse of the hermitian dielectric matrix 1811 !has been correctly generated 1812 !do ipw1=1,npwdiel 1813 !do ipw2=1,npwdiel 1814 !elementr=0.0_dp 1815 !elementi=0.0_dp 1816 !do ipw3=1,npwdiel 1817 !elementr=elementr+dielinv(1,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1)& 1818 !& -dielinv(2,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1) 1819 !elementi=elementi+dielinv(1,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1)& 1820 !& +dielinv(2,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1) 1821 !end do 1822 !if(elementr**2+elementi**2 > 1.0d-12)then 1823 !if( ipw1 /= ipw2 .or. & 1824 !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 1825 !write(std_out,*)' dielmt : the inversion procedure is not correct ' 1826 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 1827 !write(std_out,*)' elementr,elementi=',elementr,elementi 1828 !stop 1829 !end if 1830 !end if 1831 !end do 1832 !end do 1833 !write(std_out,*)'dielmt : matrix has been inverted successfully ' 1834 !stop 1835 !ENDDEBUG 1836 1837 !Then get the inverse of the asymmetric 1838 !dielectric matrix, as required for the preconditioning. 1839 1840 !Inverse of the dielectric matrix : ( 1 - 4pi (1/G^2) chi0(G,Gp) )^(-1) 1841 !In dielinv there is now (1 - 4pi (1/G) chi0(G,Gp) (1/Gp) )^(-1) 1842 !So, evaluate dielinv_after(G,Gp) = 1843 !(4pi/G^2)^(1/2) dielinv_before(G,Gp) (4pi/Gp^2)^(-1/2) 1844 !In RPA, can focus on the spin-averaged quantities 1845 do ipw1=1,npwdiel 1846 kg_red1=dble(kg_diel(1,ipw1)) 1847 kg_red2=dble(kg_diel(2,ipw1)) 1848 kg_red3=dble(kg_diel(3,ipw1)) 1849 gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 & 1850 & +two*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + & 1851 & gmet(2,3)*kg_red2*kg_red3) ) 1852 ! Distinguish G=0 from other elements 1853 if(gsquar>tol12)then 1854 gfact=sqrt(four_pi/gsquar) 1855 gfactinv=one/gfact 1856 do ipw2=1,npwdiel 1857 ! Must multiply both rows and columns 1858 dielinv(1,ipw2,1,ipw1,1)=dielinv(1,ipw2,1,ipw1,1)*gfactinv 1859 dielinv(2,ipw2,1,ipw1,1)=dielinv(2,ipw2,1,ipw1,1)*gfactinv 1860 dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)*gfact 1861 dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)*gfact 1862 end do 1863 else 1864 ! Zero the G=0 elements, head 1865 do ipw2=1,npwdiel 1866 if (ipw2/=ipw1) dielinv(1:2,ipw1,1,ipw2,1)=zero 1867 end do 1868 end if 1869 end do 1870 1871 ABI_FREE(dielmat) 1872 1873 call timab(90,2,tsec) 1874 1875 end subroutine dielmt
ABINIT/dieltcel [ 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.
SOURCE
1915 subroutine dieltcel(dielinv,gmet,kg_diel,kxc,nfft,ngfft,nkxc,npwdiel,nspden,occopt,option,prtvol,susmat) 1916 1917 !Arguments ------------------------------------ 1918 !scalars 1919 integer,intent(in) :: nfft,nkxc,npwdiel,nspden,occopt,option 1920 integer,intent(in) :: prtvol 1921 !arrays 1922 integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18) 1923 real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc) 1924 real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 1925 real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden) 1926 1927 !Local variables------------------------------- 1928 !scalars 1929 integer :: i1,i2,i3,ieig,ier,ifft,ii,index,ipw0,ipw1,ipw2,ispden,j1 1930 integer :: j2,j3,jj,k1,k2,k3,n1,n2,n3 1931 real(dp) :: ai,ai2,ar,ar2,eiginv,kg_red1,kg_red2,kg_red3,gsquar,si 1932 real(dp) :: sr,tpisq 1933 character(len=500) :: message 1934 type(MPI_type) :: mpi_enreg_seq 1935 !arrays 1936 real(dp) :: tsec(2) 1937 real(dp),allocatable :: eig_msusinvsqr(:),eig_msussqr(:) 1938 real(dp),allocatable :: eig_sus(:),eig_sym(:),invsqrsus(:,:,:,:,:) 1939 real(dp),allocatable :: khxc(:,:,:,:,:),kxcg(:,:),sqrsus(:,:,:,:,:),sush(:) 1940 real(dp),allocatable :: susvec(:,:,:),symdielmat(:,:,:,:,:),symh(:) 1941 real(dp),allocatable :: symvec(:,:,:,:,:),wkxc(:),work(:,:,:,:,:) 1942 real(dp),allocatable :: work2(:,:,:,:,:),zhpev1(:,:),zhpev2(:) 1943 !no_abirules 1944 !integer :: ipw3 1945 !real(dp) :: elementi,elementr 1946 !DEBUG 1947 !Used to moderate divergence effect near rho=0 of the Kxc 1948 !this limit value is truly empirical (exprmt on small Sr cell). 1949 !real(dp) :: kxc_min=-200.0 1950 !ENDDEBUG 1951 1952 ! ************************************************************************* 1953 1954 call timab(96,1,tsec) 1955 1956 !tpisq is (2 Pi) **2: 1957 tpisq=(two_pi)**2 1958 1959 if(nspden/=1 .and. (occopt>=3 .and. occopt<=8) )then 1960 write(message, '(a,a,a)' )& 1961 & 'In the present version of the code, one cannot produce',ch10,& 1962 & 'the dielectric matrix in the metallic, spin-polarized case.' 1963 ABI_BUG(message) 1964 end if 1965 1966 if(nspden==4)then 1967 write(message,'(a,a,a)')& 1968 & 'In the present version of the code, one cannot produce',ch10,& 1969 & 'the dielectric matrix in the non-collinear spin-polarized case.' 1970 ABI_ERROR(message) 1971 end if 1972 1973 1974 !-Diagonalize the susceptibility matrix 1975 1976 ABI_MALLOC(sush,(npwdiel*(npwdiel+1))) 1977 ABI_MALLOC(susvec,(2,npwdiel,npwdiel)) 1978 ABI_MALLOC(eig_msusinvsqr,(npwdiel)) 1979 ABI_MALLOC(eig_msussqr,(npwdiel)) 1980 ABI_MALLOC(eig_sus,(npwdiel)) 1981 ABI_MALLOC(zhpev1,(2,2*npwdiel-1)) 1982 ABI_MALLOC(zhpev2,(3*npwdiel-2)) 1983 ABI_MALLOC(work,(2,npwdiel,nspden,npwdiel,nspden)) 1984 ABI_MALLOC(work2,(2,npwdiel,nspden,npwdiel,nspden)) 1985 ABI_MALLOC(sqrsus,(2,npwdiel,nspden,npwdiel,nspden)) 1986 ABI_MALLOC(invsqrsus,(2,npwdiel,nspden,npwdiel,nspden)) 1987 1988 !At some time, should take care of different spin channels 1989 do ispden=1,nspden 1990 1991 if(nspden/=1)then 1992 ABI_ERROR('dieltcel : stop, nspden/=1') 1993 end if 1994 1995 ! Store the susceptibility matrix in proper mode before calling zhpev 1996 index=1 1997 do ii=1,npwdiel 1998 do jj=1,ii 1999 sush(index )=susmat(1,jj,1,ii,1) 2000 sush(index+1)=susmat(2,jj,1,ii,1) 2001 index=index+2 2002 end do 2003 end do 2004 2005 ier=0 2006 call ZHPEV ('V','U',npwdiel,sush,eig_sus,susvec,npwdiel,zhpev1,zhpev2,ier) 2007 2008 ! DEBUG 2009 ! write(std_out,*)' dieltcel : print eigenvalues of the susceptibility matrix' 2010 ! do ii=1,npwdiel 2011 ! write(std_out,'(i5,es16.6)' )ii,eig_sus(ii) 2012 ! end do 2013 ! ENDDEBUG 2014 2015 do ii=1,npwdiel 2016 if(-eig_sus(ii)>1.d-12)then 2017 eig_msussqr(ii)=sqrt(-eig_sus(ii)) 2018 eig_msusinvsqr(ii)=1._dp/eig_msussqr(ii) 2019 else if(-eig_sus(ii)< -1.d-12)then 2020 message = "Found positive eigenvalue of susceptibility matrix." 2021 ABI_BUG(message) 2022 else 2023 ! Set the eigenvalue corresponding to a constant potential change to 1, 2024 ! while it will be set to zero in Khx. 2025 eig_msussqr(ii)=1._dp 2026 eig_msusinvsqr(ii)=1._dp 2027 end if 2028 end do 2029 2030 ! Compute square root of minus susceptibility matrix 2031 ! and inverse square root of minus susceptibility matrix 2032 do ii=1,npwdiel 2033 work(:,:,1,ii,1)=susvec(:,:,ii)*eig_msussqr(ii) 2034 work2(:,:,1,ii,1)=susvec(:,:,ii)*eig_msusinvsqr(ii) 2035 end do 2036 do ipw2=1,npwdiel 2037 do ipw1=ipw2,npwdiel 2038 ar=0._dp ; ai=0._dp ; ar2=0._dp ; ai2=0._dp 2039 do ii=1,npwdiel 2040 sr=susvec(1,ipw2,ii) ; si=susvec(2,ipw2,ii) 2041 ar =ar +work(1,ipw1,1,ii,1)*sr +work(2,ipw1,1,ii,1)*si 2042 ai =ai +work(2,ipw1,1,ii,1)*sr -work(1,ipw1,1,ii,1)*si 2043 ar2=ar2 +work2(1,ipw1,1,ii,1)*sr +work2(2,ipw1,1,ii,1)*si 2044 ai2=ai2 +work2(2,ipw1,1,ii,1)*sr -work2(1,ipw1,1,ii,1)*si 2045 end do 2046 sqrsus(1,ipw1,1,ipw2,1)=ar 2047 sqrsus(2,ipw1,1,ipw2,1)=ai 2048 invsqrsus(1,ipw1,1,ipw2,1)=ar2 2049 invsqrsus(2,ipw1,1,ipw2,1)=ai2 2050 if(ipw1/=ipw2)then 2051 sqrsus(1,ipw2,1,ipw1,1)=ar 2052 sqrsus(2,ipw2,1,ipw1,1)=-ai 2053 invsqrsus(1,ipw2,1,ipw1,1)=ar2 2054 invsqrsus(2,ipw2,1,ipw1,1)=-ai2 2055 end if 2056 end do 2057 end do 2058 2059 ! DEBUG 2060 ! Checks whether sqrsus and invsqrsus are inverse of each other. 2061 ! do ipw1=1,npwdiel 2062 ! do ipw2=1,npwdiel 2063 ! elementr=0.0_dp 2064 ! elementi=0.0_dp 2065 ! do ipw3=1,npwdiel 2066 ! elementr=elementr+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1)& 2067 ! & -sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1) 2068 ! elementi=elementi+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1)& 2069 ! & +sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1) 2070 ! end do 2071 ! if(elementr**2+elementi**2 > 1.0d-12)then 2072 ! if( ipw1 /= ipw2 .or. & 2073 ! & ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 2074 ! write(std_out,*)' dieltcel : sqrsus and invsqrsus are not (pseudo)',& 2075 ! & 'inverse of each other' 2076 ! write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 2077 ! write(std_out,*)' elementr,elementi=',elementr,elementi 2078 ! stop 2079 ! end if 2080 ! end if 2081 ! end do 2082 ! end do 2083 ! ENDDEBUG 2084 2085 ! End loop over spins 2086 end do 2087 2088 ABI_FREE(eig_msusinvsqr) 2089 ABI_FREE(eig_msussqr) 2090 ABI_FREE(eig_sus) 2091 ABI_FREE(sush) 2092 ABI_FREE(susvec) 2093 2094 !-Compute the Hxc kernel 2095 2096 ABI_MALLOC(khxc,(2,npwdiel,nspden,npwdiel,nspden)) 2097 ABI_MALLOC(symdielmat,(2,npwdiel,nspden,npwdiel,nspden)) 2098 2099 khxc(:,:,:,:,:)=0.0_dp 2100 2101 !Compute Hartree kernel 2102 do ipw1=1,npwdiel 2103 kg_red1=dble(kg_diel(1,ipw1)) 2104 kg_red2=dble(kg_diel(2,ipw1)) 2105 kg_red3=dble(kg_diel(3,ipw1)) 2106 gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 & 2107 & +2.0_dp*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + & 2108 & gmet(2,3)*kg_red2*kg_red3) ) 2109 ! Distinguish G=0 from other elements 2110 if(gsquar>1.0d-12)then 2111 khxc(1,ipw1,1,ipw1,1)= 4.0_dp*pi/gsquar 2112 else 2113 ! G=0 2114 ipw0=ipw1 2115 end if 2116 end do 2117 2118 !Eventually add the xc part 2119 if(option>=2)then 2120 2121 ABI_MALLOC(wkxc,(nfft)) 2122 ABI_MALLOC(kxcg,(2,nfft)) 2123 wkxc(:)=kxc(:,1) 2124 ! DEBUG 2125 ! Used to moderate divergenc effect near rho=0 of the Kxc (see above). 2126 ! wkxc(:)=merge(kxc(:,1), kxc_min, kxc(:,1) > kxc_min) 2127 ! ENDDEBUG 2128 call initmpi_seq(mpi_enreg_seq) 2129 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all') 2130 call fourdp(1,kxcg,wkxc,-1,mpi_enreg_seq,nfft,1,ngfft,0) ! trsfrm R to G 2131 call destroy_mpi_enreg(mpi_enreg_seq) 2132 2133 ! Compute difference in G vectors 2134 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 2135 do ipw2=1,npwdiel 2136 if(ipw2/=ipw0)then 2137 2138 j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2) 2139 ! Fills diagonal 2140 khxc(1,ipw2,1,ipw2,1)=khxc(1,ipw2,1,ipw2,1)+kxcg(1,1) 2141 khxc(2,ipw2,1,ipw2,1)=khxc(2,ipw2,1,ipw2,1)+kxcg(2,1) 2142 2143 if(ipw2/=npwdiel)then 2144 ! Fills off-diagonal part of the matrix, except G=0 2145 do ipw1=ipw2+1,npwdiel 2146 if(ipw1/=ipw0)then 2147 i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1) 2148 ! Use of two mod calls handles both i1-j1>=ndiel1 AND i1-j1<0 2149 k1=mod(n1+mod(i1-j1,n1),n1) 2150 k2=mod(n2+mod(i2-j2,n2),n2) 2151 k3=mod(n3+mod(i3-j3,n3),n3) 2152 ifft=k1+1+n1*(k2+n2*k3) 2153 ! The signs of imaginary contributions have been checked 2154 khxc(1,ipw1,1,ipw2,1)=kxcg(1,ifft) 2155 khxc(2,ipw1,1,ipw2,1)=kxcg(2,ifft) 2156 khxc(1,ipw2,1,ipw1,1)=kxcg(1,ifft) 2157 khxc(2,ipw2,1,ipw1,1)=-kxcg(2,ifft) 2158 end if 2159 end do 2160 end if 2161 2162 end if 2163 end do 2164 2165 ABI_FREE(wkxc) 2166 ABI_FREE(kxcg) 2167 2168 ! Endif option 2 2169 end if 2170 2171 !Now, get the symmetric dielectric matrix 2172 !Premultiplication by square root of minus susceptibility matrix 2173 do ipw2=1,npwdiel 2174 do ipw1=1,npwdiel 2175 ar=0._dp ; ai=0._dp 2176 do ii=1,npwdiel 2177 ar=ar+sqrsus(1,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) & 2178 & -sqrsus(2,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1) 2179 ai=ai+sqrsus(2,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) & 2180 & +sqrsus(1,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1) 2181 end do 2182 work(1,ipw1,1,ipw2,1)=ar 2183 work(2,ipw1,1,ipw2,1)=ai 2184 end do 2185 end do 2186 !Postmultiplication by square root of minus susceptibility matrix 2187 do ipw2=1,npwdiel 2188 ! do ipw1=ipw2,npwdiel 2189 do ipw1=1,npwdiel 2190 ar=0._dp ; ai=0._dp 2191 do ii=1,npwdiel 2192 ar=ar+work(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2193 & -work(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2194 ai=ai+work(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2195 & +work(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2196 end do 2197 symdielmat(1,ipw1,1,ipw2,1)=ar 2198 symdielmat(2,ipw1,1,ipw2,1)=ai 2199 ! if(ipw1/=ipw2)then 2200 ! symdielmat(1,ipw2,1,ipw1,1)=ar 2201 ! symdielmat(2,ipw2,1,ipw1,1)=-ai 2202 ! end if 2203 end do 2204 ! Add the unity matrix 2205 symdielmat(1,ipw2,1,ipw2,1)=1._dp+symdielmat(1,ipw2,1,ipw2,1) 2206 end do 2207 2208 ABI_FREE(khxc) 2209 2210 ABI_MALLOC(symh,(npwdiel*(npwdiel+1))) 2211 ABI_MALLOC(symvec,(2,npwdiel,nspden,npwdiel,nspden)) 2212 ABI_MALLOC(eig_sym,(npwdiel)) 2213 2214 !Store the symmetrized dielectric matrix in proper mode before calling zhpev 2215 index=1 2216 do ii=1,npwdiel 2217 do jj=1,ii 2218 symh(index )=symdielmat(1,jj,1,ii,1) 2219 symh(index+1)=symdielmat(2,jj,1,ii,1) 2220 index=index+2 2221 end do 2222 end do 2223 2224 ier=0 2225 call ZHPEV ('V','U',npwdiel,symh,eig_sym,symvec,npwdiel,zhpev1,& 2226 & zhpev2,ier) 2227 2228 if(prtvol>=10)then 2229 write(message, '(a,a,a,5es12.4)' )ch10,& 2230 & ' Five largest eigenvalues of the symmetrized dielectric matrix:',& 2231 & ch10,eig_sym(npwdiel:npwdiel-4:-1) 2232 call wrtout(ab_out,message,'COLL') 2233 end if 2234 2235 write(message,'(2a)')ch10,' dieltcel : 15 largest eigenvalues of the symmetrized dielectric matrix' 2236 call wrtout(std_out,message,'COLL') 2237 write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(npwdiel:npwdiel-4:-1) 2238 call wrtout(std_out,message,'COLL') 2239 write(message, '(a,5es12.5)' )' 6-10 :',eig_sym(npwdiel-5:npwdiel-9:-1) 2240 call wrtout(std_out,message,'COLL') 2241 write(message, '(a,5es12.5)' )' 11-15:',eig_sym(npwdiel-10:npwdiel-14:-1) 2242 call wrtout(std_out,message,'COLL') 2243 write(message, '(2a)' )ch10,' dieltcel : 5 smallest eigenvalues of the symmetrized dielectric matrix' 2244 call wrtout(std_out,message,'COLL') 2245 write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(1:5) 2246 call wrtout(std_out,message,'COLL') 2247 2248 !Invert the hermitian dielectric matrix, 2249 work(:,:,:,:,:)=0.0_dp 2250 do ieig=1,npwdiel 2251 eiginv=1.0_dp/eig_sym(ieig) 2252 do ipw2=1,npwdiel 2253 ! do ipw1=ipw2,npwdiel 2254 do ipw1=1,npwdiel 2255 work(1,ipw1,1,ipw2,1)=work(1,ipw1,1,ipw2,1)+& 2256 & (symvec(1,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)+ & 2257 & symvec(2,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv 2258 work(2,ipw1,1,ipw2,1)=work(2,ipw1,1,ipw2,1)+& 2259 & (symvec(2,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)- & 2260 & symvec(1,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv 2261 end do 2262 end do 2263 end do 2264 !if(npwdiel>1)then 2265 !do ipw2=2,npwdiel 2266 !do ipw1=1,ipw2-1 2267 !work(1,ipw1,1,ipw2,1)= work(1,ipw2,1,ipw1,1) 2268 !work(2,ipw1,1,ipw2,1)=-work(2,ipw2,1,ipw1,1) 2269 !end do 2270 !end do 2271 !end if 2272 2273 ABI_FREE(eig_sym) 2274 ABI_FREE(symh) 2275 ABI_FREE(symvec) 2276 2277 !DEBUG 2278 !Checks whether the inverse of the symmetric dielectric matrix 2279 !has been correctly generated 2280 !do ipw1=1,npwdiel 2281 !do ipw2=1,npwdiel 2282 !elementr=0.0_dp 2283 !elementi=0.0_dp 2284 !do ipw3=1,npwdiel 2285 !elementr=elementr+work(1,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1)& 2286 !& -work(2,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1) 2287 !elementi=elementi+work(1,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1)& 2288 !& +work(2,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1) 2289 !end do 2290 !if(elementr**2+elementi**2 > 1.0d-12)then 2291 !if( ipw1 /= ipw2 .or. & 2292 !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 2293 !write(std_out,*)' dieltcel : the inversion procedure is not correct ' 2294 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 2295 !write(std_out,*)' elementr,elementi=',elementr,elementi 2296 !stop 2297 !end if 2298 !end if 2299 !end do 2300 !end do 2301 !write(std_out,*)'dieltcel : matrix has been inverted successfully ' 2302 !ENDDEBUG 2303 2304 ABI_FREE(symdielmat) 2305 2306 !Then get the inverse of the asymmetric 2307 !dielectric matrix, as required for the preconditioning. 2308 !Premultiplication by square root of minus susceptibility matrix 2309 do ipw2=1,npwdiel 2310 do ipw1=1,npwdiel 2311 ar=0._dp ; ai=0._dp 2312 do ii=1,npwdiel 2313 ar=ar+invsqrsus(1,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) & 2314 & -invsqrsus(2,ipw1,1,ii,1)*work(2,ii,1,ipw2,1) 2315 ai=ai+invsqrsus(2,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) & 2316 & +invsqrsus(1,ipw1,1,ii,1)*work(2,ii,1,ipw2,1) 2317 end do 2318 work2(1,ipw1,1,ipw2,1)=ar 2319 work2(2,ipw1,1,ipw2,1)=ai 2320 end do 2321 end do 2322 !Postmultiplication by square root of minus susceptibility matrix 2323 do ipw2=1,npwdiel 2324 do ipw1=1,npwdiel 2325 ar=0._dp ; ai=0._dp 2326 do ii=1,npwdiel 2327 ar=ar+work2(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2328 & -work2(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2329 ai=ai+work2(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2330 & +work2(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2331 end do 2332 dielinv(1,ipw1,1,ipw2,1)=ar 2333 dielinv(2,ipw1,1,ipw2,1)=ai 2334 end do 2335 end do 2336 2337 ABI_FREE(invsqrsus) 2338 ABI_FREE(sqrsus) 2339 ABI_FREE(work) 2340 ABI_FREE(work2) 2341 ABI_FREE(zhpev1) 2342 ABI_FREE(zhpev2) 2343 2344 call timab(96,2,tsec) 2345 2346 end subroutine dieltcel
ABINIT/linmin [ Functions ]
NAME
linmin
FUNCTION
minimizes a function along a gradient line: first bracket the minimum then perform the minimization
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, MT) This file is distributed under the terms of the 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
SOURCE
2977 subroutine linmin(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,v,grad,fmin) 2978 2979 !Arguments ------------------------------------ 2980 include "dummy_functions.inc" 2981 !scalars 2982 integer,intent(in) :: nv1,nv2 2983 real(dp),intent(out) :: fmin 2984 !arrays 2985 real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2) 2986 2987 !Local variables------------------------------- 2988 !scalars 2989 real(dp),parameter :: maglimit=10000.0_dp,tol=tol8*tol8*tol8 2990 real(dp) :: a,b,fa,fb,fx,x,xmin 2991 !no_abirules 2992 2993 !************************************************************************ 2994 a=zero 2995 x=ninth*real(1e-4,dp) 2996 call bracketing (nv1,nv2,dp_dum_v2dp,v,grad,a,x,b,fa,fx,fb) 2997 !DEBUG 2998 !write(std_out,*) 'linmin (01cg) : linmin after bracketing' 2999 !write(std_out,*) 'linmin (01cg) : point',a,x,b,'value',fa,fx,fb 3000 !ENDDEBUG 3001 fmin =brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,6,v,grad,a,x,b,tol,xmin) 3002 3003 end subroutine linmin
ABINIT/m_prcref [ Modules ]
NAME
m_prcref
FUNCTION
Routines to precondition residual potential (or density) and forces.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, MT, PMA) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_prcref 23 24 use defs_basis 25 use defs_wvltypes 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_xcdata 30 use m_frskerker1 31 use m_frskerker2 32 use mod_prc_memory 33 use m_dtset 34 35 use defs_datatypes, only : pseudopotential_type 36 use defs_abitypes, only : MPI_type 37 use m_time, only : timab 38 use m_numeric_tools, only : dotproduct 39 use m_geometry, only : xcart2xred, metric 40 use m_cgtools, only : dotprod_vn, mean_fftr 41 use m_mpinfo, only : ptabs_fourdp, destroy_mpi_enreg, initmpi_seq 42 use m_pawtab, only : pawtab_type 43 use m_pawrhoij, only : pawrhoij_type 44 use m_fftcore, only : kgindex 45 use m_fft, only : zerosym, indirect_parallel_fourier, fourdp 46 use m_kg, only : getph 47 use m_spacepar, only : hartre, laplacian 48 use m_distribfft, only : init_distribfft_seq 49 use m_forces, only : fresid 50 use m_atm2fft, only : atm2fft 51 use m_rhotoxc, only : rhotoxc 52 use m_mklocl, only : mklocl 53 use m_mkcore, only : mkcore 54 55 implicit none 56 57 private
ABINIT/moddiel [ 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.
mpi_enreg=information about MPI parallelization
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
SOURCE
1330 subroutine moddiel(cplex,dielar,mpi_enreg,nfft,ngfft,nspden,optreal,optres,qphon,rprimd,vresid,vrespc) 1331 1332 !Arguments------------------------------- 1333 !scalars 1334 integer,intent(in) :: cplex,nfft,nspden,optreal,optres 1335 type(MPI_type),intent(in) :: mpi_enreg 1336 !arrays 1337 integer,intent(in) :: ngfft(18) 1338 real(dp),intent(in) :: dielar(7),qphon(3),rprimd(3,3) 1339 real(dp),intent(in) :: vresid(cplex*nfft,nspden) 1340 real(dp),intent(out) :: vrespc(cplex*nfft,nspden) 1341 1342 !Local variables------------------------------- 1343 !scalars 1344 integer,parameter :: im=2,re=1 1345 integer :: i1,i2,i23,i3,ifft,ig,ii,ii1,ing,ispden,me_fft,mg,n1,n2,n3,nproc_fft 1346 integer :: nspden_eff,qeq0 1347 logical :: magn_precon 1348 real(dp) :: dielng,diemac,diemac_inv,diemix,diemixmag,diemix_eff,factor,gqg2p3,gqgm12,gqgm13 1349 real(dp) :: gqgm23,gs,gs2,gs3,l2g2,length2,ucvol 1350 character(len=500) :: message 1351 !arrays 1352 integer :: id(3) 1353 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1354 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1355 real(dp) :: gmet(3,3),gprimd(3,3),potg0(4),rmet(3,3) 1356 real(dp),allocatable :: gq(:,:),work1(:,:),work2(:) 1357 1358 ! ************************************************************************* 1359 1360 !Check that cplex has an allowed value 1361 if(cplex/=1 .and. cplex/=2)then 1362 write(message,'(a,i0,a,a)')& 1363 & ' From the calling routine, cplex=',cplex,ch10,& 1364 & ' but the only value allowed are 1 and 2.' 1365 ABI_BUG(message) 1366 end if 1367 1368 if(cplex==1.and.optreal==2)then 1369 ABI_BUG('When optreal=2, cplex must be 2.') 1370 end if 1371 1372 !This is to allow q=0 1373 qeq0=0 1374 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1 1375 1376 !If cplex=1 then qphon should be 0 0 0 1377 if (cplex==1.and. qeq0/=1) then 1378 write(message,'(a,3e12.4,a)' )' cplex=1 but qphon=',qphon,' qphon should be 0 0 0.' 1379 ABI_BUG(message) 1380 end if 1381 1382 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1383 me_fft=ngfft(11) 1384 nproc_fft=ngfft(10) 1385 1386 !Get the distrib associated with this fft_grid 1387 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1388 1389 !Compute different geometric tensor, as well as ucvol, from rprimd 1390 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1391 1392 dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ; diemixmag=dielar(7) 1393 1394 magn_precon=(diemixmag>=zero) ! Set to true if magnetization has to be preconditionned 1395 diemixmag=abs(diemixmag) 1396 1397 !write(std_out,*)' moddiel : diemac, diemix, diemixmag =',diemac,diemix,diemixmag 1398 1399 if(abs(diemac-1.0_dp)<1.0d-6)then 1400 1401 ! Here, simple mixing is required, through macroscopic 1402 ! dielectric constant set to 1.0_dp . 1403 vrespc(:,1)=diemix*vresid(:,1) 1404 if (nspden/=1) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden) 1405 else 1406 1407 ! Magnetization is not preconditionned 1408 if (optres==1.and.nspden>1.and.(.not.magn_precon)) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden) 1409 1410 ! Here, model dielectric function (G-diagonal operator) 1411 1412 length2=(two_pi*dielng)**2 1413 diemac_inv=1.0_dp/diemac 1414 ABI_MALLOC(work1,(2,nfft)) 1415 if (optreal==1) then 1416 ABI_MALLOC(work2,(cplex*nfft)) 1417 end if 1418 1419 ! In order to speed the routine, precompute the components of g 1420 mg=maxval(ngfft) 1421 ABI_MALLOC(gq,(3,mg)) 1422 do ii=1,3 1423 id(ii)=ngfft(ii)/2+2 1424 do ing=1,ngfft(ii) 1425 ig=ing-(ing/id(ii))*ngfft(ii)-1 1426 gq(ii,ing)=ig+qphon(ii) 1427 end do 1428 end do 1429 1430 ! Do-loop on spins 1431 ! Note XG 010922 : I doubt the preconditioner is OK for the magnetization 1432 nspden_eff=nspden;if (optres==1.and.(.not.magn_precon)) nspden_eff=1 1433 do ispden=1,nspden_eff 1434 1435 diemix_eff=diemix;if (ispden>1) diemix_eff=diemixmag 1436 1437 ! Do fft from real space (work2) to G space (work1) 1438 if (optreal==1) then 1439 work2(:)=vresid(:,ispden) 1440 call fourdp(cplex,work1,work2,-1,mpi_enreg,nfft,1,ngfft,0) 1441 else 1442 ! work1(:,:)=reshape(vresid(:,ispden),(/2,nfft/)) 1443 ! Reshape function does not work with big arrays for some compilers 1444 do ifft=1,nfft 1445 work1(1,ifft)=vresid(2*ifft-1,ispden) 1446 work1(2,ifft)=vresid(2*ifft ,ispden) 1447 end do 1448 end if 1449 1450 ! Store G=0 value 1451 potg0(ispden)=work1(re,1) 1452 1453 ! Triple loop, for the three dimensions 1454 do i3=1,n3 1455 ! Precompute some products that do not depend on i2 and i1 1456 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 1457 gqgm23=gq(3,i3)*gmet(2,3)*2 1458 gqgm13=gq(3,i3)*gmet(1,3)*2 1459 do i2=1,n2 1460 if (fftn2_distrib(i2)==me_fft) then 1461 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 1462 gqgm12=gq(2,i2)*gmet(1,2)*2 1463 gqg2p3=gqgm13+gqgm12 1464 i23=n1*( ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1)) 1465 1466 ! Do the test that eliminates the Gamma point outside 1467 ! of the inner loop 1468 ii1=1 1469 if(i2 == 1 .and. i3 == 1 .and. qeq0==1)then 1470 ! if(i23==0 .and. qeq0==1)then: this changes with the number of fft procs... 1471 ! and seems to be wrong.Pls check 1472 ii1=2 1473 end if 1474 1475 ! Here, unlike in hartre.f, the G=0 term is not eliminated, albeit 1476 ! not changed. 1477 do i1=ii1,n1 1478 1479 ! One obtains the square of the norm of q+G (defined by i1,i2,i3) 1480 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 1481 ifft=i1+i23 1482 1483 l2g2=length2*gs 1484 ! The model dielectric function is now computed 1485 factor = (l2g2+diemac_inv)/(l2g2+1.0_dp) * diemix_eff 1486 work1(re,ifft)=work1(re,ifft)*factor 1487 work1(im,ifft)=work1(im,ifft)*factor 1488 1489 end do 1490 end if 1491 end do 1492 end do 1493 1494 ! Might get rid of the G=0 term 1495 ! if(qeq0==1)then 1496 ! work1(re,1)=0.0_dp 1497 ! work1(im,1)=0.0_dp 1498 ! end if 1499 1500 ! Fourier transform 1501 if (optreal==1) then 1502 call fourdp(cplex,work1,work2,1,mpi_enreg,nfft,1,ngfft,0) 1503 vrespc(:,ispden)=work2(:) 1504 else 1505 ! vrespc(:,ispden)=reshape(work1(:,:),(/nfft*2/)) 1506 ! Reshape function does not work with big arrays for some compilers 1507 do ifft=1,nfft 1508 vrespc(2*ifft-1,ispden)=work1(1,ifft) 1509 vrespc(2*ifft ,ispden)=work1(2,ifft) 1510 end do 1511 end if 1512 1513 ! End of loop on spin polarizations 1514 end do 1515 1516 ABI_FREE(gq) 1517 ABI_FREE(work1) 1518 if (optreal==1) then 1519 ABI_FREE(work2) 1520 end if 1521 1522 ! End condition diemac/=1.0 1523 end if 1524 1525 end subroutine moddiel
ABINIT/prcref [ Functions ]
NAME
prcref
FUNCTION
Compute preconditioned residual potential (or density) and forces. iprcel, densfor_pred and iprcfc govern the choice of the preconditioner. Three tasks are done: 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 | intxc=control xc quadrature | 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. mpi_enreg=information about MPI parallelization 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
SOURCE
179 subroutine prcref(atindx,dielar,dielinv,& 180 & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,& 181 & istep,kg_diel,kxc,& 182 & mgfft,moved_atm_inside,mpi_enreg,my_natom,& 183 & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 184 & optreal,optres,pawrhoij,pawtab,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,& 185 & susmat,vhartr,vpsp,vresid,vrespc,vxc,wvl,wvl_den,xred) 186 187 !Arguments------------------------------- 188 !scalars 189 integer,intent(in) :: dielstrt,istep,my_natom,mgfft,moved_atm_inside,n1xccc 190 integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat,optreal,optres 191 real(dp),intent(in) :: etotal,gsqcut 192 type(MPI_type),intent(in) :: mpi_enreg 193 type(dataset_type),intent(in) :: dtset 194 type(pseudopotential_type),intent(in) :: psps 195 type(wvl_internal_type), intent(in) :: wvl 196 type(wvl_denspot_type), intent(inout) :: wvl_den 197 !arrays 198 integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft)) 199 integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18) 200 real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),rhog(2,nfft) 201 real(dp),intent(in) :: rhor(nfft,dtset%nspden) 202 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 203 real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden) 204 real(dp),intent(in) :: vxc(nfft,dtset%nspden) 205 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 206 real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc) 207 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft) 208 real(dp),intent(inout) :: xred(3,dtset%natom) 209 real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix),rprimd(3,3) 210 real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden) 211 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 212 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 213 214 !Local variables------------------------------- 215 !scalars 216 integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1 217 integer :: ipw2,iq,iq0,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm 218 integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method 219 real(dp) :: ai,ar,diemix,diemixmag,eei,enxc 220 real(dp) :: mixfac 221 real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg 222 logical :: computediel,non_magnetic_xc 223 character(len=500) :: message 224 type(xcdata_type) :: xcdata 225 !arrays 226 integer :: qprtrb(3) 227 integer,allocatable :: indpw_prc(:) 228 real(dp) :: dummy6(6),dummy7(6),gprimd(3,3),qphon(3),rmet(3,3),strsxc(6) 229 real(dp) :: vmean(dtset%nspden),vprtrb(2) 230 real(dp),allocatable :: dummy(:),dummy1(:),dummy2(:),dummy3(:),dummy4(:),dummy5(:),dummy8(:),dummy9(:) 231 real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:) 232 real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:) 233 real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:),rhor_new(:,:) 234 real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:) 235 real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:) 236 real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:) 237 logical,allocatable :: mask(:) 238 239 ! ************************************************************************* 240 241 !Compute different geometric tensor, as well as ucvol, from rprimd 242 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 243 244 !1) Eventually take care of the forces 245 246 if(moved_atm_inside==1)then 247 ABI_MALLOC(fcart_pc,(3,dtset%natom)) 248 249 if(dtset%iprcfc==0)then 250 fcart_pc(:,:)=fcart(:,:) 251 else 252 fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:) 253 end if 254 255 ! Compute preconditioned delta xred from preconditioned fcart and rprimd 256 call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc) 257 258 ABI_FREE(fcart_pc) 259 end if 260 261 !####################################################################### 262 263 !2) Take care of the potential residual 264 265 !Compute the residuals corresponding to the solution 266 !of an approximate realspace dielectric function according 267 !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]] 268 if(dtset%iprcel>=71.and.dtset%iprcel<=79) then 269 if (nfft==nfftprc) then 270 if (dtset%iprcel<=78) then 271 call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, & 272 & gprimd,vresid,vrespc,rhor(:,1)) 273 else 274 call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, & 275 & vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol) 276 end if 277 else 278 ! If preconditionning has to be done on a coarse grid, 279 ! has to transfer several arrays 280 ABI_MALLOC(work1,(nfftprc,dtset%nspden)) 281 ABI_MALLOC(work3,(nfftprc,dtset%nspden)) 282 ABI_MALLOC(work,(2*nfftprc)) 283 do ispden=1,dtset%nspden 284 work(:)=vresid(:,ispden) 285 call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,1,ngfftprc,0) 286 end do 287 ABI_FREE(work) 288 if (dtset%iprcel<=78) then 289 ABI_MALLOC(rhog_wk,(2,nfftprc)) 290 rhog_wk(:,:)=zero 291 if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then 292 nfftot=PRODUCT(ngfft(1:3)) 293 call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,& 294 & ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot) 295 else 296 do ii=1,nfft 297 if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii) 298 end do 299 end if 300 call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),& 301 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 302 ABI_MALLOC(work,(nfftprc)) 303 call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,1,ngfftprc,0) 304 call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, & 305 & gprimd,work1,work3,work) 306 ABI_FREE(work) 307 else 308 call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, & 309 & work1,work3,dtset%natom,xred,mpi_enreg,ucvol) 310 end if 311 do ispden=1,dtset%nspden 312 call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,1,ngfftprc,0) 313 end do 314 ABI_FREE(work1) 315 ABI_FREE(work3) 316 end if 317 318 else 319 320 if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then 321 cplex=optreal 322 qphon(:)=zero 323 ! Simple scalar multiplication, or model dielectric function 324 call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,qphon,rprimd,vresid,vrespc) 325 326 ! Use the inverse dielectric matrix in a small G sphere 327 else if( (istep>=dielstrt .and. dtset%iprcel>=21) .or. modulo(dtset%iprcel,100)>=41 )then 328 329 ! With dielop=1, the matrices will be computed when istep=dielstrt 330 ! With dielop=2, the matrices will be computed when istep=dielstrt and 1 331 dielop=1 332 if(modulo(dtset%iprcel,100)>=41)dielop=2 333 computediel = dtset%testsusmat(dielop, dielstrt, istep) !test if the matrix is to be computed 334 if(computediel) then 335 ! Compute the inverse dielectric matrix from the susceptibility matrix 336 ! There are two routines for the RPA matrix, while for the electronic 337 ! dielectric matrix, only dieltcel will do the work 338 if(modulo(dtset%iprcel,100)<=49)then 339 call dielmt(dielinv,gmet,kg_diel,& 340 & npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat) 341 else 342 option=1 343 if(modulo(dtset%iprcel,100)>=61)option=2 344 call dieltcel(dielinv,gmet,kg_diel,kxc,& 345 & nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%prtvol,susmat) 346 end if 347 end if 348 349 ABI_MALLOC(work1,(2,nfftprc)) 350 ABI_MALLOC(work2,(optreal*nfftprc)) 351 352 ! Presently, one uses the inverse of the RPA dielectric matrix, 353 ! for which spin must be averaged. 354 355 ! Do fft from real space (work2) to G space (work1) 356 if (optreal==1) then 357 work2(:)=vresid(:,1) 358 ! Must average over spins in the case of a potential residual 359 if(dtset%nspden/=1.and.optres==0)work2(:)=(work2(:)+vresid(:,2))*half 360 call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,1,ngfftprc,0) 361 else 362 work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/)) 363 if(dtset%nspden/=1.and.optres==0)work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half 364 end if 365 366 ! Multiply by restricted inverse of dielectric matrix. 367 ! Must first copy relevant elements of work1 to a npwdiel-dimensioned array, 368 ! then zero work1, operate with the dielinv matrix, and store in work1. 369 370 ABI_MALLOC(vres_diel,(2,npwdiel)) 371 ABI_MALLOC(indpw_prc,(npwdiel)) 372 ABI_MALLOC(mask,(npwdiel)) 373 mask(:)=.true. 374 call kgindex(indpw_prc,kg_diel,mask,mpi_enreg,ngfftprc,npwdiel) 375 376 do ipw1=1,npwdiel 377 if(mask(ipw1)) then 378 vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1)) 379 vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1)) 380 end if 381 end do 382 work1(:,:)=zero 383 do ipw1=1,npwdiel 384 ar=zero ; ai=zero 385 386 ! Use inverse of dielectric matrix (potential mixing) 387 if (optres==0) then 388 do ipw2=1,npwdiel 389 if(mask(ipw2))then 390 ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 391 & -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 392 ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 393 & +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 394 end if 395 end do 396 else 397 ! Use symetric of inverse of dielectric matrix (density mixing) 398 do ipw2=1,npwdiel 399 if(mask(ipw2))then 400 ar=ar+dielinv(1,ipw2,1,ipw1,1)*vres_diel(1,ipw2) & 401 & +dielinv(2,ipw2,1,ipw1,1)*vres_diel(2,ipw2) 402 ai=ai-dielinv(2,ipw2,1,ipw1,1)*vres_diel(1,ipw2) & 403 & +dielinv(1,ipw2,1,ipw1,1)*vres_diel(2,ipw2) 404 end if 405 end do 406 end if 407 ! Must be careful not to count the diagonal 1 twice : it is added later, 408 ! so must be subtracted now. 409 call xmpi_sum(ar,mpi_enreg%comm_fft,ier) 410 call xmpi_sum(ai,mpi_enreg%comm_fft,ier) 411 if(mask(ipw1)) then 412 work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1) 413 work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1) 414 end if !mask(ipw1) 415 end do ! ipw1 416 ABI_FREE(vres_diel) 417 ABI_FREE(indpw_prc) 418 ABI_FREE(mask) 419 ! Fourier transform 420 if (optreal==1) then 421 call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,1,ngfftprc,0) 422 else 423 work2(:)=reshape(work1(:,:),(/nfftprc*2/)) 424 end if 425 426 ! Add to get the preconditioned vresid, must be careful about spins. 427 if(dtset%iprcel>=30)then 428 diemix=dielar(4);diemixmag=abs(dielar(7)) 429 vrespc(:,1)=diemix*(vresid(:,1)+work2(:)) 430 if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:)) 431 if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=diemixmag*vresid(:,3:4) 432 if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=diemixmag*vresid(:,2:dtset%nspden) 433 else 434 vrespc(:,1)=vresid(:,1)+work2(:) 435 if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=vresid(:,2)+work2(:) 436 if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=vresid(:,3:4) 437 if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=vresid(:,2:dtset%nspden) 438 end if 439 440 ABI_FREE(work1) 441 ABI_FREE(work2) 442 443 ! Other choice ? 444 else 445 write(message, '(a,i3,a,a,a,a)' )& 446 & 'From the calling routine, iprcel=',dtset%iprcel,ch10,& 447 & 'The only allowed values are 0 or larger than 20.',ch10,& 448 & 'Action: correct your input file.' 449 ABI_ERROR(message) 450 end if 451 end if 452 !####################################################################### 453 454 !3) PAW only : precondition the rhoij quantities (augmentation 455 !occupancies) residuals. Use a simple preconditionning 456 !with the same mixing factor as the model dielectric function. 457 458 if (psps%usepaw==1.and.my_natom>0) then 459 if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then 460 mixfac=one;mixfacmag=one 461 else 462 mixfac=dielar(4);mixfacmag=abs(dielar(7)) 463 end if 464 if (pawrhoij(1)%cplex_rhoij==1) then 465 index=0 466 do iatom=1,my_natom 467 do iq=1,pawrhoij(iatom)%qphase 468 iq0=merge(0,pawrhoij(iatom)%lmn2_size,iq==1) 469 do ispden=1,pawrhoij(iatom)%nspden 470 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 471 do kmix=1,pawrhoij(iatom)%lmnmix_sz 472 index=index+1;klmn=iq0+pawrhoij(iatom)%kpawmix(kmix) 473 rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 474 end do 475 end do 476 end do 477 end do 478 else 479 index=-1 480 do iatom=1,my_natom 481 do iq=1,pawrhoij(iatom)%qphase 482 iq0=merge(0,2*pawrhoij(iatom)%lmn2_size,iq==1) 483 do ispden=1,pawrhoij(iatom)%nspden 484 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 485 do kmix=1,pawrhoij(iatom)%lmnmix_sz 486 index=index+2;klmn=iq0+2*pawrhoij(iatom)%kpawmix(kmix)-1 487 rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 488 end do 489 end do 490 end do 491 end do 492 end if 493 end if 494 495 !####################################################################### 496 497 !4) Take care of the change of atomic positions 498 !Note : this part is very demanding on memory... 499 !however, since this algorithm is still in development, 500 !it was NOT included in the estimation provided by memory.f 501 if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then 502 503 ! Not yet compatible with resid given in reciprocal space 504 if (optreal/=1) then 505 write(message, '(5a)' )& 506 & 'From the calling routine, densfor_pred=3',ch10,& 507 & 'You cannot use residuals in reciprocal space.',ch10,& 508 & 'Action: correct your input file.' 509 ABI_ERROR(message) 510 end if 511 ! Not compatible with non-collinear magnetism 512 if(dtset%nspden==4)then 513 ABI_ERROR('densfor_pred=3 does not work for nspden=4 !') 514 end if 515 516 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 517 nfftot=PRODUCT(ngfft(1:3)) 518 519 if (optres==0) then ! Array vresid contains a potential residual 520 ! ----------------------------------------------------------------- 521 522 ! First subtract the current local, hartree and exchange correlation potentials 523 do ispden=1,min(dtset%nspden,2) 524 vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden) 525 end do 526 if (dtset%nspden==4) then 527 do ispden=3,4 528 vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden) 529 end do 530 end if 531 532 ! Compute the modified density, in rhor_wk 533 option=2 534 ABI_MALLOC(gresid,(3,dtset%natom)) 535 ABI_MALLOC(grxc,(3,dtset%natom)) 536 ABI_MALLOC(rhor_wk,(nfft,dtset%nspden)) 537 ABI_MALLOC(rhor_wk0,(nfft,dtset%nspden)) 538 ABI_MALLOC(xred_wk,(3,dtset%natom)) 539 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 540 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 541 & ntypat,option,pawtab,rhor,rprimd,& 542 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 543 544 ! Compute up+down rhog_wk(G) by fft 545 ABI_MALLOC(work,(nfft)) 546 ABI_MALLOC(rhog_wk,(2,nfft)) 547 work(:)=rhor_wk(:,1) 548 call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,1,ngfft,0) 549 ABI_FREE(work) 550 551 ! Compute structure factor phases for new atomic pos: 552 call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk) 553 554 ! Compute local ionic pseudopotential vpsp: 555 ! and core electron density xccc3d, if needed. 556 n3xccc=0;if (n1xccc/=0) n3xccc=nfft 557 ABI_MALLOC(xccc3d,(n3xccc)) 558 ABI_MALLOC(vpsp_wk,(nfft)) 559 vprtrb(1:2)=zero 560 561 ! Determine by which method the local ionic potential and/or 562 ! the pseudo core charge density contributions have to be computed 563 ! Local ionic potential: 564 ! Method 1: PAW 565 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 566 vloc_method=1;if (psps%usepaw==0) vloc_method=2 567 if (dtset%icoulomb>0) vloc_method=2 568 if (psps%usewvl==1) vloc_method=2 569 ! Pseudo core charge density: 570 ! Method 1: PAW, nc_xccc_gspace 571 ! Method 2: Norm-conserving PP, wavelets 572 coredens_method=1;if (psps%usepaw==0) coredens_method=2 573 if (psps%nc_xccc_gspace==1) coredens_method=1 574 if (psps%nc_xccc_gspace==0) coredens_method=2 575 if (psps%usewvl==1) coredens_method=2 576 577 ! Local ionic potential and/or pseudo core charge by method 1 578 if (vloc_method==1.or.coredens_method==1) then 579 optv=0;if (vloc_method==1) optv=1 580 optn=0;if (coredens_method==1) optn=n3xccc/nfft 581 optatm=1;optdyfr=0;optgr=0;optstr=0;optn2=1;opteltfr=0 582 ! Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused... 583 call atm2fft(atindx,xccc3d,vpsp,dummy,dummy2,dummy9,dummy1,gmet,gprimd,dummy3,dummy4,gsqcut,& 584 & mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,& 585 & optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dtset%rcut,dummy5,rprimd,dummy6,dummy7,& 586 & ucvol,psps%usepaw,dummy8,dummy8,dummy8,vprtrb,psps%vlspl,& 587 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 588 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 589 end if 590 591 ! Local ionic potential by method 2 592 if (vloc_method==2) then 593 option=1 594 ABI_MALLOC(dyfrlo_indx,(3,3,dtset%natom)) 595 ABI_MALLOC(grtn_indx,(3,dtset%natom)) 596 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,& 597 & mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,& 598 & ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,& 599 & ucvol,vprtrb,vpsp_wk,wvl,wvl_den,xred) 600 ABI_FREE(dyfrlo_indx) 601 ABI_FREE(grtn_indx) 602 end if 603 604 ! Pseudo core electron density by method 2 605 if (coredens_method==2.and.n1xccc/=0) then 606 option=1 607 ABI_MALLOC(dyfrx2,(3,3,dtset%natom)) 608 ABI_MALLOC(grxc_indx,(3,dtset%natom)) 609 call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 610 & n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,& 611 & psps%xccc1d,xccc3d,xred_wk) 612 ABI_FREE(dyfrx2) 613 ABI_FREE(grxc_indx) 614 end if 615 616 ! Compute Hartree+xc potentials 617 ABI_MALLOC(vxc_wk,(nfft,dtset%nspden)) 618 ABI_MALLOC(vhartr_wk,(nfft)) 619 option=1 620 621 call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfft,ngfft,& 622 &dtset%nkpt,dtset%rcut,rhog_wk,rprimd,dtset%vcutgeo,vhartr_wk) 623 624 ! Prepare the call to rhotoxc 625 call xcdata_init(xcdata,dtset=dtset) 626 nk3xc=1 ; non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 627 call rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft,& 628 & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor_wk,rprimd,strsxc,1,& 629 & vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk) 630 ABI_FREE(xccc3d) 631 632 ! Sum all contributions 633 do ispden=1,min(dtset%nspden,2) 634 do ifft=1,nfft 635 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden) 636 end do 637 end do 638 if (dtset%nspden==4) then 639 do ispden=3,4 640 do ifft=1,nfft 641 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden) 642 end do 643 end do 644 end if 645 call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft) 646 if(dtset%nspden==2) then 647 vmean(1)=half*(vmean(1)+vmean(2)) 648 vmean(2)=vmean(1) 649 end if 650 do ispden=1,dtset%nspden 651 vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden) 652 end do 653 ABI_FREE(gresid) 654 ABI_FREE(grxc) 655 ABI_FREE(rhog_wk) 656 ABI_FREE(rhor_wk) 657 ABI_FREE(rhor_wk0) 658 ABI_FREE(xred_wk) 659 ABI_FREE(vhartr_wk) 660 ABI_FREE(vpsp_wk) 661 ABI_FREE(vxc_wk) 662 663 else ! Array vresid contains a density residual 664 ! ----------------------------------------------------------------- 665 666 ! Only have to compute the modified preconditionned density residual 667 option=2 668 ABI_MALLOC(gresid,(3,dtset%natom)) 669 ABI_MALLOC(grxc,(3,dtset%natom)) 670 ABI_MALLOC(rhor_new,(nfft,dtset%nspden)) 671 ABI_MALLOC(rhor_wk,(nfft,dtset%nspden)) 672 ABI_MALLOC(rhor_wk0,(nfft,dtset%nspden)) 673 ABI_MALLOC(xred_wk,(3,dtset%natom)) 674 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 675 rhor_new(:,1)=rhor(:,1)+vrespc(:,1) 676 if (dtset%nspden==2) then 677 rhor_new(:,1)=rhor_new(:,1)+vrespc(:,2) 678 rhor_new(:,2)=rhor(:,2)+vrespc(:,1) 679 end if 680 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 681 & ntypat,option,pawtab,rhor,rprimd,& 682 & ucvol,rhor_wk0,xred_wk,xred,psps%znuclpsp) 683 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 684 & ntypat,option,pawtab,rhor_new,rprimd,& 685 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 686 vrespc(:,1)=rhor_wk(:,dtset%nspden)-rhor_wk0(:,dtset%nspden) 687 if (dtset%nspden==2) vrespc(:,2)=rhor_wk(:,1)-rhor_wk0(:,1)-vrespc(:,1) 688 ABI_FREE(gresid) 689 ABI_FREE(grxc) 690 ABI_FREE(rhor_new) 691 ABI_FREE(rhor_wk) 692 ABI_FREE(rhor_wk0) 693 ABI_FREE(xred_wk) 694 end if 695 696 end if 697 698 end subroutine prcref
ABINIT/prcref_PMA [ Functions ]
NAME
prcref_PMA
FUNCTION
Compute preconditioned residual potential (or density) and forces. iprcel, densfor_pred and iprcfc govern the choice of the preconditioner. Three tasks are done: 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 | intxc=control xc quadrature | 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. mpi_enreg=information about MPI parallelization 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
SOURCE
813 subroutine prcref_PMA(atindx,dielar,dielinv,& 814 & dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,& 815 & istep,kg_diel,kxc,& 816 & mgfft,moved_atm_inside,mpi_enreg,my_natom,& 817 & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 818 & optreal,optres,pawrhoij,ph1d,psps,rhog, rhoijrespc,rhor,rprimd,& 819 & susmat,vhartr,vpsp,vresid,vrespc,vxc,xred,& 820 & etotal,pawtab,wvl) 821 822 !Arguments------------------------------- 823 !variables used for tfvw 824 !scalars 825 integer,intent(in) :: dielstrt,istep,mgfft,moved_atm_inside,my_natom,n1xccc 826 integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat 827 integer,intent(in) :: optreal,optres 828 real(dp),intent(in) :: etotal,gsqcut 829 type(MPI_type),intent(in) :: mpi_enreg 830 type(dataset_type),intent(in) :: dtset 831 type(pseudopotential_type),intent(in) :: psps 832 type(wvl_data), intent(inout) :: wvl 833 !arrays 834 integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft)) 835 integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18) 836 real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom) 837 real(dp),intent(in) :: rhog(2,nfft) 838 real(dp),intent(in) :: rhor(nfft,dtset%nspden) 839 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 840 real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden) 841 real(dp),intent(in) :: vxc(nfft,dtset%nspden) 842 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 843 real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc) 844 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft) 845 real(dp),intent(inout) :: xred(3,dtset%natom),rprimd(3,3) 846 real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix) 847 real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden) 848 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 849 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 850 851 !Local variables------------------------------- 852 !scalars 853 integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1 854 integer :: ipw2,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm 855 integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method 856 real(dp) :: ai,ar,diemix,diemixmag,eei,enxc 857 real(dp) :: mixfac 858 real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg 859 logical :: computediel,non_magnetic_xc 860 character(len=500) :: message 861 type(xcdata_type) :: xcdata 862 !arrays 863 integer :: qprtrb(3) 864 integer,allocatable :: indpw_prc(:) 865 real(dp) :: dummy6(6),gprimd(3,3),qphon(3),rmet(3,3),strsxc(6) 866 real(dp) :: vmean(dtset%nspden),vprtrb(2) 867 real(dp),allocatable :: dummy_in(:) 868 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0) 869 real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:) 870 real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:) 871 real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:) 872 real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:) 873 real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:) 874 real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:) 875 logical,allocatable :: mask(:) 876 877 ! ************************************************************************* 878 879 if(optres==1)then 880 ABI_ERROR('density mixing (optres=1) not admitted!') 881 end if 882 883 !Compute different geometric tensor, as well as ucvol, from rprimd 884 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 885 886 !1) Eventually take care of the forces 887 888 if(moved_atm_inside==1)then 889 ABI_MALLOC(fcart_pc,(3,dtset%natom)) 890 891 if(dtset%iprcfc==0)then 892 fcart_pc(:,:)=fcart(:,:) 893 else 894 fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:) 895 end if 896 897 ! Compute preconditioned delta xred from preconditioned fcart and rprimd 898 call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc) 899 900 ABI_FREE(fcart_pc) 901 end if 902 903 !####################################################################### 904 905 !2) Take care of the potential residual 906 907 !Compute the residuals corresponding to the solution 908 !of an approximate realspace dielectric function according 909 !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]] 910 if(dtset%iprcel>=71.and.dtset%iprcel<=79) then 911 if (nfft==nfftprc) then 912 if (dtset%iprcel<=78) then 913 call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, & 914 & gprimd,vresid,vrespc,rhor(:,1)) 915 else 916 call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, & 917 & vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol) 918 end if 919 else 920 ! If preconditionning has to be done on a coarse grid, 921 ! has to transfer several arrays 922 ABI_MALLOC(work1,(nfftprc,dtset%nspden)) 923 ABI_MALLOC(work3,(nfftprc,dtset%nspden)) 924 ABI_MALLOC(work,(2*nfftprc)) 925 do ispden=1,dtset%nspden 926 work(:)=vresid(:,ispden) 927 call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,1,ngfftprc,0) 928 end do 929 ABI_FREE(work) 930 if (dtset%iprcel<=78) then 931 ABI_MALLOC(rhog_wk,(2,nfftprc)) 932 rhog_wk(:,:)=zero 933 if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then 934 nfftot=PRODUCT(ngfft(1:3)) 935 call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,& 936 & ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot) 937 else 938 do ii=1,nfft 939 if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii) 940 end do 941 end if 942 call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),& 943 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 944 ABI_MALLOC(work,(nfftprc)) 945 call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,1,ngfftprc,0) 946 call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, & 947 & gprimd,work1,work3,work) 948 ABI_FREE(work) 949 else 950 call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, & 951 & work1,work3,dtset%natom,xred,mpi_enreg,ucvol) 952 end if 953 do ispden=1,dtset%nspden 954 call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,1,ngfftprc,0) 955 end do 956 ABI_FREE(work1) 957 ABI_FREE(work3) 958 end if 959 960 else 961 962 if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then 963 cplex=optreal 964 qphon(:)=zero 965 ! Simple scalar multiplication, or model dielectric function 966 call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,qphon,rprimd,vresid,vrespc) 967 968 ! Use the inverse dielectric matrix in a small G sphere 969 else if( (istep>=dielstrt .and. dtset%iprcel>=21) .or. modulo(dtset%iprcel,100)>=41 )then 970 971 ! Wnith dielop=1, the matrices will be computed when istep=dielstrt 972 ! With dielop=2, the matrices will be computed when istep=dielstrt and 1 973 dielop=1 974 if(modulo(dtset%iprcel,100)>=41)dielop=2 975 computediel = dtset%testsusmat(dielop, dielstrt, istep) !test if the matrix is to be computed 976 if(computediel) then 977 ! Compute the inverse dielectric matrix from the susceptibility matrix 978 ! There are two routines for the RPA matrix, while for the electronic 979 ! dielectric matrix, only dieltcel will do the work 980 if(modulo(dtset%iprcel,100)<=49)then 981 call dielmt(dielinv,gmet,kg_diel,& 982 & npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat) 983 else 984 option=1 985 if(modulo(dtset%iprcel,100)>=61)option=2 986 call dieltcel(dielinv,gmet,kg_diel,kxc,& 987 & nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%prtvol,susmat) 988 end if 989 end if 990 991 ABI_MALLOC(work1,(2,nfftprc)) 992 ABI_MALLOC(work2,(optreal*nfftprc)) 993 994 ! Presently, one uses the inverse of the RPA dielectric matrix, 995 ! for which spin must be averaged. 996 997 ! Do fft from real space (work2) to G space (work1) 998 if (optreal==1) then 999 work2(:)=vresid(:,1) 1000 ! Must average over spins if needed. 1001 if(dtset%nspden/=1)work2(:)=(work2(:)+vresid(:,2))*half 1002 call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,1,ngfftprc,0) 1003 else 1004 work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/)) 1005 if (dtset%nspden/=1) work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half 1006 end if 1007 1008 ! Multiply by restricted inverse of dielectric matrix. 1009 ! Must first copy relevant elements of work1 to a npwdiel-dimensioned array, 1010 ! then zero work1, operate with the dielinv matrix, and store in work1. 1011 1012 ABI_MALLOC(vres_diel,(2,npwdiel)) 1013 ABI_MALLOC(indpw_prc,(npwdiel)) 1014 ABI_MALLOC(mask,(npwdiel)) 1015 mask(:)=.true. 1016 call kgindex(indpw_prc,kg_diel,mask,mpi_enreg,ngfftprc,npwdiel) 1017 do ipw1=1,npwdiel 1018 if(mask(ipw1)) then 1019 vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1)) 1020 vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1)) 1021 end if 1022 end do 1023 1024 work1(:,:)=zero 1025 do ipw1=1,npwdiel 1026 ar=zero ; ai=zero 1027 1028 ! Use inverse of dielectric matrix (potential mixing) 1029 do ipw2=1,npwdiel 1030 if(mask(ipw2))then 1031 ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 1032 & -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 1033 ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 1034 & +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 1035 end if 1036 end do 1037 ! Must be careful not to count the diagonal 1 twice : it is added later, 1038 ! so must be subtracted now. 1039 call xmpi_sum(ar,mpi_enreg%comm_fft,ier) 1040 call xmpi_sum(ai,mpi_enreg%comm_fft,ier) 1041 if(mask(ipw1)) then 1042 work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1) 1043 work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1) 1044 end if !mask(ipw1) 1045 end do ! ipw1 1046 ABI_FREE(vres_diel) 1047 ABI_FREE(indpw_prc) 1048 ABI_FREE(mask) 1049 1050 ! Fourier transform 1051 if (optreal==1) then 1052 call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,1,ngfftprc,0) 1053 else 1054 work2(:)=reshape(work1(:,:),(/nfftprc*2/)) 1055 end if 1056 1057 ! Add to get the preconditioned vresid, must be careful about spins. 1058 if(dtset%iprcel>=30)then 1059 diemix=dielar(4);diemixmag=abs(dielar(7)) 1060 vrespc(:,1)=diemix*(vresid(:,1)+work2(:)) 1061 if(dtset%nspden/=1)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:)) 1062 if(dtset%nspden==4)vrespc(:,3:4)=diemixmag*vresid(:,3:4) 1063 else 1064 vrespc(:,1)=vresid(:,1)+work2(:) 1065 if(dtset%nspden/=1)vrespc(:,2)=vresid(:,2)+work2(:) 1066 if(dtset%nspden==4)vrespc(:,3:4)=vresid(:,3:4) 1067 end if 1068 1069 ABI_FREE(work1) 1070 ABI_FREE(work2) 1071 1072 ! Other choice ? 1073 else 1074 write(message, '(a,i0,a,a,a,a)' )& 1075 & 'From the calling routine, iprcel= ',dtset%iprcel,ch10,& 1076 & 'The only allowed values are 0 or larger than 20.',ch10,& 1077 & 'Action: correct your input file.' 1078 ABI_ERROR(message) 1079 end if 1080 end if 1081 !####################################################################### 1082 1083 !3) PAW only : precondition the rhoij quantities (augmentation 1084 !occupancies) residuals. Use a simple preconditionning 1085 !with the same mixing factor as the model dielectric function. 1086 1087 if (psps%usepaw==1.and.my_natom>0) then 1088 if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then 1089 mixfac=one;mixfacmag=one 1090 else 1091 mixfac=dielar(4);mixfacmag=abs(dielar(7)) 1092 end if 1093 if (pawrhoij(1)%cplex_rhoij==1) then 1094 index=0 1095 do iatom=1,my_natom 1096 do ispden=1,pawrhoij(iatom)%nspden 1097 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 1098 do kmix=1,pawrhoij(iatom)%lmnmix_sz 1099 index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix) 1100 rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 1101 end do 1102 end do 1103 end do 1104 else 1105 index=-1 1106 do iatom=1,my_natom 1107 do ispden=1,pawrhoij(iatom)%nspden 1108 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 1109 do kmix=1,pawrhoij(iatom)%lmnmix_sz 1110 index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 1111 rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 1112 end do 1113 end do 1114 end do 1115 end if 1116 end if 1117 !####################################################################### 1118 1119 !4) Take care of the change of atomic positions 1120 !Note : this part is very demanding on memory... 1121 !however, since this algorithm is still in development, 1122 !it was NOT included in the estimation provided by memory.f 1123 if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then 1124 1125 ! Not yet compatible with resid given in reciprocal space 1126 if (optreal/=1) then 1127 write(message, '(5a)' )& 1128 & 'From the calling routine, densfor_pred=3',ch10,& 1129 & 'You cannot use residuals in reciprocal space.',ch10,& 1130 & 'Action: correct your input file.' 1131 ABI_ERROR(message) 1132 end if 1133 1134 ! Not compatible with non-collinear magnetism 1135 if(dtset%nspden==4)then 1136 ABI_ERROR('densfor_pred=3 does not work for nspden=4!') 1137 end if 1138 1139 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1140 nfftot=PRODUCT(ngfft(1:3)) 1141 1142 ! First subtract the current local, hartree and exchange correlation potentials 1143 do ispden=1,min(dtset%nspden,2) 1144 vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden) 1145 end do 1146 if (dtset%nspden==4) then 1147 do ispden=3,4 1148 vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden) 1149 end do 1150 end if 1151 1152 ! Compute the modified density, in rhor_wk 1153 option=2 1154 ABI_MALLOC(gresid,(3,dtset%natom)) 1155 ABI_MALLOC(grxc,(3,dtset%natom)) 1156 ABI_MALLOC(rhor_wk,(nfft,dtset%nspden)) 1157 ABI_MALLOC(rhor_wk0,(nfft,dtset%nspden)) 1158 ABI_MALLOC(xred_wk,(3,dtset%natom)) 1159 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 1160 1161 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 1162 & ntypat,option,pawtab,rhor,rprimd,& 1163 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 1164 1165 ! Compute up+down rhog_wk(G) by fft 1166 ABI_MALLOC(work,(nfft)) 1167 ABI_MALLOC(rhog_wk,(2,nfft)) 1168 work(:)=rhor_wk(:,1) 1169 call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,1,ngfft,0) 1170 ABI_FREE(work) 1171 1172 ! Compute structure factor phases for new atomic pos: 1173 call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk) 1174 1175 ! Compute local ionic pseudopotential vpsp: 1176 ! and core electron density xccc3d, if needed. 1177 n3xccc=0;if (n1xccc/=0) n3xccc=nfft 1178 ABI_MALLOC(xccc3d,(n3xccc)) 1179 ABI_MALLOC(vpsp_wk,(nfft)) 1180 vprtrb(1:2)=zero 1181 1182 ! Determine by which method the local ionic potential and/or 1183 ! the pseudo core charge density have to be computed 1184 ! Local ionic potential: 1185 ! Method 1: PAW 1186 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 1187 vloc_method=1;if (psps%usepaw==0) vloc_method=2 1188 if (dtset%icoulomb>0) vloc_method=2 1189 if (psps%usewvl==1) vloc_method=2 1190 ! Pseudo core charge density: 1191 ! Method 1: PAW, nc_xccc_gspace 1192 ! Method 2: Norm-conserving PP, wavelets 1193 coredens_method=1;if (psps%usepaw==0) coredens_method=2 1194 if (psps%nc_xccc_gspace==1) coredens_method=1 1195 if (psps%nc_xccc_gspace==0) coredens_method=2 1196 if (psps%usewvl==1) coredens_method=2 1197 1198 ! Local ionic potential and/or pseudo core charge by method 1 1199 if (vloc_method==1.or.coredens_method==1) then 1200 optv=0;if (vloc_method==1) optv=1 1201 optn=0;if (coredens_method==1) optn=n3xccc/nfft 1202 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1 1203 ! Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused... 1204 call atm2fft(atindx,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,gmet,& 1205 & gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,& 1206 & nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,& 1207 & psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dtset%rcut,dummy_in,rprimd,dummy_out6,dummy_out7,ucvol,& 1208 & psps%usepaw,dummy_in,dummy_in,dummy_in,vprtrb,psps%vlspl,& 1209 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 1210 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 1211 end if 1212 1213 ! Local ionic potential by method 2 1214 if (vloc_method==2) then 1215 option=1 1216 ABI_MALLOC(dyfrlo_indx,(3,3,dtset%natom)) 1217 ABI_MALLOC(grtn_indx,(3,dtset%natom)) 1218 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,& 1219 & mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,& 1220 & ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,& 1221 & ucvol,vprtrb,vpsp_wk,wvl%descr,wvl%den,xred) 1222 ABI_FREE(dyfrlo_indx) 1223 ABI_FREE(grtn_indx) 1224 end if 1225 1226 ! Pseudo core electron density by method 2 1227 if (coredens_method==2.and.n1xccc/=0) then 1228 option=1 1229 ABI_MALLOC(dyfrx2,(3,3,dtset%natom)) 1230 ABI_MALLOC(grxc_indx,(3,dtset%natom)) 1231 call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 1232 & n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,& 1233 & psps%xccc1d,xccc3d,xred_wk) 1234 ABI_FREE(dyfrx2) 1235 ABI_FREE(grxc_indx) 1236 end if 1237 1238 ! Compute Hartree+xc potentials 1239 ABI_MALLOC(vxc_wk,(nfft,dtset%nspden)) 1240 ABI_MALLOC(vhartr_wk,(nfft)) 1241 option=1 1242 1243 call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfft,ngfft,& 1244 &dtset%nkpt,dtset%rcut,rhog_wk,rprimd,dtset%vcutgeo,vhartr_wk) 1245 1246 ! Prepare the call to rhotoxc 1247 call xcdata_init(xcdata,dtset=dtset) 1248 nk3xc=1 ; non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 1249 ABI_MALLOC(work,(0)) 1250 call rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft,& 1251 & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor_wk,rprimd,strsxc,1,& 1252 & vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk) 1253 ABI_FREE(work) 1254 ABI_FREE(xccc3d) 1255 1256 ! Sum all contributions 1257 do ispden=1,min(dtset%nspden,2) 1258 do ifft=1,nfft 1259 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden) 1260 end do 1261 end do 1262 if (dtset%nspden==4) then 1263 do ispden=3,4 1264 do ifft=1,nfft 1265 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden) 1266 end do 1267 end do 1268 end if 1269 call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft) 1270 if(dtset%nspden==2) then 1271 vmean(1)=half*(vmean(1)+vmean(2)) 1272 vmean(2)=vmean(1) 1273 end if 1274 do ispden=1,dtset%nspden 1275 vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden) 1276 end do 1277 ABI_FREE(gresid) 1278 ABI_FREE(grxc) 1279 ABI_FREE(rhog_wk) 1280 ABI_FREE(rhor_wk) 1281 ABI_FREE(rhor_wk0) 1282 ABI_FREE(xred_wk) 1283 ABI_FREE(vhartr_wk) 1284 ABI_FREE(vpsp_wk) 1285 ABI_FREE(vxc_wk) 1286 1287 end if 1288 1289 end subroutine prcref_PMA
ABINIT/prcrskerker1 [ 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...
SOURCE
2380 subroutine prcrskerker1(dtset,mpi_enreg,nfft,nspden,ngfft,dielar,etotal,gprimd,vresid,vrespc,base) 2381 2382 !Arguments ------------------------------------ 2383 !scalars 2384 integer,intent(in) :: nfft,nspden 2385 real(dp) :: etotal 2386 type(MPI_type),intent(in) :: mpi_enreg 2387 type(dataset_type),intent(in) :: dtset 2388 !arrays 2389 integer,intent(in) :: ngfft(18) 2390 real(dp),intent(in) :: base(nfft),dielar(7),gprimd(3,3) 2391 real(dp),intent(in) :: vresid(nfft,nspden) 2392 real(dp),intent(out) :: vrespc(nfft,nspden) 2393 2394 !Local variables------------------------------- 2395 !scalars 2396 integer :: ifft,ispden,n1,n2,n3 2397 real(dp) :: base_delta,base_max,base_min,dielng,diemac,diemix 2398 real(dp) :: diemixmag 2399 real(dp) :: rdummy1,rdummy2 2400 logical :: new_prc_func 2401 !arrays 2402 real(dp) :: deltaW(nfft,nspden) 2403 real(dp) :: g2cart(nfft) 2404 real(dp) :: mat(nfft,nspden) 2405 2406 ! ************************************************************************* 2407 2408 !DEBUG 2409 !write(std_out,*)' prckerker1 : enter ' 2410 !ENDDEBUG 2411 !if(cycle==0) then 2412 call prc_mem_init(nfft) 2413 2414 if(cycle==0) then 2415 new_prc_func=.TRUE. 2416 energy_min=etotal 2417 else if(etotal < energy_min) then 2418 new_prc_func=.TRUE. 2419 energy_min=etotal 2420 else 2421 new_prc_func=.FALSE. 2422 end if 2423 2424 2425 dielng=dielar(2) 2426 diemac=dielar(3) 2427 diemix=dielar(4) 2428 diemixmag=dielar(7) 2429 !****************************************************************** 2430 !compute the diemac(r) ** 2431 !****************************************************************** 2432 !this task will be devoted to a general function later 2433 n1=ngfft(1) 2434 n2=ngfft(2) 2435 n3=ngfft(3) 2436 !base_cp=base 2437 if(new_prc_func) then 2438 base_min=base(1) 2439 base_max=base(1) 2440 do ifft=1,nfft 2441 base_min = min(base_min,base(ifft)) 2442 base_max = max(base_max,base(ifft)) 2443 end do 2444 base_delta = base_max - base_min 2445 ! if(cycle.lt.2) then 2446 rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one) 2447 ! else 2448 ! rdiemac(:) = rdiemac(:)*0.5_dp+0.5_dp*(((base(:)-base_min) / (base_delta) ) *(diemac-one) + one) 2449 ! end if 2450 ! if(cycle==0) rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one) 2451 ! rdiemac(:) = exp(((base(:)-base_min) / (base_delta) *log(diemac))) 2452 end if 2453 cycle=cycle+1 2454 !if(cycle==5) cycle=0 2455 !end if 2456 !****************************************************************** 2457 !compute deltaW ** 2458 !****************************************************************** 2459 vrespc=vresid !starting point 2460 ! put the laplacian of the residuals into deltaW 2461 call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,rdfuncr=vrespc,laplacerdfuncr=deltaW,g2cart_out=g2cart) 2462 2463 !call laplacian(vrespc,buffer,ngfft,gprimd) ! put the laplacian of the residuals into deltaW 2464 !do ifft=1,nfft 2465 !if (buffer(ifft,1)/=deltaW(ifft,1)) then 2466 !stop 2467 !end if 2468 !end do 2469 deltaW(:,1)= diemix*(((one/rdiemac(:))*vresid(:,1))-(((dielng)**2)*deltaW(:,1))) 2470 if (nspden>1.and.(diemixmag>=zero)) then 2471 do ispden=2,nspden 2472 deltaW(:,ispden)= abs(diemixmag)*(((one/rdiemac(:))*vresid(:,ispden))-(((dielng)**2)*deltaW(:,ispden))) 2473 end do 2474 end if 2475 !call random_number(deltaW) 2476 !call random_number(vrespc) 2477 !****************************************************************** 2478 !Finding the preconditionned residuals which minimizes ** 2479 !half*(vrespc*(1-dielng2/4pi2 nabla2) vrespc) - vrespc * deltaW ** 2480 !*********************************************************************** 2481 vrespc(:,1)=diemix*vrespc(:,1) 2482 if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*vrespc(:,2:nspden) 2483 !buffer=vrespc 2484 2485 2486 !============================================================================== 2487 !============================================================================== 2488 !! Original loop 2489 !============================================================================== 2490 !============================================================================== 2491 2492 call frskerker1__init(dtset,mpi_enreg,nfft,ngfft,nspden,dielng,deltaW,gprimd,mat,g2cart) 2493 2494 !call cgpr(pf_rscgres,dpf_rscgres,newvres,real(1e-40,dp),700,vrespc,rdummy1,rdummy2) 2495 !rdummy1 = pf_rscgres(nfft,nspden,vrespc) 2496 call cgpr(nfft,nspden,frskerker1__pf,frskerker1__dpf,frskerker1__newvres,& 2497 & real(1e-10,dp),700,vrespc,rdummy1,rdummy2) 2498 call frskerker1__end() 2499 2500 !============================================================================== 2501 !============================================================================== 2502 !! Original loop end 2503 !============================================================================== 2504 !============================================================================== 2505 2506 2507 !cplex=1 2508 !qphon(:)=zero 2509 !call moddiel(cplex,dielar,nfft,ngfft,nspden,1,0,qphon,rprimd,vresid,buffer) 2510 !c1=0 2511 !do ifft=1,nfft,1 2512 !if((abs(buffer(ifft,1)-vrespc(ifft,1))/(abs(buffer(ifft,1)+vrespc(ifft,1))*half)) > 5e-3) then 2513 !c1=c1+1 2514 !end if 2515 !end do 2516 !call laplacian(vrespc,buffer,ngfft,gprimd) 2517 !buffer=vrespc(:,:)-buffer(:,:)*dielng**2 2518 !c2=0 2519 !do ifft=1,nfft,1 2520 !if((abs(buffer(ifft,1)-deltaW(ifft,1))/(abs(buffer(ifft,1)+deltaW(ifft,1))*half)) > 5e-3) then 2521 !c2=c2+1 2522 !end if 2523 !end do 2524 !!! !stop 2525 !call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,& 2526 !& g2cart_out=g2cart) 2527 2528 !vrespc=vresid 2529 !do ispden=1,nspden 2530 !call fourdp(1, gvrespc(:,:,ispden), vrespc(:,ispden),-1,mpi_enreg,nfft,ngfft,0) 2531 !end do 2532 !filtering 2533 !do ispden=1,nspden 2534 !do ifft=1,nfft 2535 !! gvrespc(:,ifft,ispden)=(one-exp(-g2cart(ifft)*15.0_dp))*gvrespc(:,ifft,ispden) 2536 !! gvrespc(:,ifft,ispden)=(exp(-g2cart(ifft)*10.0_dp))*gvrespc(:,ifft,ispden) 2537 !! gvrespc(:,ifft,ispden)=(one-one/(exp(-0.002_dp/g2cart(ifft)**2)+one))*gvrespc(:,ifft,ispden) 2538 !gvrespc(:,ifft,ispden)=(two-2_dp/(exp(-0.008_dp/(g2cart(ifft)+0.0012_dp))+one))*gvrespc(:,ifft,ispden) 2539 !gvrespc(:,ifft,ispden)=min(one,(sqrt(g2cart(ifft)/0.006_dp))**(one))*gvrespc(:,ifft,ispden) 2540 !end do 2541 !end do 2542 !change resulting potential to real space 2543 !do ispden=1,nspden 2544 !call fourdp(1,gvrespc(:,:,ispden),vrespc(:,ispden),1,mpi_enreg,nfft,ngfft,0) 2545 !end do 2546 !vrespc=vrespc*diemix 2547 !maxg2=g2cart(1) 2548 !ming2=g2cart(5) 2549 !do ifft=1,nfft 2550 !maxg2=max(g2cart(ifft),maxg2) 2551 !if(g2cart(ifft) .gt. zero) ming2=min(g2cart(ifft),ming2) 2552 !end do 2553 !stop 2554 2555 !DEBUG 2556 !write(std_out,*)' prckerker1 : exit ' 2557 !ENDDEBUG 2558 2559 end subroutine prcrskerker1
ABINIT/prcrskerker2 [ 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
SOURCE
2592 subroutine prcrskerker2(dtset,nfft,nspden,ngfft,dielar,gprimd,rprimd,vresid,vrespc,natom,xred,mpi_enreg,ucvol) 2593 2594 !Arguments ------------------------------------ 2595 !scalars 2596 integer,intent(in) :: natom,nfft,nspden 2597 real(dp),intent(in) :: ucvol 2598 type(MPI_type),intent(in) :: mpi_enreg 2599 type(dataset_type),intent(in) :: dtset 2600 !arrays 2601 integer,intent(in) :: ngfft(18) 2602 real(dp),intent(in) :: dielar(7),gprimd(3,3),rprimd(3,3),vresid(nfft,nspden) 2603 real(dp),intent(in) :: xred(3,natom) 2604 real(dp),intent(out) :: vrespc(nfft,nspden) 2605 2606 !Local variables------------------------------- 2607 !logical,save ::ok=.FALSE. 2608 !scalars 2609 integer :: cplex,i1,i2,i3,iatom,iatom27,ifft,ispden,n1,n2,n3,natom27,nfftotf 2610 integer :: option 2611 real(dp),save :: lastp1=one,lastp2=one 2612 real(dp) :: C1,C2,DE,core,dielng,diemac,diemix,diemixmag,doti,dr,l1,l2,l3,l4,r 2613 real(dp) :: rdummy1,rdummy2,rmin,xr,y,yr,zr 2614 !arrays 2615 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 2616 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 2617 real(dp) :: V1(nfft,nspden),V2(nfft,nspden),buffer(nfft,nspden) 2618 real(dp) :: deltaW(nfft,nspden) 2619 real(dp) :: mat(nfft,nspden) 2620 real(dp) :: rdielng(nfft),rdiemac(nfft),xcart(3,natom) 2621 real(dp) :: xcart27(3,natom*27) 2622 2623 ! ************************************************************************* 2624 2625 dielng=dielar(2) 2626 diemac=dielar(3) 2627 diemix=dielar(4) 2628 diemixmag=dielar(7) 2629 !****************************************************************** 2630 !compute the diemac(r) ** 2631 !****************************************************************** 2632 !this task will be devoted to a general function later 2633 n1=ngfft(1) 2634 n2=ngfft(2) 2635 n3=ngfft(3) 2636 nfftotf=n1*n2*n3 2637 !if(.not.ok) then 2638 xcart(1,:)=xred(1,:)*rprimd(1,1)+xred(2,:)*rprimd(1,2)+xred(3,:)*rprimd(1,3) 2639 xcart(2,:)=xred(1,:)*rprimd(2,1)+xred(2,:)*rprimd(2,2)+xred(3,:)*rprimd(2,3) 2640 xcart(3,:)=xred(1,:)*rprimd(3,1)+xred(2,:)*rprimd(3,2)+xred(3,:)*rprimd(3,3) 2641 2642 iatom27=0 2643 do i1=-1,1 2644 do i2=-1,1 2645 do i3=-1,1 2646 do iatom=1,natom 2647 iatom27=iatom27+1 2648 xcart27(:,iatom27)=xcart(:,iatom)+rprimd(:,1)*i1+rprimd(:,2)*i2+rprimd(:,3)*i3 2649 end do 2650 end do 2651 end do 2652 end do 2653 2654 !stop 2655 natom27=27*natom 2656 2657 l1=0.34580850339844665 2658 !l2=0.5123510203906797 !0.41242551019533985 2659 !l3=0.8001489796093203 !0.90007448980466009 2660 2661 l2=0.41242551019533985 2662 l3=0.90007448980466009 2663 l4=0.9666914966015534 2664 2665 2666 l1=0.31387233559896449 2667 l2=0.35828367346355994 2668 l3=0.9333829932031068 2669 l4=0.9777943310677023 2670 2671 l1=3.5 2672 l2=11.5 2673 l3=2.5 2674 l4=6.5 2675 !l1=30. !cellules pleines 2676 2677 rdielng=zero 2678 core=1. !(value of Er at the core of atoms) 2679 dr=2.65 ! radius of atoms=2.65165 2680 y=1. ! value of Er in the empty region 2681 2682 !Get the distrib associated with this fft_grid 2683 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2684 2685 do i3=1,n3 2686 ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft) 2687 do i2=1,n2 2688 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 2689 do i1=1,n1 2690 ifft=ifft+1 2691 ! !!!!!!!!!!!!!!!!!!!!!!!!! 2692 ! ! calculation of the simplest part void/metal 2693 ! !! x=real(real(i3,dp)/real(n3,dp),dp) 2694 ! !! !x=i3/n3 2695 ! !! if(x < l1) then 2696 ! !! rdiemac(ifft)=diemac 2697 ! !! rdielng(ifft)=dielng 2698 ! !! else if(x < l2) then 2699 ! !! xp=(l2-x)/(l2-l1) 2700 ! !! rdiemac(ifft)=y+(diemac-y)& 2701 ! !! & * (1.-(1.-xp)**4)**4 2702 ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4 2703 ! !! else if(x < l3) then 2704 ! !! rdiemac(ifft)=y 2705 ! !! rdielng(ifft)=zero 2706 ! !! else if(x < l4) then 2707 ! !! xp=(l3-x)/(l3-l4) 2708 ! !! rdiemac(ifft)=y+(diemac-y)& 2709 ! !! & * (1.-(1.-xp)**4)**4 2710 ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4 2711 ! !! else 2712 ! !! rdiemac(ifft)=diemac 2713 ! !! rdielng(ifft)=dielng 2714 ! !! end if 2715 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2716 ! !!!! calculation of atomic core dielectric 2717 ! !! rmin=1e16 2718 ! !! xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)& 2719 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(1,3) 2720 ! !! yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)& 2721 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(2,3) 2722 ! !! zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)& 2723 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(3,3) 2724 ! !! do iatom=1,natom27 2725 ! !! r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2 2726 ! !! if (r<rmin) then 2727 ! !! rmin=r 2728 ! !! end if 2729 ! !! end do 2730 ! !! if(rmin < dr**2) then 2731 ! !! rdiemac(ifft)=min(rdiemac(ifft),core+(diemac-core)*(1.-(1.-sqrt(rmin)/dr)**2)**2) 2732 ! !! rdielng(ifft)=dielng-dielng*(1.-(1.-sqrt(rmin)/dr)**4)**4 2733 ! !! else 2734 ! !! rdiemac(ifft)=min(rdiemac(ifft),diemac) 2735 ! !! end if 2736 rmin=1e16 2737 xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)& 2738 & +real((i3-1),dp)/real(n3,dp)*rprimd(1,3) 2739 yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)& 2740 & +real((i3-1),dp)/real(n3,dp)*rprimd(2,3) 2741 zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)& 2742 & +real((i3-1),dp)/real(n3,dp)*rprimd(3,3) 2743 2744 rdiemac(ifft)=y 2745 rdielng(ifft)=zero 2746 do iatom=1,natom27 2747 r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2 2748 if (r<rmin) then 2749 rmin=r 2750 end if 2751 if(r < l1) then 2752 rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y) 2753 else if(r < l2) then 2754 rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y)*(one-((sqrt(r)-l1)/(l2-l1))**2)**2 2755 else 2756 rdiemac(ifft)=rdiemac(ifft) 2757 end if 2758 if(r < l3) then 2759 rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng) 2760 else if(r < l4) then 2761 rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng) *(one-((sqrt(r)-l3)/(l4-l3))**2)**2 2762 end if 2763 end do 2764 2765 rdielng(ifft)=min(rdielng(ifft),dielng) 2766 ! rdielng(ifft)=dielng 2767 rdiemac(ifft)=min(rdiemac(ifft),diemac) 2768 rdiemac(ifft)=diemac 2769 end do 2770 end if 2771 end do 2772 end do 2773 !rdielng(:)=dielng 2774 2775 !**************************************************************************************** 2776 !**************************************************************************************** 2777 !**************************************************************************************** 2778 !**************************************************************************************** 2779 !****************************************************************** 2780 !compute V1 2781 !****************************************************************** 2782 V1=vresid 2783 call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,rdfuncr=V1,laplacerdfuncr=deltaW) 2784 deltaW(:,1)= (((one/rdiemac(:))*V1(:,1))-(((rdielng(:))**2)*deltaW(:,1))) 2785 !deltaW(:,1)= -diemix*(((rdielng(:))**2)*deltaW(:,ispden)) 2786 if (nspden>1) then 2787 do ispden=2,nspden 2788 deltaW(:,ispden)= (((one/rdiemac(:))*V1(:,ispden))-(((rdielng(:))**2)*deltaW(:,ispden))) 2789 ! deltaW(:,ispden)= -abs(diemixmag)*(((rdielng(:))**2)*deltaW(:,ispden)) 2790 end do 2791 end if 2792 call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat) 2793 call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,& 2794 & frskerker2__newvres2,lastp1*real(1e-6 ,dp),700,V1,rdummy1,rdummy2) 2795 lastp1=min(abs(rdummy1),1e-6_dp) 2796 call frskerker2__end() 2797 2798 !****************************************************************** 2799 !compute V2 2800 !****************************************************************** 2801 V2=vresid 2802 do ispden=1,nspden 2803 deltaW(:,ispden)= (rdielng(:)**2) 2804 end do 2805 call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat) 2806 call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,& 2807 & frskerker2__newvres2,lastp2*real(1e-6,dp),700,V2,rdummy1,rdummy2) 2808 lastp2=min(abs(rdummy1),1e-6_dp) 2809 call frskerker2__end() 2810 2811 2812 !****************************************************************** 2813 !compute C1, C2 & DE 2814 !****************************************************************** 2815 cplex=1; 2816 option=1; 2817 call dotprod_vn(cplex,& !complex density/pot 2818 &rdielng,& !the density 2819 &DE,& !resulting dorproduct integrated over r ! here DE is used has a buffer 2820 &doti,& !imaginary part of the integral 2821 &size(rdielng,1),& !number of localy(cpu) attributed grid point 2822 &nfftotf,& !real total number of grid point 2823 &nspden,& !nspden 2824 &option,& !1=compute only the real part 2=compute also the imaginary part 2825 &rdielng,& !the potential 2826 &ucvol,& !cell volume 2827 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 2828 do ispden=1,nspden 2829 buffer(:,ispden)=rdielng(:)*V1(:,ispden) 2830 end do 2831 call dotprod_vn(cplex,& !complex density/pot 2832 &rdielng,& !the density 2833 &C1,& !resulting dorproduct integrated over r ! here DE is used has a buffer 2834 &doti,& !imaginary part of the integral 2835 &size(rdielng,1),& !number of localy(cpu) attributed grid point 2836 &nfftotf,& !real total number of grid point 2837 &nspden,& !nspden 2838 &option,& !1=compute only the real part 2=compute also the imaginary part 2839 &buffer,& !the potential 2840 &ucvol,& !cell volume 2841 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 2842 do ispden=1,nspden 2843 buffer(:,ispden)=rdielng(:)*V2(:,ispden) 2844 end do 2845 call dotprod_vn(cplex,& !complex density/pot 2846 &rdielng,& !the density 2847 &C2,& !resulting dorproduct integrated over r ! here DE is used has a buffer 2848 &doti,& !imaginary part of the integral 2849 &size(rdielng,1),& !number of localy(cpu) attributed grid point 2850 &nfftotf,& !real total number of grid point 2851 &nspden,& !nspden 2852 &option,& !1=compute only the real part 2=compute also the imaginary part 2853 &buffer,& !the potential 2854 &ucvol,& !cell volume 2855 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 2856 C1=C1/DE 2857 C2=C2/DE 2858 DE=C1/(one-C2) 2859 2860 !****************************************************************** 2861 !compute the new preconditionned residuals 2862 !****************************************************************** 2863 vrespc(:,1)=diemix*(V1(:,1)+DE*V2(:,1)) 2864 if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*(V1(:,2:nspden)+DE*V2(:,2:nspden)) 2865 2866 end subroutine prcrskerker2