TABLE OF CONTENTS


m_mkrho/atomden [ Functions ]

[ Top ] [ m_mkrho ] [ Functions ]

NAME

 atomden

FUNCTION

 Construct atomic proto-bulk density (i.e. the superposed density
 from neutral, isolated atoms at the bulk atomic positions).
 This is useful if one wants to construct the bonding density:

 rho^{bnd} = rho^{bulk}(r)
                 - \sum_{\alpha}\rho^{atm}_{\alpha}(r-R_{\alpha})

 Where rho^{bulk} is the bulk density, rho^{atm} the atomic density
 and the index \alpha sums over all atoms. the R_{\alpha} are the
 atomic positions in the bulk. This routine calculates the sum over
 rho^{atm}_{\alpha}(r-R_{\alpha}) on a grid.

 Units are atomic.

COPYRIGHT

 Copyright (C) 2005-2022 ABINIT group (SM,VR,FJ,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 .

INPUTS

 calctype : type of calculation
          'replace' zero the input/output density array
          'add'     add to the input/output density array
 natom : number of atoms in cell
 ntypat : number of different types of atoms in cell
 typat(natom) : type of each atom
 ngrid : number of gridpoints
 r_vec_grid(3,ngrid) : real (non-reduced) coordinates for grid points
 rho(ngrid) : input/output density array
 a(3),b(3),c(3) : real-space basis vectors
 atom_pos(3,natom) : reduced coordinates for atomic positions
 natomgr(ntypat) : number of gridpoints for each atomic density grid
 natomgrmax : max(natomgr(ntypat))
 atomrgrid(natomgrmax,ntypat)
 density(natomgrmax,ntypat)

OUTPUT

 rho(ngrid) : input/output density array

SIDE EFFECTS

NOTES

 There are two ways to compile the proto density in real space
 for a solid. One alternative is that the density is calculated
 for an extended grid encompassing the sphere of points around
 one atom, and the results are folded back into the unit cell.
 On the other hand one can, around each grid point, identify the
 number of atoms in a sphere equivalent to the length of the radial
 grid for each type of atom.
 The second approach, with some modification, is taken here. The
 numer of atoms in a supercell cell are listed such that the supercell
 encompasses the atoms which could contribute to any point in the grid.
 That list is kept and cycled through, to avoid recalculating it at
 each point.
 Note that the density calculated from the atom is the spherical
 average, since there is no preferred direction without any
 external field (and it's simpler)

SOURCE

2124 subroutine atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,atom_pos, &
2125 &                  natomgr,natomgrmax,atomrgrid,density,prtvol,calctype)
2126 
2127 !Arguments ------------------------------------
2128 !scalars
2129  integer,intent(in) :: natom,ntypat,ngrid,natomgrmax,prtvol
2130  character(len=7),intent(in) :: calctype
2131 !arrays
2132  type(MPI_type),intent(in) :: MPI_enreg
2133  integer,intent(in) :: typat(natom),natomgr(ntypat)
2134  real(dp),intent(in) :: r_vec_grid(3,ngrid),a(3),b(3),c(3)
2135  real(dp),intent(in) :: atom_pos(3,natom),atomrgrid(natomgrmax,ntypat)
2136  real(dp),intent(in) :: density(natomgrmax,ntypat)
2137  real(dp),intent(inout) :: rho(ngrid)
2138 
2139 !Local variables-------------------------------
2140 !scalars
2141  character(len=500) :: message
2142  integer :: cnt,delta,i,l,m,n,iatom,itypat,igrid,ncells,n_grid_p
2143  integer :: ierr,spaceComm,nprocs,master,rank,remainder
2144  real(dp) :: a_norm,b_norm,c_norm
2145  real(dp) :: r_max,R_sphere_max,dp_dummy,ybcbeg,ybcend
2146 !arrays
2147  integer :: n_equiv_atoms(ntypat),grid_index(ngrid)
2148  integer :: my_start_equiv_atoms(ntypat)
2149  integer :: my_end_equiv_atoms(ntypat)
2150  integer :: l_min(ntypat),m_min(ntypat),n_min(ntypat)
2151  integer :: l_max(ntypat),m_max(ntypat),n_max(ntypat)
2152  real(dp) :: center(3),dp_vec_dummy(3),delta_a(3),delta_b(3),delta_c(3)
2153  real(dp) :: r_atom(3),grid_distances(ngrid)
2154  integer, allocatable :: new_index(:),i_1d_dummy(:)
2155  real(dp),allocatable :: equiv_atom_dist(:,:),equiv_atom_pos(:,:,:),rho_temp(:,:)
2156  real(dp),allocatable :: dp_1d_dummy(:),dp_2d_dummy(:,:),ypp(:)
2157  real(dp),allocatable :: x_fit(:),y_fit(:)
2158 
2159 
2160 ! ************************************************************************
2161 
2162 !initialise and check parallel execution
2163  spaceComm=MPI_enreg%comm_cell
2164  nprocs=xmpi_comm_size(spaceComm)
2165  rank=MPI_enreg%me_kpt
2166 
2167  master=0
2168 
2169 !initialise variables and vectors
2170  a_norm = sqrt(dot_product(a,a))
2171  b_norm = sqrt(dot_product(b,b))
2172  c_norm = sqrt(dot_product(c,c))
2173  center = (a+b+c)*half
2174  dp_dummy = dot_product(a,b)/(b_norm*b_norm)
2175  dp_vec_dummy = dp_dummy*b
2176  delta_a = a - dp_vec_dummy
2177  dp_dummy = dot_product(b,a)/(a_norm*a_norm)
2178  dp_vec_dummy = dp_dummy*a
2179  delta_b = b - dp_vec_dummy
2180  dp_dummy = dot_product(c,(a+b))/(dot_product((a+b),(a+b)))
2181  dp_vec_dummy = dp_dummy*(a+b)
2182  delta_c = c - dp_vec_dummy
2183  ABI_MALLOC(rho_temp,(ngrid,ntypat))
2184  rho_temp = zero
2185 
2186 !write(std_out,*) '*** --- In atomden --- ***'
2187 !write(std_out,*) ' a_norm:',a_norm,' b_norm:',b_norm,' c_norm:',c_norm
2188 !write(std_out,*) 'delta_a:',delta_a,'delta_b:',delta_b,'delta_c:',delta_c
2189 !write(std_out,*) ' center:',center
2190 
2191 !Find supercell which will contain all possible contributions
2192 !for all atoms, and enumerate positions for all atoms
2193 !TODO list of atoms can be "pruned", i.e identify all atoms
2194 !that can't possibly contribute and remove from list.
2195 !Should be most important for very oblique cells
2196  do itypat=1,ntypat
2197    R_sphere_max = atomrgrid(natomgr(itypat),itypat)
2198    l_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_a,delta_a)))
2199    l_max(itypat) = -l_min(itypat)
2200    m_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_b,delta_b)))
2201    m_max(itypat) = -m_min(itypat)
2202    n_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_c,delta_c)))
2203    n_max(itypat) = -n_min(itypat)
2204    ncells = (l_max(itypat)-l_min(itypat)+1) &
2205 &   *(m_max(itypat)-m_min(itypat)+1) &
2206 &   *(n_max(itypat)-n_min(itypat)+1)
2207    n_equiv_atoms(itypat) = 0
2208    do iatom=1,natom
2209      if (typat(iatom)==itypat) then
2210        n_equiv_atoms(itypat) = n_equiv_atoms(itypat) + ncells
2211      end if ! if type=itypat
2212    end do ! number of atoms per cell
2213    if ((rank==master).and.(prtvol>9)) then
2214      write(message,'(a)') '*** --- In atomden --- find box ***'
2215      call wrtout(std_out,message,'COLL')
2216      write(message,'(a,I4)') ' itypat:',itypat
2217      call wrtout(std_out,message,'COLL')
2218      write(message,'(2(a,I4))') ' l_min:',l_min(itypat),' l_max:',l_max(itypat)
2219      call wrtout(std_out,message,'COLL')
2220      write(message,'(2(a,I4))') ' m_min:',m_min(itypat),' m_max:',m_max(itypat)
2221      call wrtout(std_out,message,'COLL')
2222      write(message,'(2(a,I4))') ' n_min:',n_min(itypat),' n_max:',n_max(itypat)
2223      call wrtout(std_out,message,'COLL')
2224      write(message,'(2(a,I4))') ' n_equiv_atoms:',n_equiv_atoms(itypat)
2225      call wrtout(std_out,message,'COLL')
2226    end if
2227  end do !atom type
2228 
2229 !allocate arrays
2230  n = maxval(n_equiv_atoms)
2231  ABI_MALLOC(equiv_atom_pos,(3,n,ntypat))
2232  ABI_MALLOC(equiv_atom_dist,(n,ntypat))
2233  equiv_atom_pos = zero
2234  equiv_atom_dist = zero
2235 
2236 !Find positions and distance of atoms from center of cell
2237  do itypat=1,ntypat
2238    i = 1
2239    do l=l_min(itypat),l_max(itypat)
2240      do m=m_min(itypat),m_max(itypat)
2241        do n=n_min(itypat),n_max(itypat)
2242          do iatom=1,natom
2243            if (typat(iatom)==itypat) then
2244              if (i>n_equiv_atoms(itypat)) then
2245                ABI_ERROR('atomden: i>n_equiv_atoms')
2246              end if
2247              equiv_atom_pos(:,i,itypat) = (atom_pos(1,iatom)+dble(l))*a &
2248 &             + (atom_pos(2,iatom)+dble(m))*b &
2249 &             + (atom_pos(3,iatom)+dble(n))*c
2250              dp_vec_dummy = equiv_atom_pos(:,i,itypat)-center
2251              equiv_atom_dist(i,itypat) = &
2252 &             sqrt(dot_product(dp_vec_dummy,dp_vec_dummy))
2253              i = i + 1
2254            end if
2255          end do
2256        end do !n
2257      end do !m
2258    end do !l
2259 !  write(std_out,*) '*** --- In atomden --- find equiv ***'
2260 !  write(std_out,*) ' itypat:',itypat
2261 !  write(std_out,*) ' equiv_atom_pos:'
2262 !  write(std_out,*) equiv_atom_pos(:,:,itypat)
2263 !  write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat)
2264  end do !atom type
2265 
2266 !Sort the atoms after distance so that the density from the ones
2267 !furthest away can be added first. This is to prevent truncation error.
2268  do itypat=1,ntypat
2269    n = n_equiv_atoms(itypat)
2270    ABI_MALLOC(dp_1d_dummy,(n))
2271    ABI_MALLOC(new_index,(n))
2272    ABI_MALLOC(dp_2d_dummy,(3,n))
2273    dp_1d_dummy = equiv_atom_dist(1:n,itypat)
2274    dp_2d_dummy = equiv_atom_pos(1:3,1:n,itypat)
2275    do i=1,n
2276      new_index(i) = i
2277    end do
2278    call sort_dp(n,dp_1d_dummy,new_index,tol14)
2279    do i=1,n
2280 !    write(std_out,*) i,' -> ',new_index(i)
2281      equiv_atom_pos(1:3,n+1-i,itypat) = dp_2d_dummy(1:3,new_index(i))
2282      equiv_atom_dist(1:n,itypat) = dp_1d_dummy
2283    end do
2284    ABI_FREE(dp_1d_dummy)
2285    ABI_FREE(new_index)
2286    ABI_FREE(dp_2d_dummy)
2287 !  write(std_out,*) '*** --- In atomden ---  sorting atoms ***'
2288 !  write(std_out,*) ' itypat:',itypat
2289 !  write(std_out,*) ' equiv_atom_pos:'
2290 !  write(std_out,*) equiv_atom_pos(:,:,itypat)
2291 !  write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat)
2292  end do ! atom type
2293 
2294 !Divide the work in case of parallel execution
2295  if (nprocs==1) then ! Make sure everything runs with one proc
2296    if (prtvol>9) then
2297      write(message,'(a)') '  In atomden - number of processors:     1'
2298      call wrtout(std_out,message,'COLL')
2299      write(message,'(a)') '  Calculation of proto-atomic density done in serial'
2300      call wrtout(std_out,message,'COLL')
2301    end if
2302    do itypat=1,ntypat
2303      if (prtvol>9) then
2304        write(message,'(a,I6)') '  Number of equivalent atoms:',n_equiv_atoms(itypat)
2305        call wrtout(std_out,message,'COLL')
2306      end if
2307      my_start_equiv_atoms(itypat) = 1
2308      my_end_equiv_atoms(itypat) = n_equiv_atoms(itypat)
2309    end do
2310  else
2311    if (rank==master.and.prtvol>9) then
2312      write(message,'(a,I5)') '  In atomden - number of processors:',nprocs
2313      call wrtout(std_out,message,'COLL')
2314      write(message,'(a)') '  Calculation of proto-atomic density done in parallel'
2315      call wrtout(std_out,message,'COLL')
2316    end if
2317    do itypat=1,ntypat
2318      if (rank==master.and.prtvol>9) then
2319        write(message,'(a,I6)') '  Number of equivalent atoms:',n_equiv_atoms(itypat)
2320        call wrtout(std_out,message,'COLL')
2321      end if
2322 !    Divide the atoms among the processors by shuffling indices
2323      delta = int(floor(real(n_equiv_atoms(itypat))/real(nprocs)))
2324      remainder = n_equiv_atoms(itypat)-nprocs*delta
2325      my_start_equiv_atoms(itypat) = 1+rank*delta
2326      my_end_equiv_atoms(itypat) = (rank+1)*delta
2327 !    Divide the remainder points among the processors
2328 !    by shuffling indices
2329      if ((rank+1)>remainder) then
2330        my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + remainder
2331        my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + remainder
2332      else
2333        my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + rank
2334        my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + rank + 1
2335      end if
2336      if (prtvol>9) then
2337        write(message,'(a,I3)') '          For atom type: ',itypat
2338        call wrtout(std_out,message,'PERS')
2339 !      write(message,'(a,I6)') '  I''ll take atoms from: ',my_start_equiv_atoms(itypat)
2340 !      call wrtout(std_out,message,'PERS')
2341 !      write(message,'(a,I6)') '           total for me: ',my_end_equiv_atoms(itypat)
2342 !      call wrtout(std_out,message,'PERS')
2343        write(message,'(a,I6)') '            total for me: ', &
2344 &       my_end_equiv_atoms(itypat)+1-my_start_equiv_atoms(itypat)
2345        call wrtout(std_out,message,'PERS')
2346      end if
2347    end do
2348  end if
2349 
2350 !Loop over types of atoms and equivalent atoms and
2351 !interpolate density onto grid
2352  do itypat=1,ntypat
2353 !  do iatom=my_start_equiv_atoms(itypat),my_end_equiv_atoms(itypat)
2354 
2355    cnt = 0
2356    iatom = rank+1 - nprocs
2357 !  Parallel execution of loop
2358    do
2359      cnt = cnt + 1
2360      iatom = iatom + nprocs
2361      if (iatom>n_equiv_atoms(itypat)) exit ! Exit if index is too large
2362 
2363      if (mod(cnt,100)==0.and.prtvol>0) then
2364        write(message,'(2(a,I6))') ' atoms so far',cnt,' of: ',n_equiv_atoms(itypat)/nprocs
2365        call wrtout(std_out,message,'PERS')
2366      end if
2367 
2368      r_max = atomrgrid(natomgr(itypat),itypat)
2369      r_atom = equiv_atom_pos(:,iatom,itypat)
2370 
2371 !    Set up an array with the gridpoint distances
2372      i = 1
2373      grid_distances = zero
2374      grid_index = 0
2375      do igrid=1,ngrid
2376        dp_vec_dummy(:) = r_vec_grid(:,igrid) - r_atom(:)
2377        dp_dummy = sqrt(dot_product(dp_vec_dummy,dp_vec_dummy))
2378        if (dp_dummy <= r_max) then
2379          grid_distances(i) = dp_dummy
2380          grid_index(i) = igrid
2381          i = i + 1
2382        else
2383          cycle ! cycle if point is too far away
2384        end if
2385      end do
2386      n_grid_p = i - 1
2387 
2388      if (n_grid_p==0) cycle ! Cycle if no point needs
2389 !    to be interpolated
2390 
2391 !    Sort points to be interpolated in ascending order
2392      ABI_MALLOC(dp_1d_dummy,(n_grid_p))
2393      ABI_MALLOC(new_index,(n_grid_p))
2394      ABI_MALLOC(i_1d_dummy,(n_grid_p))
2395      dp_1d_dummy = grid_distances(1:n_grid_p)
2396      do i=1,n_grid_p
2397        new_index(i) = i
2398      end do
2399      call sort_dp(n_grid_p,dp_1d_dummy,new_index,tol16)
2400      grid_distances(1:n_grid_p) = dp_1d_dummy
2401      i_1d_dummy = grid_index(1:n_grid_p)
2402      do i=1,n_grid_p
2403 !      write(std_out,*) i_1d_dummy(i),' -> ',i_1d_dummy(new_index(i))
2404        grid_index(i) = i_1d_dummy(new_index(i))
2405      end do
2406      ABI_FREE(dp_1d_dummy)
2407      ABI_FREE(new_index)
2408      ABI_FREE(i_1d_dummy)
2409 
2410 !    Interpolate density onto all grid points
2411      ABI_MALLOC(ypp,(natomgr(itypat)))
2412      ABI_MALLOC(x_fit,(n_grid_p))
2413      ABI_MALLOC(y_fit,(n_grid_p))
2414      ypp = zero; y_fit = zero
2415      ybcbeg = zero; ybcend = zero
2416      x_fit = grid_distances(1:n_grid_p)
2417      call spline(atomrgrid(1:natomgr(itypat),itypat), &
2418 &     density(1:natomgr(itypat),itypat), &
2419 &     natomgr(itypat),ybcbeg,ybcend,ypp)
2420      call splint(natomgr(itypat),atomrgrid(1:natomgr(itypat),itypat), &
2421 &     density(1:natomgr(itypat),itypat),ypp,n_grid_p, &
2422 &     x_fit,y_fit)
2423 
2424 !    Save the interpolated points to grid
2425      do i=1,n_grid_p
2426        rho_temp(grid_index(i),itypat) = rho_temp(grid_index(i),itypat) + y_fit(i)
2427      end do
2428      ABI_FREE(ypp)
2429      ABI_FREE(x_fit)
2430      ABI_FREE(y_fit)
2431 
2432    end do ! n equiv atoms
2433  end do ! type of atom
2434 
2435 !Collect all contributions to rho_temp if
2436 !we are running in parallel
2437  if (nprocs>1) then
2438    call xmpi_barrier(spaceComm)
2439    call xmpi_sum_master(rho_temp,master,spaceComm,ierr)
2440    call xmpi_barrier(spaceComm)
2441    if (prtvol>9) then
2442      write(message,'(a)') '  In atomden - contributions to rho_temp collected'
2443      call wrtout(std_out,message,'PERS')
2444    end if
2445  end if
2446 
2447 !Now rho_temp contains the atomic protodensity for each atom.
2448 !Check whether this is to replace or be added to the input/output array
2449 !and sum up contributions
2450  if (trim(calctype)=='replace') rho = zero
2451  do itypat=1,ntypat
2452    rho(:) = rho(:) + rho_temp(:,itypat)
2453  end do
2454 
2455 !deallocations
2456  if (allocated(rho_temp))  then
2457    ABI_FREE(rho_temp)
2458  end if
2459  if (allocated(equiv_atom_pos))  then
2460    ABI_FREE(equiv_atom_pos)
2461  end if
2462  if (allocated(equiv_atom_dist))  then
2463    ABI_FREE(equiv_atom_dist)
2464  end if
2465 !if (allocated()) deallocate()
2466 
2467  return
2468 
2469  end subroutine atomden

m_mkrho/initro [ Functions ]

[ Top ] [ m_mkrho ] [ Functions ]

NAME

 initro

FUNCTION

 Initialize the density using either:
  - a gaussian of adjustable decay length (norm-conserving psp)
  - PS atomic valence density from psp file (PAW or NC psps with valence change in the pp file)

INPUTS

 atindx(natom)=index table for atoms (see gstate.f)
 densty(ntypat,4)=parameters for initialisation of the density of each atom type
 gmet(3,3)=reciprocal space metric (Bohr**-2)
 gsqcut=cutoff G**2 for included G s in fft box (larger sphere).
 izero=if 1, unbalanced components of rho(g) have to be set to zero
 mgfft=maximum size of 1D FFTs
 mpi_enreg=information about mpi parallelization
 mqgrid=number of grid pts in q array for n^AT(q) spline.
 natom=number of atoms in cell.
 nattyp(ntypat)=number of atoms of each type in cell.
 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
 ntypat=number of types of atoms in cell.
 nspden=number of spin-density components
 psps<type(pseudopotential_type)>=variables related to pseudopotentials
 pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
 ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase information for given atom coordinates.
 qgrid(mqgrid)=q grid for spline atomic valence density n^AT(q) from 0 to qmax.
 spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
 ucvol=unit cell volume (Bohr**3).
 usepaw= 0 for non paw calculation; =1 for paw calculation
 zion(ntypat)=charge on each type of atom (real number)
 znucl(ntypat)=atomic number, for each type of atom

OUTPUT

 rhog(2,nfft)=initialized total density in reciprocal space
 rhor(nfft,nspden)=initialized total density in real space.
         as well as spin-up part if spin-polarized

SOURCE

 828 subroutine initro(atindx,densty,gmet,gsqcut,izero,mgfft,mpi_enreg,mqgrid,natom,nattyp,&
 829 &  nfft,ngfft,nspden,ntypat,psps,pawtab,ph1d,qgrid,rhog,rhor,spinat,ucvol,usepaw,zion,znucl)
 830 
 831 !Arguments ------------------------------------
 832 !scalars
 833  integer,intent(in) :: izero,mgfft,mqgrid,natom,nfft,nspden,ntypat
 834  integer,intent(in) :: usepaw
 835  real(dp),intent(in) :: gsqcut,ucvol
 836  type(mpi_type),intent(in) :: mpi_enreg
 837  type(pseudopotential_type),intent(in) :: psps
 838 !arrays
 839  integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
 840  real(dp),intent(in) :: densty(ntypat,4),gmet(3,3),ph1d(2,3*(2*mgfft+1)*natom)
 841  real(dp),intent(in) :: qgrid(mqgrid),spinat(3,natom),zion(ntypat)
 842  real(dp),intent(in) :: znucl(ntypat)
 843  real(dp),intent(out) :: rhog(2,nfft),rhor(nfft,nspden)
 844  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
 845 
 846 !Local variables-------------------------------
 847 !The decay lengths should be optimized element by element, and even pseudopotential by pseudopotential.
 848 !scalars
 849  integer,parameter :: im=2,re=1
 850  integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,ispden
 851  integer :: itypat,jj,jtemp,me_fft,n1,n2,n3,nproc_fft
 852  real(dp),parameter :: tolfix=1.000000001_dp
 853  real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqm1,fact,fact0,gmag
 854  real(dp) :: gsquar,rhoat,sfi,sfr
 855  real(dp) :: xnorm
 856  character(len=500) :: message
 857 !arrays
 858  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
 859  real(dp),allocatable :: length(:),spinat_indx(:,:),work(:)
 860  logical,allocatable :: use_gaussian(:)
 861 
 862 ! *************************************************************************
 863 
 864  if(nspden==4)then
 865    write(std_out,*)' initro : might work yet for nspden=4 (not checked)'
 866    write(std_out,*)' spinat',spinat(1:3,1:natom)
 867 !  stop
 868  end if
 869 
 870  n1=ngfft(1)
 871  n2=ngfft(2)
 872  n3=ngfft(3)
 873  me_fft=ngfft(11)
 874  nproc_fft=ngfft(10)
 875  ABI_MALLOC(work,(nfft))
 876  ABI_MALLOC(spinat_indx,(3,natom))
 877 
 878  ! Get the distrib associated with this fft_grid
 879  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
 880 
 881 !Transfer the spinat array to an array in which the atoms have the proper order, type by type.
 882  do ia=1,natom
 883    spinat_indx(:,atindx(ia))=spinat(:,ia)
 884  end do
 885 
 886 !Check whether the values of spinat are acceptable
 887  if(nspden==2)then
 888    ia1=1
 889    do itypat=1,ntypat
 890 !    ia1,ia2 sets range of loop over atoms:
 891      ia2=ia1+nattyp(itypat)-1
 892      do ia=ia1,ia2
 893        if( sqrt(spinat_indx(1,ia)**2+spinat_indx(2,ia)**2+spinat_indx(3,ia)**2) &
 894 &       > abs(zion(itypat))*(1.0_dp + epsilon(0.0_dp)) ) then
 895          write(message, '(a,a,a,a,i4,a,a,3es11.4,a,a,a,es11.4)' ) ch10,&
 896 &         ' initro : WARNING - ',ch10,&
 897 &         '  For type-ordered atom number ',ia,ch10,&
 898 &         '  input spinat=',spinat_indx(:,ia),'  is larger, in magnitude,',ch10,&
 899 &         '  than zion(ia)=',zion(itypat)
 900          call wrtout(std_out,message,'COLL')
 901          call wrtout(ab_out,message,'COLL')
 902        end if
 903      end do
 904      ia1=ia2+1
 905    end do
 906  end if
 907 
 908 !Compute the decay length of each type of atom
 909  ABI_MALLOC(length,(ntypat))
 910  ABI_MALLOC(use_gaussian,(ntypat))
 911  jtemp=0
 912  do itypat=1,ntypat
 913 
 914    use_gaussian(itypat)=.true.
 915    if (usepaw==0) use_gaussian(itypat) = .not. psps%nctab(itypat)%has_tvale
 916    if (usepaw==1) use_gaussian(itypat)=(pawtab(itypat)%has_tvale==0)
 917    if (.not.use_gaussian(itypat)) jtemp=jtemp+1
 918 
 919    if (use_gaussian(itypat)) then
 920      length(itypat) = atom_length(densty(itypat,1),zion(itypat),znucl(itypat))
 921      write(message,'(a,i3,a,f12.4,a,a,a,f12.4,a,i3,a,es12.4,a)' )&
 922 &     ' initro: for itypat=',itypat,', take decay length=',length(itypat),',',ch10,&
 923 &     ' initro: indeed, coreel=',znucl(itypat)-zion(itypat),', nval=',int(zion(itypat)),' and densty=',densty(itypat,1),'.'
 924      call wrtout(std_out,message,'COLL')
 925    else
 926      write(message,"(a,i3,a)")' initro: for itypat=',itypat,", take pseudo charge density from pp file"
 927      call wrtout(std_out,message,"COLL")
 928    end if
 929 
 930  end do
 931 
 932  if (jtemp>0) then
 933    dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
 934    dqm1=1.0_dp/dq
 935    dq2div6=dq**2/6.0_dp
 936  end if
 937 
 938  cutoff=gsqcut*tolfix
 939  xnorm=1.0_dp/ucvol
 940 
 941  id1=n1/2+2
 942  id2=n2/2+2
 943  id3=n3/2+2
 944 
 945  if(nspden /= 4) then
 946 
 947    do ispden=nspden,1,-1
 948 !    This loop overs spins will actually be as follows :
 949 !    ispden=2 for spin up
 950 !    ispden=1 for total spin (also valid for non-spin-polarized calculations)
 951 !    The reverse ispden order is chosen, in order to end up with
 952 !    rhog containing the proper total density.
 953 
 954      rhog(:,:)=zero
 955 
 956      ia1=1
 957      do itypat=1,ntypat
 958 
 959        if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2
 960 
 961 !      ia1,ia2 sets range of loop over atoms:
 962        ia2=ia1+nattyp(itypat)-1
 963        ii=0
 964        jtemp=0
 965 
 966        do i3=1,n3
 967          ig3=i3-(i3/id3)*n3-1
 968          do i2=1,n2
 969            ig2=i2-(i2/id2)*n2-1
 970            if (fftn2_distrib(i2)==me_fft) then
 971              do i1=1,n1
 972 
 973                ig1=i1-(i1/id1)*n1-1
 974                ii=ii+1
 975 !              gsquar=gsq_ini(ig1,ig2,ig3)
 976                gsquar=dble(ig1*ig1)*gmet(1,1)+dble(ig2*ig2)*gmet(2,2)+&
 977 &               dble(ig3*ig3)*gmet(3,3)+dble(2*ig1*ig2)*gmet(1,2)+&
 978 &               dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1)
 979 
 980 !              Skip G**2 outside cutoff:
 981                if (gsquar<=cutoff) then
 982 
 983 !                Assemble structure factor over all atoms of given type,
 984 !                also taking into account the spin-charge on each atom:
 985                  sfr=zero;sfi=zero
 986                  if(ispden==1)then
 987                    do ia=ia1,ia2
 988                      sfr=sfr+phre_ini(ig1,ig2,ig3,ia)
 989                      sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)
 990                    end do
 991                    if (use_gaussian(itypat)) then
 992                      sfr=sfr*zion(itypat)
 993                      sfi=sfi*zion(itypat)
 994                    end if
 995                  else
 996                    fact0=half;if (.not.use_gaussian(itypat)) fact0=half/zion(itypat)
 997                    do ia=ia1,ia2
 998 !                    Here, take care only of the z component
 999                      fact=fact0*(zion(itypat)+spinat_indx(3,ia))
1000                      sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact
1001                      sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact
1002                    end do
1003                  end if
1004 
1005 !                Charge density integrating to one
1006                  if (use_gaussian(itypat)) then
1007                    rhoat=xnorm*exp(-gsquar*alf2pi2)
1008 !                  Multiply structure factor times rhoat (atomic density in reciprocal space)
1009                    rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1010                    rhog(im,ii)=rhog(im,ii)+sfi*rhoat
1011                  else
1012                    gmag=sqrt(gsquar)
1013                    jj=1+int(gmag*dqm1)
1014                    diff=gmag-qgrid(jj)
1015                    bb = diff*dqm1
1016                    aa = one-bb
1017                    cc = aa*(aa**2-one)*dq2div6
1018                    dd = bb*(bb**2-one)*dq2div6
1019                    if (usepaw == 1) then
1020                      rhoat=(aa*pawtab(itypat)%tvalespl(jj,1)+bb*pawtab(itypat)%tvalespl(jj+1,1)+&
1021 &                     cc*pawtab(itypat)%tvalespl(jj,2)+dd*pawtab(itypat)%tvalespl(jj+1,2)) *xnorm
1022                    else if (usepaw == 0) then
1023                      rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+&
1024                      cc*psps%nctab(itypat)%tvalespl(jj,2)+dd*psps%nctab(itypat)%tvalespl(jj+1,2))*xnorm
1025                    else
1026                      ABI_BUG('Initialization of density is non consistent.')
1027                    end if
1028 !                  Multiply structure factor times rhoat (atomic density in reciprocal space)
1029                    rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1030                    rhog(im,ii)=rhog(im,ii)+sfi*rhoat
1031                  end if
1032 
1033                else
1034                  jtemp=jtemp+1
1035                end if
1036 
1037              end do ! End loop on i1
1038            end if
1039          end do ! End loop on i2
1040        end do ! End loop on i3
1041        ia1=ia2+1
1042 
1043      end do ! End loop on type of atoms
1044 
1045 !    Set contribution of unbalanced components to zero
1046      if (izero==1) then
1047        call zerosym(rhog,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
1048      end if
1049      !write(std_out,*)"initro: ispden, ucvol * rhog(:2,1)",ispden, ucvol * rhog(:2,1)
1050 
1051 !    Note, we end with ispden=1, so that rhog contains the total density
1052      call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
1053      rhor(:,ispden)=work(:)
1054    end do ! End loop on spins
1055 
1056  else if(nspden==4) then
1057    do ispden=nspden,1,-1
1058 !    This loop overs spins will actually be as follows :
1059 !    ispden=2,3,4 for mx,my,mz
1060 !    ispden=1 for total spin (also valid for non-spin-polarized calculations)
1061 !    The reverse ispden order is chosen, in order to end up with
1062 !    rhog containing the proper total density.
1063 
1064      rhog(:,:)=zero
1065 
1066      ia1=1
1067      do itypat=1,ntypat
1068 
1069        if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2
1070 
1071 !      ia1,ia2 sets range of loop over atoms:
1072        ia2=ia1+nattyp(itypat)-1
1073        ii=0
1074        jtemp=0
1075        do i3=1,n3
1076          ig3=i3-(i3/id3)*n3-1
1077          do i2=1,n2
1078            ig2=i2-(i2/id2)*n2-1
1079            if (fftn2_distrib(i2)==me_fft) then
1080              do i1=1,n1
1081 
1082                ig1=i1-(i1/id1)*n1-1
1083                ii=ii+1
1084 !              gsquar=gsq_ini(ig1,ig2,ig3)
1085                gsquar=dble(ig1*ig1)*gmet(1,1)+dble(ig2*ig2)*gmet(2,2)+&
1086 &               dble(ig3*ig3)*gmet(3,3)+dble(2*ig1*ig2)*gmet(1,2)+&
1087 &               dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1)
1088 
1089 !              Skip G**2 outside cutoff:
1090                if (gsquar<=cutoff) then
1091 
1092 !                Assemble structure factor over all atoms of given type,
1093 !                also taking into account the spin-charge on each atom:
1094                  sfr=zero;sfi=zero
1095                  if(ispden==1)then
1096                    do ia=ia1,ia2
1097                      sfr=sfr+phre_ini(ig1,ig2,ig3,ia)
1098                      sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)
1099                    end do
1100                    if (use_gaussian(itypat)) then
1101                      sfr=sfr*zion(itypat)
1102                      sfi=sfi*zion(itypat)
1103                    end if
1104                  else
1105                    fact0=one;if (.not.use_gaussian(itypat)) fact0=one/zion(itypat)
1106                    do ia=ia1,ia2
1107 !                    Here, take care of the components of m
1108                      fact=fact0*spinat_indx(ispden-1,ia)
1109                      sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact
1110                      sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact
1111                    end do
1112                  end if
1113 
1114 !                Charge density integrating to one
1115                  if (use_gaussian(itypat)) then
1116                    rhoat=xnorm*exp(-gsquar*alf2pi2)
1117                  else
1118                    gmag=sqrt(gsquar)
1119                    jj=1+int(gmag*dqm1)
1120                    diff=gmag-qgrid(jj)
1121                    bb = diff*dqm1
1122                    aa = one-bb
1123                    cc = aa*(aa**2-one)*dq2div6
1124                    dd = bb*(bb**2-one)*dq2div6
1125                    if (usepaw == 1) then
1126                      rhoat=(aa*pawtab(itypat)%tvalespl(jj,1)+bb*pawtab(itypat)%tvalespl(jj+1,1)+&
1127 &                     cc*pawtab(itypat)%tvalespl(jj,2)+dd*pawtab(itypat)%tvalespl(jj+1,2)) *xnorm
1128                    else if (usepaw == 0) then
1129                      rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+&
1130                      cc*psps%nctab(itypat)%tvalespl(jj,2)+dd*psps%nctab(itypat)%tvalespl(jj+1,2))*xnorm
1131                    else
1132                      ABI_BUG('Initialization of density is non consistent.')
1133                    end if
1134                  end if
1135 
1136 !                Multiply structure factor times rhoat (atomic density in reciprocal space)
1137                  rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1138                  rhog(im,ii)=rhog(im,ii)+sfi*rhoat
1139                else
1140                  jtemp=jtemp+1
1141                end if
1142 
1143              end do ! End loop on i1
1144            end if
1145          end do ! End loop on i2
1146        end do ! End loop on i3
1147        ia1=ia2+1
1148      end do ! End loop on type of atoms
1149 
1150 !    Set contribution of unbalanced components to zero
1151      if (izero==1) then
1152        call zerosym(rhog,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
1153      end if
1154      !write(std_out,*)"initro: ispden, ucvol * rhog(:2,1)",ispden, ucvol * rhog(:2,1)
1155 
1156 !    Note, we end with ispden=1, so that rhog contains the total density
1157      call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
1158      rhor(:,ispden)=work(:)
1159 
1160    end do ! End loop on spins
1161 
1162 !  Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities
1163 !    Add a small real to the magnetization
1164    if (all(abs(spinat(:,:))<tol10)) rhor(:,4)=rhor(:,4)+tol14
1165 
1166  end if ! nspden==4
1167 
1168  ABI_FREE(length)
1169  ABI_FREE(use_gaussian)
1170  ABI_FREE(spinat_indx)
1171  ABI_FREE(work)
1172 
1173  contains
1174 
1175 !Real and imaginary parts of phase.
1176    function phr_ini(x1,y1,x2,y2,x3,y3)
1177 
1178    real(dp) :: phr_ini
1179    real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1180    phr_ini=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1181  end function phr_ini
1182 
1183    function phi_ini(x1,y1,x2,y2,x3,y3)
1184 
1185    real(dp) :: phi_ini
1186    real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1187    phi_ini=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1188  end function phi_ini
1189 
1190    function ph1_ini(nri,ig1,ia)
1191 
1192    real(dp) :: ph1_ini
1193    integer,intent(in) :: nri,ig1,ia
1194    ph1_ini=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1195  end function ph1_ini
1196 
1197    function ph2_ini(nri,ig2,ia)
1198 
1199    real(dp) :: ph2_ini
1200    integer,intent(in) :: nri,ig2,ia
1201    ph2_ini=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1202  end function ph2_ini
1203 
1204    function ph3_ini(nri,ig3,ia)
1205 
1206    real(dp) :: ph3_ini
1207    integer,intent(in) :: nri,ig3,ia
1208    ph3_ini=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1209  end function ph3_ini
1210 
1211    function phre_ini(ig1,ig2,ig3,ia)
1212 
1213    real(dp) :: phre_ini
1214    integer,intent(in) :: ig1,ig2,ig3,ia
1215    phre_ini=phr_ini(ph1_ini(re,ig1,ia),ph1_ini(im,ig1,ia),&
1216 &   ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia))
1217  end function phre_ini
1218 
1219    function phimag_ini(ig1,ig2,ig3,ia)
1220 
1221    real(dp) :: phimag_ini
1222    integer,intent(in) :: ig1,ig2,ig3,ia
1223    phimag_ini=phi_ini(ph1_ini(re,ig1,ia),ph1_ini(im,ig1,ia),&
1224 &   ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia))
1225  end function phimag_ini
1226 
1227 end subroutine initro

m_mkrho/m_mkrho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mkrho

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, LSI, AR, MB, 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_mkrho
23 
24  use defs_basis
25  use defs_wvltypes
26  use m_abicore
27  use m_xmpi
28  use m_errors
29  use m_dtset
30  use m_extfpmd
31 
32  use defs_abitypes,  only : MPI_type
33  use m_time,         only : timab
34  use m_fftcore,      only : sphereboundary
35  use m_fft,          only : fftpac, zerosym, fourwf, fourdp
36  use m_hamiltonian,  only : gs_hamiltonian_type
37  use m_bandfft_kpt,  only : bandfft_kpt_set_ikpt
38  use m_paw_dmft,     only : paw_dmft_type
39  use m_spacepar,     only : symrhg
40  use defs_datatypes, only : pseudopotential_type
41  use m_atomdata,     only : atom_length
42  use m_mpinfo,       only : ptabs_fourdp, proc_distrb_cycle
43  use m_pawtab,       only : pawtab_type
44  use m_io_tools,     only : open_file
45  use m_splines,      only : spline,splint
46  use m_sort,         only : sort_dp
47  use m_prep_kgb,     only : prep_fourwf
48  use m_wvl_rho,      only : wvl_mkrho
49  use m_rot_cg,       only : rot_cg
50 
51  implicit none
52 
53  private

m_mkrho/mkrho [ Functions ]

[ Top ] [ m_mkrho ] [ Functions ]

NAME

 mkrho

FUNCTION

 Depending on option argument value:
 --Compute charge density rho(r) and rho(G) in electrons/bohr**3
   from input wavefunctions, band occupations, and k point wts.
 --Compute kinetic energy density tau(r) and tau(G) in bohr**-5
   from input wavefunctions, band occupations, and k point wts.
 --Compute a given element of the kinetic energy density tensor
   tau_{alpha,beta}(r) and tau_{alpha,beta}(G) in bohr**-5
   from input wavefunctions, band occupations, and k point wts.

INPUTS

  cg(2,mcg)=wf in G space
  dtset <type(dataset_type)>=all input variables for this dataset
   | istwfk(nkpt)=input option parameter that describes the storage of wfs
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs
   | mkmem=Number of k points treated by this node
   | mpw=maximum allowed value for npw
   | nband(nkpt*nsppol)=number of bands to be included in summation
   |  at each k point for each spin channel
   | 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
   | nkpt=number of k points
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | nsym=number of symmetry elements in group (at least 1 for identity)
   | symafm(nsym)=(anti)ferromagnetic part of symmetry operations
   | symrel(3,3,nsym)=symmetry matrices in real space (integers)
   | wtk(nkpt)=k point weights (they sum to 1.0)
  gprimd(3,3)=dimensional reciprocal space primitive translations
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  kg(3,mpw*mkmem)=reduced planewave coordinates
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mpi_enreg=information about MPI parallelization
  npwarr(nkpt)=number of planewaves and boundary planewaves at each k
  occ(mband*nkpt*nsppol)=
          occupation numbers for each band (usually 2.0) at each k point
  option if 0: compute rhor (electron density)
         if 1: compute taur (kinetic energy density)
               (i.e. Trace over the kinetic energy density tensor)
         if 2: compute taur_{alpha,beta} !!NOT YET IMPLEMENTED
               (a given element of the kinetic energy density tensor)
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  rprimd(3,3)=dimensional real space primitive translations
  tim_mkrho=timing code of the calling routine(can be set to 0 if not attributed)
  ucvol=unit cell volume (Bohr**3)
  wvl_den <type(wvl_denspot_type)>=density information for wavelets
  wvl_wfs <type(wvl_projector_type)>=wavefunctions information for wavelets

OUTPUT

 rhog(2,nfft)=total electron density in G space
 rhor(nfft,nspden)=electron density in r space
   (if spin polarized, array contains total density in first half and spin-up density in second half)
   (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz in units of hbar/2)

SOURCE

128 subroutine mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
129 &                rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wvl_wfs,&
130 &                extfpmd,option) !optional
131 
132 !Arguments ------------------------------------
133 !scalars
134  integer,intent(in) :: mcg,tim_mkrho
135  integer,intent(in),optional :: option
136  real(dp),intent(in) :: ucvol
137  type(extfpmd_type),intent(in),pointer,optional :: extfpmd
138  type(MPI_type),intent(inout) :: mpi_enreg
139  type(dataset_type),intent(in) :: dtset
140  type(paw_dmft_type), intent(in)  :: paw_dmft
141  type(wvl_wf_type),intent(inout) :: wvl_wfs
142  type(wvl_denspot_type), intent(inout) :: wvl_den
143 !no_abirules
144 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
145  integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,  &
146 &               (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
147  integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
148  real(dp), intent(in) :: gprimd(3,3)
149  real(dp), intent(in) :: cg(2,mcg)
150  real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
151 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
152  real(dp), intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym),  &
153 &                                 (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
154  real(dp), intent(in) :: rprimd(3,3)
155  real(dp), intent(out) :: rhor(dtset%nfft,dtset%nspden),rhog(2,dtset%nfft)
156 
157 !Local variables-------------------------------
158 !scalars
159  integer,save :: nskip=0
160  integer :: alpha,use_nondiag_occup_dmft,bdtot_index,beta,blocksize,iband,iband1,ibandc1,ib,iblock,icg,ierr
161  integer :: ifft,ikg,ikpt,ioption,ipw,ipwsp,ishf,ispden,ispinor,ispinor_index
162  integer :: isppol,istwf_k,jspinor_index
163  integer :: me,my_nspinor,n1,n2,n3,n4,n5,n6,nalpha,nband_k,nbandc1,nbdblock,nbeta
164  integer :: ndat,nfftot,npw_k,spaceComm,tim_fourwf
165  integer :: iband_me
166  integer :: mband_mem
167  real(dp) :: kpt_cart,kg_k_cart,gp2pi1,gp2pi2,gp2pi3,cwftmp
168  real(dp) :: weight,weight_i
169  !character(len=500) :: message
170 !arrays
171  integer,allocatable :: gbound(:,:),kg_k(:,:)
172  logical :: locc_test,nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
173  real(dp) :: dummy(2,1),tsec(2)
174  real(dp),allocatable :: cwavef(:,:,:),cwavefb(:,:,:),cwavef_x(:,:)
175  real(dp),allocatable :: cwavef_y(:,:),cwavefb_x(:,:),cwavefb_y(:,:),kg_k_cart_block(:)
176  real(dp),allocatable :: occ_k(:),rhoaug(:,:,:),rhoaug_down(:,:,:)
177  real(dp),allocatable :: rhoaug_mx(:,:,:),rhoaug_my(:,:,:),rhoaug_up(:,:,:)
178  real(dp),allocatable :: taur_alphabeta(:,:,:,:),wfraug(:,:,:,:)
179  real(dp),allocatable :: occ_diag(:)
180 ! real(dp),allocatable :: occ_nd(2, :, :)
181  real(dp),allocatable :: cwavef_rot(:,:,:,:)
182 
183 ! *************************************************************************
184 
185  DBG_ENTER("COLL")
186 
187  call timab(790+tim_mkrho,1,tsec)
188  call timab(799,1,tsec)
189 
190  if(.not.(present(option))) then
191    ioption=0
192  else
193    ioption=option
194  end if
195 
196  if(ioption/=0.and.paw_dmft%use_sc_dmft==1) then
197    ABI_ERROR('option argument value of this routines should be 0 if usedmft=1.')
198  end if
199  if(paw_dmft%use_sc_dmft/=0) then
200    nbandc1=(paw_dmft%mbandc-1)*paw_dmft%use_sc_dmft+1
201  else
202    nbandc1=1
203  end if
204  use_nondiag_occup_dmft=0
205 
206 
207 !if(dtset%nspinor==2.and.paw_dmft%use_sc_dmft==1) then
208 !write(message, '(a,a,a,a)' )ch10,&
209 !&   ' mkrho : ERROR -',ch10,&
210 !&   '  nspinor argument value of this routines should be 1 if usedmft=1. '
211 !call wrtout(std_out,message,'COLL')
212 !end if
213 
214  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
215  if (mpi_enreg%paral_spinor==0) then
216    ispinor_index=1;jspinor_index=1
217    nspinor1TreatedByThisProc=.true.
218    nspinor2TreatedByThisProc=(dtset%nspinor==2)
219  else
220    ispinor_index=mpi_enreg%me_spinor+1;jspinor_index=3-ispinor_index
221    nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
222    nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
223  end if
224 
225 !Set local variable which depend on option argument
226 
227 !nalpha*nbeta is the number of element of the kinetic energy density tensor
228 !to be computed in the irreducible Brillouin Zone (BZ) to get the result in the full BZ.
229 !In case of electron density calculation, nalpha=nbeta=1
230  select case (ioption)
231  case (0)
232    nalpha = 1
233    nbeta = 1
234  case (1)
235    nalpha = 3
236    nbeta = 1
237    ABI_MALLOC(taur_alphabeta,(dtset%nfft,dtset%nspden,3,1))
238  case (2)
239    nalpha = 3
240    nbeta = 3
241    ABI_MALLOC(taur_alphabeta,(dtset%nfft,dtset%nspden,3,3))
242  case default
243    ABI_BUG('ioption argument value should be 0,1 or 2.')
244  end select
245 
246 !Init me
247  me=mpi_enreg%me_kpt
248 !zero the charge density array in real space
249 !$OMP PARALLEL DO COLLAPSE(2)
250  do ispden=1,dtset%nspden
251    do ifft=1,dtset%nfft
252      rhor(ifft,ispden)=zero
253    end do
254  end do
255 
256 !WVL - Branching with a separate mkrho procedure in wavelet.
257  if (dtset%usewvl == 1) then
258    select case(ioption)
259    case (0)
260      call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den)
261      return
262    case (1)
263      !call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den)
264      ABI_ERROR("kinetic energy density (taur) is not yet implemented in wavelet formalism.")
265    case (2)
266      !call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den)
267      ABI_BUG('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented in wavelet formalism.')
268    end select
269  end if
270 !WVL - Following is done in plane waves.
271 
272 !start loop over alpha and beta
273 
274  do alpha=1,nalpha
275    do beta=1,nbeta
276 
277 !    start loop over spin and k points
278      bdtot_index=0
279      icg=0
280 
281 !    n4,n5,n6 are FFT dimensions, modified to avoir cache trashing
282      n1 = dtset%ngfft(1) ; n2 = dtset%ngfft(2) ; n3 = dtset%ngfft(3)
283      n4 = dtset%ngfft(4) ; n5 = dtset%ngfft(5) ; n6 = dtset%ngfft(6)
284      ndat = 1 ; if (mpi_enreg%paral_kgb==1) ndat = mpi_enreg%bandpp
285      ABI_MALLOC(cwavef,(2,dtset%mpw,my_nspinor))
286      ABI_MALLOC(rhoaug,(n4,n5,n6))
287      ABI_MALLOC(wfraug,(2,n4,n5,n6*ndat))
288      ABI_MALLOC(cwavefb,(2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor))
289      if(dtset%nspden==4) then
290        ABI_MALLOC(rhoaug_up,(n4,n5,n6))
291        ABI_MALLOC(rhoaug_down,(n4,n5,n6))
292        ABI_MALLOC(rhoaug_mx,(n4,n5,n6))
293        ABI_MALLOC(rhoaug_my,(n4,n5,n6))
294        rhoaug_up(:,:,:)=zero
295        rhoaug_down(:,:,:)=zero
296        rhoaug_mx(:,:,:)=zero
297        rhoaug_my(:,:,:)=zero
298      end if
299 
300      do isppol=1,dtset%nsppol
301        ikg=0
302 
303        rhoaug(:,:,:)=zero
304        do ikpt=1,dtset%nkpt
305 
306          nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
307          mband_mem = nband_k
308          npw_k=npwarr(ikpt)
309          istwf_k = dtset%istwfk(ikpt)
310 
311          if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
312            bdtot_index=bdtot_index+nband_k
313            cycle
314          end if
315 
316          ABI_MALLOC(gbound,(2*dtset%mgfft+8,2))
317          ABI_MALLOC(kg_k,(3,npw_k))
318 
319          kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
320          call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k)
321 
322 !        Loop over bands to fft and square for rho(r)
323 !        Shoulb be changed to treat bands by batch always
324 
325          if(mpi_enreg%paral_kgb /= 1) then  ! Not yet parallelized on spinors
326            mband_mem = nband_k/mpi_enreg%nproc_band
327            iband_me = 0
328            do iband=1,nband_k
329 !            if(paw_dmft%use_sc_dmft==1) then
330 !            write(std_out,*) 'iband  ',iband,occ(iband+bdtot_index),paw_dmft%occnd(iband,iband,ikpt,isppol)
331 !            else
332 !            write(std_out,*) 'iband  ',iband,occ(iband+bdtot_index)
333 !            endif
334              if(mpi_enreg%paralbd==1)then
335                if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle
336              end if
337              iband_me = iband_me + 1
338              do ibandc1=1,nbandc1 ! in case of DMFT
339 !              Check if DMFT and only treat occupied states (check on occ.)
340                if(paw_dmft%use_sc_dmft == 1) then
341                  iband1 = paw_dmft%include_bands(ibandc1)
342                  if(paw_dmft%band_in(iband)) then
343                    if(.not. paw_dmft%band_in(iband1))  stop
344                    use_nondiag_occup_dmft = 1
345                    locc_test = abs(paw_dmft%occnd(1,iband,iband1,ikpt,isppol)) +&
346 &                   abs(paw_dmft%occnd(2,iband,iband1,ikpt,isppol))>tol8
347 !                  write(std_out,*) "mkrho,ikpt,iband,use_occnd",ikpt,iband
348                  else
349                    use_nondiag_occup_dmft = 0
350                    locc_test = abs(occ(iband+bdtot_index))>tol8
351                    if(ibandc1 /=1 .and. .not. paw_dmft%band_in(iband)) cycle
352                  end if
353                else
354                  use_nondiag_occup_dmft = 0
355                  locc_test = abs(occ(iband+bdtot_index))>tol8
356                end if
357 
358                if (locc_test) then
359 !                Obtain Fourier transform in fft box and accumulate the density or the kinetic energy density
360 !                Not yet parallise on nspinor if paral_kgb non equal to 1
361                  ipwsp=(iband_me-1)*npw_k*my_nspinor +icg
362                  cwavef(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k)
363                  if (my_nspinor==2) cwavef(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k)
364 
365                  if(ioption==1)then
366 !                  Multiplication by 2pi i (k+G)_alpha
367                    gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi
368                    kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt)
369                    do ispinor=1,my_nspinor
370                      do ipw=1,npw_k
371                        kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart
372                        cwftmp=-cwavef(2,ipw,ispinor)*kg_k_cart
373                        cwavef(2,ipw,ispinor)=cwavef(1,ipw,ispinor)*kg_k_cart
374                        cwavef(1,ipw,ispinor)=cwftmp
375                      end do
376                    end do
377                  else if(ioption==2)then
378                    ABI_ERROR('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.')
379                  end if
380 
381 !                Non diag occupation in DMFT.
382 ! TODO : this will break in full distrib of band memory
383                  if(use_nondiag_occup_dmft==1) then
384                    ipwsp=(iband1-1)*npw_k*my_nspinor +icg
385                    cwavefb(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k)
386                    if (my_nspinor==2) cwavefb(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k)
387                    weight  =paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol
388                    weight_i=paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol
389 
390                  else
391                    weight=occ(iband+bdtot_index)*dtset%wtk(ikpt)/ucvol
392                    weight_i=weight
393                  end if
394 
395                  if(mpi_enreg%paralbd==0) tim_fourwf=3
396                  if(mpi_enreg%paralbd==1) tim_fourwf=6
397 
398 !                The same section of code is also found in vtowfk.F90 : should be rationalized !
399 
400                  call fourwf(1,rhoaug,cwavef(:,:,1),dummy,wfraug,gbound,gbound,&
401 &                 istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
402 &                 npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
403 &                 use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,1),&
404 &                 use_gpu_cuda=dtset%use_gpu_cuda)
405 
406 
407                  if(dtset%nspinor==2)then
408 !                  DEBUG GZ !To obtain a x-directed magnetization(test)
409 !                  cwavef1(1,1:npw_k)=-cwavef(2,1:npw_k)
410 !                  cwavef1(2,1:npw_k)= cwavef(1,1:npw_k)
411 !                  ENDDEBUG
412 
413                    if(dtset%nspden==1) then
414 !                    We need only the total density : accumulation continues on top of rhoaug
415 
416                      call fourwf(1,rhoaug,cwavef(:,:,2),dummy,wfraug,gbound,gbound,&
417 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
418 &                     npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
419 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),use_gpu_cuda=dtset%use_gpu_cuda)
420 
421 
422                    else if(dtset%nspden==4) then
423                      ! Build the four components of rho. We use only norm quantities and, so fourwf.
424                      ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
425                      ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
426                      ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
427                      ABI_MALLOC(cwavef_x,(2,npw_k))
428                      ABI_MALLOC(cwavef_y,(2,npw_k))
429                      ABI_MALLOC(cwavefb_x,(2,npw_k*paw_dmft%use_sc_dmft))
430                      ABI_MALLOC(cwavefb_y,(2,npw_k*paw_dmft%use_sc_dmft))
431                      ! $(\Psi^{1}+\Psi^{2})$
432                      cwavef_x(:,:)=cwavef(:,1:npw_k,1)+cwavef(:,1:npw_k,2)
433                      ! $(\Psi^{1}-i \Psi^{2})$
434                      cwavef_y(1,:)=cwavef(1,1:npw_k,1)+cwavef(2,1:npw_k,2)
435                      cwavef_y(2,:)=cwavef(2,1:npw_k,1)-cwavef(1,1:npw_k,2)
436                      if(use_nondiag_occup_dmft==1) then
437                        cwavefb_x(:,:)=cwavefb(:,1:npw_k,1)+cwavefb(:,1:npw_k,2)
438                        cwavefb_y(1,:)=cwavefb(1,1:npw_k,1)+cwavefb(2,1:npw_k,2)
439                        cwavefb_y(2,:)=cwavefb(2,1:npw_k,1)-cwavefb(1,1:npw_k,2)
440                      end if
441                      rhoaug_up(:,:,:)=rhoaug(:,:,:) !Already computed
442 
443                      call fourwf(1,rhoaug_down,cwavef(:,:,2),dummy,wfraug,gbound,gbound,&
444 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
445 &                     npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
446 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),use_gpu_cuda=dtset%use_gpu_cuda)
447 
448                      call fourwf(1,rhoaug_mx,cwavef_x,dummy,wfraug,gbound,gbound,&
449 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
450 &                     npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
451 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_x,use_gpu_cuda=dtset%use_gpu_cuda)
452 
453                      call fourwf(1,rhoaug_my,cwavef_y,dummy,wfraug,gbound,gbound,&
454 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
455 &                     npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
456 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_y,use_gpu_cuda=dtset%use_gpu_cuda)
457 
458                      ABI_FREE(cwavef_x)
459                      ABI_FREE(cwavef_y)
460                      ABI_FREE(cwavefb_x)
461                      ABI_FREE(cwavefb_y)
462 
463                    end if ! dtset%nspden/=4
464                  end if
465 
466                else
467 !                Accumulate the number of one-way 3D ffts skipped
468                  nskip=nskip+1
469                end if ! abs(occ(iband+bdtot_index))>tol8
470 !              End loop on iband
471              end do ! iband1=1,(nband_k-1)*paw_dmft%use_sc_dmft+1
472            end do ! iband=1,nband_k
473 
474          else !paral_kgb==1
475            if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
476 
477            call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
478 
479            nbdblock=nband_k/(mpi_enreg%nproc_band * mpi_enreg%bandpp)
480            blocksize=nband_k/nbdblock
481            if(allocated(cwavef))  then
482              ABI_FREE(cwavef)
483            end if
484            ABI_MALLOC(cwavef,(2,npw_k*blocksize,dtset%nspinor))
485            if(ioption==1)  then
486              ABI_MALLOC(kg_k_cart_block,(npw_k))
487            end if
488            ABI_MALLOC(occ_k,(nband_k))
489            occ_k(:)=occ(bdtot_index+1:bdtot_index+nband_k)
490 
491 ! ---------- DMFT
492            if(allocated(cwavef_rot))  then
493              ABI_FREE(cwavef_rot)
494              ABI_FREE(occ_diag)
495              ! ABI_FREE(occ_nd)
496            end if
497            if(paw_dmft%use_sc_dmft==1) then
498              ! Allocation of DMFT temporaries arrays
499              ABI_MALLOC(cwavef_rot,(2,npw_k,blocksize,dtset%nspinor))
500              ABI_MALLOC(occ_diag,(blocksize))
501              ! ABI_MALLOC(occ_nd,(2, blocksize, blocksize, dtset%nspinor))
502            end if
503 ! ---------- END DMFT
504 
505            do iblock=1,nbdblock
506              if (dtset%nspinor==1) then
507                cwavef(:,1:npw_k*blocksize,1)=cg(:,1+(iblock-1)*npw_k*blocksize+icg:iblock*npw_k*blocksize+icg)
508              else
509                if (mpi_enreg%paral_spinor==0) then
510                  ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg
511                  do ib=1,blocksize
512                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,1)=cg(:,1+(2*ib-2)*npw_k+ishf:(2*ib-1)*npw_k+ishf)
513                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,2)=cg(:,1+(2*ib-1)*npw_k+ishf:ib*2*npw_k+ishf)
514                  end do
515                else
516                  ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg
517                  do ib=1,blocksize
518                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,ispinor_index)=&
519 &                   cg(:,1+(ib-1)*npw_k+ishf:ib*npw_k+ishf)
520                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,jspinor_index)=zero
521                  end do
522                  call xmpi_sum(cwavef,mpi_enreg%comm_spinor,ierr)
523                end if
524              end if
525 
526              if(ioption==1)then
527 !              Multiplication by 2pi i (k+G)_alpha
528                gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi
529                kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt)
530                kg_k_cart_block(1:npw_k)=gp2pi1*kg_k(1,1:npw_k)+gp2pi2*kg_k(2,1:npw_k)+gp2pi3*kg_k(3,1:npw_k)+kpt_cart
531                do ib=1,blocksize
532                  do ipw=1,npw_k
533                    cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
534                    cwavef(2,ipw,1)=cwavef(1,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
535                    cwavef(1,ipw,1)=cwftmp
536                    if (my_nspinor==2) then
537                      cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
538                      cwavef(2,ipw,2)=cwavef(1,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
539                      cwavef(1,ipw,2)=cwftmp
540                    end if
541                  end do
542                end do
543              else if(ioption==2)then
544                ABI_ERROR("kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.")
545              end if
546 
547 ! ---------- DMFT
548              if(paw_dmft%use_sc_dmft==1) then
549                ! initialisation of DMFT arrays
550                cwavef_rot(:,:,:,:) = zero
551                occ_diag(:) = zero
552                ! occ_nd(:,:,:,:) = paw_dmft%occnd(:,:,:,ikpt,:)
553 
554                do ib=1,blocksize
555                  cwavef_rot(:, :, ib, :) = cwavef(:, 1+(ib-1)*npw_k:ib*npw_k, :)
556                end do
557 
558                call rot_cg(paw_dmft%occnd(:,:,:,ikpt,isppol), cwavef_rot, npw_k, nband_k, blocksize,&
559 &                          dtset%nspinor, paw_dmft%include_bands(1), paw_dmft%mbandc, occ_diag)
560                do ib=1,blocksize
561                  cwavef(:, 1+(ib-1)*npw_k:ib*npw_k, :) = cwavef_rot(:, :, ib, :)
562                end do
563 
564                occ_k(:) = occ_diag(:)
565              end if
566 ! ---------- END DMFT
567 
568              call timab(538,1,tsec)
569              if (nspinor1TreatedByThisProc) then
570                call prep_fourwf(rhoaug,blocksize,cwavef(:,:,1),wfraug,iblock,istwf_k,dtset%mgfft,mpi_enreg,&
571 &               nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
572 &               dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
573              end if
574              call timab(538,2,tsec)
575              if(dtset%nspinor==2)then
576                if (dtset%nspden==1) then
577                  if (nspinor2TreatedByThisProc) then
578                    call prep_fourwf(rhoaug,blocksize,cwavef(:,:,2),wfraug,&
579 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
580 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
581 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
582                  end if
583                else if(dtset%nspden==4 ) then
584                  ABI_MALLOC(cwavef_x,(2,npw_k*blocksize))
585                  ABI_MALLOC(cwavef_y,(2,npw_k*blocksize))
586                  cwavef_x(:,:)=cwavef(:,:,1)+cwavef(:,:,2)
587                  cwavef_y(1,:)=cwavef(1,:,1)+cwavef(2,:,2)
588                  cwavef_y(2,:)=cwavef(2,:,1)-cwavef(1,:,2)
589                  call timab(538,1,tsec)
590                  if (nspinor1TreatedByThisProc) then
591                    call prep_fourwf(rhoaug_down,blocksize,cwavef(:,:,2),wfraug,&
592 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
593 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
594 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
595                  end if
596                  if (nspinor2TreatedByThisProc) then
597                    call prep_fourwf(rhoaug_mx,blocksize,cwavef_x,wfraug,&
598 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
599 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
600 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
601                    call prep_fourwf(rhoaug_my,blocksize,cwavef_y,wfraug,&
602 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
603 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
604 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
605                  end if
606                  call timab(538,2,tsec)
607                  ABI_FREE(cwavef_x)
608                  ABI_FREE(cwavef_y)
609                end if
610              end if
611            end do !iblock
612            if(ioption==1)  then
613              ABI_FREE(kg_k_cart_block)
614            end if
615            if (allocated(cwavef))  then
616              ABI_FREE(cwavef)
617            end if
618            ABI_FREE(occ_k)
619          end if
620 
621          ABI_FREE(gbound)
622          ABI_FREE(kg_k)
623 
624          bdtot_index=bdtot_index+nband_k
625 
626          if (dtset%mkmem/=0) then
627            icg=icg+npw_k*my_nspinor*mband_mem !iband_me
628            ikg=ikg+npw_k
629          end if
630 
631        end do ! ikpt
632 
633        if(mpi_enreg%paral_kgb == 1) then
634          call bandfft_kpt_set_ikpt(-1,mpi_enreg)
635          if (dtset%nspden==4) then
636 !          Sum the contribution of the band and of the FFT
637            call xmpi_sum(rhoaug     ,mpi_enreg%comm_bandspinorfft, ierr)
638            call xmpi_sum(rhoaug_down,mpi_enreg%comm_bandspinorfft, ierr)
639            call xmpi_sum(rhoaug_mx ,mpi_enreg%comm_bandspinorfft, ierr)
640            call xmpi_sum(rhoaug_my ,mpi_enreg%comm_bandspinorfft, ierr)
641            rhoaug_up(:,:,:) = rhoaug(:,:,:)
642          else
643            call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr)
644          end if
645        end if
646 
647 !      Transfer density on augmented fft grid to normal fft grid in real space
648 !      Take also into account the spin, to place it correctly in rhor.
649        if(dtset%nspden==1 .or. dtset%nspden==2) then
650          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug,1)
651        else if(dtset%nspden==4) then
652          ispden=1
653          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_up,1)
654          ispden=2
655          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_mx,1)
656          ispden=3
657          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_my,1)
658          ispden=4
659          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_down,1)
660          ABI_FREE(rhoaug_up)
661          ABI_FREE(rhoaug_down)
662          ABI_FREE(rhoaug_mx)
663          ABI_FREE(rhoaug_my)
664        end if
665 
666      end do ! isppol
667 
668      if(allocated(cwavef))  then
669        ABI_FREE(cwavef)
670      end if
671      ABI_FREE(cwavefb)
672      ABI_FREE(rhoaug)
673      ABI_FREE(wfraug)
674 
675      if(allocated(cwavef_rot))  then
676        ABI_FREE(cwavef_rot)
677        ABI_FREE(occ_diag)
678        ! ABI_FREE(occ_nd)
679      end if
680 
681 !    Recreate full rhor on all proc.
682      call timab(48,1,tsec)
683      call timab(71,1,tsec)
684      spaceComm=mpi_enreg%comm_cell
685      if (mpi_enreg%paral_hf==1)spaceComm=mpi_enreg%comm_kpt
686      if(mpi_enreg%paral_kgb==1)spaceComm=mpi_enreg%comm_kpt
687      call xmpi_sum(rhor,spaceComm,ierr)
688      call timab(71,2,tsec)
689      call timab(48,2,tsec)
690 
691      call timab(799,2,tsec)
692      call timab(549,1,tsec)
693 
694      if(ioption==1 .or. ioption==2) then
695 !$OMP PARALLEL DO COLLAPSE(2)
696        do ispden=1,dtset%nspden
697          do ifft=1,dtset%nfft
698            taur_alphabeta(ifft,ispden,alpha,beta) = rhor(ifft,ispden)
699          end do
700        end do
701      end if
702 
703    end do !  beta=1,nbeta
704  end do !  alpha=1,nalpha
705 
706 !Compute the trace over the kinetic energy density tensor. i.e. Sum of the 3 diagonal elements.
707  if(ioption==1)then
708 !  zero rhor array in real space
709    do ispden=1,dtset%nspden
710 !$OMP PARALLEL DO
711      do ifft=1,dtset%nfft
712        rhor(ifft,ispden)=zero
713      end do
714    end do
715    do alpha = 1, nalpha
716 !$OMP PARALLEL DO COLLAPSE(2)
717      do ispden=1,dtset%nspden
718        do ifft=1,dtset%nfft
719          rhor(ifft,ispden) = rhor(ifft,ispden) + taur_alphabeta(ifft,ispden,alpha,1)
720        end do
721      end do
722    end do
723  end if
724 
725  nfftot=dtset%ngfft(1) * dtset%ngfft(2) * dtset%ngfft(3)
726 
727 !Add extfpmd free electrons contribution to density
728  if(present(extfpmd)) then
729    if(associated(extfpmd)) then
730      if(extfpmd%version==1.or.extfpmd%version==2.or.extfpmd%version==3.or.extfpmd%version==4) then
731        rhor(:,:)=rhor(:,:)+extfpmd%nelect/ucvol/dtset%nspden
732      else if(extfpmd%version==10) then
733        do ispden=1,dtset%nspden
734          do ifft=1,dtset%nfft
735            rhor(ifft,ispden)=rhor(ifft,ispden)+extfpmd%nelectarr(ifft,ispden)/ucvol/dtset%nspden
736          end do
737        end do
738      end if
739      rhog(1,1)=rhog(1,1)+extfpmd%nelect/ucvol/dtset%nspden
740    end if
741  end if
742 
743  select case (ioption)
744  case(0, 1)
745    call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
746                phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
747    if(ioption==1)then
748 !$OMP PARALLEL DO
749      do ifft=1,dtset%nfft
750        do ispden=1,dtset%nspden
751          rhor(ifft,ispden)=1.0d0/2.0d0*rhor(ifft,ispden)
752        end do
753        rhog(:,ifft)=1.0d0/2.0d0*rhog(:,ifft)
754      end do
755    end if
756  case(2)
757    ABI_BUG('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.')
758    !call symtaug(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
759    !dtset%paral_kgb,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel)
760  end select
761 
762  call timab(549,2,tsec)
763 
764 !We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2
765 !we also have the spin-up density, symmetrized, in rhor(:,2).
766 !In case of non collinear magnetism, we have rho,mx,my,mz. No symmetry is applied
767 
768  call timab(799,1,tsec)
769 
770  if(ioption==1 .or. ioption==2)  then
771    ABI_FREE(taur_alphabeta)
772  end if
773 
774 !Find and print minimum and maximum total electron density
775 !(or total kinetic energy density, or total element of kinetic energy density tensor) and locations
776  call wrtout(std_out,' mkrho: echo density (plane-wave part only)','COLL')
777  call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,optrhor=ioption,ucvol=ucvol)
778 
779  call timab(799,2,tsec)
780  call timab(790+tim_mkrho,2,tsec)
781 
782  DBG_EXIT("COLL")
783 
784 end subroutine mkrho

m_mkrho/prtrhomxmn [ Functions ]

[ Top ] [ m_mkrho ] [ Functions ]

NAME

 prtrhomxmn

FUNCTION

 If option==1, compute the maximum and minimum of the density (and spin-polarization if nspden==2), and print it.
 If option==2, also compute and print the second maximum or minimum

INPUTS

  iout=unit for output file
  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
  option, see above
  optrhor=option for rhor (If optrhor==0, rhor is expected to be the electron density)
                          (If optrhor==1, rhor is expected to be the kinetic energy density (taur))
                          (If optrhor==2, rhor is expected to be the gradient of the electron density (grhor))
                          (If optrhor==3, rhor is expected to be the laplacian of the electron density (lrhor))
                          (If optrhor==4, rhor is expected to be the ELF (elfr))
  rhor(nfft,nspden)=electron density (electrons/bohr^3)

OUTPUT

NOTES

  The tolerance tol12 aims at giving a machine-independent ordering.
  (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)

SOURCE

1260 subroutine prtrhomxmn(iout,mpi_enreg,nfft,ngfft,nspden,option,rhor,optrhor,ucvol)
1261 
1262 !Arguments ------------------------------------
1263 !scalars
1264  integer,intent(in) :: iout,nfft,nspden,option
1265  integer,intent(in),optional :: optrhor
1266  real(dp),intent(in),optional :: ucvol
1267  type(MPI_type),intent(in) :: mpi_enreg
1268 !arrays
1269  integer,intent(in) :: ngfft(18)
1270  real(dp),intent(in) :: rhor(nfft,nspden)
1271 
1272 !Local variables-------------------------------
1273 !scalars
1274  integer :: i1,i2,i3,ierr,ifft,ii,iisign,iitems,index1,ioptrhor
1275  integer :: index2,indsign,iproc,istart,me,n1,n2,n3,nitems
1276  integer :: nfft_,nfftot,nproc,spaceComm
1277  real(dp) :: temp,value1,value2
1278  character(len=500) :: message,txt1_in_mssg,txt2_in_mssg,txt3_in_mssg
1279  logical :: reduce=.false.
1280 !arrays
1281  integer,allocatable :: iindex(:,:,:),index_fft(:,:,:,:)
1282  real(dp) :: rhomn1(4),rhomn2(4),rhomx1(4),rhomx2(4),ri_rhomn1(3,4)
1283  real(dp) :: ri_rhomn2(3,4),ri_rhomx1(3,4),ri_rhomx2(3,4),ri_zetmn1(3,2)
1284  real(dp) :: ri_zetmn2(3,2),ri_zetmx1(3,2),ri_zetmx2(3,2),zetmn1(2)
1285  real(dp) :: zetmn2(2),zetmx1(2),zetmx2(2)
1286  real(dp),allocatable :: array(:),coord(:,:,:,:),value(:,:,:),integrated(:)
1287  real(dp),allocatable :: value_fft(:,:,:)
1288 
1289 ! *************************************************************************
1290 
1291  if(.not.(present(optrhor))) then
1292    ioptrhor=0
1293  else
1294    ioptrhor=optrhor
1295  end if
1296 
1297  if(option/=1 .and. option/=2)then
1298    write(message, '(a,i0)' )' Option must be 1 or 2, while it is ',option
1299    ABI_BUG(message)
1300  end if
1301 
1302  if (mpi_enreg%nproc_wvl>1) then
1303 !  nfft is always the potential size (in GGA, the density has buffers).
1304    nfft_ = ngfft(1) * ngfft(2) * mpi_enreg%nscatterarr(mpi_enreg%me_wvl, 2)
1305    n1 = ngfft(1)
1306    n2 = ngfft(2)
1307    n3 = sum(mpi_enreg%nscatterarr(:, 2))
1308    istart = mpi_enreg%nscatterarr(mpi_enreg%me_wvl, 4)
1309  else
1310    nfft_ = nfft
1311    n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1312    istart = 0
1313  end if
1314 
1315 !--------------------------------------------------------------------------
1316 !One has to determine the maximum and minimum (etc...) values
1317 !over all space, and then output it, as well as to identify
1318 !the point at which it occurs ...
1319 !This will require a bit of data exchange, and correct indirect indexing ...
1320 
1321 !For the local processor, find different items :
1322 !maximum and minimum total electron density and locations
1323 !and also spin-polarisation and magnetization
1324 !also keep the second maximal or minimal value
1325  if(nspden==1)nitems=1   ! Simply the total density
1326  if(nspden==2)nitems=5   ! Total density, spin up, spin down, magnetization, zeta
1327  if(nspden==4)nitems=6   ! Total density, x, y, z, magnetization, zeta
1328 
1329  ABI_MALLOC(value,(2,2,nitems))
1330  ABI_MALLOC(iindex,(2,2,nitems))
1331  ABI_MALLOC(array,(nfft))
1332  ABI_MALLOC(integrated,(nitems))
1333 
1334  do iitems=1,nitems
1335 
1336 !  Copy the correct values into the array
1337 !  First set of items : the density, for each spin component
1338    if(iitems<=nspden)then
1339      array(:)=rhor(:,iitems)
1340    end if
1341 !  Case nspden==2, some computation to be done
1342    if(nspden==2)then
1343      if(iitems==3)then ! Spin down
1344        array(:)=rhor(:,1)-rhor(:,2)
1345      else if(iitems==4)then  ! Magnetization
1346        array(:)=2*rhor(:,2)-rhor(:,1)
1347      else if(iitems==5)then  ! zeta = relative magnetization
1348        ! Avoid 0/0: the limit of (x - y) / (x+ y) depends on the direction.
1349        array(:)=zero
1350        where (abs(rhor(:,1)) > tol12) array(:)=(2*rhor(:,2)-rhor(:,1))/rhor(:,1)
1351      end if
1352 !    Case nspden==4, some other computation to be done
1353    else if(nspden==4)then
1354      if(iitems==5)then ! Magnetization
1355        array(:)=sqrt(rhor(:,2)**2+rhor(:,3)**2+rhor(:,4)**2)
1356      else if(iitems==6)then ! zeta = relative magnetization
1357        array(:)=(sqrt(rhor(:,2)**2+rhor(:,3)**2+rhor(:,4)**2))/rhor(:,1)
1358      end if
1359    end if
1360 
1361 !  Zero all the absolute values that are lower than tol8, for portability reasons.
1362    do ifft = 1, nfft_
1363      if(abs(array(ifft))<tol8)array(ifft)=zero
1364    end do
1365 
1366 !  DEBUG
1367 !  write(std_out,*) ' iitems,array(1:2)=',iitems,array(1:2)
1368 !  ENDDEBUG
1369 
1370    do indsign=1,2 ! Find alternatively the maximum and the minimum
1371      iisign=3-2*indsign
1372 
1373      if (nfft_ > 1) then
1374 !      Initialize the two first values
1375        value1=array(istart + 1) ; value2=array(istart + 2)
1376        index1=1 ; index2=2
1377 
1378 !      Ordering, if needed
1379        if( iisign*(value2+tol12) > iisign*(value1)) then
1380          temp=value2 ; value2=value1 ; value1=temp
1381          index1=2 ; index2=1
1382        end if
1383 
1384 !      Integration, if relevant
1385        if(present(ucvol).and. indsign==1)then
1386          integrated(iitems) = array(istart + 1)+array(istart + 2)
1387        end if
1388      else
1389        value1 = zero; value2 = zero
1390        index1 = 0;    index2 = 0
1391      end if
1392 
1393 !    DEBUG
1394 !    write(std_out,*) ' value1,value2,index1,index2=',value1,value2,index1,index2
1395 !    ENDDEBUG
1396 
1397 !    Loop over all points
1398      do ifft = 3, nfft_
1399 
1400        temp=array(istart + ifft)
1401        if(present(ucvol).and. indsign==1)integrated(iitems) = integrated(iitems)+temp
1402 !      Compares it to the second value
1403        if( iisign*(temp+tol12) > iisign*value2 ) then
1404 !        Compare it to the first value
1405          if( iisign*(temp+tol12) > iisign*value1 ) then
1406            value2=value1 ; index2=index1
1407            value1=temp   ; index1=ifft
1408          else
1409            value2=temp   ; index2=ifft
1410          end if
1411        end if
1412 
1413      end do ! ifft
1414 
1415      value(1,indsign,iitems)=value1
1416      value(2,indsign,iitems)=value2
1417      iindex(1,indsign,iitems)=index1
1418      iindex(2,indsign,iitems)=index2
1419 
1420 !    DEBUG
1421 !    write(std_out,*) ' it,v1,i1=',iitems, value1,index1
1422 !    write(std_out,*) ' it,v2,i2=',iitems, value2,index2
1423 !    ENDDEBUG
1424 
1425    end do ! indsign
1426 
1427    if(present(ucvol))then
1428      nfftot=ngfft(1) * ngfft(2) * ngfft(3)
1429      integrated(iitems)=integrated(iitems)*ucvol/nfftot
1430    end if
1431 
1432 !  Integrate the array
1433 !  integrated(iitems)=zero
1434 !  do ifft=1,nfft_
1435 !  integrated(iitems) = integrated(iitems) + array(istart + ifft)
1436 !  enddo
1437 !  if(present(ucvol))integrated(iitems) = integrated(iitems)*ucvol/nfft_
1438 !  write(std_err,*)present(ucvol)
1439 !  if(present(ucvol))then
1440 !  write(std_err,*)ucvol
1441 !  endif
1442 
1443  end do ! iitems
1444 
1445  ABI_FREE(array)
1446 
1447 !-------------------------------------------------------------------
1448 !Enter section for FFT parallel case
1449 !if(mpi_enreg%paral_kgb>1) spaceComm=mpi_enreg%comm_fft; reduce=.true.
1450  spaceComm=mpi_enreg%comm_fft; reduce=.false.
1451  if(mpi_enreg%nproc_fft>1) then
1452    spaceComm=mpi_enreg%comm_fft; reduce=.true.
1453  else if(mpi_enreg%nproc_wvl>1) then
1454    spaceComm=mpi_enreg%comm_wvl; reduce=.true.
1455  end if
1456  nproc=xmpi_comm_size(spaceComm)
1457  me=xmpi_comm_rank(spaceComm)
1458 
1459  if (reduce) then
1460 
1461 !  Communicate all data to all processors with only two global communications
1462    ABI_MALLOC(value_fft,(5,nitems,nproc))
1463    ABI_MALLOC(index_fft,(2,2,nitems,nproc))
1464    value_fft(:,:,:)=zero
1465    index_fft(:,:,:,:)=0
1466    value_fft(1,:,me + 1)=value(1,1,:)
1467    value_fft(2,:,me + 1)=value(2,1,:)
1468    value_fft(3,:,me + 1)=value(1,2,:)
1469    value_fft(4,:,me + 1)=value(2,2,:)
1470    if(present(ucvol))value_fft(5,:,me + 1)=integrated(:)
1471    index_fft(:,:,:,me + 1)=iindex(:,:,:)
1472    call xmpi_sum(value_fft,spaceComm,ierr)
1473    call xmpi_sum(index_fft,spaceComm,ierr)
1474 
1475 !  Determine the global optimum and second optimum for each item
1476 !  Also, the integrated quantities, if relevant.
1477    do iitems=1,nitems
1478 
1479      if(present(ucvol))integrated(iitems)=sum(value_fft(5,iitems,1:nproc))
1480 
1481      do indsign=1,2 ! Find alternatively the maximum and the minimum
1482        iisign=3-2*indsign
1483 
1484 !      Initialisation
1485        value1=value_fft(2*indsign-1,iitems,1)
1486        value2=value_fft(2*indsign  ,iitems,1)
1487        index1=index_fft(1,indsign,iitems,1)
1488        index2=index_fft(2,indsign,iitems,1)
1489 
1490 !      Loop
1491        do iproc=1, nproc, 1
1492          do ii=1,2
1493            if(iproc>1 .or. ii==2)then
1494 
1495              temp=value_fft(ii+2*(indsign-1),iitems,iproc)
1496 !            Compares it to the second value
1497              if( iisign*(temp+tol12) > iisign*value2 ) then
1498 !              Compare it to the first value
1499                if( iisign*(temp+tol12) > iisign*value1 ) then
1500                  value2=value1 ; index2=index1
1501                  value1=temp   ; index1=index_fft(ii,indsign,iitems,iproc)
1502                else
1503                  value2=temp   ; index2=index_fft(ii,indsign,iitems,iproc)
1504                end if
1505              end if
1506 
1507            end if ! if(iproc>1 .or. ii==2)
1508          end do ! ii
1509        end do ! iproc
1510 
1511        value(1,indsign,iitems)=value1
1512        value(2,indsign,iitems)=value2
1513        iindex(1,indsign,iitems)=index1
1514        iindex(2,indsign,iitems)=index2
1515 
1516      end do ! iisign
1517    end do ! iitems
1518 
1519    ABI_FREE(value_fft)
1520    ABI_FREE(index_fft)
1521 
1522  end if !if(reduce)
1523 
1524 !-------------------------------------------------------------------
1525 
1526 !Determines the reduced coordinates of the min and max for each item
1527  ABI_MALLOC(coord,(3,2,2,nitems))
1528  do iitems=1,nitems
1529    do indsign=1,2
1530      do ii=1,2
1531        index1=iindex(ii,indsign,iitems)
1532        i3=(index1-1)/n1/n2
1533        i2=(index1-1-i3*n1*n2)/n1
1534        i1=index1-1-i3*n1*n2-i2*n1
1535        coord(1,ii,indsign,iitems)=dble(i1)/dble(n1)+tol12
1536        coord(2,ii,indsign,iitems)=dble(i2)/dble(n2)+tol12
1537        coord(3,ii,indsign,iitems)=dble(i3)/dble(n3)+tol12
1538 !      DEBUG
1539 !      write(std_out,*)' ii,indsign,iitems,coord(1:3)=',ii,indsign,iitems,coord(:,ii,indsign,iitems)
1540 !      write(std_out,*)' value ', value(ii, indsign, iitems)
1541 !      ENDDEBUG
1542      end do
1543    end do
1544  end do
1545 
1546 !-------------------------------------------------------------------------
1547 !Output
1548  if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then
1549    if(.true.)then
1550      do iitems=1,nitems
1551 
1552        if(ioptrhor==4 .and. iitems>2)exit
1553 
1554        select case (ioptrhor)
1555        case(0)
1556 
1557          if(iitems==1) write(message,'(a)')' Total charge density [el/Bohr^3]'
1558          if(nspden==2)then
1559            if(iitems==2) write(message,'(a)')' Spin up density      [el/Bohr^3]'
1560            if(iitems==3) write(message,'(a)')' Spin down density    [el/Bohr^3]'
1561            if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^3]'
1562            if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1563          else if(nspden==4)then
1564            if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^3]'
1565            if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^3]'
1566            if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^3]'
1567            if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^3]'
1568            if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1569          end if
1570 
1571        case(1)
1572 
1573          if(iitems==1) write(message,'(a)')' Total kinetic energy density [Ha/Bohr^3]'
1574          if(nspden==2)then
1575            if(iitems==2) write(message,'(a)')' Spin up density      [Ha/Bohr^3]'
1576            if(iitems==3) write(message,'(a)')' Spin down density    [Ha/Bohr^3]'
1577            if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [Ha/Bohr^3]'
1578            if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1579          else if(nspden==4)then
1580            if(iitems==2) write(message,'(a)')' x component of magnetization [Ha/Bohr^3]'
1581            if(iitems==3) write(message,'(a)')' y component of magnetization [Ha/Bohr^3]'
1582            if(iitems==4) write(message,'(a)')' z component of magnetization [Ha/Bohr^3]'
1583            if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [Ha/Bohr^3]'
1584            if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1585          end if
1586 
1587        case(2)
1588 
1589          if(iitems==1) write(message,'(a)')' Gradient of the electronic density [el/Bohr^4]'
1590          if(nspden==2)then
1591            if(iitems==2) write(message,'(a)')' Spin up density      [el/Bohr^4]'
1592            if(iitems==3) write(message,'(a)')' Spin down density    [el/Bohr^4]'
1593            if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]'
1594            if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1595          else if(nspden==4)then
1596            if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^4]'
1597            if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^4]'
1598            if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^4]'
1599            if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^4]'
1600            if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1601          end if
1602 
1603        case(3)
1604 
1605          if(iitems==1) write(message,'(a)')' Laplacian of the electronic density [el/Bohr^5]'
1606          if(nspden==2)then
1607            if(iitems==2) write(message,'(a)')' Spin up density      [el/Bohr^5]'
1608            if(iitems==3) write(message,'(a)')' Spin down density    [el/Bohr^5]'
1609            if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^5]'
1610            if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1611          else if(nspden==4)then
1612            if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^5]'
1613            if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^5]'
1614            if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^5]'
1615            if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^5]'
1616            if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1617          end if
1618 
1619        case(4)
1620 
1621          if(iitems==1) write(message,'(a)')' Electron Localization Function (ELF) [min:0;max:1]'
1622          if(nspden==2)then
1623            if(iitems==2) write(message,'(a)')' Spin up ELF      [min:0;max:1]'
1624 !            if(iitems==3) write(message,'(a)')' Spin down ELF    [min:0;max:1]'
1625 !            if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]'
1626 !            if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1627          else if(nspden==4)then
1628 !            if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^4]'
1629 !            if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^4]'
1630 !            if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^4]'
1631 !            if(iitems==5) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]'
1632 !            if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1)   '
1633          end if
1634        end select
1635 
1636        call wrtout(iout,message,'COLL')
1637 
1638        write(message,'(a,es13.4,a,3f10.4)')   ')     Maximum= ',&
1639 &       value(1,1,iitems),'  at reduced coord.',coord(:,1,1,iitems)
1640        call wrtout(iout,message,'COLL')
1641        if(option==2)then
1642          write(message,'(a,es13.4,a,3f10.4)') ')Next maximum= ',&
1643 &         value(2,1,iitems),'  at reduced coord.',coord(:,2,1,iitems)
1644          call wrtout(iout,message,'COLL')
1645        end if
1646        write(message,'(a,es13.4,a,3f10.4)')   ')     Minimum= ',&
1647 &       value(1,2,iitems),'  at reduced coord.',coord(:,1,2,iitems)
1648        call wrtout(iout,message,'COLL')
1649        if(option==2)then
1650          write(message,'(a,es13.4,a,3f10.4)') ')Next minimum= ',&
1651 &         value(2,2,iitems),'  at reduced coord.',coord(:,2,2,iitems)
1652          call wrtout(iout,message,'COLL')
1653        end if
1654        if(present(ucvol))then
1655          if(.not.(nspden==2.and.iitems==5) .and. .not.(nspden==4.and.iitems==6))then
1656            if(abs(integrated(iitems))<tol10)integrated(iitems)=zero
1657            write(message,'(a,es13.4)')'   Integrated= ',integrated(iitems)
1658            call wrtout(iout,message,'COLL')
1659          end if
1660        end if
1661 
1662      end do ! iitems
1663    end if
1664 
1665    if(.false.)then
1666 
1667      select case(optrhor)
1668      case(0)
1669        write(txt1_in_mssg, '(a)')" Min el dens="
1670        write(txt2_in_mssg, '(a)')" el/bohr^3 at reduced coord."
1671        write(txt3_in_mssg, '(a)')" Max el dens="
1672      case(1)
1673        write(txt1_in_mssg, '(a)')" Min kin energy dens="
1674        write(txt2_in_mssg, '(a)')" bohr^(-5) at reduced coord."
1675        write(txt3_in_mssg, '(a)')" Max kin energy dens="
1676      end select
1677 
1678      write(message, '(a,a,1p,e12.4,a,0p,3f8.4)' ) ch10,&
1679 &     trim(txt1_in_mssg),value(1,2,1),&
1680 &     trim(txt2_in_mssg),coord(:,1,2,1)
1681      call wrtout(iout,message,'COLL')
1682      if(option==2)then
1683        write(message, '(a,1p,e12.4,a,0p,3f8.4)' ) &
1684 &       ',   next min=',value(2,2,1),&
1685 &       trim(txt2_in_mssg),coord(:,2,2,1)
1686        call wrtout(iout,message,'COLL')
1687      end if
1688      write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
1689 &     trim(txt3_in_mssg),value(1,1,1),&
1690 &     trim(txt2_in_mssg),coord(:,1,1,1)
1691      call wrtout(iout,message,'COLL')
1692      if(option==2)then
1693        write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
1694 &       ',   next max=',value(2,1,1),&
1695 &       trim(txt2_in_mssg),coord(:,2,1,1)
1696        call wrtout(iout,message,'COLL')
1697      end if
1698 
1699      if(nspden>=2)then
1700        write(message, '(a,a,1p,e12.4,a,0p,3f8.4)' ) ch10,&
1701 &       ',Min spin pol zeta=',value(1,2,4+nspden/2),&
1702 &       ' at reduced coord.',coord(:,1,2,4+nspden/2)
1703        call wrtout(iout,message,'COLL')
1704        if(option==2)then
1705          write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
1706 &         ',         next min=',value(2,2,4+nspden/2),&
1707 &         ' at reduced coord.',coord(:,2,2,4+nspden/2)
1708          call wrtout(iout,message,'COLL')
1709        end if
1710        write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
1711 &       ',Max spin pol zeta=',value(1,1,4+nspden/2),&
1712 &       ' at reduced coord.',coord(:,1,1,4+nspden/2)
1713        call wrtout(iout,message,'COLL')
1714        if(option==2)then
1715          write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
1716 &         ',         next max=',value(2,1,4+nspden/2),&
1717 &         ' at reduced coord.',coord(:,2,1,4+nspden/2)
1718          call wrtout(iout,message,'COLL')
1719        end if
1720      end if ! nspden
1721 
1722    end if ! second section always true
1723 
1724    if(nspden==2 .and. .false.)then
1725      write(message,'(a)')&
1726 &     '                               Position in reduced coord.       (  x         y         z )'
1727      call wrtout(iout,message,'COLL')
1728      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Total  el-den) : [el/Bohr^3]',&
1729 &     rhomn1(1),'  at',ri_rhomn1(1,1),ri_rhomn1(2,1),ri_rhomn1(3,1)
1730      call wrtout(iout,message,'COLL')
1731      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Spin-up   den) : [el/Bohr^3]',&
1732 &     rhomn1(2),'  at',ri_rhomn1(1,2),ri_rhomn1(2,2),ri_rhomn1(3,2)
1733      call wrtout(iout,message,'COLL')
1734      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Spin-down den) : [el/Bohr^3]',&
1735 &     zetmn1(1),'  at',ri_zetmn1(1,1),ri_zetmn1(2,1),ri_zetmn1(3,1)
1736      call wrtout(iout,message,'COLL')
1737      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Spin pol zeta) :   [m/|m|]  ',&
1738 &     zetmn1(2),'  at',ri_zetmn1(1,2),ri_zetmn1(2,2),ri_zetmn1(3,2)
1739      call wrtout(iout,message,'COLL')
1740      if(option==2)then
1741        write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Total  el-den) : [el/Bohr^3]',&
1742 &       rhomn2(1),'  at',ri_rhomn2(1,1),ri_rhomn2(2,1),ri_rhomn2(3,1)
1743        call wrtout(iout,message,'COLL')
1744        write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin-up   den) : [el/Bohr^3]',&
1745 &       rhomn2(2),'  at',ri_rhomn2(1,2),ri_rhomn2(2,2),ri_rhomn2(3,2)
1746        call wrtout(iout,message,'COLL')
1747        write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin-down den) : [el/Bohr^3]',&
1748 &       zetmn2(1),'  at',ri_zetmn2(1,1),ri_zetmn2(2,1),ri_zetmn2(3,1)
1749        call wrtout(iout,message,'COLL')
1750        write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin pol zeta) :   [m/|m|]  ',&
1751 &       zetmn2(2),'  at',ri_zetmn2(1,2),ri_zetmn2(2,2),ri_zetmn2(3,2)
1752        call wrtout(iout,message,'COLL')
1753      end if
1754      write(message,*)' '
1755      call wrtout(iout,message,'COLL')
1756      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Total  el-den) : [el/Bohr^3]',&
1757 &     rhomx1(1),'  at',ri_rhomx1(1,1),ri_rhomx1(2,1),ri_rhomx1(3,1)
1758      call wrtout(iout,message,'COLL')
1759      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Spin-up   den) : [el/Bohr^3]',&
1760 &     rhomx1(2),'  at',ri_rhomx1(1,2),ri_rhomx1(2,2),ri_rhomx1(3,2)
1761      call wrtout(iout,message,'COLL')
1762      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Spin-down den) : [el/Bohr^3]',&
1763 &     zetmx1(1),'  at',ri_zetmx1(1,1),ri_zetmx1(2,1),ri_zetmx1(3,1)
1764      call wrtout(iout,message,'COLL')
1765      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Spin pol zeta) :   [m/|m|]  ',&
1766 &     zetmx1(2),'  at',ri_zetmx1(1,2),ri_zetmx1(2,2),ri_zetmx1(3,2)
1767      call wrtout(iout,message,'COLL')
1768      if(option==2)then
1769        write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Total  el-den) : [el/Bohr^3]',&
1770 &       rhomx2(1),'  at',ri_rhomx2(1,1),ri_rhomx2(2,1),ri_rhomx2(3,1)
1771        call wrtout(iout,message,'COLL')
1772        write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin-up   den) : [el/Bohr^3]',&
1773 &       rhomx2(2),'  at',ri_rhomx2(1,2),ri_rhomx2(2,2),ri_rhomx2(3,2)
1774        call wrtout(iout,message,'COLL')
1775        write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin-down den) : [el/Bohr^3]',&
1776 &       zetmx2(1),'  at',ri_zetmx2(1,1),ri_zetmx2(2,1),ri_zetmx2(3,1)
1777        call wrtout(iout,message,'COLL')
1778        write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin pol zeta) :   [m/|m|]  ',&
1779 &       zetmx2(2),'  at',ri_zetmx2(1,2),ri_zetmx2(2,2),ri_zetmx2(3,2)
1780        call wrtout(iout,message,'COLL')
1781      end if
1782    end if
1783 
1784    if(nspden==4 .and. .false.)then
1785      write(message,'(a)')&
1786 &     '                               Position in reduced coord.       (  x         y         z )'
1787      call wrtout(iout,message,'COLL')
1788      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Total  el-den) : [el/Bohr^3]',&
1789 &     rhomn1(1),'  at',ri_rhomn1(1,1),ri_rhomn1(2,1),ri_rhomn1(3,1)
1790      call wrtout(iout,message,'COLL')
1791      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Magnetizat.-x) :   [m/|m|]  ',&
1792 &     rhomn1(2),'  at',ri_rhomn1(1,2),ri_rhomn1(2,2),ri_rhomn1(3,2)
1793      call wrtout(iout,message,'COLL')
1794      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Magnetizat.-y) :   [m/|m|]  ',&
1795 &     rhomn1(3),'  at',ri_rhomn1(1,3),ri_rhomn1(2,3),ri_rhomn1(3,3)
1796      call wrtout(iout,message,'COLL')
1797      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Magnetizat.-z) :   [m/|m|]  ',&
1798 &     rhomn1(4),'  at',ri_rhomn1(1,4),ri_rhomn1(2,4),ri_rhomn1(3,4)
1799      call wrtout(iout,message,'COLL')
1800      write(message,'(a,es13.4,a,3f10.4)')'      Minimum (Spin pol zeta) :   [m/|m|]  ',&
1801 &     zetmn1(1),'  at',ri_zetmn1(1,1),ri_zetmn1(2,1),ri_zetmn1(3,1)
1802      call wrtout(iout,message,'COLL')
1803      if(option==2)then
1804        write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Total  el-den) : [el/Bohr^3]',&
1805 &       rhomn2(1),'  at',ri_rhomn2(1,1),ri_rhomn2(2,1),ri_rhomn2(3,1)
1806        call wrtout(iout,message,'COLL')
1807        write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-x) :   [m/|m|]  ',&
1808 &       rhomn2(2),'  at',ri_rhomn2(1,2),ri_rhomn2(2,2),ri_rhomn2(3,2)
1809        call wrtout(iout,message,'COLL')
1810        write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-y) :   [m/|m|]  ',&
1811 &       rhomn2(3),'  at',ri_rhomn2(1,3),ri_rhomn2(2,3),ri_rhomn2(3,3)
1812        call wrtout(iout,message,'COLL')
1813        write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-z) :   [m/|m|]  ',&
1814 &       rhomn2(4),'  at',ri_rhomn2(1,4),ri_rhomn2(2,4),ri_rhomn2(3,4)
1815        call wrtout(iout,message,'COLL')
1816        write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Spin pol zeta) :   [m/|m|]  ',&
1817 &       zetmn2(1),'  at',ri_zetmn2(1,1),ri_zetmn2(2,1),ri_zetmn2(3,1)
1818        call wrtout(iout,message,'COLL')
1819      end if
1820      write(message,*)' '
1821      call wrtout(iout,message,'COLL')
1822      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Total  el-den) : [el/Bohr^3]',&
1823 &     rhomx1(1),'  at',ri_rhomx1(1,1),ri_rhomx1(2,1),ri_rhomx1(3,1)
1824      call wrtout(iout,message,'COLL')
1825      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Magnetizat.-x) :   [m/|m|]  ',&
1826 &     rhomx1(2),'  at',ri_rhomx1(1,2),ri_rhomx1(2,2),ri_rhomx1(3,2)
1827      call wrtout(iout,message,'COLL')
1828      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Magnetizat.-y) :   [m/|m|]  ',&
1829 &     rhomx1(3),'  at',ri_rhomx1(1,3),ri_rhomx1(2,3),ri_rhomx1(3,3)
1830      call wrtout(iout,message,'COLL')
1831      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Magnetizat.-z) :   [m/|m|]  ',&
1832 &     rhomx1(4),'  at',ri_rhomx1(1,4),ri_rhomx1(2,4),ri_rhomx1(3,4)
1833      call wrtout(iout,message,'COLL')
1834      write(message,'(a,es13.4,a,3f10.4)')'      Maximum (Spin pol zeta) :   [m/|m|]  ',&
1835 &     zetmx1(1),'  at',ri_zetmx1(1,1),ri_zetmx1(2,1),ri_zetmx1(3,1)
1836      call wrtout(iout,message,'COLL')
1837      if(option==2)then
1838        write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Total  el-den) : [el/Bohr^3]',&
1839 &       rhomx2(1),'  at',ri_rhomx2(1,1),ri_rhomx2(2,1),ri_rhomx2(3,1)
1840        call wrtout(iout,message,'COLL')
1841        write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-x) :   [m/|m|]  ',&
1842 &       rhomx2(2),'  at',ri_rhomx2(1,2),ri_rhomx2(2,2),ri_rhomx2(3,2)
1843        call wrtout(iout,message,'COLL')
1844        write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-y) :   [m/|m|]  ',&
1845 &       rhomx2(3),'  at',ri_rhomx2(1,3),ri_rhomx2(2,3),ri_rhomx2(3,3)
1846        call wrtout(iout,message,'COLL')
1847        write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-z) :   [m/|m|]  ',&
1848 &       rhomx2(4),'  at',ri_rhomx2(1,4),ri_rhomx2(2,4),ri_rhomx2(3,4)
1849        call wrtout(iout,message,'COLL')
1850        write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Spin pol zeta) :   [m/|m|]  ',&
1851 &       zetmx2(1),'  at',ri_zetmx2(1,1),ri_zetmx2(2,1),ri_zetmx2(3,1)
1852        call wrtout(iout,message,'COLL')
1853      end if
1854    end if
1855  end if
1856 
1857  ABI_FREE(coord)
1858  ABI_FREE(value)
1859  ABI_FREE(iindex)
1860  ABI_FREE(integrated)
1861 
1862 end subroutine prtrhomxmn

m_mkrho/read_atomden [ Functions ]

[ Top ] [ m_mkrho ] [ Functions ]

NAME

 read_atomden

FUNCTION

COPYRIGHT

 Copyright (C) 2005-2022 ABINIT group (SM,VR,FJ,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 .

INPUTS

 natom : number of atoms in cell
 nfft=(effective) number of FFT grid points (for this processor) - fine grid
 ngfft(18)=contain all needed information about 3D FFT,
 nspden : number of spin densities
 ntypat : number of types of atoms in the cell
 typat(natom) : list of atom types

OUTPUT

 rhor_atm(nfft,nspden) : full electron density on the (fine) grid

SIDE EFFECTS

NOTES

SOURCE

1894 subroutine read_atomden(MPI_enreg,natom,nfft,ngfft,nspden,ntypat, &
1895 &                       rhor_atm,typat,rprimd,xred,prtvol,file_prefix)
1896 
1897 !Arguments ------------------------------------
1898 !scalars
1899  integer,intent(in) :: natom,nfft,nspden,ntypat,prtvol
1900 !arrays
1901  type(MPI_type),intent(in) :: MPI_enreg
1902  integer,intent(in) :: ngfft(18),typat(natom)
1903  real(dp), intent(in) :: rprimd(3,3),xred(3,natom)
1904  real(dp),intent(inout) :: rhor_atm(nfft,nspden)
1905  character(len=7), intent(in) :: file_prefix
1906 
1907 !Local variables-------------------------------
1908 !scalars
1909  character(len=500) :: message
1910  character(len=120) :: filename
1911  character(len=7) :: calctype='replace'
1912  integer :: igrid,i,i1,i2,i3,io_err,itypat,unt
1913  integer :: natomgrmax,nlines,ngrid,n1,n2,n3
1914  real(dp) :: difx,dify,difz,ucvol!,norm
1915 !arrays
1916  integer :: natomgr(ntypat)
1917  real(dp) :: a(3),b(3),c(3)
1918  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1919  real(dp),allocatable :: atomrgrid(:,:),r_vec_grid(:,:),density(:,:)
1920  real(dp),allocatable :: rho(:)
1921 
1922 ! ************************************************************************
1923 
1924 !Initialise various variables
1925  ngrid = nfft
1926  a(:) = rprimd(:,1)
1927  b(:) = rprimd(:,2)
1928  c(:) = rprimd(:,3)
1929  ABI_MALLOC(rho,(ngrid))
1930  if (nspden/=1) then
1931    ABI_ERROR('read_atomden: Only nspden=1 allowed.')
1932  end if
1933  rho = rhor_atm(:,1)
1934  gmet=zero;gprimd=zero;rmet=zero;ucvol=zero
1935 
1936 
1937 !Calculate the r vector (reduced coord.) of the fine gridpoints
1938  ABI_MALLOC(r_vec_grid,(3,ngrid))
1939  igrid = 0
1940  n1 = ngfft(1)
1941  n2 = ngfft(2)
1942  n3 = ngfft(3)
1943  do i3=0,n3-1
1944    difz=dble(i3)/dble(n3)
1945    do i2=0,n2-1
1946      dify=dble(i2)/dble(n2)
1947      do i1=0,n1-1
1948        difx=dble(i1)/dble(n1)
1949        igrid = igrid + 1
1950        r_vec_grid(1,igrid)=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
1951        r_vec_grid(2,igrid)=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1952        r_vec_grid(3,igrid)=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1953      end do
1954    end do
1955  end do
1956  if (igrid/=ngrid) then
1957    ABI_ERROR('read_atomden: igrid not equal to ngrid')
1958  end if
1959 
1960 !Read in atomic density data for each atom type
1961 !first check how many datapoints are in each file
1962  do itypat=1,ntypat
1963    filename='';io_err=0;
1964    if (itypat>0)  write(filename,'(a,a,i1,a)') trim(file_prefix), &
1965 &   '_density_atom_type',itypat,'.dat'
1966    if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), &
1967 &   '_density_atom_type',itypat,'.dat'
1968    if (open_file(filename, message, newunit=unt, status='old',action='read') /= 0) then
1969      write(std_out,*) 'ERROR in read_atomden: Could not open file: ',filename
1970      write(std_out,*) ' Current implementation requires this file to be present'
1971      write(std_out,*) ' for each type of atom.'
1972      write(std_out,*)trim(message)
1973      ABI_ERROR("Cannot continue")
1974    end if
1975 !  Check number of lines in file
1976    nlines = 1;io_err=0;
1977    do
1978      read(unt,*,iostat=io_err)
1979      if (io_err<0) exit
1980      nlines = nlines + 1
1981    end do
1982    close(unt)
1983    natomgr(itypat) = nlines - 2
1984  end do ! Atom type
1985 !Allocate arrays and read in data
1986  natomgrmax = maxval(natomgr)
1987  ABI_MALLOC(atomrgrid,(natomgrmax,ntypat))
1988  ABI_MALLOC(density,(natomgrmax,ntypat))
1989  atomrgrid = zero ; density = zero
1990  do itypat=1,ntypat
1991    filename='';io_err=0;
1992    if (itypat>0)  write(filename,'(a,a,i1,a)') trim(file_prefix), &
1993 &   '_density_atom_type',itypat,'.dat'
1994    if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), &
1995 &   '_density_atom_type',itypat,'.dat'
1996    if (open_file(filename,message,newunit=unt,status='old',action='read') /= 0) then
1997      ABI_ERROR(message)
1998    end if
1999    read(unt,*) ! Skip comment line
2000    do i=1,natomgr(itypat)
2001      read(unt,*) atomrgrid(i,itypat),density(i,itypat)
2002    end do
2003    close(unt)
2004    if (atomrgrid(1,itypat)/=zero) then
2005      write(std_out,*) 'ERROR in read_atomden, in file: ',filename
2006      write(std_out,*) ' First gridpoint has to be the origin.'
2007      ABI_ERROR("Cannot continue")
2008    end if
2009  end do ! Atom type
2010 
2011 !write(std_out,*) '*** --- In read_atomden before call--- ***'
2012 !write(std_out,*) '  calctype:',calctype,' natom:',natom
2013 !write(std_out,*) '    ntypat:',ntypat,' typat:',typat
2014 !write(std_out,*) '     ngrid:',ngrid
2015 !write(std_out,*) '         a:',a
2016 !write(std_out,*) '         b:',b
2017 !write(std_out,*) '         c:',c
2018 !write(std_out,*) '      xred:',xred
2019 !write(std_out,*) '   natomgr:',natomgr
2020 !write(std_out,*) 'natomgrmax:',natomgrmax
2021 !write(std_out,*) ' atomrgrid:',atomrgrid
2022 !write(std_out,*) '   density:',density
2023 !write(std_out,*) 'r_vec_grid:'
2024 !write(std_out,*) r_vec_grid
2025 
2026 !Call atomden
2027  call atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,xred, &
2028 & natomgr,natomgrmax,atomrgrid,density,prtvol,calctype)
2029 
2030 !if (prtvol>9) then ! calculate norm
2031 !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2032 !norm = SUM(rho(:))*ucvol/dble(n1*n2*n3)
2033 !write(message,'(a,F8.4)') '  In read_atomden - NORM OF DENSITY: ',norm
2034 !call wrtout(std_out,message,'COLL')
2035 !end if
2036 
2037  rhor_atm(:,1) = rho
2038 
2039  if (allocated(atomrgrid))  then
2040    ABI_FREE(atomrgrid)
2041  end if
2042  if (allocated(density))  then
2043    ABI_FREE(density)
2044  end if
2045  if (allocated(r_vec_grid))  then
2046    ABI_FREE(r_vec_grid)
2047  end if
2048  if (allocated(rho))  then
2049    ABI_FREE(rho)
2050  end if
2051 
2052  return
2053 
2054 end subroutine read_atomden