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)
PARENTS
linmin
CHILDREN
SOURCE
3186 subroutine bracketing (nv1,nv2,dp_dum_v2dp,v,grad,a,x,b,fa,fx,fb) 3187 3188 3189 !This section has been created automatically by the script Abilint (TD). 3190 !Do not modify the following lines by hand. 3191 #undef ABI_FUNC 3192 #define ABI_FUNC 'bracketing' 3193 !End of the abilint section 3194 3195 implicit none 3196 3197 !Arguments ------------------------------------ 3198 include "dummy_functions.inc" 3199 !scalars 3200 integer,intent(in) :: nv1,nv2 3201 real(dp),intent(inout) :: a,x 3202 real(dp),intent(out) :: b,fa,fb,fx 3203 !arrays 3204 real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2) 3205 3206 !Local variables------------------------------- 3207 !scalars 3208 real(dp),parameter :: maglimit=10000.0_dp 3209 real(dp) :: c,fu,q,r,u,ulim 3210 3211 ! ************************************************************************* 3212 3213 fa=dp_dum_v2dp(nv1,nv2,v(:,:)+(a*grad(:,:))) 3214 fx=dp_dum_v2dp(nv1,nv2,(x*grad(:,:))+v(:,:)) 3215 if(fx > fa) then 3216 c=a 3217 a=x 3218 x=c 3219 c=fa 3220 fa=fx 3221 fx=c 3222 end if 3223 b=x+gold*(x-a) 3224 fb=dp_dum_v2dp(nv1,nv2,(b*grad(:,:))+v(:,:)) 3225 do 3226 if (fx <= fb) return 3227 r=(x-a)*(fx-fb) 3228 q=(x-b)*(fx-fa) 3229 u=x-((x-b)*q-(x-a)*r)/(two*sign(max(abs(q-r),smallest_real),q-r)) 3230 ulim=x+maglimit*(b-x) 3231 if((x-u)*(u-b) > zero) then 3232 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3233 if(fu < fb) then 3234 a=x 3235 fa=fx 3236 x=u 3237 fx=fu 3238 return 3239 else if (fx < fu) then 3240 b=u 3241 fb=fu 3242 return 3243 end if 3244 u=b+gold*(b-x) 3245 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3246 else if((b-u)*(u-ulim) > zero) then 3247 fu=dp_dum_v2dp(nv1,nv2,u*grad(:,:)+v(:,:)) 3248 if(fu<fb) then 3249 x=b 3250 b=u 3251 u=b+gold*(b-x) 3252 fx=fb 3253 fb=fu 3254 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3255 end if 3256 else if((u-ulim)*(ulim-b) >= zero) then 3257 u=ulim 3258 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3259 else 3260 u=b+gold*(b-x) 3261 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3262 end if 3263 a=x 3264 x=b 3265 b=u 3266 fa=fx 3267 fx=fb 3268 fb=fu 3269 end do 3270 3271 end subroutine bracketing
ABINIT/brent [ Functions ]
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
PARENTS
linmin
CHILDREN
dotproduct
SOURCE
3304 function brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,itmax,v,grad,ax,xx,bx,tol,xmin) 3305 3306 3307 !This section has been created automatically by the script Abilint (TD). 3308 !Do not modify the following lines by hand. 3309 #undef ABI_FUNC 3310 #define ABI_FUNC 'brent' 3311 !End of the abilint section 3312 3313 implicit none 3314 3315 !Arguments ------------------------------------ 3316 include "dummy_functions.inc" 3317 !scalars 3318 integer,intent(in) :: itmax,nv1,nv2 3319 real(dp) :: brent 3320 real(dp),intent(in) :: ax,bx,tol,xx 3321 real(dp),intent(out) :: xmin 3322 !arrays 3323 real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2) 3324 3325 !Local variables------------------------------- 3326 !scalars 3327 integer :: iter 3328 real(dp) :: a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,vv,w 3329 real(dp) :: x,xm,zeps 3330 logical :: ok1,ok2,ok3,ok4 3331 3332 !************************************************************************ 3333 zeps=epsilon(ax*real(1e-2,dp)) 3334 a=min(ax,bx) 3335 b=max(ax,bx) 3336 vv=xx 3337 w=xx 3338 x=xx 3339 e=zero 3340 fx=dp_dum_v2dp(nv1,nv2,x*grad(:,:)+v(:,:)) 3341 fv=fx 3342 fw=fx 3343 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of 3344 !v(:,:)=v(:,:)+(grad(:,:)*x) 3345 !but for instance renormilizing the density if brent is used on a density... 3346 !vp(:,:) = v(:,:) 3347 !sub_dum_dp_v2dp_v2dp(x,grad(:,:),vp(:,:) 3348 !dx=dotproduct(v2dp_dum_v2dp(vp(:,:)),grad(:,:)) 3349 dx=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,v(:,:)+x*grad(:,:)),grad(:,:)) 3350 dv=dx 3351 dw=dx 3352 do iter=1,itmax 3353 xm=half*(a+b) 3354 tol1=tol*abs(x)+zeps 3355 tol2=two*tol1 3356 if(abs(x-xm) <= (tol2-half*(b-a))) then 3357 exit 3358 end if 3359 if(abs(e) > tol1) then 3360 d1=two*(b-a) 3361 d2=d1 3362 if(dw /= dx) d1=(w-x)*dx/(dx-dw) 3363 if(dv /= dx) d2=(vv-x)*dx/(dx-dv) 3364 u1=x+d1 3365 u2=x+d2 3366 ok1=((a-u1)*(u1-b)>zero).and.(dx*d1<=zero) 3367 ok2=((a-u2)*(u2-b)>zero).and.(dx*d2<=zero) 3368 olde=e 3369 e=d 3370 if(ok1.or.ok2) then 3371 if(ok1.and.ok2) then 3372 d=merge(d1,d2,abs(d1)<abs(d2)) 3373 else 3374 d=merge(d1,d2,ok1) 3375 end if 3376 if(abs(d)<=abs(half*olde)) then 3377 u=x+d 3378 if(((u-a)<tol2).or.((b-u)<tol2)) d=sign(tol1,xm-x) 3379 else 3380 e=merge(a,b,dx>=zero)-x 3381 d=half*e 3382 end if 3383 else 3384 e=merge(a,b,dx>=zero)-x 3385 d=half*e 3386 end if 3387 else 3388 e=merge(a,b,dx>=zero)-x 3389 d=half*e 3390 end if 3391 3392 if(abs(d) >=tol1)then 3393 u=x+d 3394 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3395 else 3396 u=x+sign(tol1,d) 3397 fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)) 3398 if(fu>fx) then 3399 exit 3400 end if 3401 end if 3402 du=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)),grad(:,:)) 3403 if(fu<=fx)then 3404 if(u>=x)then 3405 a=x 3406 else 3407 b=x 3408 end if 3409 vv=w 3410 fv=fw 3411 dv=dw 3412 w=x 3413 fw=fx 3414 dw=dx 3415 x=u 3416 dx=du 3417 fx=fu 3418 else 3419 if(u<x) then 3420 a=u 3421 else 3422 b=u 3423 end if 3424 ok3=(w==x).or.(fu.le.fw) 3425 ok4=(vv==w).or.(vv==x).or.(fu.lt.fv) 3426 if(ok3) then 3427 vv=w 3428 fv=fw 3429 dv=dw 3430 w=u 3431 fw=fu 3432 dw=du 3433 else if( ok4 ) then 3434 vv=u 3435 fv=fu 3436 dv=du 3437 end if 3438 end if 3439 end do 3440 xmin=x 3441 !the function sub_dum_dp_v2dp_v2dp must do the equivalent of 3442 !v(:,:)=v(:,:)+(grad(:,:)*x) 3443 !but for instance renormilizing the density if brent is used on a density... 3444 call sub_dum_dp_v2dp_v2dp(nv1,nv2,x,grad(:,:),v(:,:)) 3445 brent=fx 3446 3447 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
PARENTS
prcrskerker1,prcrskerker2
CHILDREN
linmin
SOURCE
3021 subroutine cgpr(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,dtol,itmax,v,fmin,delta) 3022 3023 3024 !This section has been created automatically by the script Abilint (TD). 3025 !Do not modify the following lines by hand. 3026 #undef ABI_FUNC 3027 #define ABI_FUNC 'cgpr' 3028 !End of the abilint section 3029 3030 implicit none 3031 3032 !Arguments ------------------------------------ 3033 include "dummy_functions.inc" 3034 !scalars 3035 integer,intent(in) :: itmax,nv1,nv2 3036 real(dp),intent(in) :: dtol 3037 real(dp),intent(out) :: delta,fmin 3038 !arrays 3039 real(dp),intent(inout) :: v(nv1,nv2) 3040 3041 !Local variables------------------------------- 3042 !scalars 3043 integer :: iiter 3044 real(dp) :: fv,gam,gscal,gscal2,sto 3045 !arrays 3046 real(dp) :: grad0(nv1,nv2),grad1(nv1,nv2),grad2(nv1,nv2),grad3(nv1,nv2) 3047 !no_abirules 3048 3049 !************************************************************************ 3050 fv = dp_dum_v2dp(nv1,nv2,v(:,:)) 3051 grad0(:,:) = -v2dp_dum_v2dp(nv1,nv2,v(:,:)) 3052 grad1(:,:) = grad0(:,:) 3053 grad2(:,:) = grad0(:,:) 3054 do iiter=1,itmax 3055 call linmin(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,v,grad0,fmin) 3056 ! return if the min is reached 3057 sto=dtol*(abs(fmin)+abs(fv)+tol14) 3058 delta=abs(fv-fmin) 3059 delta=abs(delta) 3060 if((delta.lt.sto).or.(iiter==itmax)) then 3061 ! DEBUG 3062 ! write(std_out,*) 'cgpr (01cg) : stop cond for cgpr:',sto,'delta:',delta,'fv:',fv,'fmin:',fmin 3063 ! ENDDEBUG 3064 return 3065 end if 3066 ! a new step 3067 fv=fmin 3068 grad0(:,:)=v2dp_dum_v2dp(nv1,nv2,v(:,:)) 3069 gscal=dotproduct(nv1,nv2,grad1(:,:),grad1(:,:)) 3070 grad3(:,:)=grad0(:,:)+grad1(:,:) 3071 gscal2=dotproduct(nv1,nv2,grad3(:,:),grad0(:,:)) 3072 gam=gscal2/gscal 3073 grad1(:,:)=-grad0(:,:) 3074 grad2(:,:)=grad1(:,:)+gam*grad2(:,:) 3075 grad0(:,:)=grad2(:,:) 3076 ! DEBUG 3077 ! write(std_out,*) 'cgpr (01cg) :=================================================================================' 3078 ! write(std_out,*) 'cgpr (01cg) : step',iiter,'delta:',delta ,'fv',fv,'fmin',fmin 3079 ! write(std_out,*) 'cgpr (01cg) :=================================================================================' 3080 ! ENDDEBUG 3081 end do 3082 3083 end subroutine cgpr
ABINIT/dielmt [ Functions ]
NAME
dielmt
FUNCTION
Compute dielectric matrix from susceptibility matrix Diagonalize it, then invert it.
INPUTS
gmet(3,3)=reciprocal space metric tensor in bohr**-2. kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. npwdiel=size of the dielinv and susmat arrays. nspden=number of spin-density components occopt=option for occupancies prtvol=control print volume and debugging output susmat(2,npwdiel,nspden,npwdiel,nspden)= the susceptibility (or density-density response) matrix in reciprocal space
OUTPUT
dielinv(2,npwdiel,(nspden+4)/3,npwdiel,(nspden+4)/3)=inverse of the (non-hermitian)
TC dielectric matrix in reciprocal space.
NOTES
Warning : will not work in the spin-polarized, metallic case. Output (not cleaned) !!! Spin behaviour is not obvious !!!
TODO
Write equation below (hermitian matrix)
PARENTS
prcref,prcref_PMA
CHILDREN
timab,wrtout,zhpev
SOURCE
1626 subroutine dielmt(dielinv,gmet,kg_diel,npwdiel,nspden,occopt,prtvol,susmat) 1627 1628 1629 !This section has been created automatically by the script Abilint (TD). 1630 !Do not modify the following lines by hand. 1631 #undef ABI_FUNC 1632 #define ABI_FUNC 'dielmt' 1633 !End of the abilint section 1634 1635 implicit none 1636 1637 !Arguments ------------------------------------ 1638 !scalars 1639 integer,intent(in) :: npwdiel,nspden,occopt,prtvol 1640 !arrays 1641 integer,intent(in) :: kg_diel(3,npwdiel) 1642 real(dp),intent(in) :: gmet(3,3) 1643 real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 1644 real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden) 1645 1646 !Local variables------------------------------- 1647 !scalars 1648 integer :: ieig,ier,ii,index,ipw,ipw1,ipw2,isp,jj,npwsp 1649 real(dp) :: ai1,ai2,ar1,ar2,eiginv,gfact,gfactinv,gred1,gred2,gred3,gsquar 1650 real(dp) :: tpisq 1651 character(len=500) :: message 1652 !arrays 1653 real(dp) :: tsec(2) 1654 real(dp),allocatable :: dielh(:),dielmat(:,:,:,:,:),dielvec(:,:,:) 1655 real(dp),allocatable :: eig_diel(:),zhpev1(:,:),zhpev2(:) 1656 !no_abirules 1657 !integer :: ipw3 1658 !real(dp) :: elementi,elementr 1659 1660 ! ************************************************************************* 1661 1662 !DEBUG 1663 !write(std_out,*)' dielmt : enter ' 1664 !if(.true.)stop 1665 !ENDDEBUG 1666 1667 !tpisq is (2 Pi) **2: 1668 tpisq=(two_pi)**2 1669 1670 call timab(90,1,tsec) 1671 1672 !-Compute now the hermitian dielectric matrix------------------------------ 1673 !Following remarks are only valid within RPA approximation (Kxc=0): 1674 1675 !for the spin-unpolarized case, 1 - 4pi (1/G) chi0(G,Gp) (1/Gp) 1676 1677 !for the spin-polarized case, 1678 !( 1 0 ) - 4pi ( 1/G 1/G ) ( chi0 upup chi0 updn ) ( 1/Gp 1/Gp ) 1679 !( 0 1 ) ( 1/G 1/G ) ( chi0 dnup chi0 dndn ) ( 1/Gp 1/Gp ) 1680 !which is equal to 1681 !( 1 0 ) - 4pi (1/G 0 ) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) (1/Gp 0 ) 1682 !( 0 1 ) ( 0 1/G) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) ( 0 1/Gp) 1683 !So, if spin-polarized, sum all spin contributions 1684 !Note: chi0 updn = chi0 dnup = zero for non-metallic systems 1685 1686 !In the case of non-collinear magnetism, within RPA, this is the same because: 1687 !chi0_(s1,s2),(s3,s4) = delta_s1,s2 * delta_s3,s4 * chi0_(s1,s1),(s3,s3) 1688 !Only chi_upup,upup, chi_dndn,dndn, chi_upup,dndn and chi_dndn,upup 1689 !have to be taken into account (stored, susmat(:,ipw1,1:2,ipw2,1:2) 1690 1691 ABI_ALLOCATE(dielmat,(2,npwdiel,min(nspden,2),npwdiel,min(nspden,2))) 1692 1693 if(nspden/=1)then 1694 if (occopt<3) then 1695 do ipw2=1,npwdiel 1696 do ipw1=1,npwdiel 1697 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2) 1698 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2) 1699 end do 1700 end do 1701 else 1702 do ipw2=1,npwdiel 1703 do ipw1=1,npwdiel 1704 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2)+susmat(1,ipw1,1,ipw2,2)+susmat(1,ipw1,2,ipw2,1) 1705 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2)+susmat(2,ipw1,1,ipw2,2)+susmat(2,ipw1,2,ipw2,1) 1706 end do 1707 end do 1708 end if 1709 else 1710 do ipw2=1,npwdiel 1711 do ipw1=1,npwdiel 1712 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1) 1713 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1) 1714 end do 1715 end do 1716 end if 1717 !Compute 1/G factors and include them in the dielectric matrix 1718 do ipw1=1,npwdiel 1719 gred1=dble(kg_diel(1,ipw1)) 1720 gred2=dble(kg_diel(2,ipw1)) 1721 gred3=dble(kg_diel(3,ipw1)) 1722 gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 & 1723 & +two*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 + & 1724 & gmet(2,3)*gred2*gred3) ) 1725 ! Distinguish G=0 from other elements 1726 if(gsquar>tol12)then 1727 ! !$ gfact=\sqrt (4.0_dp \pi/gsquar/dble(nspden))$ 1728 gfact=sqrt(four_pi/gsquar) 1729 do ipw2=1,npwdiel 1730 ! Must multiply both rows and columns, and also changes the sign 1731 dielmat(1,ipw2,1,ipw1,1)=-dielmat(1,ipw2,1,ipw1,1)*gfact 1732 dielmat(2,ipw2,1,ipw1,1)=-dielmat(2,ipw2,1,ipw1,1)*gfact 1733 dielmat(1,ipw1,1,ipw2,1)= dielmat(1,ipw1,1,ipw2,1)*gfact 1734 dielmat(2,ipw1,1,ipw2,1)= dielmat(2,ipw1,1,ipw2,1)*gfact 1735 end do 1736 else 1737 ! Zero the G=0 elements, head and wings 1738 do ipw2=1,npwdiel 1739 dielmat(1,ipw2,1,ipw1,1)=zero 1740 dielmat(2,ipw2,1,ipw1,1)=zero 1741 dielmat(1,ipw1,1,ipw2,1)=zero 1742 dielmat(2,ipw1,1,ipw2,1)=zero 1743 end do 1744 end if 1745 end do 1746 1747 !Complete the matrix in the spin-polarized case 1748 !should this be nspden==2?? 1749 if(nspden/=1)then 1750 do ipw1=1,npwdiel 1751 do ipw2=1,npwdiel 1752 dielmat(1,ipw1,1,ipw2,2)=dielmat(1,ipw1,1,ipw2,1) 1753 dielmat(2,ipw1,1,ipw2,2)=dielmat(2,ipw1,1,ipw2,1) 1754 dielmat(1,ipw1,2,ipw2,1)=dielmat(1,ipw1,1,ipw2,1) 1755 dielmat(2,ipw1,2,ipw2,1)=dielmat(2,ipw1,1,ipw2,1) 1756 dielmat(1,ipw1,2,ipw2,2)=dielmat(1,ipw1,1,ipw2,1) 1757 dielmat(2,ipw1,2,ipw2,2)=dielmat(2,ipw1,1,ipw2,1) 1758 end do 1759 end do 1760 end if 1761 1762 !DEBUG 1763 !write(std_out,*)' dielmt : make dielmat equal to identity matrix ' 1764 !do ipw1=1,npwdiel 1765 !do ipw2=1,npwdiel 1766 !dielmat(1,ipw1,1,ipw2,1)=0.0_dp 1767 !dielmat(2,ipw1,1,ipw2,1)=0.0_dp 1768 !end do 1769 !end do 1770 !ENDDEBUG 1771 1772 !Add the diagonal part 1773 do isp=1,min(nspden,2) 1774 do ipw=1,npwdiel 1775 dielmat(1,ipw,isp,ipw,isp)=one+dielmat(1,ipw,isp,ipw,isp) 1776 end do 1777 end do 1778 1779 !-The hermitian dielectric matrix is computed ------------------------------ 1780 !-Now, diagonalize it ------------------------------------------------------ 1781 1782 !In RPA, everything is projected on the spin-symmetrized 1783 !space. This was coded here (for the time being). 1784 1785 !Diagonalize the hermitian dielectric matrix 1786 1787 !npwsp=npwdiel*nspden 1788 npwsp=npwdiel 1789 1790 ABI_ALLOCATE(dielh,(npwsp*(npwsp+1))) 1791 ABI_ALLOCATE(dielvec,(2,npwsp,npwsp)) 1792 ABI_ALLOCATE(eig_diel,(npwsp)) 1793 ABI_ALLOCATE(zhpev1,(2,2*npwsp-1)) 1794 ABI_ALLOCATE(zhpev2,(3*npwsp-2)) 1795 ier=0 1796 !Store the dielectric matrix in proper mode before calling zhpev 1797 index=1 1798 do ii=1,npwdiel 1799 do jj=1,ii 1800 dielh(index )=dielmat(1,jj,1,ii,1) 1801 dielh(index+1)=dielmat(2,jj,1,ii,1) 1802 index=index+2 1803 end do 1804 end do 1805 !If spin-polarized and non RPA, need to store other parts of the matrix 1806 !if(nspden/=1)then 1807 !do ii=1,npwdiel 1808 !Here, spin-flip contribution 1809 !do jj=1,npwdiel 1810 !dielh(index )=dielmat(1,jj,1,ii,2) 1811 !dielh(index+1)=dielmat(2,jj,1,ii,2) 1812 !index=index+2 1813 !end do 1814 !Here spin down-spin down upper matrix 1815 !do jj=1,ii 1816 !dielh(index )=dielmat(1,jj,2,ii,2) 1817 !dielh(index+1)=dielmat(2,jj,2,ii,2) 1818 !index=index+2 1819 !end do 1820 !end do 1821 !end if 1822 1823 call ZHPEV ('V','U',npwsp,dielh,eig_diel,dielvec,npwdiel,zhpev1,& 1824 & zhpev2,ier) 1825 ABI_DEALLOCATE(zhpev1) 1826 ABI_DEALLOCATE(zhpev2) 1827 1828 if(prtvol>=10)then 1829 write(message, '(a,a,a,5es12.4)' )ch10,& 1830 & ' Five largest eigenvalues of the hermitian RPA dielectric matrix:',& 1831 & ch10,eig_diel(npwdiel:npwdiel-4:-1) 1832 call wrtout(ab_out,message,'COLL') 1833 end if 1834 1835 write(message, '(a,a)' )ch10,& 1836 & ' dielmt : 15 largest eigenvalues of the hermitian RPA dielectric matrix' 1837 call wrtout(std_out,message,'COLL') 1838 write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(npwdiel:npwdiel-4:-1) 1839 call wrtout(std_out,message,'COLL') 1840 write(message, '(a,5es12.5)' )' 6-10 :',eig_diel(npwdiel-5:npwdiel-9:-1) 1841 call wrtout(std_out,message,'COLL') 1842 write(message, '(a,5es12.5)' )' 11-15:',eig_diel(npwdiel-10:npwdiel-14:-1) 1843 call wrtout(std_out,message,'COLL') 1844 write(message, '(a,a)' )ch10,& 1845 & ' dielmt : 5 smallest eigenvalues of the hermitian RPA dielectric matrix' 1846 call wrtout(std_out,message,'COLL') 1847 write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(1:5) 1848 call wrtout(std_out,message,'COLL') 1849 1850 !Invert the hermitian dielectric matrix, 1851 !Should use a BLAS call ! 1852 do ipw2=1,npwdiel 1853 do ipw1=ipw2,npwdiel 1854 dielinv(1,ipw1,1,ipw2,1)=zero 1855 dielinv(2,ipw1,1,ipw2,1)=zero 1856 end do 1857 end do 1858 do ieig=1,npwdiel 1859 eiginv=one/eig_diel(ieig) 1860 do ipw2=1,npwdiel 1861 do ipw1=ipw2,npwdiel 1862 ar1=dielvec(1,ipw1,ieig) 1863 ai1=dielvec(2,ipw1,ieig) 1864 ar2=dielvec(1,ipw2,ieig) 1865 ai2=dielvec(2,ipw2,ieig) 1866 dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)+& 1867 & (ar1*ar2+ai1*ai2)*eiginv 1868 dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)+& 1869 & (ai1*ar2-ar1*ai2)*eiginv 1870 end do 1871 end do 1872 end do 1873 do ipw2=1,npwdiel-1 1874 do ipw1=ipw2+1,npwdiel 1875 dielinv(1,ipw2,1,ipw1,1)= dielinv(1,ipw1,1,ipw2,1) 1876 dielinv(2,ipw2,1,ipw1,1)=-dielinv(2,ipw1,1,ipw2,1) 1877 end do 1878 end do 1879 1880 ABI_DEALLOCATE(dielh) 1881 ABI_DEALLOCATE(dielvec) 1882 ABI_DEALLOCATE(eig_diel) 1883 1884 !DEBUG 1885 !Checks whether the inverse of the hermitian dielectric matrix 1886 !has been correctly generated 1887 !do ipw1=1,npwdiel 1888 !do ipw2=1,npwdiel 1889 !elementr=0.0_dp 1890 !elementi=0.0_dp 1891 !do ipw3=1,npwdiel 1892 !elementr=elementr+dielinv(1,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1)& 1893 !& -dielinv(2,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1) 1894 !elementi=elementi+dielinv(1,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1)& 1895 !& +dielinv(2,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1) 1896 !end do 1897 !if(elementr**2+elementi**2 > 1.0d-12)then 1898 !if( ipw1 /= ipw2 .or. & 1899 !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 1900 !write(std_out,*)' dielmt : the inversion procedure is not correct ' 1901 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 1902 !write(std_out,*)' elementr,elementi=',elementr,elementi 1903 !stop 1904 !end if 1905 !end if 1906 !end do 1907 !end do 1908 !write(std_out,*)'dielmt : matrix has been inverted successfully ' 1909 !stop 1910 !ENDDEBUG 1911 1912 !Then get the inverse of the asymmetric 1913 !dielectric matrix, as required for the preconditioning. 1914 1915 !Inverse of the dielectric matrix : ( 1 - 4pi (1/G^2) chi0(G,Gp) )^(-1) 1916 !In dielinv there is now (1 - 4pi (1/G) chi0(G,Gp) (1/Gp) )^(-1) 1917 !So, evaluate dielinv_after(G,Gp) = 1918 !(4pi/G^2)^(1/2) dielinv_before(G,Gp) (4pi/Gp^2)^(-1/2) 1919 !In RPA, can focus on the spin-averaged quantities 1920 do ipw1=1,npwdiel 1921 gred1=dble(kg_diel(1,ipw1)) 1922 gred2=dble(kg_diel(2,ipw1)) 1923 gred3=dble(kg_diel(3,ipw1)) 1924 gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 & 1925 & +two*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 + & 1926 & gmet(2,3)*gred2*gred3) ) 1927 ! Distinguish G=0 from other elements 1928 if(gsquar>tol12)then 1929 gfact=sqrt(four_pi/gsquar) 1930 gfactinv=one/gfact 1931 do ipw2=1,npwdiel 1932 ! Must multiply both rows and columns 1933 dielinv(1,ipw2,1,ipw1,1)=dielinv(1,ipw2,1,ipw1,1)*gfactinv 1934 dielinv(2,ipw2,1,ipw1,1)=dielinv(2,ipw2,1,ipw1,1)*gfactinv 1935 dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)*gfact 1936 dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)*gfact 1937 end do 1938 else 1939 ! Zero the G=0 elements, head 1940 do ipw2=1,npwdiel 1941 if (ipw2/=ipw1) dielinv(1:2,ipw1,1,ipw2,1)=zero 1942 end do 1943 end if 1944 end do 1945 1946 ABI_DEALLOCATE(dielmat) 1947 1948 call timab(90,2,tsec) 1949 1950 end subroutine dielmt
ABINIT/dieltcel [ Functions ]
NAME
dieltcel
FUNCTION
Compute either test charge or electronic dielectric matrices from susceptibility matrix Diagonalize it, then invert it.
INPUTS
gmet(3,3)=reciprocal space metric tensor in bohr**-2. kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. kxc(nfft,nkxc)=exchange-correlation kernel, needed if the electronic dielectric matrix is computed nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkxc=second dimension of the array kxc, see rhotoxc.f for a description npwdiel=size of the dielinv and susmat arrays. nspden=number of spin-density components occopt=option for occupancies option=1 for Test Charge dielectric matrix, 2 for electronic dielectric matrix prtvol=control print volume and debugging output susmat(2,npwdiel,nspden,npwdiel,nspden)= the susceptibility (or density-density response) matrix in reciprocal space
OUTPUT
dielinv(2,npwdiel,nspden,npwdiel,nspden)=inverse of the (non-hermitian)
TC dielectric matrix in reciprocal space.
NOTES
Output (not cleaned) !!! Spin behaviour is not obvious !!! Will not work in the spin-polarized, metallic case.
PARENTS
prcref,prcref_PMA
CHILDREN
destroy_mpi_enreg,fourdp,init_distribfft_seq,initmpi_seq,timab,wrtout zhpev
SOURCE
1997 subroutine dieltcel(dielinv,gmet,kg_diel,kxc,& 1998 & nfft,ngfft,nkxc,npwdiel,nspden,occopt,option,paral_kgb,prtvol,susmat) 1999 2000 2001 !This section has been created automatically by the script Abilint (TD). 2002 !Do not modify the following lines by hand. 2003 #undef ABI_FUNC 2004 #define ABI_FUNC 'dieltcel' 2005 !End of the abilint section 2006 2007 implicit none 2008 2009 !Arguments ------------------------------------ 2010 !scalars 2011 integer,intent(in) :: nfft,nkxc,npwdiel,nspden,occopt,option,paral_kgb 2012 integer,intent(in) :: prtvol 2013 !arrays 2014 integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18) 2015 real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc) 2016 real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 2017 real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden) 2018 2019 !Local variables------------------------------- 2020 !scalars 2021 integer :: i1,i2,i3,ieig,ier,ifft,ii,index,ipw0,ipw1,ipw2,ispden,j1 2022 integer :: j2,j3,jj,k1,k2,k3,n1,n2,n3 2023 real(dp) :: ai,ai2,ar,ar2,eiginv,gred1,gred2,gred3,gsquar,si 2024 real(dp) :: sr,tpisq 2025 character(len=500) :: message 2026 type(MPI_type) :: mpi_enreg_seq 2027 !arrays 2028 real(dp) :: tsec(2) 2029 real(dp),allocatable :: eig_msusinvsqr(:),eig_msussqr(:) 2030 real(dp),allocatable :: eig_sus(:),eig_sym(:),invsqrsus(:,:,:,:,:) 2031 real(dp),allocatable :: khxc(:,:,:,:,:),kxcg(:,:),sqrsus(:,:,:,:,:),sush(:) 2032 real(dp),allocatable :: susvec(:,:,:),symdielmat(:,:,:,:,:),symh(:) 2033 real(dp),allocatable :: symvec(:,:,:,:,:),wkxc(:),work(:,:,:,:,:) 2034 real(dp),allocatable :: work2(:,:,:,:,:),zhpev1(:,:),zhpev2(:) 2035 !no_abirules 2036 !integer :: ipw3 2037 !real(dp) :: elementi,elementr 2038 !DEBUG 2039 !Used to moderate divergence effect near rho=0 of the Kxc 2040 !this limit value is truly empirical (exprmt on small Sr cell). 2041 !real(dp) :: kxc_min=-200.0 2042 !ENDDEBUG 2043 2044 ! ************************************************************************* 2045 2046 call timab(96,1,tsec) 2047 2048 !tpisq is (2 Pi) **2: 2049 tpisq=(two_pi)**2 2050 2051 if(nspden/=1 .and. (occopt>=3 .and. occopt<=8) )then 2052 write(message, '(a,a,a)' )& 2053 & 'In the present version of the code, one cannot produce',ch10,& 2054 & 'the dielectric matrix in the metallic, spin-polarized case.' 2055 MSG_BUG(message) 2056 end if 2057 2058 if(nspden==4)then 2059 write(message,'(a,a,a)')& 2060 & 'In the present version of the code, one cannot produce',ch10,& 2061 & 'the dielectric matrix in the non-collinear spin-polarized case.' 2062 MSG_ERROR(message) 2063 end if 2064 2065 2066 !-Diagonalize the susceptibility matrix 2067 2068 ABI_ALLOCATE(sush,(npwdiel*(npwdiel+1))) 2069 ABI_ALLOCATE(susvec,(2,npwdiel,npwdiel)) 2070 ABI_ALLOCATE(eig_msusinvsqr,(npwdiel)) 2071 ABI_ALLOCATE(eig_msussqr,(npwdiel)) 2072 ABI_ALLOCATE(eig_sus,(npwdiel)) 2073 ABI_ALLOCATE(zhpev1,(2,2*npwdiel-1)) 2074 ABI_ALLOCATE(zhpev2,(3*npwdiel-2)) 2075 ABI_ALLOCATE(work,(2,npwdiel,nspden,npwdiel,nspden)) 2076 ABI_ALLOCATE(work2,(2,npwdiel,nspden,npwdiel,nspden)) 2077 ABI_ALLOCATE(sqrsus,(2,npwdiel,nspden,npwdiel,nspden)) 2078 ABI_ALLOCATE(invsqrsus,(2,npwdiel,nspden,npwdiel,nspden)) 2079 2080 !At some time, should take care of different spin channels 2081 do ispden=1,nspden 2082 2083 if(nspden/=1)then 2084 MSG_ERROR('dieltcel : stop, nspden/=1') 2085 end if 2086 2087 ! Store the susceptibility matrix in proper mode before calling zhpev 2088 index=1 2089 do ii=1,npwdiel 2090 do jj=1,ii 2091 sush(index )=susmat(1,jj,1,ii,1) 2092 sush(index+1)=susmat(2,jj,1,ii,1) 2093 index=index+2 2094 end do 2095 end do 2096 2097 ier=0 2098 call ZHPEV ('V','U',npwdiel,sush,eig_sus,susvec,npwdiel,zhpev1,zhpev2,ier) 2099 2100 ! DEBUG 2101 ! write(std_out,*)' dieltcel : print eigenvalues of the susceptibility matrix' 2102 ! do ii=1,npwdiel 2103 ! write(std_out,'(i5,es16.6)' )ii,eig_sus(ii) 2104 ! end do 2105 ! ENDDEBUG 2106 2107 do ii=1,npwdiel 2108 if(-eig_sus(ii)>1.d-12)then 2109 eig_msussqr(ii)=sqrt(-eig_sus(ii)) 2110 eig_msusinvsqr(ii)=1._dp/eig_msussqr(ii) 2111 else if(-eig_sus(ii)< -1.d-12)then 2112 message = "Found positive eigenvalue of susceptibility matrix." 2113 MSG_BUG(message) 2114 else 2115 ! Set the eigenvalue corresponding to a constant potential change to 1, 2116 ! while it will be set to zero in Khx. 2117 eig_msussqr(ii)=1._dp 2118 eig_msusinvsqr(ii)=1._dp 2119 end if 2120 end do 2121 2122 ! Compute square root of minus susceptibility matrix 2123 ! and inverse square root of minus susceptibility matrix 2124 do ii=1,npwdiel 2125 work(:,:,1,ii,1)=susvec(:,:,ii)*eig_msussqr(ii) 2126 work2(:,:,1,ii,1)=susvec(:,:,ii)*eig_msusinvsqr(ii) 2127 end do 2128 do ipw2=1,npwdiel 2129 do ipw1=ipw2,npwdiel 2130 ar=0._dp ; ai=0._dp ; ar2=0._dp ; ai2=0._dp 2131 do ii=1,npwdiel 2132 sr=susvec(1,ipw2,ii) ; si=susvec(2,ipw2,ii) 2133 ar =ar +work(1,ipw1,1,ii,1)*sr +work(2,ipw1,1,ii,1)*si 2134 ai =ai +work(2,ipw1,1,ii,1)*sr -work(1,ipw1,1,ii,1)*si 2135 ar2=ar2 +work2(1,ipw1,1,ii,1)*sr +work2(2,ipw1,1,ii,1)*si 2136 ai2=ai2 +work2(2,ipw1,1,ii,1)*sr -work2(1,ipw1,1,ii,1)*si 2137 end do 2138 sqrsus(1,ipw1,1,ipw2,1)=ar 2139 sqrsus(2,ipw1,1,ipw2,1)=ai 2140 invsqrsus(1,ipw1,1,ipw2,1)=ar2 2141 invsqrsus(2,ipw1,1,ipw2,1)=ai2 2142 if(ipw1/=ipw2)then 2143 sqrsus(1,ipw2,1,ipw1,1)=ar 2144 sqrsus(2,ipw2,1,ipw1,1)=-ai 2145 invsqrsus(1,ipw2,1,ipw1,1)=ar2 2146 invsqrsus(2,ipw2,1,ipw1,1)=-ai2 2147 end if 2148 end do 2149 end do 2150 2151 ! DEBUG 2152 ! Checks whether sqrsus and invsqrsus are inverse of each other. 2153 ! do ipw1=1,npwdiel 2154 ! do ipw2=1,npwdiel 2155 ! elementr=0.0_dp 2156 ! elementi=0.0_dp 2157 ! do ipw3=1,npwdiel 2158 ! elementr=elementr+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1)& 2159 ! & -sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1) 2160 ! elementi=elementi+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1)& 2161 ! & +sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1) 2162 ! end do 2163 ! if(elementr**2+elementi**2 > 1.0d-12)then 2164 ! if( ipw1 /= ipw2 .or. & 2165 ! & ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 2166 ! write(std_out,*)' dieltcel : sqrsus and invsqrsus are not (pseudo)',& 2167 ! & 'inverse of each other' 2168 ! write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 2169 ! write(std_out,*)' elementr,elementi=',elementr,elementi 2170 ! stop 2171 ! end if 2172 ! end if 2173 ! end do 2174 ! end do 2175 ! ENDDEBUG 2176 2177 ! End loop over spins 2178 end do 2179 2180 ABI_DEALLOCATE(eig_msusinvsqr) 2181 ABI_DEALLOCATE(eig_msussqr) 2182 ABI_DEALLOCATE(eig_sus) 2183 ABI_DEALLOCATE(sush) 2184 ABI_DEALLOCATE(susvec) 2185 2186 !-Compute the Hxc kernel 2187 2188 ABI_ALLOCATE(khxc,(2,npwdiel,nspden,npwdiel,nspden)) 2189 ABI_ALLOCATE(symdielmat,(2,npwdiel,nspden,npwdiel,nspden)) 2190 2191 khxc(:,:,:,:,:)=0.0_dp 2192 2193 !Compute Hartree kernel 2194 do ipw1=1,npwdiel 2195 gred1=dble(kg_diel(1,ipw1)) 2196 gred2=dble(kg_diel(2,ipw1)) 2197 gred3=dble(kg_diel(3,ipw1)) 2198 gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 & 2199 & +2.0_dp*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 + & 2200 & gmet(2,3)*gred2*gred3) ) 2201 ! Distinguish G=0 from other elements 2202 if(gsquar>1.0d-12)then 2203 khxc(1,ipw1,1,ipw1,1)= 4.0_dp*pi/gsquar 2204 else 2205 ! G=0 2206 ipw0=ipw1 2207 end if 2208 end do 2209 2210 !Eventually add the xc part 2211 if(option>=2)then 2212 2213 ABI_ALLOCATE(wkxc,(nfft)) 2214 ABI_ALLOCATE(kxcg,(2,nfft)) 2215 wkxc(:)=kxc(:,1) 2216 ! DEBUG 2217 ! Used to moderate divergenc effect near rho=0 of the Kxc (see above). 2218 ! wkxc(:)=merge(kxc(:,1), kxc_min, kxc(:,1) > kxc_min) 2219 ! ENDDEBUG 2220 call initmpi_seq(mpi_enreg_seq) 2221 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all') 2222 call fourdp(1,kxcg,wkxc,-1,mpi_enreg_seq,nfft,ngfft,paral_kgb,0) ! trsfrm R to G 2223 call destroy_mpi_enreg(mpi_enreg_seq) 2224 2225 ! Compute difference in G vectors 2226 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 2227 do ipw2=1,npwdiel 2228 if(ipw2/=ipw0)then 2229 2230 j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2) 2231 ! Fills diagonal 2232 khxc(1,ipw2,1,ipw2,1)=khxc(1,ipw2,1,ipw2,1)+kxcg(1,1) 2233 khxc(2,ipw2,1,ipw2,1)=khxc(2,ipw2,1,ipw2,1)+kxcg(2,1) 2234 2235 if(ipw2/=npwdiel)then 2236 ! Fills off-diagonal part of the matrix, except G=0 2237 do ipw1=ipw2+1,npwdiel 2238 if(ipw1/=ipw0)then 2239 i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1) 2240 ! Use of two mod calls handles both i1-j1>=ndiel1 AND i1-j1<0 2241 k1=mod(n1+mod(i1-j1,n1),n1) 2242 k2=mod(n2+mod(i2-j2,n2),n2) 2243 k3=mod(n3+mod(i3-j3,n3),n3) 2244 ifft=k1+1+n1*(k2+n2*k3) 2245 ! The signs of imaginary contributions have been checked 2246 khxc(1,ipw1,1,ipw2,1)=kxcg(1,ifft) 2247 khxc(2,ipw1,1,ipw2,1)=kxcg(2,ifft) 2248 khxc(1,ipw2,1,ipw1,1)=kxcg(1,ifft) 2249 khxc(2,ipw2,1,ipw1,1)=-kxcg(2,ifft) 2250 end if 2251 end do 2252 end if 2253 2254 end if 2255 end do 2256 2257 ABI_DEALLOCATE(wkxc) 2258 ABI_DEALLOCATE(kxcg) 2259 2260 ! Endif option 2 2261 end if 2262 2263 !Now, get the symmetric dielectric matrix 2264 !Premultiplication by square root of minus susceptibility matrix 2265 do ipw2=1,npwdiel 2266 do ipw1=1,npwdiel 2267 ar=0._dp ; ai=0._dp 2268 do ii=1,npwdiel 2269 ar=ar+sqrsus(1,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) & 2270 & -sqrsus(2,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1) 2271 ai=ai+sqrsus(2,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) & 2272 & +sqrsus(1,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1) 2273 end do 2274 work(1,ipw1,1,ipw2,1)=ar 2275 work(2,ipw1,1,ipw2,1)=ai 2276 end do 2277 end do 2278 !Postmultiplication by square root of minus susceptibility matrix 2279 do ipw2=1,npwdiel 2280 ! do ipw1=ipw2,npwdiel 2281 do ipw1=1,npwdiel 2282 ar=0._dp ; ai=0._dp 2283 do ii=1,npwdiel 2284 ar=ar+work(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2285 & -work(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2286 ai=ai+work(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2287 & +work(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2288 end do 2289 symdielmat(1,ipw1,1,ipw2,1)=ar 2290 symdielmat(2,ipw1,1,ipw2,1)=ai 2291 ! if(ipw1/=ipw2)then 2292 ! symdielmat(1,ipw2,1,ipw1,1)=ar 2293 ! symdielmat(2,ipw2,1,ipw1,1)=-ai 2294 ! end if 2295 end do 2296 ! Add the unity matrix 2297 symdielmat(1,ipw2,1,ipw2,1)=1._dp+symdielmat(1,ipw2,1,ipw2,1) 2298 end do 2299 2300 ABI_DEALLOCATE(khxc) 2301 2302 ABI_ALLOCATE(symh,(npwdiel*(npwdiel+1))) 2303 ABI_ALLOCATE(symvec,(2,npwdiel,nspden,npwdiel,nspden)) 2304 ABI_ALLOCATE(eig_sym,(npwdiel)) 2305 2306 !Store the symmetrized dielectric matrix in proper mode before calling zhpev 2307 index=1 2308 do ii=1,npwdiel 2309 do jj=1,ii 2310 symh(index )=symdielmat(1,jj,1,ii,1) 2311 symh(index+1)=symdielmat(2,jj,1,ii,1) 2312 index=index+2 2313 end do 2314 end do 2315 2316 ier=0 2317 call ZHPEV ('V','U',npwdiel,symh,eig_sym,symvec,npwdiel,zhpev1,& 2318 & zhpev2,ier) 2319 2320 if(prtvol>=10)then 2321 write(message, '(a,a,a,5es12.4)' )ch10,& 2322 & ' Five largest eigenvalues of the symmetrized dielectric matrix:',& 2323 & ch10,eig_sym(npwdiel:npwdiel-4:-1) 2324 call wrtout(ab_out,message,'COLL') 2325 end if 2326 2327 write(message, '(a,a)' )ch10,& 2328 & ' dieltcel : 15 largest eigenvalues of the symmetrized dielectric matrix' 2329 call wrtout(std_out,message,'COLL') 2330 write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(npwdiel:npwdiel-4:-1) 2331 call wrtout(std_out,message,'COLL') 2332 write(message, '(a,5es12.5)' )' 6-10 :',eig_sym(npwdiel-5:npwdiel-9:-1) 2333 call wrtout(std_out,message,'COLL') 2334 write(message, '(a,5es12.5)' )' 11-15:',eig_sym(npwdiel-10:npwdiel-14:-1) 2335 call wrtout(std_out,message,'COLL') 2336 write(message, '(a,a)' )ch10,& 2337 & ' dieltcel : 5 smallest eigenvalues of the symmetrized dielectric matrix' 2338 call wrtout(std_out,message,'COLL') 2339 write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(1:5) 2340 call wrtout(std_out,message,'COLL') 2341 2342 !Invert the hermitian dielectric matrix, 2343 work(:,:,:,:,:)=0.0_dp 2344 do ieig=1,npwdiel 2345 eiginv=1.0_dp/eig_sym(ieig) 2346 do ipw2=1,npwdiel 2347 ! do ipw1=ipw2,npwdiel 2348 do ipw1=1,npwdiel 2349 work(1,ipw1,1,ipw2,1)=work(1,ipw1,1,ipw2,1)+& 2350 & (symvec(1,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)+ & 2351 & symvec(2,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv 2352 work(2,ipw1,1,ipw2,1)=work(2,ipw1,1,ipw2,1)+& 2353 & (symvec(2,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)- & 2354 & symvec(1,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv 2355 end do 2356 end do 2357 end do 2358 !if(npwdiel>1)then 2359 !do ipw2=2,npwdiel 2360 !do ipw1=1,ipw2-1 2361 !work(1,ipw1,1,ipw2,1)= work(1,ipw2,1,ipw1,1) 2362 !work(2,ipw1,1,ipw2,1)=-work(2,ipw2,1,ipw1,1) 2363 !end do 2364 !end do 2365 !end if 2366 2367 ABI_DEALLOCATE(eig_sym) 2368 ABI_DEALLOCATE(symh) 2369 ABI_DEALLOCATE(symvec) 2370 2371 !DEBUG 2372 !Checks whether the inverse of the symmetric dielectric matrix 2373 !has been correctly generated 2374 !do ipw1=1,npwdiel 2375 !do ipw2=1,npwdiel 2376 !elementr=0.0_dp 2377 !elementi=0.0_dp 2378 !do ipw3=1,npwdiel 2379 !elementr=elementr+work(1,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1)& 2380 !& -work(2,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1) 2381 !elementi=elementi+work(1,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1)& 2382 !& +work(2,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1) 2383 !end do 2384 !if(elementr**2+elementi**2 > 1.0d-12)then 2385 !if( ipw1 /= ipw2 .or. & 2386 !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 2387 !write(std_out,*)' dieltcel : the inversion procedure is not correct ' 2388 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 2389 !write(std_out,*)' elementr,elementi=',elementr,elementi 2390 !stop 2391 !end if 2392 !end if 2393 !end do 2394 !end do 2395 !write(std_out,*)'dieltcel : matrix has been inverted successfully ' 2396 !ENDDEBUG 2397 2398 ABI_DEALLOCATE(symdielmat) 2399 2400 !Then get the inverse of the asymmetric 2401 !dielectric matrix, as required for the preconditioning. 2402 !Premultiplication by square root of minus susceptibility matrix 2403 do ipw2=1,npwdiel 2404 do ipw1=1,npwdiel 2405 ar=0._dp ; ai=0._dp 2406 do ii=1,npwdiel 2407 ar=ar+invsqrsus(1,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) & 2408 & -invsqrsus(2,ipw1,1,ii,1)*work(2,ii,1,ipw2,1) 2409 ai=ai+invsqrsus(2,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) & 2410 & +invsqrsus(1,ipw1,1,ii,1)*work(2,ii,1,ipw2,1) 2411 end do 2412 work2(1,ipw1,1,ipw2,1)=ar 2413 work2(2,ipw1,1,ipw2,1)=ai 2414 end do 2415 end do 2416 !Postmultiplication by square root of minus susceptibility matrix 2417 do ipw2=1,npwdiel 2418 do ipw1=1,npwdiel 2419 ar=0._dp ; ai=0._dp 2420 do ii=1,npwdiel 2421 ar=ar+work2(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2422 & -work2(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2423 ai=ai+work2(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 2424 & +work2(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 2425 end do 2426 dielinv(1,ipw1,1,ipw2,1)=ar 2427 dielinv(2,ipw1,1,ipw2,1)=ai 2428 end do 2429 end do 2430 2431 ABI_DEALLOCATE(invsqrsus) 2432 ABI_DEALLOCATE(sqrsus) 2433 ABI_DEALLOCATE(work) 2434 ABI_DEALLOCATE(work2) 2435 ABI_DEALLOCATE(zhpev1) 2436 ABI_DEALLOCATE(zhpev2) 2437 2438 call timab(96,2,tsec) 2439 2440 end subroutine dieltcel
ABINIT/linmin [ Functions ]
NAME
linmin
FUNCTION
minimizes a function along a gradient line: first bracket the minimum then perform the minimization
COPYRIGHT
Copyright (C) 1998-2018 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
PARENTS
cgpr
CHILDREN
bracketing
SOURCE
3120 subroutine linmin(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,v,grad,fmin) 3121 3122 3123 !This section has been created automatically by the script Abilint (TD). 3124 !Do not modify the following lines by hand. 3125 #undef ABI_FUNC 3126 #define ABI_FUNC 'linmin' 3127 !End of the abilint section 3128 3129 implicit none 3130 3131 !Arguments ------------------------------------ 3132 include "dummy_functions.inc" 3133 !scalars 3134 integer,intent(in) :: nv1,nv2 3135 real(dp),intent(out) :: fmin 3136 !arrays 3137 real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2) 3138 3139 !Local variables------------------------------- 3140 !scalars 3141 real(dp),parameter :: maglimit=10000.0_dp,tol=tol8*tol8*tol8 3142 real(dp) :: a,b,fa,fb,fx,x,xmin 3143 !no_abirules 3144 3145 !************************************************************************ 3146 a=zero 3147 x=ninth*real(1e-4,dp) 3148 call bracketing (nv1,nv2,dp_dum_v2dp,v,grad,a,x,b,fa,fx,fb) 3149 !DEBUG 3150 !write(std_out,*) 'linmin (01cg) : linmin after bracketing' 3151 !write(std_out,*) 'linmin (01cg) : point',a,x,b,'value',fa,fx,fb 3152 !ENDDEBUG 3153 fmin =brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,6,v,grad,a,x,b,tol,xmin) 3154 3155 end subroutine linmin
ABINIT/m_prcref [ Modules ]
NAME
m_prcref
FUNCTION
Routines to precondition residual potential (or density) and forces.
COPYRIGHT
Copyright (C) 1998-2018 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 .
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 module m_prcref 28 29 use defs_basis 30 use defs_datatypes 31 use defs_abitypes 32 use defs_wvltypes 33 use m_errors 34 use m_abicore 35 use m_xmpi 36 use m_xcdata 37 38 use m_frskerker1 39 use m_frskerker2 40 use mod_prc_memory 41 42 use m_time, only : timab 43 use m_numeric_tools, only : dotproduct 44 use m_geometry, only : xcart2xred, metric 45 use m_cgtools, only : dotprod_vn, mean_fftr 46 use m_mpinfo, only : ptabs_fourdp, destroy_mpi_enreg, initmpi_seq 47 use m_pawtab, only : pawtab_type 48 use m_pawrhoij, only : pawrhoij_type 49 use m_fftcore, only : kgindex 50 use m_fft, only : zerosym, indirect_parallel_fourier, fourdp 51 use m_kg, only : getph 52 use m_dtset, only : testsusmat 53 use m_spacepar, only : hartre, laplacian 54 use m_distribfft, only : init_distribfft_seq 55 use m_forces, only : fresid 56 use m_atm2fft, only : atm2fft 57 use m_rhotoxc, only : rhotoxc 58 use m_mklocl, only : mklocl 59 use m_mkcore, only : mkcore 60 61 62 implicit none 63 64 private
ABINIT/moddiel [ Functions ]
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
PARENTS
dfpt_newvtr,prcref,prcref_PMA
CHILDREN
fourdp,metric,ptabs_fourdp
SOURCE
1379 subroutine moddiel(cplex,dielar,mpi_enreg,nfft,ngfft,nspden,optreal,optres,paral_kgb,qphon,rprimd,vresid,vrespc) 1380 1381 1382 !This section has been created automatically by the script Abilint (TD). 1383 !Do not modify the following lines by hand. 1384 #undef ABI_FUNC 1385 #define ABI_FUNC 'moddiel' 1386 !End of the abilint section 1387 1388 implicit none 1389 1390 !Arguments------------------------------- 1391 !scalars 1392 integer,intent(in) :: cplex,nfft,nspden,optreal,optres,paral_kgb 1393 type(MPI_type),intent(in) :: mpi_enreg 1394 !arrays 1395 integer,intent(in) :: ngfft(18) 1396 real(dp),intent(in) :: dielar(7),qphon(3),rprimd(3,3) 1397 real(dp),intent(in) :: vresid(cplex*nfft,nspden) 1398 real(dp),intent(out) :: vrespc(cplex*nfft,nspden) 1399 1400 !Local variables------------------------------- 1401 !scalars 1402 integer,parameter :: im=2,re=1 1403 integer :: i1,i2,i23,i3,ifft,ig,ii,ii1,ing,ispden,me_fft,mg,n1,n2,n3,nproc_fft 1404 integer :: nspden_eff,qeq0 1405 logical :: magn_precon 1406 real(dp) :: dielng,diemac,diemac_inv,diemix,diemixmag,diemix_eff,factor,gqg2p3,gqgm12,gqgm13 1407 real(dp) :: gqgm23,gs,gs2,gs3,l2g2,length2,ucvol 1408 character(len=500) :: message 1409 !arrays 1410 integer :: id(3) 1411 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1412 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1413 real(dp) :: gmet(3,3),gprimd(3,3),potg0(4),rmet(3,3) 1414 real(dp),allocatable :: gq(:,:),work1(:,:),work2(:) 1415 1416 ! ************************************************************************* 1417 1418 !Check that cplex has an allowed value 1419 if(cplex/=1 .and. cplex/=2)then 1420 write(message,'(a,i0,a,a)')& 1421 & ' From the calling routine, cplex=',cplex,ch10,& 1422 & ' but the only value allowed are 1 and 2.' 1423 MSG_BUG(message) 1424 end if 1425 1426 if(cplex==1.and.optreal==2)then 1427 MSG_BUG('When optreal=2, cplex must be 2.') 1428 end if 1429 1430 !This is to allow q=0 1431 qeq0=0 1432 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1 1433 1434 !If cplex=1 then qphon should be 0 0 0 1435 if (cplex==1.and. qeq0/=1) then 1436 write(message,'(a,3e12.4,a)' )' cplex=1 but qphon=',qphon,' qphon should be 0 0 0.' 1437 MSG_BUG(message) 1438 end if 1439 1440 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1441 me_fft=ngfft(11) 1442 nproc_fft=ngfft(10) 1443 1444 !Get the distrib associated with this fft_grid 1445 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1446 1447 !Compute different geometric tensor, as well as ucvol, from rprimd 1448 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1449 1450 dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ; diemixmag=dielar(7) 1451 1452 magn_precon=(diemixmag>=zero) ! Set to true if magnetization has to be preconditionned 1453 diemixmag=abs(diemixmag) 1454 1455 !DEBUG 1456 !write(std_out,*)' moddiel : diemac, diemix, diemixmag =',diemac,diemix,diemixmag 1457 !ENDDEBUG 1458 1459 if(abs(diemac-1.0_dp)<1.0d-6)then 1460 1461 ! Here, simple mixing is required, through macroscopic 1462 ! dielectric constant set to 1.0_dp . 1463 vrespc(:,1)=diemix*vresid(:,1) 1464 if (nspden/=1) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden) 1465 else 1466 1467 ! Magnetization is not preconditionned 1468 if (optres==1.and.nspden>1.and.(.not.magn_precon)) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden) 1469 1470 ! Here, model dielectric function (G-diagonal operator) 1471 1472 length2=(two_pi*dielng)**2 1473 diemac_inv=1.0_dp/diemac 1474 ABI_ALLOCATE(work1,(2,nfft)) 1475 if (optreal==1) then 1476 ABI_ALLOCATE(work2,(cplex*nfft)) 1477 end if 1478 1479 ! In order to speed the routine, precompute the components of g 1480 mg=maxval(ngfft) 1481 ABI_ALLOCATE(gq,(3,mg)) 1482 do ii=1,3 1483 id(ii)=ngfft(ii)/2+2 1484 do ing=1,ngfft(ii) 1485 ig=ing-(ing/id(ii))*ngfft(ii)-1 1486 gq(ii,ing)=ig+qphon(ii) 1487 end do 1488 end do 1489 1490 ! Do-loop on spins 1491 ! Note XG 010922 : I doubt the preconditioner is OK for the magnetization 1492 nspden_eff=nspden;if (optres==1.and.(.not.magn_precon)) nspden_eff=1 1493 do ispden=1,nspden_eff 1494 1495 diemix_eff=diemix;if (ispden>1) diemix_eff=diemixmag 1496 1497 ! Do fft from real space (work2) to G space (work1) 1498 if (optreal==1) then 1499 work2(:)=vresid(:,ispden) 1500 call fourdp(cplex,work1,work2,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1501 else 1502 ! work1(:,:)=reshape(vresid(:,ispden),(/2,nfft/)) 1503 ! Reshape function does not work with big arrays for some compilers 1504 do ifft=1,nfft 1505 work1(1,ifft)=vresid(2*ifft-1,ispden) 1506 work1(2,ifft)=vresid(2*ifft ,ispden) 1507 end do 1508 end if 1509 1510 ! Store G=0 value 1511 potg0(ispden)=work1(re,1) 1512 1513 ! Triple loop, for the three dimensions 1514 do i3=1,n3 1515 ! Precompute some products that do not depend on i2 and i1 1516 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 1517 gqgm23=gq(3,i3)*gmet(2,3)*2 1518 gqgm13=gq(3,i3)*gmet(1,3)*2 1519 do i2=1,n2 1520 if (fftn2_distrib(i2)==me_fft) then 1521 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 1522 gqgm12=gq(2,i2)*gmet(1,2)*2 1523 gqg2p3=gqgm13+gqgm12 1524 i23=n1*( ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1)) 1525 1526 ! Do the test that eliminates the Gamma point outside 1527 ! of the inner loop 1528 ii1=1 1529 if(i2 == 1 .and. i3 == 1 .and. qeq0==1)then 1530 ! if(i23==0 .and. qeq0==1)then: this changes with the number of fft procs... 1531 ! and seems to be wrong.Pls check 1532 ii1=2 1533 end if 1534 1535 ! Here, unlike in hartre.f, the G=0 term is not eliminated, albeit 1536 ! not changed. 1537 do i1=ii1,n1 1538 1539 ! One obtains the square of the norm of q+G (defined by i1,i2,i3) 1540 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 1541 ifft=i1+i23 1542 1543 l2g2=length2*gs 1544 ! The model dielectric function is now computed 1545 factor = (l2g2+diemac_inv)/(l2g2+1.0_dp) * diemix_eff 1546 work1(re,ifft)=work1(re,ifft)*factor 1547 work1(im,ifft)=work1(im,ifft)*factor 1548 1549 end do 1550 end if 1551 end do 1552 end do 1553 1554 ! Might get rid of the G=0 term 1555 ! if(qeq0==1)then 1556 ! work1(re,1)=0.0_dp 1557 ! work1(im,1)=0.0_dp 1558 ! end if 1559 1560 ! Fourier transform 1561 if (optreal==1) then 1562 call fourdp(cplex,work1,work2,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1563 vrespc(:,ispden)=work2(:) 1564 else 1565 ! vrespc(:,ispden)=reshape(work1(:,:),(/nfft*2/)) 1566 ! Reshape function does not work with big arrays for some compilers 1567 do ifft=1,nfft 1568 vrespc(2*ifft-1,ispden)=work1(1,ifft) 1569 vrespc(2*ifft ,ispden)=work1(2,ifft) 1570 end do 1571 end if 1572 1573 ! End of loop on spin polarizations 1574 end do 1575 1576 ABI_DEALLOCATE(gq) 1577 ABI_DEALLOCATE(work1) 1578 if (optreal==1) then 1579 ABI_DEALLOCATE(work2) 1580 end if 1581 1582 ! End condition diemac/=1.0 1583 end if 1584 1585 end subroutine moddiel
ABINIT/prcref [ Functions ]
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
PARENTS
newrho
CHILDREN
atm2fft,dielmt,dieltcel,fourdp,fresid,getph,hartre indirect_parallel_fourier,kgindex,mean_fftr,metric,mkcore,mklocl moddiel,prcrskerker1,prcrskerker2,rhotoxc,testsusmat,xcart2xred xcdata_init,xmpi_sum,zerosym
SOURCE
195 subroutine prcref(atindx,dielar,dielinv,& 196 & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,& 197 & istep,kg_diel,kxc,& 198 & mgfft,moved_atm_inside,mpi_enreg,my_natom,& 199 & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 200 & optreal,optres,pawrhoij,pawtab,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,& 201 & susmat,vhartr,vpsp,vresid,vrespc,vxc,wvl,wvl_den,xred) 202 203 204 !This section has been created automatically by the script Abilint (TD). 205 !Do not modify the following lines by hand. 206 #undef ABI_FUNC 207 #define ABI_FUNC 'prcref' 208 !End of the abilint section 209 210 implicit none 211 212 !Arguments------------------------------- 213 !scalars 214 integer,intent(in) :: dielstrt,istep,my_natom,mgfft,moved_atm_inside,n1xccc 215 integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat,optreal,optres 216 real(dp),intent(in) :: etotal,gsqcut 217 type(MPI_type),intent(in) :: mpi_enreg 218 type(dataset_type),intent(in) :: dtset 219 type(pseudopotential_type),intent(in) :: psps 220 type(wvl_internal_type), intent(in) :: wvl 221 type(wvl_denspot_type), intent(inout) :: wvl_den 222 !arrays 223 integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft)) 224 integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18) 225 real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),rhog(2,nfft) 226 real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3) 227 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 228 real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden) 229 real(dp),intent(in) :: vxc(nfft,dtset%nspden) 230 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 231 real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc) 232 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft) 233 real(dp),intent(inout) :: xred(3,dtset%natom) 234 real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix) 235 real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden) 236 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 237 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 238 239 !Local variables------------------------------- 240 !scalars 241 integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1 242 integer :: ipw2,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm 243 integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method 244 real(dp) :: ai,ar,diemix,diemixmag,eei,enxc 245 real(dp) :: mixfac 246 real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg 247 logical :: computediel 248 logical :: non_magnetic_xc 249 character(len=500) :: message 250 type(xcdata_type) :: xcdata 251 !arrays 252 integer :: qprtrb(3) 253 integer,allocatable :: indpw_prc(:) 254 real(dp) :: dummy6(6),dummy7(6),gprimd(3,3),qphon(3),rmet(3,3),strsxc(6) 255 real(dp) :: vmean(dtset%nspden),vprtrb(2) 256 real(dp),allocatable :: dummy(:),dummy1(:),dummy2(:),dummy3(:),dummy4(:),dummy5(:),dummy8(:),dummy9(:) 257 real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:) 258 real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:) 259 real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:),rhor_new(:,:) 260 real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:) 261 real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:) 262 real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:) 263 logical,allocatable :: mask(:) 264 265 ! ************************************************************************* 266 267 !Compute different geometric tensor, as well as ucvol, from rprimd 268 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 269 270 ! Initialise non_magnetic_xc for rhohxc 271 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 272 273 !1) Eventually take care of the forces 274 275 if(moved_atm_inside==1)then 276 ABI_ALLOCATE(fcart_pc,(3,dtset%natom)) 277 278 if(dtset%iprcfc==0)then 279 fcart_pc(:,:)=fcart(:,:) 280 else 281 fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:) 282 end if 283 284 ! Compute preconditioned delta xred from preconditioned fcart and rprimd 285 call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc) 286 287 ABI_DEALLOCATE(fcart_pc) 288 end if 289 290 !####################################################################### 291 292 !2) Take care of the potential residual 293 294 !Compute the residuals corresponding to the solution 295 !of an approximate realspace dielectric function according 296 !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]] 297 if(dtset%iprcel>=71.and.dtset%iprcel<=79) then 298 if (nfft==nfftprc) then 299 if (dtset%iprcel<=78) then 300 call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, & 301 & gprimd,vresid,vrespc,rhor(:,1)) 302 else 303 call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, & 304 & vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol) 305 end if 306 else 307 ! If preconditionning has to be done on a coarse grid, 308 ! has to transfer several arrays 309 ABI_ALLOCATE(work1,(nfftprc,dtset%nspden)) 310 ABI_ALLOCATE(work3,(nfftprc,dtset%nspden)) 311 ABI_ALLOCATE(work,(2*nfftprc)) 312 do ispden=1,dtset%nspden 313 work(:)=vresid(:,ispden) 314 call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 315 end do 316 ABI_DEALLOCATE(work) 317 if (dtset%iprcel<=78) then 318 ABI_ALLOCATE(rhog_wk,(2,nfftprc)) 319 rhog_wk(:,:)=zero 320 if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then 321 nfftot=PRODUCT(ngfft(1:3)) 322 call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,& 323 & ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot) 324 else 325 do ii=1,nfft 326 if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii) 327 end do 328 end if 329 call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),& 330 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 331 ABI_ALLOCATE(work,(nfftprc)) 332 call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 333 call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, & 334 & gprimd,work1,work3,work) 335 ABI_DEALLOCATE(work) 336 else 337 call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, & 338 & work1,work3,dtset%natom,xred,mpi_enreg,ucvol) 339 end if 340 do ispden=1,dtset%nspden 341 call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 342 end do 343 ABI_DEALLOCATE(work1) 344 ABI_DEALLOCATE(work3) 345 end if 346 347 else 348 349 if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then 350 cplex=optreal 351 qphon(:)=zero 352 ! Simple scalar multiplication, or model dielectric function 353 call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,dtset%paral_kgb,qphon,rprimd,vresid,vrespc) 354 355 ! Use the inverse dielectric matrix in a small G sphere 356 else if( (istep>=dielstrt .and. dtset%iprcel>=21) .or. modulo(dtset%iprcel,100)>=41 )then 357 358 ! With dielop=1, the matrices will be computed when istep=dielstrt 359 ! With dielop=2, the matrices will be computed when istep=dielstrt and 1 360 dielop=1 361 if(modulo(dtset%iprcel,100)>=41)dielop=2 362 call testsusmat(computediel,dielop,dielstrt,dtset,istep) !test if the matrix is to be computed 363 if(computediel) then 364 ! Compute the inverse dielectric matrix from the susceptibility matrix 365 ! There are two routines for the RPA matrix, while for the electronic 366 ! dielectric matrix, only dieltcel will do the work 367 if(modulo(dtset%iprcel,100)<=49)then 368 call dielmt(dielinv,gmet,kg_diel,& 369 & npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat) 370 else 371 option=1 372 if(modulo(dtset%iprcel,100)>=61)option=2 373 call dieltcel(dielinv,gmet,kg_diel,kxc,& 374 & nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%paral_kgb,dtset%prtvol,susmat) 375 end if 376 end if 377 378 ABI_ALLOCATE(work1,(2,nfftprc)) 379 ABI_ALLOCATE(work2,(optreal*nfftprc)) 380 381 ! Presently, one uses the inverse of the RPA dielectric matrix, 382 ! for which spin must be averaged. 383 384 ! Do fft from real space (work2) to G space (work1) 385 if (optreal==1) then 386 work2(:)=vresid(:,1) 387 ! Must average over spins in the case of a potential residual 388 if(dtset%nspden/=1.and.optres==0)work2(:)=(work2(:)+vresid(:,2))*half 389 call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 390 else 391 work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/)) 392 if(dtset%nspden/=1.and.optres==0)work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half 393 end if 394 395 ! Multiply by restricted inverse of dielectric matrix. 396 ! Must first copy relevant elements of work1 to a npwdiel-dimensioned array, 397 ! then zero work1, operate with the dielinv matrix, and store in work1. 398 399 ABI_ALLOCATE(vres_diel,(2,npwdiel)) 400 ABI_ALLOCATE(indpw_prc,(npwdiel)) 401 ABI_ALLOCATE(mask,(npwdiel)) 402 mask(:)=.true. 403 call kgindex(indpw_prc,kg_diel,mask,mpi_enreg,ngfftprc,npwdiel) 404 405 do ipw1=1,npwdiel 406 if(mask(ipw1)) then 407 vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1)) 408 vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1)) 409 end if 410 end do 411 work1(:,:)=zero 412 do ipw1=1,npwdiel 413 ar=zero ; ai=zero 414 415 ! Use inverse of dielectric matrix (potential mixing) 416 if (optres==0) then 417 do ipw2=1,npwdiel 418 if(mask(ipw2))then 419 ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 420 & -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 421 ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 422 & +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 423 end if 424 end do 425 else 426 ! Use symetric of inverse of dielectric matrix (density mixing) 427 do ipw2=1,npwdiel 428 if(mask(ipw2))then 429 ar=ar+dielinv(1,ipw2,1,ipw1,1)*vres_diel(1,ipw2) & 430 & +dielinv(2,ipw2,1,ipw1,1)*vres_diel(2,ipw2) 431 ai=ai-dielinv(2,ipw2,1,ipw1,1)*vres_diel(1,ipw2) & 432 & +dielinv(1,ipw2,1,ipw1,1)*vres_diel(2,ipw2) 433 end if 434 end do 435 end if 436 ! Must be careful not to count the diagonal 1 twice : it is added later, 437 ! so must be subtracted now. 438 call xmpi_sum(ar,mpi_enreg%comm_fft,ier) 439 call xmpi_sum(ai,mpi_enreg%comm_fft,ier) 440 if(mask(ipw1)) then 441 work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1) 442 work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1) 443 end if !mask(ipw1) 444 end do ! ipw1 445 ABI_DEALLOCATE(vres_diel) 446 ABI_DEALLOCATE(indpw_prc) 447 ABI_DEALLOCATE(mask) 448 ! Fourier transform 449 if (optreal==1) then 450 call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 451 else 452 work2(:)=reshape(work1(:,:),(/nfftprc*2/)) 453 end if 454 455 ! Add to get the preconditioned vresid, must be careful about spins. 456 if(dtset%iprcel>=30)then 457 diemix=dielar(4);diemixmag=abs(dielar(7)) 458 vrespc(:,1)=diemix*(vresid(:,1)+work2(:)) 459 if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:)) 460 if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=diemixmag*vresid(:,3:4) 461 if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=diemixmag*vresid(:,2:dtset%nspden) 462 else 463 vrespc(:,1)=vresid(:,1)+work2(:) 464 if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=vresid(:,2)+work2(:) 465 if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=vresid(:,3:4) 466 if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=vresid(:,2:dtset%nspden) 467 end if 468 469 ABI_DEALLOCATE(work1) 470 ABI_DEALLOCATE(work2) 471 472 ! Other choice ? 473 else 474 write(message, '(a,i3,a,a,a,a)' )& 475 & 'From the calling routine, iprcel=',dtset%iprcel,ch10,& 476 & 'The only allowed values are 0 or larger than 20.',ch10,& 477 & 'Action: correct your input file.' 478 MSG_ERROR(message) 479 end if 480 end if 481 !####################################################################### 482 483 !3) PAW only : precondition the rhoij quantities (augmentation 484 !occupancies) residuals. Use a simple preconditionning 485 !with the same mixing factor as the model dielectric function. 486 487 if (psps%usepaw==1.and.my_natom>0) then 488 if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then 489 mixfac=one;mixfacmag=one 490 else 491 mixfac=dielar(4);mixfacmag=abs(dielar(7)) 492 end if 493 if (pawrhoij(1)%cplex==1) then 494 index=0 495 do iatom=1,my_natom 496 do ispden=1,pawrhoij(iatom)%nspden 497 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 498 do kmix=1,pawrhoij(iatom)%lmnmix_sz 499 index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix) 500 rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 501 end do 502 end do 503 end do 504 else 505 index=-1 506 do iatom=1,my_natom 507 do ispden=1,pawrhoij(iatom)%nspden 508 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 509 do kmix=1,pawrhoij(iatom)%lmnmix_sz 510 index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 511 rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 512 end do 513 end do 514 end do 515 end if 516 end if 517 518 !####################################################################### 519 520 !4) Take care of the change of atomic positions 521 !Note : this part is very demanding on memory... 522 !however, since this algorithm is still in development, 523 !it was NOT included in the estimation provided by memory.f 524 if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then 525 526 ! Not yet compatible with resid given in reciprocal space 527 if (optreal/=1) then 528 write(message, '(5a)' )& 529 & 'From the calling routine, densfor_pred=3',ch10,& 530 & 'You cannot use residuals in reciprocal space.',ch10,& 531 & 'Action: correct your input file.' 532 MSG_ERROR(message) 533 end if 534 ! Not compatible with non-collinear magnetism 535 if(dtset%nspden==4)then 536 MSG_ERROR('densfor_pred=3 does not work for nspden=4 !') 537 end if 538 539 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 540 nfftot=PRODUCT(ngfft(1:3)) 541 542 if (optres==0) then ! Array vresid contains a potential residual 543 ! ----------------------------------------------------------------- 544 545 ! First subtract the current local, hartree and exchange correlation potentials 546 do ispden=1,min(dtset%nspden,2) 547 vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden) 548 end do 549 if (dtset%nspden==4) then 550 do ispden=3,4 551 vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden) 552 end do 553 end if 554 555 ! Compute the modified density, in rhor_wk 556 option=2 557 ABI_ALLOCATE(gresid,(3,dtset%natom)) 558 ABI_ALLOCATE(grxc,(3,dtset%natom)) 559 ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden)) 560 ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden)) 561 ABI_ALLOCATE(xred_wk,(3,dtset%natom)) 562 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 563 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 564 & ntypat,option,pawtab,rhor,rprimd,& 565 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 566 567 ! Compute up+down rhog_wk(G) by fft 568 ABI_ALLOCATE(work,(nfft)) 569 ABI_ALLOCATE(rhog_wk,(2,nfft)) 570 work(:)=rhor_wk(:,1) 571 call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 572 ABI_DEALLOCATE(work) 573 574 ! Compute structure factor phases for new atomic pos: 575 call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk) 576 577 ! Compute local ionic pseudopotential vpsp: 578 ! and core electron density xccc3d, if needed. 579 n3xccc=0;if (n1xccc/=0) n3xccc=nfft 580 ABI_ALLOCATE(xccc3d,(n3xccc)) 581 ABI_ALLOCATE(vpsp_wk,(nfft)) 582 vprtrb(1:2)=zero 583 584 ! Determine by which method the local ionic potential and/or 585 ! the pseudo core charge density contributions have to be computed 586 ! Local ionic potential: 587 ! Method 1: PAW 588 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 589 vloc_method=1;if (psps%usepaw==0) vloc_method=2 590 if (dtset%icoulomb>0) vloc_method=2 591 if (psps%usewvl==1) vloc_method=2 592 ! Pseudo core charge density: 593 ! Method 1: PAW, nc_xccc_gspace 594 ! Method 2: Norm-conserving PP, wavelets 595 coredens_method=1;if (psps%usepaw==0) coredens_method=2 596 if (psps%nc_xccc_gspace==1) coredens_method=1 597 if (psps%nc_xccc_gspace==0) coredens_method=2 598 if (psps%usewvl==1) coredens_method=2 599 600 ! Local ionic potential and/or pseudo core charge by method 1 601 if (vloc_method==1.or.coredens_method==1) then 602 optv=0;if (vloc_method==1) optv=1 603 optn=0;if (coredens_method==1) optn=n3xccc/nfft 604 optatm=1;optdyfr=0;optgr=0;optstr=0;optn2=1;opteltfr=0 605 ! Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused... 606 call atm2fft(atindx,xccc3d,vpsp,dummy,dummy2,dummy9,dummy1,gmet,gprimd,dummy3,dummy4,gsqcut,& 607 & mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,& 608 & optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dummy5,dummy6,dummy7,& 609 & ucvol,psps%usepaw,dummy8,dummy8,dummy8,vprtrb,psps%vlspl,& 610 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 611 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 612 end if 613 614 ! Local ionic potential by method 2 615 if (vloc_method==2) then 616 option=1 617 ABI_ALLOCATE(dyfrlo_indx,(3,3,dtset%natom)) 618 ABI_ALLOCATE(grtn_indx,(3,dtset%natom)) 619 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,& 620 & mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,& 621 & ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,& 622 & ucvol,vprtrb,vpsp_wk,wvl,wvl_den,xred) 623 ABI_DEALLOCATE(dyfrlo_indx) 624 ABI_DEALLOCATE(grtn_indx) 625 end if 626 627 ! Pseudo core electron density by method 2 628 if (coredens_method==2.and.n1xccc/=0) then 629 option=1 630 ABI_ALLOCATE(dyfrx2,(3,3,dtset%natom)) 631 ABI_ALLOCATE(grxc_indx,(3,dtset%natom)) 632 call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 633 & n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,& 634 & psps%xccc1d,xccc3d,xred_wk) 635 ABI_DEALLOCATE(dyfrx2) 636 ABI_DEALLOCATE(grxc_indx) 637 end if 638 639 ! Compute Hartree+xc potentials 640 ABI_ALLOCATE(vxc_wk,(nfft,dtset%nspden)) 641 ABI_ALLOCATE(vhartr_wk,(nfft)) 642 option=1 643 644 call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog_wk,rprimd,vhartr_wk) 645 646 ! Prepare the call to rhotoxc 647 call xcdata_init(xcdata,dtset=dtset) 648 nk3xc=1 649 call rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft,& 650 & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor_wk,rprimd,strsxc,1,& 651 & vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk) 652 ABI_DEALLOCATE(xccc3d) 653 654 ! Sum all contributions 655 do ispden=1,min(dtset%nspden,2) 656 do ifft=1,nfft 657 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden) 658 end do 659 end do 660 if (dtset%nspden==4) then 661 do ispden=3,4 662 do ifft=1,nfft 663 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden) 664 end do 665 end do 666 end if 667 call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft) 668 if(dtset%nspden==2) then 669 vmean(1)=half*(vmean(1)+vmean(2)) 670 vmean(2)=vmean(1) 671 end if 672 do ispden=1,dtset%nspden 673 vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden) 674 end do 675 ABI_DEALLOCATE(gresid) 676 ABI_DEALLOCATE(grxc) 677 ABI_DEALLOCATE(rhog_wk) 678 ABI_DEALLOCATE(rhor_wk) 679 ABI_DEALLOCATE(rhor_wk0) 680 ABI_DEALLOCATE(xred_wk) 681 ABI_DEALLOCATE(vhartr_wk) 682 ABI_DEALLOCATE(vpsp_wk) 683 ABI_DEALLOCATE(vxc_wk) 684 685 else ! Array vresid contains a density residual 686 ! ----------------------------------------------------------------- 687 688 ! Only have to compute the modified preconditionned density residual 689 option=2 690 ABI_ALLOCATE(gresid,(3,dtset%natom)) 691 ABI_ALLOCATE(grxc,(3,dtset%natom)) 692 ABI_ALLOCATE(rhor_new,(nfft,dtset%nspden)) 693 ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden)) 694 ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden)) 695 ABI_ALLOCATE(xred_wk,(3,dtset%natom)) 696 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 697 rhor_new(:,1)=rhor(:,1)+vrespc(:,1) 698 if (dtset%nspden==2) then 699 rhor_new(:,1)=rhor_new(:,1)+vrespc(:,2) 700 rhor_new(:,2)=rhor(:,2)+vrespc(:,1) 701 end if 702 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 703 & ntypat,option,pawtab,rhor,rprimd,& 704 & ucvol,rhor_wk0,xred_wk,xred,psps%znuclpsp) 705 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 706 & ntypat,option,pawtab,rhor_new,rprimd,& 707 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 708 vrespc(:,1)=rhor_wk(:,dtset%nspden)-rhor_wk0(:,dtset%nspden) 709 if (dtset%nspden==2) vrespc(:,2)=rhor_wk(:,1)-rhor_wk0(:,1)-vrespc(:,1) 710 ABI_DEALLOCATE(gresid) 711 ABI_DEALLOCATE(grxc) 712 ABI_DEALLOCATE(rhor_new) 713 ABI_DEALLOCATE(rhor_wk) 714 ABI_DEALLOCATE(rhor_wk0) 715 ABI_DEALLOCATE(xred_wk) 716 end if 717 718 end if 719 720 end subroutine prcref
ABINIT/prcref_PMA [ Functions ]
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
PARENTS
newvtr
CHILDREN
atm2fft,dielmt,dieltcel,fourdp,fresid,getph,hartre indirect_parallel_fourier,kgindex,mean_fftr,metric,mkcore,mklocl moddiel,prcrskerker1,prcrskerker2,rhotoxc,testsusmat,xcart2xred xcdata_init,xmpi_sum,zerosym
SOURCE
844 subroutine prcref_PMA(atindx,dielar,dielinv,& 845 & dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,& 846 & istep,kg_diel,kxc,& 847 & mgfft,moved_atm_inside,mpi_enreg,my_natom,& 848 & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 849 & optreal,optres,pawrhoij,ph1d,psps,rhog, rhoijrespc,rhor,rprimd,& 850 & susmat,vhartr,vpsp,vresid,vrespc,vxc,xred,& 851 & etotal,pawtab,wvl) 852 853 854 !This section has been created automatically by the script Abilint (TD). 855 !Do not modify the following lines by hand. 856 #undef ABI_FUNC 857 #define ABI_FUNC 'prcref_PMA' 858 !End of the abilint section 859 860 implicit none 861 862 !Arguments------------------------------- 863 !variables used for tfvw 864 !scalars 865 integer,intent(in) :: dielstrt,istep,mgfft,moved_atm_inside,my_natom,n1xccc 866 integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat 867 integer,intent(in) :: optreal,optres 868 real(dp),intent(in) :: etotal,gsqcut 869 type(MPI_type),intent(in) :: mpi_enreg 870 type(dataset_type),intent(in) :: dtset 871 type(pseudopotential_type),intent(in) :: psps 872 type(wvl_data), intent(inout) :: wvl 873 !arrays 874 integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft)) 875 integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18) 876 real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom) 877 real(dp),intent(in) :: rhog(2,nfft) 878 real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3) 879 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 880 real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden) 881 real(dp),intent(in) :: vxc(nfft,dtset%nspden) 882 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 883 real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc) 884 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft) 885 real(dp),intent(inout) :: xred(3,dtset%natom) 886 real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix) 887 real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden) 888 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 889 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 890 891 !Local variables------------------------------- 892 !scalars 893 integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1 894 integer :: ipw2,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm 895 integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method 896 real(dp) :: ai,ar,diemix,diemixmag,eei,enxc 897 real(dp) :: mixfac 898 real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg 899 logical :: computediel 900 logical :: non_magnetic_xc 901 character(len=500) :: message 902 type(xcdata_type) :: xcdata 903 !arrays 904 integer :: qprtrb(3) 905 integer,allocatable :: indpw_prc(:) 906 real(dp) :: dummy6(6),gprimd(3,3),qphon(3),rmet(3,3),strsxc(6) 907 real(dp) :: vmean(dtset%nspden),vprtrb(2) 908 real(dp),allocatable :: dummy_in(:) 909 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0) 910 real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:) 911 real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:) 912 real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:) 913 real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:) 914 real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:) 915 real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:) 916 logical,allocatable :: mask(:) 917 918 ! ************************************************************************* 919 920 if(optres==1)then 921 MSG_ERROR('density mixing (optres=1) not admitted!') 922 end if 923 924 !Compute different geometric tensor, as well as ucvol, from rprimd 925 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 926 927 ! Initialise non_magnetic_xc for rhohxc 928 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 929 930 !1) Eventually take care of the forces 931 932 if(moved_atm_inside==1)then 933 ABI_ALLOCATE(fcart_pc,(3,dtset%natom)) 934 935 if(dtset%iprcfc==0)then 936 fcart_pc(:,:)=fcart(:,:) 937 else 938 fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:) 939 end if 940 941 ! Compute preconditioned delta xred from preconditioned fcart and rprimd 942 call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc) 943 944 ABI_DEALLOCATE(fcart_pc) 945 end if 946 947 !####################################################################### 948 949 !2) Take care of the potential residual 950 951 !Compute the residuals corresponding to the solution 952 !of an approximate realspace dielectric function according 953 !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]] 954 if(dtset%iprcel>=71.and.dtset%iprcel<=79) then 955 if (nfft==nfftprc) then 956 if (dtset%iprcel<=78) then 957 call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, & 958 & gprimd,vresid,vrespc,rhor(:,1)) 959 else 960 call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, & 961 & vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol) 962 end if 963 else 964 ! If preconditionning has to be done on a coarse grid, 965 ! has to transfer several arrays 966 ABI_ALLOCATE(work1,(nfftprc,dtset%nspden)) 967 ABI_ALLOCATE(work3,(nfftprc,dtset%nspden)) 968 ABI_ALLOCATE(work,(2*nfftprc)) 969 do ispden=1,dtset%nspden 970 work(:)=vresid(:,ispden) 971 call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 972 end do 973 ABI_DEALLOCATE(work) 974 if (dtset%iprcel<=78) then 975 ABI_ALLOCATE(rhog_wk,(2,nfftprc)) 976 rhog_wk(:,:)=zero 977 if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then 978 nfftot=PRODUCT(ngfft(1:3)) 979 call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,& 980 & ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot) 981 else 982 do ii=1,nfft 983 if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii) 984 end do 985 end if 986 call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),& 987 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 988 ABI_ALLOCATE(work,(nfftprc)) 989 call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 990 call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, & 991 & gprimd,work1,work3,work) 992 ABI_DEALLOCATE(work) 993 else 994 call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, & 995 & work1,work3,dtset%natom,xred,mpi_enreg,ucvol) 996 end if 997 do ispden=1,dtset%nspden 998 call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 999 end do 1000 ABI_DEALLOCATE(work1) 1001 ABI_DEALLOCATE(work3) 1002 end if 1003 1004 else 1005 1006 if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then 1007 cplex=optreal 1008 qphon(:)=zero 1009 ! Simple scalar multiplication, or model dielectric function 1010 call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,dtset%paral_kgb,qphon,rprimd,vresid,vrespc) 1011 1012 ! Use the inverse dielectric matrix in a small G sphere 1013 else if( (istep>=dielstrt .and. dtset%iprcel>=21) .or. modulo(dtset%iprcel,100)>=41 )then 1014 1015 ! Wnith dielop=1, the matrices will be computed when istep=dielstrt 1016 ! With dielop=2, the matrices will be computed when istep=dielstrt and 1 1017 dielop=1 1018 if(modulo(dtset%iprcel,100)>=41)dielop=2 1019 call testsusmat(computediel,dielop,dielstrt,dtset,istep) !test if the matrix is to be computed 1020 if(computediel) then 1021 ! Compute the inverse dielectric matrix from the susceptibility matrix 1022 ! There are two routines for the RPA matrix, while for the electronic 1023 ! dielectric matrix, only dieltcel will do the work 1024 if(modulo(dtset%iprcel,100)<=49)then 1025 call dielmt(dielinv,gmet,kg_diel,& 1026 & npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat) 1027 else 1028 option=1 1029 if(modulo(dtset%iprcel,100)>=61)option=2 1030 call dieltcel(dielinv,gmet,kg_diel,kxc,& 1031 & nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%paral_kgb,dtset%prtvol,susmat) 1032 end if 1033 end if 1034 1035 ABI_ALLOCATE(work1,(2,nfftprc)) 1036 ABI_ALLOCATE(work2,(optreal*nfftprc)) 1037 1038 ! Presently, one uses the inverse of the RPA dielectric matrix, 1039 ! for which spin must be averaged. 1040 1041 ! Do fft from real space (work2) to G space (work1) 1042 if (optreal==1) then 1043 work2(:)=vresid(:,1) 1044 ! Must average over spins if needed. 1045 if(dtset%nspden/=1)work2(:)=(work2(:)+vresid(:,2))*half 1046 call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 1047 else 1048 work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/)) 1049 if (dtset%nspden/=1) work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half 1050 end if 1051 1052 ! Multiply by restricted inverse of dielectric matrix. 1053 ! Must first copy relevant elements of work1 to a npwdiel-dimensioned array, 1054 ! then zero work1, operate with the dielinv matrix, and store in work1. 1055 1056 ABI_ALLOCATE(vres_diel,(2,npwdiel)) 1057 ABI_ALLOCATE(indpw_prc,(npwdiel)) 1058 ABI_ALLOCATE(mask,(npwdiel)) 1059 mask(:)=.true. 1060 call kgindex(indpw_prc,kg_diel,mask,mpi_enreg,ngfftprc,npwdiel) 1061 do ipw1=1,npwdiel 1062 if(mask(ipw1)) then 1063 vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1)) 1064 vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1)) 1065 end if 1066 end do 1067 1068 work1(:,:)=zero 1069 do ipw1=1,npwdiel 1070 ar=zero ; ai=zero 1071 1072 ! Use inverse of dielectric matrix (potential mixing) 1073 do ipw2=1,npwdiel 1074 if(mask(ipw2))then 1075 ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 1076 & -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 1077 ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 1078 & +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 1079 end if 1080 end do 1081 ! Must be careful not to count the diagonal 1 twice : it is added later, 1082 ! so must be subtracted now. 1083 call xmpi_sum(ar,mpi_enreg%comm_fft,ier) 1084 call xmpi_sum(ai,mpi_enreg%comm_fft,ier) 1085 if(mask(ipw1)) then 1086 work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1) 1087 work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1) 1088 end if !mask(ipw1) 1089 end do ! ipw1 1090 ABI_DEALLOCATE(vres_diel) 1091 ABI_DEALLOCATE(indpw_prc) 1092 ABI_DEALLOCATE(mask) 1093 1094 ! Fourier transform 1095 if (optreal==1) then 1096 call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 1097 else 1098 work2(:)=reshape(work1(:,:),(/nfftprc*2/)) 1099 end if 1100 1101 ! Add to get the preconditioned vresid, must be careful about spins. 1102 if(dtset%iprcel>=30)then 1103 diemix=dielar(4);diemixmag=abs(dielar(7)) 1104 vrespc(:,1)=diemix*(vresid(:,1)+work2(:)) 1105 if(dtset%nspden/=1)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:)) 1106 if(dtset%nspden==4)vrespc(:,3:4)=diemixmag*vresid(:,3:4) 1107 else 1108 vrespc(:,1)=vresid(:,1)+work2(:) 1109 if(dtset%nspden/=1)vrespc(:,2)=vresid(:,2)+work2(:) 1110 if(dtset%nspden==4)vrespc(:,3:4)=vresid(:,3:4) 1111 end if 1112 1113 ABI_DEALLOCATE(work1) 1114 ABI_DEALLOCATE(work2) 1115 1116 ! Other choice ? 1117 else 1118 write(message, '(a,i0,a,a,a,a)' )& 1119 & 'From the calling routine, iprcel= ',dtset%iprcel,ch10,& 1120 & 'The only allowed values are 0 or larger than 20.',ch10,& 1121 & 'Action: correct your input file.' 1122 MSG_ERROR(message) 1123 end if 1124 end if 1125 !####################################################################### 1126 1127 !3) PAW only : precondition the rhoij quantities (augmentation 1128 !occupancies) residuals. Use a simple preconditionning 1129 !with the same mixing factor as the model dielectric function. 1130 1131 if (psps%usepaw==1.and.my_natom>0) then 1132 if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then 1133 mixfac=one;mixfacmag=one 1134 else 1135 mixfac=dielar(4);mixfacmag=abs(dielar(7)) 1136 end if 1137 if (pawrhoij(1)%cplex==1) then 1138 index=0 1139 do iatom=1,my_natom 1140 do ispden=1,pawrhoij(iatom)%nspden 1141 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 1142 do kmix=1,pawrhoij(iatom)%lmnmix_sz 1143 index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix) 1144 rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 1145 end do 1146 end do 1147 end do 1148 else 1149 index=-1 1150 do iatom=1,my_natom 1151 do ispden=1,pawrhoij(iatom)%nspden 1152 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 1153 do kmix=1,pawrhoij(iatom)%lmnmix_sz 1154 index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 1155 rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 1156 end do 1157 end do 1158 end do 1159 end if 1160 end if 1161 !####################################################################### 1162 1163 !4) Take care of the change of atomic positions 1164 !Note : this part is very demanding on memory... 1165 !however, since this algorithm is still in development, 1166 !it was NOT included in the estimation provided by memory.f 1167 if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then 1168 1169 ! Not yet compatible with resid given in reciprocal space 1170 if (optreal/=1) then 1171 write(message, '(5a)' )& 1172 & 'From the calling routine, densfor_pred=3',ch10,& 1173 & 'You cannot use residuals in reciprocal space.',ch10,& 1174 & 'Action: correct your input file.' 1175 MSG_ERROR(message) 1176 end if 1177 1178 ! Not compatible with non-collinear magnetism 1179 if(dtset%nspden==4)then 1180 MSG_ERROR('densfor_pred=3 does not work for nspden=4!') 1181 end if 1182 1183 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1184 nfftot=PRODUCT(ngfft(1:3)) 1185 1186 ! First subtract the current local, hartree and exchange correlation potentials 1187 do ispden=1,min(dtset%nspden,2) 1188 vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden) 1189 end do 1190 if (dtset%nspden==4) then 1191 do ispden=3,4 1192 vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden) 1193 end do 1194 end if 1195 1196 ! Compute the modified density, in rhor_wk 1197 option=2 1198 ABI_ALLOCATE(gresid,(3,dtset%natom)) 1199 ABI_ALLOCATE(grxc,(3,dtset%natom)) 1200 ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden)) 1201 ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden)) 1202 ABI_ALLOCATE(xred_wk,(3,dtset%natom)) 1203 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 1204 1205 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 1206 & ntypat,option,pawtab,rhor,rprimd,& 1207 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 1208 1209 ! Compute up+down rhog_wk(G) by fft 1210 ABI_ALLOCATE(work,(nfft)) 1211 ABI_ALLOCATE(rhog_wk,(2,nfft)) 1212 work(:)=rhor_wk(:,1) 1213 call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 1214 ABI_DEALLOCATE(work) 1215 1216 ! Compute structure factor phases for new atomic pos: 1217 call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk) 1218 1219 ! Compute local ionic pseudopotential vpsp: 1220 ! and core electron density xccc3d, if needed. 1221 n3xccc=0;if (n1xccc/=0) n3xccc=nfft 1222 ABI_ALLOCATE(xccc3d,(n3xccc)) 1223 ABI_ALLOCATE(vpsp_wk,(nfft)) 1224 vprtrb(1:2)=zero 1225 1226 ! Determine by which method the local ionic potential and/or 1227 ! the pseudo core charge density have to be computed 1228 ! Local ionic potential: 1229 ! Method 1: PAW 1230 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 1231 vloc_method=1;if (psps%usepaw==0) vloc_method=2 1232 if (dtset%icoulomb>0) vloc_method=2 1233 if (psps%usewvl==1) vloc_method=2 1234 ! Pseudo core charge density: 1235 ! Method 1: PAW, nc_xccc_gspace 1236 ! Method 2: Norm-conserving PP, wavelets 1237 coredens_method=1;if (psps%usepaw==0) coredens_method=2 1238 if (psps%nc_xccc_gspace==1) coredens_method=1 1239 if (psps%nc_xccc_gspace==0) coredens_method=2 1240 if (psps%usewvl==1) coredens_method=2 1241 1242 ! Local ionic potential and/or pseudo core charge by method 1 1243 if (vloc_method==1.or.coredens_method==1) then 1244 optv=0;if (vloc_method==1) optv=1 1245 optn=0;if (coredens_method==1) optn=n3xccc/nfft 1246 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1 1247 ! Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused... 1248 call atm2fft(atindx,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,gmet,& 1249 & gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,& 1250 & nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,& 1251 & psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dummy_in,dummy_out6,dummy_out7,ucvol,& 1252 & psps%usepaw,dummy_in,dummy_in,dummy_in,vprtrb,psps%vlspl,& 1253 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 1254 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 1255 end if 1256 1257 ! Local ionic potential by method 2 1258 if (vloc_method==2) then 1259 option=1 1260 ABI_ALLOCATE(dyfrlo_indx,(3,3,dtset%natom)) 1261 ABI_ALLOCATE(grtn_indx,(3,dtset%natom)) 1262 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,& 1263 & mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,& 1264 & ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,& 1265 & ucvol,vprtrb,vpsp_wk,wvl%descr,wvl%den,xred) 1266 ABI_DEALLOCATE(dyfrlo_indx) 1267 ABI_DEALLOCATE(grtn_indx) 1268 end if 1269 1270 ! Pseudo core electron density by method 2 1271 if (coredens_method==2.and.n1xccc/=0) then 1272 option=1 1273 ABI_ALLOCATE(dyfrx2,(3,3,dtset%natom)) 1274 ABI_ALLOCATE(grxc_indx,(3,dtset%natom)) 1275 call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 1276 & n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,& 1277 & psps%xccc1d,xccc3d,xred_wk) 1278 ABI_DEALLOCATE(dyfrx2) 1279 ABI_DEALLOCATE(grxc_indx) 1280 end if 1281 1282 ! Compute Hartree+xc potentials 1283 ABI_ALLOCATE(vxc_wk,(nfft,dtset%nspden)) 1284 ABI_ALLOCATE(vhartr_wk,(nfft)) 1285 option=1 1286 1287 call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog_wk,rprimd,vhartr_wk) 1288 1289 ! Prepare the call to rhotoxc 1290 call xcdata_init(xcdata,dtset=dtset) 1291 nk3xc=1 1292 ABI_ALLOCATE(work,(0)) 1293 call rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft,& 1294 & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor_wk,rprimd,strsxc,1,& 1295 & vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk) 1296 ABI_DEALLOCATE(work) 1297 ABI_DEALLOCATE(xccc3d) 1298 1299 ! Sum all contributions 1300 do ispden=1,min(dtset%nspden,2) 1301 do ifft=1,nfft 1302 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden) 1303 end do 1304 end do 1305 if (dtset%nspden==4) then 1306 do ispden=3,4 1307 do ifft=1,nfft 1308 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden) 1309 end do 1310 end do 1311 end if 1312 call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft) 1313 if(dtset%nspden==2) then 1314 vmean(1)=half*(vmean(1)+vmean(2)) 1315 vmean(2)=vmean(1) 1316 end if 1317 do ispden=1,dtset%nspden 1318 vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden) 1319 end do 1320 ABI_DEALLOCATE(gresid) 1321 ABI_DEALLOCATE(grxc) 1322 ABI_DEALLOCATE(rhog_wk) 1323 ABI_DEALLOCATE(rhor_wk) 1324 ABI_DEALLOCATE(rhor_wk0) 1325 ABI_DEALLOCATE(xred_wk) 1326 ABI_DEALLOCATE(vhartr_wk) 1327 ABI_DEALLOCATE(vpsp_wk) 1328 ABI_DEALLOCATE(vxc_wk) 1329 1330 end if 1331 1332 end subroutine prcref_PMA
ABINIT/prcrskerker1 [ Functions ]
NAME
prcrskerker1
FUNCTION
preconditionning by a real-space conjugate gradient on residual using a model dielectric function in real space
INPUTS
nfft=number of fft grid points nspden=number of spin-density components ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft dielar(7)=input parameters for dielectric matrix: diecut,dielng,diemac,diemix,diegap,dielam,diemixmag. gprimd(3,3)=dimensional primitive translations in fourier space (bohr**-1) rprimd(3,3)=dimensional primitive translations in real space (bohr) vresid(nfft,nspden)=residual potential base(nfft) = real space function used as a basis to guess a fine dielectric funtion see the calling routine to know the content
OUTPUT
vrespc(nfft,nspden)=preconditioned residual of the potential
WARNINGS
This is experimental code : input, ouptput, results and any other feature may vary greatly.
NOTES
needs severe cleaning and this is abuse of modules as common blocks...
PARENTS
prcref,prcref_PMA
CHILDREN
cgpr,frskerker1__end,frskerker1__init,laplacian,prc_mem_init
SOURCE
2480 subroutine prcrskerker1(dtset,mpi_enreg,nfft,nspden,ngfft,dielar,etotal,gprimd,vresid,vrespc,base) 2481 2482 2483 !This section has been created automatically by the script Abilint (TD). 2484 !Do not modify the following lines by hand. 2485 #undef ABI_FUNC 2486 #define ABI_FUNC 'prcrskerker1' 2487 !End of the abilint section 2488 2489 implicit none 2490 2491 !Arguments ------------------------------------ 2492 !scalars 2493 integer,intent(in) :: nfft,nspden 2494 real(dp) :: etotal 2495 type(MPI_type),intent(in) :: mpi_enreg 2496 type(dataset_type),intent(in) :: dtset 2497 !arrays 2498 integer,intent(in) :: ngfft(18) 2499 real(dp),intent(in) :: base(nfft),dielar(7),gprimd(3,3) 2500 real(dp),intent(in) :: vresid(nfft,nspden) 2501 real(dp),intent(out) :: vrespc(nfft,nspden) 2502 2503 !Local variables------------------------------- 2504 !scalars 2505 integer :: ifft,ispden,n1,n2,n3 2506 real(dp) :: base_delta,base_max,base_min,dielng,diemac,diemix 2507 real(dp) :: diemixmag 2508 real(dp) :: rdummy1,rdummy2 2509 logical :: new_prc_func 2510 !arrays 2511 real(dp) :: deltaW(nfft,nspden) 2512 real(dp) :: g2cart(nfft) 2513 real(dp) :: mat(nfft,nspden) 2514 2515 ! ************************************************************************* 2516 2517 !DEBUG 2518 !write(std_out,*)' prckerker1 : enter ' 2519 !ENDDEBUG 2520 !if(cycle==0) then 2521 call prc_mem_init(nfft) 2522 2523 if(cycle==0) then 2524 new_prc_func=.TRUE. 2525 energy_min=etotal 2526 else if(etotal < energy_min) then 2527 new_prc_func=.TRUE. 2528 energy_min=etotal 2529 else 2530 new_prc_func=.FALSE. 2531 end if 2532 2533 2534 dielng=dielar(2) 2535 diemac=dielar(3) 2536 diemix=dielar(4) 2537 diemixmag=dielar(7) 2538 !****************************************************************** 2539 !compute the diemac(r) ** 2540 !****************************************************************** 2541 !this task will be devoted to a general function later 2542 n1=ngfft(1) 2543 n2=ngfft(2) 2544 n3=ngfft(3) 2545 !base_cp=base 2546 if(new_prc_func) then 2547 base_min=base(1) 2548 base_max=base(1) 2549 do ifft=1,nfft 2550 base_min = min(base_min,base(ifft)) 2551 base_max = max(base_max,base(ifft)) 2552 end do 2553 base_delta = base_max - base_min 2554 ! if(cycle.lt.2) then 2555 rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one) 2556 ! else 2557 ! rdiemac(:) = rdiemac(:)*0.5_dp+0.5_dp*(((base(:)-base_min) / (base_delta) ) *(diemac-one) + one) 2558 ! end if 2559 ! if(cycle==0) rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one) 2560 ! rdiemac(:) = exp(((base(:)-base_min) / (base_delta) *log(diemac))) 2561 end if 2562 cycle=cycle+1 2563 !if(cycle==5) cycle=0 2564 !end if 2565 !****************************************************************** 2566 !compute deltaW ** 2567 !****************************************************************** 2568 vrespc=vresid !starting point 2569 call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,dtset%paral_kgb,rdfuncr=vrespc,laplacerdfuncr=deltaW,g2cart_out=g2cart) ! put the laplacian of the residuals into deltaW 2570 !call laplacian(vrespc,buffer,ngfft,gprimd) ! put the laplacian of the residuals into deltaW 2571 !do ifft=1,nfft 2572 !if (buffer(ifft,1)/=deltaW(ifft,1)) then 2573 !stop 2574 !end if 2575 !end do 2576 deltaW(:,1)= diemix*(((one/rdiemac(:))*vresid(:,1))-(((dielng)**2)*deltaW(:,1))) 2577 if (nspden>1.and.(diemixmag>=zero)) then 2578 do ispden=2,nspden 2579 deltaW(:,ispden)= abs(diemixmag)*(((one/rdiemac(:))*vresid(:,ispden))-(((dielng)**2)*deltaW(:,ispden))) 2580 end do 2581 end if 2582 !call random_number(deltaW) 2583 !call random_number(vrespc) 2584 !****************************************************************** 2585 !Finding the preconditionned residuals which minimizes ** 2586 !half*(vrespc*(1-dielng2/4pi2 nabla2) vrespc) - vrespc * deltaW ** 2587 !*********************************************************************** 2588 vrespc(:,1)=diemix*vrespc(:,1) 2589 if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*vrespc(:,2:nspden) 2590 !buffer=vrespc 2591 2592 2593 !============================================================================== 2594 !============================================================================== 2595 !! Original loop 2596 !============================================================================== 2597 !============================================================================== 2598 2599 call frskerker1__init(dtset,mpi_enreg,nfft,ngfft,nspden,dielng,deltaW,gprimd,mat,g2cart) 2600 2601 !call cgpr(pf_rscgres,dpf_rscgres,newvres,real(1e-40,dp),700,vrespc,rdummy1,rdummy2) 2602 !rdummy1 = pf_rscgres(nfft,nspden,vrespc) 2603 call cgpr(nfft,nspden,frskerker1__pf,frskerker1__dpf,frskerker1__newvres,& 2604 & real(1e-10,dp),700,vrespc,rdummy1,rdummy2) 2605 call frskerker1__end() 2606 2607 !============================================================================== 2608 !============================================================================== 2609 !! Original loop end 2610 !============================================================================== 2611 !============================================================================== 2612 2613 2614 !cplex=1 2615 !qphon(:)=zero 2616 !call moddiel(cplex,dielar,nfft,ngfft,nspden,1,0,qphon,rprimd,vresid,buffer) 2617 !c1=0 2618 !do ifft=1,nfft,1 2619 !if((abs(buffer(ifft,1)-vrespc(ifft,1))/(abs(buffer(ifft,1)+vrespc(ifft,1))*half)) > 5e-3) then 2620 !c1=c1+1 2621 !end if 2622 !end do 2623 !call laplacian(vrespc,buffer,ngfft,gprimd) 2624 !buffer=vrespc(:,:)-buffer(:,:)*dielng**2 2625 !c2=0 2626 !do ifft=1,nfft,1 2627 !if((abs(buffer(ifft,1)-deltaW(ifft,1))/(abs(buffer(ifft,1)+deltaW(ifft,1))*half)) > 5e-3) then 2628 !c2=c2+1 2629 !end if 2630 !end do 2631 !!! !stop 2632 !call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,& 2633 !& g2cart_out=g2cart) 2634 2635 !vrespc=vresid 2636 !do ispden=1,nspden 2637 !call fourdp(1, gvrespc(:,:,ispden), vrespc(:,ispden),-1,mpi_enreg,nfft,ngfft,0) 2638 !end do 2639 !filtering 2640 !do ispden=1,nspden 2641 !do ifft=1,nfft 2642 !! gvrespc(:,ifft,ispden)=(one-exp(-g2cart(ifft)*15.0_dp))*gvrespc(:,ifft,ispden) 2643 !! gvrespc(:,ifft,ispden)=(exp(-g2cart(ifft)*10.0_dp))*gvrespc(:,ifft,ispden) 2644 !! gvrespc(:,ifft,ispden)=(one-one/(exp(-0.002_dp/g2cart(ifft)**2)+one))*gvrespc(:,ifft,ispden) 2645 !gvrespc(:,ifft,ispden)=(two-2_dp/(exp(-0.008_dp/(g2cart(ifft)+0.0012_dp))+one))*gvrespc(:,ifft,ispden) 2646 !gvrespc(:,ifft,ispden)=min(one,(sqrt(g2cart(ifft)/0.006_dp))**(one))*gvrespc(:,ifft,ispden) 2647 !end do 2648 !end do 2649 !change resulting potential to real space 2650 !do ispden=1,nspden 2651 !call fourdp(1,gvrespc(:,:,ispden),vrespc(:,ispden),1,mpi_enreg,nfft,ngfft,0) 2652 !end do 2653 !vrespc=vrespc*diemix 2654 !maxg2=g2cart(1) 2655 !ming2=g2cart(5) 2656 !do ifft=1,nfft 2657 !maxg2=max(g2cart(ifft),maxg2) 2658 !if(g2cart(ifft) .gt. zero) ming2=min(g2cart(ifft),ming2) 2659 !end do 2660 !stop 2661 2662 !DEBUG 2663 !write(std_out,*)' prckerker1 : exit ' 2664 !ENDDEBUG 2665 2666 end subroutine prcrskerker1
ABINIT/prcrskerker2 [ Functions ]
NAME
prcrskerker2
FUNCTION
preconditionning by a real-space conjugate gradient on residual using a model dielectric function in real space differing from prcrskerker1 by the use of a linear response approach
INPUTS
nfft=number of fft grid points nspden=number of spin-density components ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft dielar(7)=input parameters for dielectric matrix: diecut,dielng,diemac,diemix,diegap,dielam,diemixmag. gprimd(3,3)=dimensional primitive translations in fourier space (bohr**-1) rprimd(3,3)=dimensional primitive translations in real space (bohr) vresid(nfft,nspden)=residual potential
OUTPUT
vrespc(nfft,nspden)=preconditioned residual of the potential
WARNINGS
This is experimental code : input, ouptput, results and any other feature may vary greatly.
NOTES
PARENTS
prcref,prcref_PMA
CHILDREN
cgpr,dotprod_vn,frskerker2__end,frskerker2__init,laplacian,ptabs_fourdp
SOURCE
2705 subroutine prcrskerker2(dtset,nfft,nspden,ngfft,dielar,gprimd,rprimd,vresid,vrespc,natom,xred,mpi_enreg,ucvol) 2706 2707 2708 !This section has been created automatically by the script Abilint (TD). 2709 !Do not modify the following lines by hand. 2710 #undef ABI_FUNC 2711 #define ABI_FUNC 'prcrskerker2' 2712 !End of the abilint section 2713 2714 implicit none 2715 2716 !Arguments ------------------------------------ 2717 !scalars 2718 integer,intent(in) :: natom,nfft,nspden 2719 real(dp),intent(in) :: ucvol 2720 type(MPI_type),intent(in) :: mpi_enreg 2721 type(dataset_type),intent(in) :: dtset 2722 !arrays 2723 integer,intent(in) :: ngfft(18) 2724 real(dp),intent(in) :: dielar(7),gprimd(3,3),rprimd(3,3),vresid(nfft,nspden) 2725 real(dp),intent(in) :: xred(3,natom) 2726 real(dp),intent(out) :: vrespc(nfft,nspden) 2727 2728 !Local variables------------------------------- 2729 !logical,save ::ok=.FALSE. 2730 !scalars 2731 integer :: cplex,i1,i2,i3,iatom,iatom27,ifft,ispden,n1,n2,n3,natom27,nfftotf 2732 integer :: option 2733 real(dp),save :: lastp1=one,lastp2=one 2734 real(dp) :: C1,C2,DE,core,dielng,diemac,diemix,diemixmag,doti,dr,l1,l2,l3,l4,r 2735 real(dp) :: rdummy1,rdummy2,rmin,xr,y,yr,zr 2736 !arrays 2737 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 2738 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 2739 real(dp) :: V1(nfft,nspden),V2(nfft,nspden),buffer(nfft,nspden) 2740 real(dp) :: deltaW(nfft,nspden) 2741 real(dp) :: mat(nfft,nspden) 2742 real(dp) :: rdielng(nfft),rdiemac(nfft),xcart(3,natom) 2743 real(dp) :: xcart27(3,natom*27) 2744 2745 ! ************************************************************************* 2746 2747 dielng=dielar(2) 2748 diemac=dielar(3) 2749 diemix=dielar(4) 2750 diemixmag=dielar(7) 2751 !****************************************************************** 2752 !compute the diemac(r) ** 2753 !****************************************************************** 2754 !this task will be devoted to a general function later 2755 n1=ngfft(1) 2756 n2=ngfft(2) 2757 n3=ngfft(3) 2758 nfftotf=n1*n2*n3 2759 !if(.not.ok) then 2760 xcart(1,:)=xred(1,:)*rprimd(1,1)+xred(2,:)*rprimd(1,2)+xred(3,:)*rprimd(1,3) 2761 xcart(2,:)=xred(1,:)*rprimd(2,1)+xred(2,:)*rprimd(2,2)+xred(3,:)*rprimd(2,3) 2762 xcart(3,:)=xred(1,:)*rprimd(3,1)+xred(2,:)*rprimd(3,2)+xred(3,:)*rprimd(3,3) 2763 2764 iatom27=0 2765 do i1=-1,1 2766 do i2=-1,1 2767 do i3=-1,1 2768 do iatom=1,natom 2769 iatom27=iatom27+1 2770 xcart27(:,iatom27)=xcart(:,iatom)+rprimd(:,1)*i1+rprimd(:,2)*i2+rprimd(:,3)*i3 2771 end do 2772 end do 2773 end do 2774 end do 2775 2776 !stop 2777 natom27=27*natom 2778 2779 l1=0.34580850339844665 2780 !l2=0.5123510203906797 !0.41242551019533985 2781 !l3=0.8001489796093203 !0.90007448980466009 2782 2783 l2=0.41242551019533985 2784 l3=0.90007448980466009 2785 l4=0.9666914966015534 2786 2787 2788 l1=0.31387233559896449 2789 l2=0.35828367346355994 2790 l3=0.9333829932031068 2791 l4=0.9777943310677023 2792 2793 l1=3.5 2794 l2=11.5 2795 l3=2.5 2796 l4=6.5 2797 !l1=30. !cellules pleines 2798 2799 rdielng=zero 2800 core=1. !(value of Er at the core of atoms) 2801 dr=2.65 ! radius of atoms=2.65165 2802 y=1. ! value of Er in the empty region 2803 2804 !Get the distrib associated with this fft_grid 2805 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2806 2807 do i3=1,n3 2808 ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft) 2809 do i2=1,n2 2810 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 2811 do i1=1,n1 2812 ifft=ifft+1 2813 ! !!!!!!!!!!!!!!!!!!!!!!!!! 2814 ! ! calculation of the simplest part void/metal 2815 ! !! x=real(real(i3,dp)/real(n3,dp),dp) 2816 ! !! !x=i3/n3 2817 ! !! if(x < l1) then 2818 ! !! rdiemac(ifft)=diemac 2819 ! !! rdielng(ifft)=dielng 2820 ! !! else if(x < l2) then 2821 ! !! xp=(l2-x)/(l2-l1) 2822 ! !! rdiemac(ifft)=y+(diemac-y)& 2823 ! !! & * (1.-(1.-xp)**4)**4 2824 ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4 2825 ! !! else if(x < l3) then 2826 ! !! rdiemac(ifft)=y 2827 ! !! rdielng(ifft)=zero 2828 ! !! else if(x < l4) then 2829 ! !! xp=(l3-x)/(l3-l4) 2830 ! !! rdiemac(ifft)=y+(diemac-y)& 2831 ! !! & * (1.-(1.-xp)**4)**4 2832 ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4 2833 ! !! else 2834 ! !! rdiemac(ifft)=diemac 2835 ! !! rdielng(ifft)=dielng 2836 ! !! end if 2837 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2838 ! !!!! calculation of atomic core dielectric 2839 ! !! rmin=1e16 2840 ! !! xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)& 2841 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(1,3) 2842 ! !! yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)& 2843 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(2,3) 2844 ! !! zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)& 2845 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(3,3) 2846 ! !! do iatom=1,natom27 2847 ! !! r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2 2848 ! !! if (r<rmin) then 2849 ! !! rmin=r 2850 ! !! end if 2851 ! !! end do 2852 ! !! if(rmin < dr**2) then 2853 ! !! rdiemac(ifft)=min(rdiemac(ifft),core+(diemac-core)*(1.-(1.-sqrt(rmin)/dr)**2)**2) 2854 ! !! rdielng(ifft)=dielng-dielng*(1.-(1.-sqrt(rmin)/dr)**4)**4 2855 ! !! else 2856 ! !! rdiemac(ifft)=min(rdiemac(ifft),diemac) 2857 ! !! end if 2858 rmin=1e16 2859 xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)& 2860 & +real((i3-1),dp)/real(n3,dp)*rprimd(1,3) 2861 yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)& 2862 & +real((i3-1),dp)/real(n3,dp)*rprimd(2,3) 2863 zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)& 2864 & +real((i3-1),dp)/real(n3,dp)*rprimd(3,3) 2865 2866 rdiemac(ifft)=y 2867 rdielng(ifft)=zero 2868 do iatom=1,natom27 2869 r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2 2870 if (r<rmin) then 2871 rmin=r 2872 end if 2873 if(r < l1) then 2874 rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y) 2875 else if(r < l2) then 2876 rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y)*(one-((sqrt(r)-l1)/(l2-l1))**2)**2 2877 else 2878 rdiemac(ifft)=rdiemac(ifft) 2879 end if 2880 if(r < l3) then 2881 rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng) 2882 else if(r < l4) then 2883 rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng) *(one-((sqrt(r)-l3)/(l4-l3))**2)**2 2884 end if 2885 end do 2886 2887 rdielng(ifft)=min(rdielng(ifft),dielng) 2888 ! rdielng(ifft)=dielng 2889 rdiemac(ifft)=min(rdiemac(ifft),diemac) 2890 rdiemac(ifft)=diemac 2891 end do 2892 end if 2893 end do 2894 end do 2895 !rdielng(:)=dielng 2896 2897 !**************************************************************************************** 2898 !**************************************************************************************** 2899 !**************************************************************************************** 2900 !**************************************************************************************** 2901 !****************************************************************** 2902 !compute V1 2903 !****************************************************************** 2904 V1=vresid 2905 call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,dtset%paral_kgb,rdfuncr=V1,laplacerdfuncr=deltaW) 2906 deltaW(:,1)= (((one/rdiemac(:))*V1(:,1))-(((rdielng(:))**2)*deltaW(:,1))) 2907 !deltaW(:,1)= -diemix*(((rdielng(:))**2)*deltaW(:,ispden)) 2908 if (nspden>1) then 2909 do ispden=2,nspden 2910 deltaW(:,ispden)= (((one/rdiemac(:))*V1(:,ispden))-(((rdielng(:))**2)*deltaW(:,ispden))) 2911 ! deltaW(:,ispden)= -abs(diemixmag)*(((rdielng(:))**2)*deltaW(:,ispden)) 2912 end do 2913 end if 2914 call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat) 2915 call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,& 2916 & frskerker2__newvres2,lastp1*real(1e-6 ,dp),700,V1,rdummy1,rdummy2) 2917 lastp1=min(abs(rdummy1),1e-6_dp) 2918 call frskerker2__end() 2919 2920 !****************************************************************** 2921 !compute V2 2922 !****************************************************************** 2923 V2=vresid 2924 do ispden=1,nspden 2925 deltaW(:,ispden)= (rdielng(:)**2) 2926 end do 2927 call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat) 2928 call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,& 2929 & frskerker2__newvres2,lastp2*real(1e-6,dp),700,V2,rdummy1,rdummy2) 2930 lastp2=min(abs(rdummy1),1e-6_dp) 2931 call frskerker2__end() 2932 2933 2934 !****************************************************************** 2935 !compute C1, C2 & DE 2936 !****************************************************************** 2937 cplex=1; 2938 option=1; 2939 call dotprod_vn(cplex,& !complex density/pot 2940 &rdielng,& !the density 2941 &DE,& !resulting dorproduct integrated over r ! here DE is used has a buffer 2942 &doti,& !imaginary part of the integral 2943 &size(rdielng,1),& !number of localy(cpu) attributed grid point 2944 &nfftotf,& !real total number of grid point 2945 &nspden,& !nspden 2946 &option,& !1=compute only the real part 2=compute also the imaginary part 2947 &rdielng,& !the potential 2948 &ucvol,& !cell volume 2949 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 2950 do ispden=1,nspden 2951 buffer(:,ispden)=rdielng(:)*V1(:,ispden) 2952 end do 2953 call dotprod_vn(cplex,& !complex density/pot 2954 &rdielng,& !the density 2955 &C1,& !resulting dorproduct integrated over r ! here DE is used has a buffer 2956 &doti,& !imaginary part of the integral 2957 &size(rdielng,1),& !number of localy(cpu) attributed grid point 2958 &nfftotf,& !real total number of grid point 2959 &nspden,& !nspden 2960 &option,& !1=compute only the real part 2=compute also the imaginary part 2961 &buffer,& !the potential 2962 &ucvol,& !cell volume 2963 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 2964 do ispden=1,nspden 2965 buffer(:,ispden)=rdielng(:)*V2(:,ispden) 2966 end do 2967 call dotprod_vn(cplex,& !complex density/pot 2968 &rdielng,& !the density 2969 &C2,& !resulting dorproduct integrated over r ! here DE is used has a buffer 2970 &doti,& !imaginary part of the integral 2971 &size(rdielng,1),& !number of localy(cpu) attributed grid point 2972 &nfftotf,& !real total number of grid point 2973 &nspden,& !nspden 2974 &option,& !1=compute only the real part 2=compute also the imaginary part 2975 &buffer,& !the potential 2976 &ucvol,& !cell volume 2977 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 2978 C1=C1/DE 2979 C2=C2/DE 2980 DE=C1/(one-C2) 2981 2982 !****************************************************************** 2983 !compute the new preconditionned residuals 2984 !****************************************************************** 2985 vrespc(:,1)=diemix*(V1(:,1)+DE*V2(:,1)) 2986 if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*(V1(:,2:nspden)+DE*V2(:,2:nspden)) 2987 2988 end subroutine prcrskerker2