TABLE OF CONTENTS
ABINIT/initaim [ Functions ]
NAME
initaim
FUNCTION
Initialization for the 3D interpolation for the AIM code: - this procedure reads the charge density of the electrons of valence on the equidistant 3D grid (*_DEN output file of ABINIT) and the core charge density of electrons from *.fc files (fhi package) - the Cholesky decomposition of the general matrix for the computation of the 1D spline coeficients in each direction is done. Warning - the procedure is modified to use periodic boundary conditions already during the decomposition - the second derivations of valence density in three directions are computed and stored in the real space grid of the density for interpolation. - the core density is stored separately in the radial grid together with th second radial derivation
COPYRIGHT
Copyright (C) 2002-2017 ABINIT group (PCasek,FF,XG) 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/doc/developers/contributors.txt .
INPUTS
aim_dtset= the structured entity containing all input variables
OUTPUT
znucl_batom= the nuclear charge of the Bader atom (see side effects)
SIDE EFFECTS
thie routine works on the data contained in the aim_fields and aim_prom modules WARNING This file does not follow the ABINIT coding rules (yet)
PARENTS
drvaim
CHILDREN
bschg1,hdr_bcast,hdr_echo,hdr_fort_read,hdr_free,hdr_ncread,inspln lubksb,ludcmp,metric,xmpi_bcast
SOURCE
49 #if defined HAVE_CONFIG_H 50 #include "config.h" 51 #endif 52 53 #include "abi_common.h" 54 55 subroutine initaim(aim_dtset,znucl_batom) 56 57 use defs_basis 58 use defs_abitypes 59 use defs_aimfields 60 use defs_aimprom 61 use defs_parameters 62 use m_profiling_abi 63 use m_errors 64 use m_xmpi 65 #ifdef HAVE_NETCDF 66 use netcdf 67 #endif 68 use m_hdr 69 70 !This section has been created automatically by the script Abilint (TD). 71 !Do not modify the following lines by hand. 72 #undef ABI_FUNC 73 #define ABI_FUNC 'initaim' 74 use interfaces_28_numeric_noabirule 75 use interfaces_41_geometry 76 use interfaces_63_bader, except_this_one => initaim 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments ------------------------------------ 82 !scalars 83 type(aim_dataset_type),intent(in) :: aim_dtset 84 85 !Local variables ------------------------------ 86 !scalars 87 integer,parameter :: master=0 88 integer :: fform0,id,ierr,ii,info,jj,kk,kod,mm,ndtmax,nn,nsa,nsb,nsc,nsym,me,nproc,npsp 89 integer :: unth,comm 90 #ifdef HAVE_NETCDF 91 integer :: den_id 92 #endif 93 real(dp) :: ss,ucvol,znucl_batom 94 real(dp) :: zz 95 type(hdr_type) :: hdr 96 !arrays 97 integer :: ipiv(3) 98 integer,allocatable :: symrel(:,:,:) 99 real(dp) :: aa(3),bb(3),gmet(3,3),gprimd(3,3),rmet(3,3),yy(3,3) 100 real(dp),allocatable :: tnons(:,:),znucl(:),zionpsp(:) 101 real(dp),pointer :: ptc(:),ptd(:),ptf(:),ptp(:),ptsd(:) 102 103 ! ********************************************************************* 104 105 !DEBUG 106 !write(std_out,*) ' initaim : enter ' 107 !ENDDEBUG 108 109 comm = xmpi_world 110 me=xmpi_comm_rank(comm) 111 nproc=xmpi_comm_size(comm) 112 113 slc=0 ! code for follow 114 115 !The use of the "hdr" routines is much better for the future 116 !maintenance of the code. Indeed, the content of the header 117 !will continue to change from time to time, and the associated 118 !changes should be done in one place only. 119 120 !Read ABINIT header ---------------------------------------------------------- 121 if(me==master)then 122 if (aim_iomode == IO_MODE_ETSF) then 123 call hdr_ncread(hdr, untad, fform0) 124 else 125 call hdr_fort_read(hdr, untad, fform0) 126 end if 127 end if 128 ABI_CHECK(fform0 /= 0, "hdr_read returned fform == 0") 129 call hdr_bcast(hdr,master,me,comm) 130 131 !Echo part of the header 132 call hdr_echo(hdr, fform0, 4, unit=std_out) 133 call hdr_echo(hdr, fform0, 4, unit=untout) 134 135 natom=hdr%natom 136 ngfft(1:3)=hdr%ngfft(:) 137 nsym=hdr%nsym 138 npsp=hdr%npsp 139 ntypat=hdr%ntypat 140 rprimd(:,:)=hdr%rprimd(:,:) 141 142 ABI_ALLOCATE(zionpsp,(npsp)) 143 ABI_ALLOCATE(znucl,(ntypat)) 144 ABI_ALLOCATE(typat,(natom)) 145 ABI_ALLOCATE(xred,(3,natom)) 146 ABI_ALLOCATE(symrel,(3,3,nsym)) 147 ABI_ALLOCATE(tnons,(3,nsym)) 148 ABI_ALLOCATE(xatm,(3,natom)) 149 150 symrel(:,:,:)=hdr%symrel(:,:,:) 151 typat(:)=hdr%typat(:) 152 tnons(:,:)=hdr%tnons(:,:) 153 znucl(:)=hdr%znucltypat(:) 154 zionpsp(:)=hdr%zionpsp(:) 155 xred(:,:)=hdr%xred(:,:) 156 157 !This is to deallocate records of hdr 158 call hdr_free(hdr) 159 160 !------------------------------------------------------------------------------- 161 162 ABI_ALLOCATE(dvl,(ngfft(1),ngfft(2),ngfft(3))) 163 164 if(me==master)then 165 if (aim_iomode == IO_MODE_ETSF) then 166 #ifdef HAVE_NETCDF 167 ! netcdf array has shape [cplex, n1, n2, n3, nspden]), here we read only the total density. 168 NCF_CHECK(nf90_inq_varid(untad, "density", den_id)) 169 NCF_CHECK(nf90_get_var(untad, den_id, dvl, start=[1,1,1,1], count=[1, ngfft(1), ngfft(2), ngfft(3), 1])) 170 #endif 171 else 172 read(untad,iostat=nn) dvl(1:ngfft(1),1:ngfft(2),1:ngfft(3)) 173 ABI_CHECK(nn==0,"error of reading !") 174 end if 175 end if 176 call xmpi_bcast(dvl, master, comm, ierr) 177 178 write(std_out,*)ch10,' initaim : the valence density has been read' ,ch10 179 180 !INITIALISATION OF SOME IMPORTANT FIELDS 181 182 !Only interpolation is computed (inside vgh_rho) in reduced 183 !coordinates. In all other routines the cart. coordinates (CC) are used. 184 185 !transformation of the atom positions to CC 186 do ii=1,natom 187 xatm(:,ii)=xred(:,ii) 188 call bschg1(xatm(:,ii),1) 189 end do 190 191 !Generation of the neighbouring cells + transf to CC 192 nn=0 193 nsa=aim_dtset%nsa ; nsb=aim_dtset%nsb ; nsc=aim_dtset%nsc 194 do ii=-nsa,nsa 195 do jj=-nsb,nsb 196 do kk=-nsc,nsc 197 nn=nn+1 198 atp(1,nn)=ii*1._dp 199 atp(2,nn)=jj*1._dp 200 atp(3,nn)=kk*1._dp 201 call bschg1(atp(:,nn),1) 202 end do 203 end do 204 end do 205 nnpos=nn 206 207 !DEBUG 208 !write(std_out,*)' initaim : nnpos=',nnpos 209 !ENDDEBUG 210 211 batcell=nsa*(2*nsb+1)*(2*nsc+1)+(2*nsc+1)*nsb+nsc+1 212 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 213 maxatdst=min(maxatdst, nsa*sqrt(rmet(1,1)), nsb*sqrt(rmet(2,2)), nsc*sqrt(rmet(3,3)) ) 214 if (maxcpdst > maxatdst) maxcpdst=0.75*maxatdst 215 216 217 !RPRIM ITS INVERSE AND TRANSPOSE 218 219 do ii=1,3 220 do jj=1,3 221 yy(ii,jj)=rprimd(ii,jj) 222 end do 223 end do 224 call ludcmp(yy,3,3,ipiv,id,info) 225 ABI_CHECK(info==0,'Error inverting rprimd') 226 227 do ii=1,3 228 do jj=1,3 229 ivrprim(ii,jj)=0._dp 230 end do 231 ivrprim(ii,ii)=1._dp 232 end do 233 do ii=1,3 234 call lubksb(yy,3,3,ipiv,ivrprim(:,ii)) 235 end do 236 do ii=1,3 237 do jj=1,3 238 trivrp(ii,jj)=ivrprim(jj,ii) 239 end do 240 end do 241 242 write(std_out,'(" INVERSE OF RPRIMD: ",/,3F16.8,/,3F16.8,/,3F16.8,/)') & 243 & ((ivrprim(ii,jj), jj=1,3), ii=1,3) 244 write(untout,'(" INVERSE OF RPRIMD: ",/,3F16.8,/,3F16.8,/,3F16.8,/)') & 245 & ((ivrprim(ii,jj), jj=1,3), ii=1,3) 246 247 write(std_out,*) "ATOMS (index,at.number,Zionic,position(xcart.))" 248 write(std_out,*) "=======================================" 249 do ii=1,natom 250 jj=typat(ii) 251 write(std_out,'(I4,2F10.6,3F16.8)') ii, znucl(jj), zionpsp(jj), (xatm(kk,ii),kk=1,3) 252 end do 253 write(untout,*) "ATOMS (index,at.number,Zionic,position(xcart.))" 254 write(untout,*) "=======================================" 255 do ii=1,natom 256 jj=typat(ii) 257 write(untout,'(I4,2F10.6,3F16.8)') ii, znucl(jj), zionpsp(jj), (xatm(kk,ii),kk=1,3) 258 end do 259 260 !STEPS IN REAL SPACE GRID (REDUCED) 261 do ii=1,3 262 dix(ii)=1._dp/ngfft(ii) 263 end do 264 265 !READING OF THE CORE DENSITY 266 write(std_out,*)ch10,' initaim : will read the core densities' ,ch10 267 268 ABI_ALLOCATE(ndat,(ntypat)) 269 ABI_ALLOCATE(rminl,(natom)) 270 ndtmax=0 271 if(me==master)then 272 do ii=1,ntypat 273 unth=unt+ii 274 ! DEBUG 275 ! write(std_out,*)' read from unit ',unth 276 ! call flush(std_out) 277 ! stop 278 ! ENDDEBUG 279 read(unth,*) ndat(ii),ss 280 if (ndat(ii)>ndtmax) ndtmax=ndat(ii) 281 end do 282 end if 283 call xmpi_bcast(ndat,master,comm,ierr) 284 call xmpi_bcast(ndtmax,master,comm,ierr) 285 call xmpi_bcast(ss,master,comm,ierr) 286 287 !FIELDS FOR STORING CORE DENSITY 288 289 ABI_ALLOCATE(rrad,(ndtmax,ntypat)) 290 ABI_ALLOCATE(crho,(ndtmax,ntypat)) 291 ABI_ALLOCATE(sp2,(ndtmax,ntypat)) 292 ABI_ALLOCATE(sp3,(ndtmax,ntypat)) 293 ABI_ALLOCATE(sp4,(ndtmax,ntypat)) 294 ABI_ALLOCATE(corlim,(ntypat)) 295 296 sp2(:,:)=zero 297 sp3(:,:)=zero 298 sp4(:,:)=zero 299 300 !Reading of the core densities 301 corlim(:)=0 302 kod=0 303 if(me==master)then 304 do ii=1,ntypat 305 unth=unt+ii 306 do jj=1,ndat(ii) 307 read(unth,*) rrad(jj,ii),crho(jj,ii),sp2(jj,ii),sp3(jj,ii) 308 ! this is the integral of the core charge read in 309 crho(jj,ii) = crho(jj,ii)/4._dp/pi 310 if ((crho(jj,ii) < aim_rhocormin) .and. (corlim(ii)==0)) corlim(ii)=jj 311 sp2(jj,ii)=sp2(jj,ii)/4._dp/pi 312 sp3(jj,ii)=sp3(jj,ii)/4._dp/pi ! ATENTION!!! in sp3 is just second derivation 313 end do 314 do jj=1,ndat(ii)-1 315 sp4(jj,ii)=(sp3(jj+1,ii)-sp3(jj,ii))/(6._dp*(rrad(jj+1,ii)-rrad(jj,ii))) 316 end do 317 ! 318 zz = crho(1,ii) * rrad(1,ii)**2 * (rrad(2,ii)-rrad(1,ii)) 319 do jj=2,ndat(ii)-1 320 zz = zz + crho(jj,ii) * rrad(jj,ii)**2 * (rrad(jj+1,ii)-rrad(jj-1,ii)) 321 end do 322 zz = zz * half * 4._dp * pi 323 if (corlim(ii)==0) corlim(ii)=ndat(ii) 324 325 ! add check on zion wrt FHI .fc file 326 ! compare zion to zionpsp(typat(aim_dtset%batom)) 327 if (abs(znucl(ii) - zz - zionpsp(ii)) > 1.e-1_dp) then 328 write (std_out,*) 'error: your core charge ', zz, ' does not correspond to the correct number' 329 write (std_out,*) ' of valence electrons', zionpsp(ii), ' and the nuclear charge ', znucl(ii) 330 write (std_out,*) ' You have probably used a pseudopotential which has more valence electrons than the' 331 write (std_out,*) ' original FHI ones. ACTION: make a .fc file with the correct core charge' 332 stop 333 end if 334 335 end do 336 end if 337 call xmpi_bcast(rrad,master,comm,ierr) 338 call xmpi_bcast(crho,master,comm,ierr) 339 call xmpi_bcast(sp2,master,comm,ierr) 340 call xmpi_bcast(sp3,master,comm,ierr) 341 call xmpi_bcast(sp4,master,comm,ierr) 342 call xmpi_bcast(corlim,master,comm,ierr) 343 344 write(std_out,*)ch10,' initaim : the core densities have been read' ,ch10 345 346 347 !CORRECTION OF THE CORE DENSITY NORMALISATION 348 crho(:,:)=1.0003*crho(:,:) 349 sp2(:,:)=1.0003*sp2(:,:) 350 sp3(:,:)=1.0003*sp3(:,:) 351 sp4(:,:)=1.0003*sp4(:,:) 352 353 !FIELDS FOR INTERPOLATIONS OF THE VALENCE DENSITY 354 355 ABI_ALLOCATE(dig1,(ngfft(1))) 356 ABI_ALLOCATE(dig2,(ngfft(2))) 357 ABI_ALLOCATE(dig3,(ngfft(3))) 358 ABI_ALLOCATE(llg1,(ngfft(1))) 359 ABI_ALLOCATE(llg2,(ngfft(2))) 360 ABI_ALLOCATE(llg3,(ngfft(3))) 361 ABI_ALLOCATE(cdig1,(ngfft(1)-1)) 362 ABI_ALLOCATE(cdig2,(ngfft(2)-1)) 363 ABI_ALLOCATE(cdig3,(ngfft(3)-1)) 364 ABI_ALLOCATE(ddx,(ngfft(1),ngfft(2),ngfft(3))) 365 ABI_ALLOCATE(ddy,(ngfft(1),ngfft(2),ngfft(3))) 366 ABI_ALLOCATE(ddz,(ngfft(1),ngfft(2),ngfft(3))) 367 368 !DECOMPOSITION OF THE MATRIX FOR THE DETERMINATION OF COEFFICIENTS 369 !FOR CUBIC SPLINE INTERPOLATION (using the periodic boundary conditions) 370 371 !MAIN DIAGONAL (aa) AND SECONDARY DIAGONAL (bb) MATRIX ELEMENTS 372 373 nmax=ngfft(1) 374 do ii=2,3 375 if (ngfft(ii) > nmax) nmax=ngfft(ii) 376 end do 377 nullify(ptf,ptsd) 378 nullify(ptd,ptc,ptp) 379 aa(:)=2.0*dix(:)**2/3.0 380 bb(:)=dix(:)**2/6.0 381 382 do ii=1,3 383 if(ii==1) then 384 ptd=>dig1;ptc=>cdig1;ptp=>llg1 385 elseif (ii==2) then 386 ptd=>dig2;ptc=>cdig2;ptp=>llg2 387 else 388 ptd=>dig3;ptc=>cdig3;ptp=>llg3 389 end if 390 ptd(1)=sqrt(aa(ii)) 391 ptc(1)=bb(ii)/ptd(1) 392 ptp(1)=ptc(1) 393 do jj=2,ngfft(ii)-1 394 ptd(jj)=aa(ii)-ptc(jj-1)**2 395 if(ptd(jj)<zero) then 396 MSG_ERROR('Matrix is not positive definite !') 397 end if 398 ptd(jj)=sqrt(ptd(jj)) 399 if (jj==ngfft(ii)-1) then 400 ptc(jj)=(bb(ii)-ptp(jj-1)*ptc(jj-1))/ptd(jj) 401 ptp(jj)=ptc(jj) 402 exit 403 end if 404 ptc(jj)=bb(ii)/ptd(jj) 405 ptp(jj)=-ptp(jj-1)*ptc(jj-1)/ptd(jj) 406 end do 407 ss=0._dp 408 do jj=1,ngfft(ii)-1 409 ss=ss+ptp(jj)**2 410 end do 411 ss=aa(ii)-ss 412 if(ss<zero) then 413 MSG_ERROR('Matrix is not positive definite !') 414 end if 415 ptd(ngfft(ii))=sqrt(ss) 416 ptp(ngfft(ii))=ptd(ngfft(ii)) 417 418 419 ! INITIALISATION OF THE SECOND DERIVATIVE FIELDS 420 421 nn=ii+1 422 if (nn>3) nn=nn-3 423 mm=ii+2 424 if (mm>3) mm=mm-3 425 do jj=1,ngfft(nn) 426 do kk=1,ngfft(mm) 427 ! The calcul of the second derivations on the grid 428 call inspln(ii,jj,kk) 429 end do 430 end do 431 nullify(ptd,ptc,ptp) 432 end do 433 nullify(ptd,ptc,ptp) 434 435 znucl_batom=znucl(typat(aim_dtset%batom)) 436 437 ABI_DEALLOCATE(znucl) 438 ABI_DEALLOCATE(zionpsp) 439 ABI_DEALLOCATE(symrel) 440 ABI_DEALLOCATE(tnons) 441 442 !the pointers are obsolete - to remove later 443 444 end subroutine initaim