TABLE OF CONTENTS
ABINIT/mkcore_paw [ Functions ]
NAME
mkcore_paw
FUNCTION
FIXME: add description.
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (TRangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
mpi_enreg=informations about MPI parallelization
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
ind_positions_,mkcore_inner,mkgrid_fft,pawrad_free,pawrad_init ptabs_fourdp,strconv,timab,wrtout,xcart2xred,xmpi_sum,xred2xcart
SOURCE
33 #if defined HAVE_CONFIG_H 34 #include "config.h" 35 #endif 36 37 #include "abi_common.h" 38 39 subroutine mkcore_paw(atindx1,corstr,dyfrx2,grxc,icoulomb,natom,mpi_enreg,& 40 & nattyp,nfft,ngfft,nspden,ntypat,n3xccc,option,pawrad,pawtab,psppar,rprimd,& 41 & ucvol,vxc,xccc3d,xred) 42 43 use defs_basis 44 use defs_abitypes 45 use m_profiling_abi 46 use m_errors 47 use m_xmpi 48 49 use m_geometry, only : xcart2xred, xred2xcart 50 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 51 use m_pawtab, only : pawtab_type 52 use m_mpinfo, only : ptabs_fourdp 53 54 !This section has been created automatically by the script Abilint (TD). 55 !Do not modify the following lines by hand. 56 #undef ABI_FUNC 57 #define ABI_FUNC 'mkcore_paw' 58 use interfaces_14_hidewrite 59 use interfaces_18_timing 60 use interfaces_41_geometry 61 use interfaces_67_common, except_this_one => mkcore_paw 62 !End of the abilint section 63 64 implicit none 65 66 !Arguments ------------------------------------ 67 !scalars 68 integer,intent(in) :: icoulomb 69 integer,intent(in) :: natom,ntypat,nfft,nspden,n3xccc,option 70 type(MPI_type),intent(in) :: mpi_enreg 71 !arrays 72 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 73 integer,intent(in) :: ngfft(18) 74 real(dp),intent(in) :: psppar(0:4,0:6,ntypat) 75 real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden) 76 real(dp),intent(in) :: ucvol 77 real(dp),intent(in) :: xred(3,natom) 78 real(dp),intent(inout) :: xccc3d(nfft) 79 real(dp),intent(out) :: corstr(6),dyfrx2(3,3,natom),grxc(3,natom) 80 type(pawtab_type),intent(in) :: pawtab(ntypat) 81 type(pawrad_type),intent(in) :: pawrad(ntypat) 82 83 !Local variables------------------------------- 84 integer :: iat,iatm,iatom,iatom_tot 85 integer :: ier,iex,iey,iez,ind 86 integer :: isx,isy,isz,itypat 87 integer :: i1,i2,i3,i3loc 88 integer :: j1,j2,j3,msz 89 integer :: ncmax,nfgd,nfgd_r0 90 integer :: me,me_fft,nfftot,nproc,nproc_fft,nu 91 integer :: n1,n2,n3,n3d 92 real(dp) :: cutoff,factor,grxc1,grxc2,grxc3 93 real(dp) :: rloc,rr2,rshp,rshpm1,rx,ry,rz 94 real(dp) :: r2shp,strdia 95 character(len=500) :: message 96 character(len=1) :: geocode 97 logical :: perx,pery,perz,gox,goy,goz 98 type(pawrad_type)::core_mesh 99 !arrays 100 integer,allocatable :: ifftsph_tmp(:) 101 real(dp) :: corfra(3,3) 102 real(dp) :: hh(3) !fine grid spacing 103 real(dp) :: rmet(3,3),tsec(2),xcart(3,natom) 104 real(dp),allocatable :: gridcart(:,:),rr(:),rred(:,:) 105 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 106 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 107 108 ! ************************************************************************* 109 110 DBG_ENTER("COLL") 111 112 if(nspden >1) then 113 write(message, '(a)')'mkcore_paw: this is not yet generalized to npsden>1' 114 MSG_ERROR(message) 115 end if 116 117 geocode='P' 118 if (icoulomb==1) geocode='F' 119 if (icoulomb==2) geocode='S' 120 121 !Compute metric tensor in real space rmet 122 do nu=1,3 123 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+& 124 & rprimd(3,:)*rprimd(3,nu) 125 end do 126 127 !MPI 128 nproc =xmpi_comm_size(mpi_enreg%comm_fft); nproc_fft= ngfft(10) 129 me =xmpi_comm_rank(mpi_enreg%comm_fft); me_fft= ngfft(11) 130 131 if(me /= me_fft .or. nproc /= nproc_fft) then 132 MSG_BUG("mkcore_paw: comm_size or comm_rank not equal to the corresponding values in ngfft") 133 end if 134 135 n1 = ngfft(1) 136 n2 = ngfft(2) 137 n3 = ngfft(3) 138 n3d = ngfft(13) 139 if(nproc==1) n3d=n3 140 141 !Get the distrib associated with this fft_grid 142 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 143 144 !Store xcart for each atom 145 call xred2xcart(natom, rprimd, xcart, xred) 146 147 !Store cartesian coordinates for each grid points 148 ABI_ALLOCATE(gridcart,(3, nfft)) 149 call mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd) 150 151 !definition of the grid spacings 152 hh(1) = rprimd(1,1)/(ngfft(1)) 153 hh(2) = rprimd(2,2)/(ngfft(2)) 154 hh(3) = rprimd(3,3)/(ngfft(3)) 155 156 if(nfft .ne. n3xccc)then 157 write(message,'(a,a,a,a,a,a,2i6)') ch10,& 158 & ' mkcore_paw: BUG -',ch10,& 159 & ' nfft and n3xccc should be equal,',ch10,& 160 & ' however, nfft and n3xccc=',nfft,n3xccc 161 MSG_BUG(message) 162 end if 163 164 if (option==1) then 165 ! Zero out array to permit accumulation over atom types below: 166 xccc3d(:)=zero 167 else if (option==2) then 168 ! Zero out gradient of Exc array 169 grxc(:,:)=zero 170 else if (option==3) then 171 ! Zero out locally defined stress array 172 corfra(:,:)=zero 173 strdia=zero 174 else if (option==4) then 175 ! Zero out fr-wf part of the dynamical matrix 176 dyfrx2(:,:,:)=zero 177 else 178 write(message, '(a,a,a,a)' ) ch10,& 179 & ' mkcore_paw: BUG -',ch10,& 180 & ' Can''t be here ! (bad option).' 181 MSG_BUG(message) 182 end if 183 184 write(message,'(a,a)') ch10,& 185 & ' mkcore_paw: Compute core density' 186 call wrtout(std_out,message,'COLL') 187 188 !conditions for periodicity in the three directions 189 perx=(geocode /= 'F') 190 pery=(geocode == 'P') 191 perz=(geocode /= 'F') 192 193 iatm=0 194 !Big loop on atom types 195 do itypat=1,ntypat 196 197 rloc=psppar(0,0,itypat) 198 cutoff=10.d0*rloc 199 200 ! Set radius size: 201 rshp=pawtab(itypat)%rcore 202 r2shp=1.0000001_dp*rshp**2 203 rshpm1=one/rshp 204 205 ! allocate arrays 206 ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3) 207 ! 1+int(1.1* factors are included just for cautioness 208 ncmax=1+int(1.1d0*((rshp/hh(1))*(rshp/hh(2))*pi)) 209 210 ABI_ALLOCATE(ifftsph_tmp,(ncmax)) 211 ABI_ALLOCATE(rr,(ncmax)) 212 if(option>1) then 213 ABI_ALLOCATE(rred,(3,ncmax)) 214 else 215 ABI_ALLOCATE(rred,(0,0)) 216 end if 217 218 ! Create mesh_core object 219 ! since core_mesh_size can be bigger than pawrad%mesh_size, 220 msz=pawtab(itypat)%core_mesh_size 221 call pawrad_init(core_mesh,mesh_size=msz,mesh_type=pawrad(itypat)%mesh_type,& 222 & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep) 223 224 ! Big loop on atoms 225 do iat=1,nattyp(itypat) 226 iatm=iatm+1;iatom=atindx1(iatm) 227 iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom) 228 229 if(option==2) then 230 grxc1=zero 231 grxc2=zero 232 grxc3=zero 233 end if 234 235 ! Define a "box" around each atom 236 rx=xcart(1,iatom) 237 ry=xcart(2,iatom) 238 rz=xcart(3,iatom) 239 240 isx=floor((rx-cutoff)/hh(1)) 241 isy=floor((ry-cutoff)/hh(2)) 242 isz=floor((rz-cutoff)/hh(3)) 243 iex=ceiling((rx+cutoff)/hh(1)) 244 iey=ceiling((ry+cutoff)/hh(2)) 245 iez=ceiling((rz+cutoff)/hh(3)) 246 247 do i3=isz,iez 248 call ind_positions_(perz,i3,n3,j3,goz) 249 250 if(fftn3_distrib(j3)==me_fft) then 251 i3loc=ffti3_local(j3) 252 253 ! Initialize counters 254 nfgd=0 255 nfgd_r0=0 256 257 do i2=isy,iey 258 call ind_positions_(pery,i2,n2,j2,goy) 259 do i1=isx,iex 260 call ind_positions_(perx,i1,n1,j1,gox) 261 ! r2=x**2+y**2+z**2 262 if (goz .and. goy .and. gox) then 263 ind=j1+(j2-1)*n1+(i3loc-1)*n1*n2 264 rr2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2 265 266 if(rr2<=r2shp) then 267 if(rr2>tol5) then 268 nfgd=nfgd+1 269 rr(nfgd)=sqrt(rr2) 270 ifftsph_tmp(nfgd)=ind 271 if(option>1) then 272 call xcart2xred(1,rprimd,gridcart(:,ind),rred(:,nfgd)) 273 end if 274 elseif (option==4) then 275 ! We save r=0 vectors only for option==4: 276 ! for other options this is ignored 277 278 ! We reuse the same variable "ifftshp_tmp", 279 ! but we start from the higher index 280 nfgd_r0=nfgd_r0+1 281 ifftsph_tmp(ncmax-nfgd_r0+1)=ind 282 283 end if !rr2>tol5 284 end if !rr2<r2shp 285 end if !gox.. 286 end do !i1 287 end do !i2 288 289 ! All of the following could be done inside or outside the loops (i2,i1,i3) 290 ! Outside the loops: the memory consuption increases. 291 ! Inside the inner loop: the time of calculation increases. 292 ! Here, I choose to do it here, somewhere in the middle. 293 if (option/=4.and.nfgd==0) cycle 294 if (option==4.and.nfgd==0.and.nfgd_r0==0) cycle 295 call mkcore_inner(corfra,core_mesh,dyfrx2,& 296 & grxc1,grxc2,grxc3,ifftsph_tmp,msz,& 297 & natom,ncmax,nfft,nfgd,nfgd_r0,nspden,n3xccc,option,pawtab(itypat),& 298 & rmet,rr,strdia,vxc,xccc3d,rred=rred) 299 300 end if !parallel fftn3 301 end do !i3 302 303 if(option==2) then 304 nfftot=product(ngfft(1:3)) 305 factor=(ucvol/real(nfftot,dp)) !/rshp 306 grxc(1,iatom)=grxc1*factor 307 grxc(2,iatom)=grxc2*factor 308 grxc(3,iatom)=grxc3*factor 309 310 if(nproc_fft > 1) then 311 call timab(539,1,tsec) 312 call xmpi_sum(grxc1,mpi_enreg%comm_fft,ier) 313 call xmpi_sum(grxc2,mpi_enreg%comm_fft,ier) 314 call xmpi_sum(grxc3,mpi_enreg%comm_fft,ier) 315 call timab(539,2,tsec) 316 end if 317 318 end if 319 end do !iatom 320 321 ! Deallocate 322 call pawrad_free(core_mesh) 323 ABI_DEALLOCATE(ifftsph_tmp) 324 ABI_DEALLOCATE(rr) 325 ABI_DEALLOCATE(rred) 326 327 end do !itypat 328 329 if (option==2) then 330 ! Apply rmet as needed to get reduced coordinate gradients 331 ! do iatom=1,natom 332 ! t1=grxc(1,iatom) 333 ! t2=grxc(2,iatom) 334 ! t3=grxc(3,iatom) 335 ! grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3 336 !! grxc(:,iatom)=rprimd(1,:)*t1+rprimd(2,:)*t2+rprimd(3,:)*t3 337 ! end do 338 339 else if (option==3) then 340 341 ! Transform stress tensor from full storage mode to symmetric storage mode 342 corstr(1)=corfra(1,1) 343 corstr(2)=corfra(2,2) 344 corstr(3)=corfra(3,3) 345 corstr(4)=corfra(3,2) 346 corstr(5)=corfra(3,1) 347 corstr(6)=corfra(2,1) 348 349 ! Transform stress tensor from reduced coordinates to cartesian coordinates 350 call strconv(corstr,rprimd,corstr) 351 352 ! Compute diagonal contribution to stress tensor (need input xccc3d) 353 ! strdia = (1/N) Sum(r) [mu_xc_avg(r) * rho_core(r)] 354 strdia=zero 355 do i3=1,n3d 356 do i2=1,n2 357 do i1=1,n1 358 ind=i1+(i2-1)*n1+(i3-1)*n1*n2 359 strdia=strdia+vxc(ind,1)*xccc3d(ind) 360 end do 361 end do 362 end do 363 strdia=strdia/real(nfft,dp) 364 ! Add diagonal term to stress tensor 365 corstr(1)=corstr(1)+strdia 366 corstr(2)=corstr(2)+strdia 367 corstr(3)=corstr(3)+strdia 368 end if 369 370 if(nproc_fft > 1) then 371 call timab(539,1,tsec) 372 if(option==3) then 373 call xmpi_sum(corstr,mpi_enreg%comm_fft,ier) 374 end if 375 if(option==2) then 376 call xmpi_sum(grxc,mpi_enreg%comm_fft,ier) 377 end if 378 call timab(539,2,tsec) 379 end if 380 381 DBG_EXIT("COLL") 382 383 end subroutine mkcore_paw