TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (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 . For the initials of contributors, see ~ABINIT/Infos/contributors .
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
SIDE EFFECTS
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
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 54 subroutine prcrskerker2(dtset,nfft,nspden,ngfft,dielar,gprimd,rprimd,vresid,vrespc,natom,xred,mpi_enreg,ucvol) 55 56 use defs_basis 57 use defs_abitypes 58 use m_errors 59 use m_profiling_abi 60 use frskerker2 61 62 use m_cgtools, only : dotprod_vn 63 use m_mpinfo, only : ptabs_fourdp 64 65 !This section has been created automatically by the script Abilint (TD). 66 !Do not modify the following lines by hand. 67 #undef ABI_FUNC 68 #define ABI_FUNC 'prcrskerker2' 69 use interfaces_56_recipspace 70 use interfaces_62_cg_noabirule 71 !End of the abilint section 72 73 implicit none 74 75 !Arguments ------------------------------------ 76 !scalars 77 integer,intent(in) :: natom,nfft,nspden 78 real(dp),intent(in) :: ucvol 79 type(MPI_type),intent(in) :: mpi_enreg 80 type(dataset_type),intent(in) :: dtset 81 !arrays 82 integer,intent(in) :: ngfft(18) 83 real(dp),intent(in) :: dielar(7),gprimd(3,3),rprimd(3,3),vresid(nfft,nspden) 84 real(dp),intent(in) :: xred(3,natom) 85 real(dp),intent(out) :: vrespc(nfft,nspden) 86 87 !Local variables------------------------------- 88 !logical,save ::ok=.FALSE. 89 !scalars 90 integer :: cplex,i1,i2,i3,iatom,iatom27,ifft,ispden,n1,n2,n3,natom27,nfftotf 91 integer :: option 92 real(dp),save :: lastp1=one,lastp2=one 93 real(dp) :: C1,C2,DE,core,dielng,diemac,diemix,diemixmag,doti,dr,l1,l2,l3,l4,r 94 real(dp) :: rdummy1,rdummy2,rmin,xr,y,yr,zr 95 !arrays 96 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 97 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 98 real(dp) :: V1(nfft,nspden),V2(nfft,nspden),buffer(nfft,nspden) 99 real(dp) :: deltaW(nfft,nspden) 100 real(dp) :: mat(nfft,nspden) 101 real(dp) :: rdielng(nfft),rdiemac(nfft),xcart(3,natom) 102 real(dp) :: xcart27(3,natom*27) 103 104 ! ************************************************************************* 105 106 dielng=dielar(2) 107 diemac=dielar(3) 108 diemix=dielar(4) 109 diemixmag=dielar(7) 110 !****************************************************************** 111 !compute the diemac(r) ** 112 !****************************************************************** 113 !this task will be devoted to a general function later 114 n1=ngfft(1) 115 n2=ngfft(2) 116 n3=ngfft(3) 117 nfftotf=n1*n2*n3 118 !if(.not.ok) then 119 xcart(1,:)=xred(1,:)*rprimd(1,1)+xred(2,:)*rprimd(1,2)+xred(3,:)*rprimd(1,3) 120 xcart(2,:)=xred(1,:)*rprimd(2,1)+xred(2,:)*rprimd(2,2)+xred(3,:)*rprimd(2,3) 121 xcart(3,:)=xred(1,:)*rprimd(3,1)+xred(2,:)*rprimd(3,2)+xred(3,:)*rprimd(3,3) 122 123 iatom27=0 124 do i1=-1,1 125 do i2=-1,1 126 do i3=-1,1 127 do iatom=1,natom 128 iatom27=iatom27+1 129 xcart27(:,iatom27)=xcart(:,iatom)+rprimd(:,1)*i1+rprimd(:,2)*i2+rprimd(:,3)*i3 130 end do 131 end do 132 end do 133 end do 134 135 !stop 136 natom27=27*natom 137 138 l1=0.34580850339844665 139 !l2=0.5123510203906797 !0.41242551019533985 140 !l3=0.8001489796093203 !0.90007448980466009 141 142 l2=0.41242551019533985 143 l3=0.90007448980466009 144 l4=0.9666914966015534 145 146 147 l1=0.31387233559896449 148 l2=0.35828367346355994 149 l3=0.9333829932031068 150 l4=0.9777943310677023 151 152 l1=3.5 153 l2=11.5 154 l3=2.5 155 l4=6.5 156 !l1=30. !cellules pleines 157 158 rdielng=zero 159 core=1. !(value of Er at the core of atoms) 160 dr=2.65 ! radius of atoms=2.65165 161 y=1. ! value of Er in the empty region 162 163 !Get the distrib associated with this fft_grid 164 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 165 166 do i3=1,n3 167 ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft) 168 do i2=1,n2 169 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 170 do i1=1,n1 171 ifft=ifft+1 172 ! !!!!!!!!!!!!!!!!!!!!!!!!! 173 ! ! calculation of the simplest part void/metal 174 ! !! x=real(real(i3,dp)/real(n3,dp),dp) 175 ! !! !x=i3/n3 176 ! !! if(x < l1) then 177 ! !! rdiemac(ifft)=diemac 178 ! !! rdielng(ifft)=dielng 179 ! !! else if(x < l2) then 180 ! !! xp=(l2-x)/(l2-l1) 181 ! !! rdiemac(ifft)=y+(diemac-y)& 182 ! !! & * (1.-(1.-xp)**4)**4 183 ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4 184 ! !! else if(x < l3) then 185 ! !! rdiemac(ifft)=y 186 ! !! rdielng(ifft)=zero 187 ! !! else if(x < l4) then 188 ! !! xp=(l3-x)/(l3-l4) 189 ! !! rdiemac(ifft)=y+(diemac-y)& 190 ! !! & * (1.-(1.-xp)**4)**4 191 ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4 192 ! !! else 193 ! !! rdiemac(ifft)=diemac 194 ! !! rdielng(ifft)=dielng 195 ! !! end if 196 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 197 ! !!!! calculation of atomic core dielectric 198 ! !! rmin=1e16 199 ! !! xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)& 200 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(1,3) 201 ! !! yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)& 202 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(2,3) 203 ! !! zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)& 204 ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(3,3) 205 ! !! do iatom=1,natom27 206 ! !! r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2 207 ! !! if (r<rmin) then 208 ! !! rmin=r 209 ! !! end if 210 ! !! end do 211 ! !! if(rmin < dr**2) then 212 ! !! rdiemac(ifft)=min(rdiemac(ifft),core+(diemac-core)*(1.-(1.-sqrt(rmin)/dr)**2)**2) 213 ! !! rdielng(ifft)=dielng-dielng*(1.-(1.-sqrt(rmin)/dr)**4)**4 214 ! !! else 215 ! !! rdiemac(ifft)=min(rdiemac(ifft),diemac) 216 ! !! end if 217 rmin=1e16 218 xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)& 219 & +real((i3-1),dp)/real(n3,dp)*rprimd(1,3) 220 yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)& 221 & +real((i3-1),dp)/real(n3,dp)*rprimd(2,3) 222 zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)& 223 & +real((i3-1),dp)/real(n3,dp)*rprimd(3,3) 224 225 rdiemac(ifft)=y 226 rdielng(ifft)=zero 227 do iatom=1,natom27 228 r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2 229 if (r<rmin) then 230 rmin=r 231 end if 232 if(r < l1) then 233 rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y) 234 else if(r < l2) then 235 rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y)*(one-((sqrt(r)-l1)/(l2-l1))**2)**2 236 else 237 rdiemac(ifft)=rdiemac(ifft) 238 end if 239 if(r < l3) then 240 rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng) 241 else if(r < l4) then 242 rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng) *(one-((sqrt(r)-l3)/(l4-l3))**2)**2 243 end if 244 end do 245 246 rdielng(ifft)=min(rdielng(ifft),dielng) 247 ! rdielng(ifft)=dielng 248 rdiemac(ifft)=min(rdiemac(ifft),diemac) 249 rdiemac(ifft)=diemac 250 end do 251 end if 252 end do 253 end do 254 !rdielng(:)=dielng 255 256 !**************************************************************************************** 257 !**************************************************************************************** 258 !**************************************************************************************** 259 !**************************************************************************************** 260 !****************************************************************** 261 !compute V1 262 !****************************************************************** 263 V1=vresid 264 call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,dtset%paral_kgb,rdfuncr=V1,laplacerdfuncr=deltaW) 265 deltaW(:,1)= (((one/rdiemac(:))*V1(:,1))-(((rdielng(:))**2)*deltaW(:,1))) 266 !deltaW(:,1)= -diemix*(((rdielng(:))**2)*deltaW(:,ispden)) 267 if (nspden>1) then 268 do ispden=2,nspden 269 deltaW(:,ispden)= (((one/rdiemac(:))*V1(:,ispden))-(((rdielng(:))**2)*deltaW(:,ispden))) 270 ! deltaW(:,ispden)= -abs(diemixmag)*(((rdielng(:))**2)*deltaW(:,ispden)) 271 end do 272 end if 273 call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat) 274 call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,& 275 & frskerker2__newvres2,lastp1*real(1e-6 ,dp),700,V1,rdummy1,rdummy2) 276 lastp1=min(abs(rdummy1),1e-6_dp) 277 call frskerker2__end() 278 279 !****************************************************************** 280 !compute V2 281 !****************************************************************** 282 V2=vresid 283 do ispden=1,nspden 284 deltaW(:,ispden)= (rdielng(:)**2) 285 end do 286 call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat) 287 call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,& 288 & frskerker2__newvres2,lastp2*real(1e-6,dp),700,V2,rdummy1,rdummy2) 289 lastp2=min(abs(rdummy1),1e-6_dp) 290 call frskerker2__end() 291 292 293 !****************************************************************** 294 !compute C1, C2 & DE 295 !****************************************************************** 296 cplex=1; 297 option=1; 298 call dotprod_vn(cplex,& !complex density/pot 299 &rdielng,& !the density 300 &DE,& !resulting dorproduct integrated over r ! here DE is used has a buffer 301 &doti,& !imaginary part of the integral 302 &size(rdielng,1),& !number of localy(cpu) attributed grid point 303 &nfftotf,& !real total number of grid point 304 &nspden,& !nspden 305 &option,& !1=compute only the real part 2=compute also the imaginary part 306 &rdielng,& !the potential 307 &ucvol,& !cell volume 308 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 309 do ispden=1,nspden 310 buffer(:,ispden)=rdielng(:)*V1(:,ispden) 311 end do 312 call dotprod_vn(cplex,& !complex density/pot 313 &rdielng,& !the density 314 &C1,& !resulting dorproduct integrated over r ! here DE is used has a buffer 315 &doti,& !imaginary part of the integral 316 &size(rdielng,1),& !number of localy(cpu) attributed grid point 317 &nfftotf,& !real total number of grid point 318 &nspden,& !nspden 319 &option,& !1=compute only the real part 2=compute also the imaginary part 320 &buffer,& !the potential 321 &ucvol,& !cell volume 322 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 323 do ispden=1,nspden 324 buffer(:,ispden)=rdielng(:)*V2(:,ispden) 325 end do 326 call dotprod_vn(cplex,& !complex density/pot 327 &rdielng,& !the density 328 &C2,& !resulting dorproduct integrated over r ! here DE is used has a buffer 329 &doti,& !imaginary part of the integral 330 &size(rdielng,1),& !number of localy(cpu) attributed grid point 331 &nfftotf,& !real total number of grid point 332 &nspden,& !nspden 333 &option,& !1=compute only the real part 2=compute also the imaginary part 334 &buffer,& !the potential 335 &ucvol,& !cell volume 336 &mpi_comm_sphgrid=mpi_enreg%comm_fft) 337 C1=C1/DE 338 C2=C2/DE 339 DE=C1/(one-C2) 340 341 !****************************************************************** 342 !compute the new preconditionned residuals 343 !****************************************************************** 344 vrespc(:,1)=diemix*(V1(:,1)+DE*V2(:,1)) 345 if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*(V1(:,2:nspden)+DE*V2(:,2:nspden)) 346 347 end subroutine prcrskerker2