TABLE OF CONTENTS


ABINIT/initaim [ Functions ]

[ Top ] [ 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