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-2018 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)

PARENTS

      read_atomden

CHILDREN

      sort_dp,spline,splint,wrtout,xmpi_barrier,xmpi_sum_master

SOURCE

2179 subroutine atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,atom_pos, &
2180 &                  natomgr,natomgrmax,atomrgrid,density,prtvol,calctype)
2181 
2182 
2183 !This section has been created automatically by the script Abilint (TD).
2184 !Do not modify the following lines by hand.
2185 #undef ABI_FUNC
2186 #define ABI_FUNC 'atomden'
2187 !End of the abilint section
2188 
2189  implicit none
2190 
2191 !Arguments ------------------------------------
2192 !scalars
2193  integer,intent(in) :: natom,ntypat,ngrid,natomgrmax,prtvol
2194  character(len=7),intent(in) :: calctype
2195 !arrays
2196  type(MPI_type),intent(in) :: MPI_enreg
2197  integer,intent(in) :: typat(natom),natomgr(ntypat)
2198  real(dp),intent(in) :: r_vec_grid(3,ngrid),a(3),b(3),c(3)
2199  real(dp),intent(in) :: atom_pos(3,natom),atomrgrid(natomgrmax,ntypat)
2200  real(dp),intent(in) :: density(natomgrmax,ntypat)
2201  real(dp),intent(inout) :: rho(ngrid)
2202 
2203 !Local variables-------------------------------
2204 !scalars
2205  character(len=500) :: message
2206  integer :: cnt,delta,i,l,m,n,iatom,itypat,igrid,ncells,n_grid_p
2207  integer :: ierr,spaceComm,nprocs,master,rank,remainder
2208  real(dp) :: a_norm,b_norm,c_norm
2209  real(dp) :: r_max,R_sphere_max,dp_dummy,ybcbeg,ybcend
2210 !arrays
2211  integer :: n_equiv_atoms(ntypat),grid_index(ngrid)
2212  integer :: my_start_equiv_atoms(ntypat)
2213  integer :: my_end_equiv_atoms(ntypat)
2214  integer :: l_min(ntypat),m_min(ntypat),n_min(ntypat)
2215  integer :: l_max(ntypat),m_max(ntypat),n_max(ntypat)
2216  real(dp) :: center(3),dp_vec_dummy(3),delta_a(3),delta_b(3),delta_c(3)
2217  real(dp) :: r_atom(3),grid_distances(ngrid)
2218  integer, allocatable :: new_index(:),i_1d_dummy(:)
2219  real(dp),allocatable :: equiv_atom_dist(:,:),equiv_atom_pos(:,:,:),rho_temp(:,:)
2220  real(dp),allocatable :: dp_1d_dummy(:),dp_2d_dummy(:,:),ypp(:)
2221  real(dp),allocatable :: x_fit(:),y_fit(:)
2222 
2223 
2224 ! ************************************************************************
2225 
2226 !initialise and check parallel execution
2227  spaceComm=MPI_enreg%comm_cell
2228  nprocs=xmpi_comm_size(spaceComm)
2229  rank=MPI_enreg%me_kpt
2230 
2231  master=0
2232 
2233 !initialise variables and vectors
2234  a_norm = sqrt(dot_product(a,a))
2235  b_norm = sqrt(dot_product(b,b))
2236  c_norm = sqrt(dot_product(c,c))
2237  center = (a+b+c)*half
2238  dp_dummy = dot_product(a,b)/(b_norm*b_norm)
2239  dp_vec_dummy = dp_dummy*b
2240  delta_a = a - dp_vec_dummy
2241  dp_dummy = dot_product(b,a)/(a_norm*a_norm)
2242  dp_vec_dummy = dp_dummy*a
2243  delta_b = b - dp_vec_dummy
2244  dp_dummy = dot_product(c,(a+b))/(dot_product((a+b),(a+b)))
2245  dp_vec_dummy = dp_dummy*(a+b)
2246  delta_c = c - dp_vec_dummy
2247  ABI_ALLOCATE(rho_temp,(ngrid,ntypat))
2248  rho_temp = zero
2249 
2250 !write(std_out,*) '*** --- In atomden --- ***'
2251 !write(std_out,*) ' a_norm:',a_norm,' b_norm:',b_norm,' c_norm:',c_norm
2252 !write(std_out,*) 'delta_a:',delta_a,'delta_b:',delta_b,'delta_c:',delta_c
2253 !write(std_out,*) ' center:',center
2254 
2255 !Find supercell which will contain all possible contributions
2256 !for all atoms, and enumerate positions for all atoms
2257 !TODO list of atoms can be "pruned", i.e identify all atoms
2258 !that can't possibly contribute and remove from list.
2259 !Should be most important for very oblique cells
2260  do itypat=1,ntypat
2261    R_sphere_max = atomrgrid(natomgr(itypat),itypat)
2262    l_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_a,delta_a)))
2263    l_max(itypat) = -l_min(itypat)
2264    m_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_b,delta_b)))
2265    m_max(itypat) = -m_min(itypat)
2266    n_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_c,delta_c)))
2267    n_max(itypat) = -n_min(itypat)
2268    ncells = (l_max(itypat)-l_min(itypat)+1) &
2269 &   *(m_max(itypat)-m_min(itypat)+1) &
2270 &   *(n_max(itypat)-n_min(itypat)+1)
2271    n_equiv_atoms(itypat) = 0
2272    do iatom=1,natom
2273      if (typat(iatom)==itypat) then
2274        n_equiv_atoms(itypat) = n_equiv_atoms(itypat) + ncells
2275      end if ! if type=itypat
2276    end do ! number of atoms per cell
2277    if ((rank==master).and.(prtvol>9)) then
2278      write(message,'(a)') '*** --- In atomden --- find box ***'
2279      call wrtout(std_out,message,'COLL')
2280      write(message,'(a,I4)') ' itypat:',itypat
2281      call wrtout(std_out,message,'COLL')
2282      write(message,'(2(a,I4))') ' l_min:',l_min(itypat),' l_max:',l_max(itypat)
2283      call wrtout(std_out,message,'COLL')
2284      write(message,'(2(a,I4))') ' m_min:',m_min(itypat),' m_max:',m_max(itypat)
2285      call wrtout(std_out,message,'COLL')
2286      write(message,'(2(a,I4))') ' n_min:',n_min(itypat),' n_max:',n_max(itypat)
2287      call wrtout(std_out,message,'COLL')
2288      write(message,'(2(a,I4))') ' n_equiv_atoms:',n_equiv_atoms(itypat)
2289      call wrtout(std_out,message,'COLL')
2290    end if
2291  end do !atom type
2292 
2293 !allocate arrays
2294  n = maxval(n_equiv_atoms)
2295  ABI_ALLOCATE(equiv_atom_pos,(3,n,ntypat))
2296  ABI_ALLOCATE(equiv_atom_dist,(n,ntypat))
2297  equiv_atom_pos = zero
2298  equiv_atom_dist = zero
2299 
2300 !Find positions and distance of atoms from center of cell
2301  do itypat=1,ntypat
2302    i = 1
2303    do l=l_min(itypat),l_max(itypat)
2304      do m=m_min(itypat),m_max(itypat)
2305        do n=n_min(itypat),n_max(itypat)
2306          do iatom=1,natom
2307            if (typat(iatom)==itypat) then
2308              if (i>n_equiv_atoms(itypat)) then
2309                MSG_ERROR('atomden: i>n_equiv_atoms')
2310              end if
2311              equiv_atom_pos(:,i,itypat) = (atom_pos(1,iatom)+dble(l))*a &
2312 &             + (atom_pos(2,iatom)+dble(m))*b &
2313 &             + (atom_pos(3,iatom)+dble(n))*c
2314              dp_vec_dummy = equiv_atom_pos(:,i,itypat)-center
2315              equiv_atom_dist(i,itypat) = &
2316 &             sqrt(dot_product(dp_vec_dummy,dp_vec_dummy))
2317              i = i + 1
2318            end if
2319          end do
2320        end do !n
2321      end do !m
2322    end do !l
2323 !  write(std_out,*) '*** --- In atomden --- find equiv ***'
2324 !  write(std_out,*) ' itypat:',itypat
2325 !  write(std_out,*) ' equiv_atom_pos:'
2326 !  write(std_out,*) equiv_atom_pos(:,:,itypat)
2327 !  write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat)
2328  end do !atom type
2329 
2330 !Sort the atoms after distance so that the density from the ones
2331 !furthest away can be added first. This is to prevent truncation error.
2332  do itypat=1,ntypat
2333    n = n_equiv_atoms(itypat)
2334    ABI_ALLOCATE(dp_1d_dummy,(n))
2335    ABI_ALLOCATE(new_index,(n))
2336    ABI_ALLOCATE(dp_2d_dummy,(3,n))
2337    dp_1d_dummy = equiv_atom_dist(1:n,itypat)
2338    dp_2d_dummy = equiv_atom_pos(1:3,1:n,itypat)
2339    do i=1,n
2340      new_index(i) = i
2341    end do
2342    call sort_dp(n,dp_1d_dummy,new_index,tol14)
2343    do i=1,n
2344 !    write(std_out,*) i,' -> ',new_index(i)
2345      equiv_atom_pos(1:3,n+1-i,itypat) = dp_2d_dummy(1:3,new_index(i))
2346      equiv_atom_dist(1:n,itypat) = dp_1d_dummy
2347    end do
2348    ABI_DEALLOCATE(dp_1d_dummy)
2349    ABI_DEALLOCATE(new_index)
2350    ABI_DEALLOCATE(dp_2d_dummy)
2351 !  write(std_out,*) '*** --- In atomden ---  sorting atoms ***'
2352 !  write(std_out,*) ' itypat:',itypat
2353 !  write(std_out,*) ' equiv_atom_pos:'
2354 !  write(std_out,*) equiv_atom_pos(:,:,itypat)
2355 !  write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat)
2356  end do ! atom type
2357 
2358 !Divide the work in case of parallel execution
2359  if (nprocs==1) then ! Make sure everything runs with one proc
2360    if (prtvol>9) then
2361      write(message,'(a)') '  In atomden - number of processors:     1'
2362      call wrtout(std_out,message,'COLL')
2363      write(message,'(a)') '  Calculation of proto-atomic density done in serial'
2364      call wrtout(std_out,message,'COLL')
2365    end if
2366    do itypat=1,ntypat
2367      if (prtvol>9) then
2368        write(message,'(a,I6)') '  Number of equivalent atoms:',n_equiv_atoms(itypat)
2369        call wrtout(std_out,message,'COLL')
2370      end if
2371      my_start_equiv_atoms(itypat) = 1
2372      my_end_equiv_atoms(itypat) = n_equiv_atoms(itypat)
2373    end do
2374  else
2375    if (rank==master.and.prtvol>9) then
2376      write(message,'(a,I5)') '  In atomden - number of processors:',nprocs
2377      call wrtout(std_out,message,'COLL')
2378      write(message,'(a)') '  Calculation of proto-atomic density done in parallel'
2379      call wrtout(std_out,message,'COLL')
2380    end if
2381    do itypat=1,ntypat
2382      if (rank==master.and.prtvol>9) then
2383        write(message,'(a,I6)') '  Number of equivalent atoms:',n_equiv_atoms(itypat)
2384        call wrtout(std_out,message,'COLL')
2385      end if
2386 !    Divide the atoms among the processors by shuffling indices
2387      delta = int(floor(real(n_equiv_atoms(itypat))/real(nprocs)))
2388      remainder = n_equiv_atoms(itypat)-nprocs*delta
2389      my_start_equiv_atoms(itypat) = 1+rank*delta
2390      my_end_equiv_atoms(itypat) = (rank+1)*delta
2391 !    Divide the remainder points among the processors
2392 !    by shuffling indices
2393      if ((rank+1)>remainder) then
2394        my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + remainder
2395        my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + remainder
2396      else
2397        my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + rank
2398        my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + rank + 1
2399      end if
2400      if (prtvol>9) then
2401        write(message,'(a,I3)') '          For atom type: ',itypat
2402        call wrtout(std_out,message,'PERS')
2403 !      write(message,'(a,I6)') '  I''ll take atoms from: ',my_start_equiv_atoms(itypat)
2404 !      call wrtout(std_out,message,'PERS')
2405 !      write(message,'(a,I6)') '           total for me: ',my_end_equiv_atoms(itypat)
2406 !      call wrtout(std_out,message,'PERS')
2407        write(message,'(a,I6)') '            total for me: ', &
2408 &       my_end_equiv_atoms(itypat)+1-my_start_equiv_atoms(itypat)
2409        call wrtout(std_out,message,'PERS')
2410      end if
2411    end do
2412  end if
2413 
2414 !Loop over types of atoms and equivalent atoms and
2415 !interpolate density onto grid
2416  do itypat=1,ntypat
2417 !  do iatom=my_start_equiv_atoms(itypat),my_end_equiv_atoms(itypat)
2418 
2419    cnt = 0
2420    iatom = rank+1 - nprocs
2421 !  Parallel execution of loop
2422    do
2423      cnt = cnt + 1
2424      iatom = iatom + nprocs
2425      if (iatom>n_equiv_atoms(itypat)) exit ! Exit if index is too large
2426 
2427      if (mod(cnt,100)==0.and.prtvol>0) then
2428        write(message,'(2(a,I6))') ' atoms so far',cnt,' of: ',n_equiv_atoms(itypat)/nprocs
2429        call wrtout(std_out,message,'PERS')
2430      end if
2431 
2432      r_max = atomrgrid(natomgr(itypat),itypat)
2433      r_atom = equiv_atom_pos(:,iatom,itypat)
2434 
2435 !    Set up an array with the gridpoint distances
2436      i = 1
2437      grid_distances = zero
2438      grid_index = 0
2439      do igrid=1,ngrid
2440        dp_vec_dummy(:) = r_vec_grid(:,igrid) - r_atom(:)
2441        dp_dummy = sqrt(dot_product(dp_vec_dummy,dp_vec_dummy))
2442        if (dp_dummy <= r_max) then
2443          grid_distances(i) = dp_dummy
2444          grid_index(i) = igrid
2445          i = i + 1
2446        else
2447          cycle ! cycle if point is too far away
2448        end if
2449      end do
2450      n_grid_p = i - 1
2451 
2452      if (n_grid_p==0) cycle ! Cycle if no point needs
2453 !    to be interpolated
2454 
2455 !    Sort points to be interpolated in ascending order
2456      ABI_ALLOCATE(dp_1d_dummy,(n_grid_p))
2457      ABI_ALLOCATE(new_index,(n_grid_p))
2458      ABI_ALLOCATE(i_1d_dummy,(n_grid_p))
2459      dp_1d_dummy = grid_distances(1:n_grid_p)
2460      do i=1,n_grid_p
2461        new_index(i) = i
2462      end do
2463      call sort_dp(n_grid_p,dp_1d_dummy,new_index,tol16)
2464      grid_distances(1:n_grid_p) = dp_1d_dummy
2465      i_1d_dummy = grid_index(1:n_grid_p)
2466      do i=1,n_grid_p
2467 !      write(std_out,*) i_1d_dummy(i),' -> ',i_1d_dummy(new_index(i))
2468        grid_index(i) = i_1d_dummy(new_index(i))
2469      end do
2470      ABI_DEALLOCATE(dp_1d_dummy)
2471      ABI_DEALLOCATE(new_index)
2472      ABI_DEALLOCATE(i_1d_dummy)
2473 
2474 !    Interpolate density onto all grid points
2475      ABI_ALLOCATE(ypp,(natomgr(itypat)))
2476      ABI_ALLOCATE(x_fit,(n_grid_p))
2477      ABI_ALLOCATE(y_fit,(n_grid_p))
2478      ypp = zero; y_fit = zero
2479      ybcbeg = zero; ybcend = zero
2480      x_fit = grid_distances(1:n_grid_p)
2481      call spline(atomrgrid(1:natomgr(itypat),itypat), &
2482 &     density(1:natomgr(itypat),itypat), &
2483 &     natomgr(itypat),ybcbeg,ybcend,ypp)
2484      call splint(natomgr(itypat),atomrgrid(1:natomgr(itypat),itypat), &
2485 &     density(1:natomgr(itypat),itypat),ypp,n_grid_p, &
2486 &     x_fit,y_fit)
2487 
2488 !    Save the interpolated points to grid
2489      do i=1,n_grid_p
2490        rho_temp(grid_index(i),itypat) = rho_temp(grid_index(i),itypat) + y_fit(i)
2491      end do
2492      ABI_DEALLOCATE(ypp)
2493      ABI_DEALLOCATE(x_fit)
2494      ABI_DEALLOCATE(y_fit)
2495 
2496    end do ! n equiv atoms
2497  end do ! type of atom
2498 
2499 !Collect all contributions to rho_temp if
2500 !we are running in parallel
2501  if (nprocs>1) then
2502    call xmpi_barrier(spaceComm)
2503    call xmpi_sum_master(rho_temp,master,spaceComm,ierr)
2504    call xmpi_barrier(spaceComm)
2505    if (prtvol>9) then
2506      write(message,'(a)') '  In atomden - contributions to rho_temp collected'
2507      call wrtout(std_out,message,'PERS')
2508    end if
2509  end if
2510 
2511 !Now rho_temp contains the atomic protodensity for each atom.
2512 !Check whether this is to replace or be added to the input/output array
2513 !and sum up contributions
2514  if (trim(calctype)=='replace') rho = zero
2515  do itypat=1,ntypat
2516    rho(:) = rho(:) + rho_temp(:,itypat)
2517  end do
2518 
2519 !deallocations
2520  if (allocated(rho_temp))  then
2521    ABI_DEALLOCATE(rho_temp)
2522  end if
2523  if (allocated(equiv_atom_pos))  then
2524    ABI_DEALLOCATE(equiv_atom_pos)
2525  end if
2526  if (allocated(equiv_atom_dist))  then
2527    ABI_DEALLOCATE(equiv_atom_dist)
2528  end if
2529 !if (allocated()) deallocate()
2530 
2531  return
2532 
2533  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=informations 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

PARENTS

      gstate,setup_positron

CHILDREN

      fourdp,ptabs_fourdp,wrtout,zerosym

SOURCE

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

m_mkrho/m_mkrho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mkrho

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_mkrho
28 
29  use defs_basis
30  use defs_abitypes
31  use defs_wvltypes
32  use m_abicore
33  use m_xmpi
34  use m_errors
35 
36  use m_time,         only : timab
37  use m_fftcore,      only : sphereboundary
38  use m_fft,          only : fftpac, zerosym, fourwf, fourdp
39  use m_hamiltonian,  only : gs_hamiltonian_type
40  use m_bandfft_kpt,  only : bandfft_kpt_set_ikpt
41  use m_paw_dmft,     only : paw_dmft_type
42  use m_spacepar,     only : symrhg
43  use defs_datatypes, only : pseudopotential_type
44  use m_atomdata,     only : atom_length
45  use m_mpinfo,       only : ptabs_fourdp, proc_distrb_cycle
46  use m_pawtab,       only : pawtab_type
47  use m_io_tools,     only : open_file
48  use m_splines,      only : spline,splint
49  use m_sort,         only : sort_dp
50  use m_prep_kgb,     only : prep_fourwf
51  use m_wvl_rho,      only : wvl_mkrho
52 
53  implicit none
54 
55  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=informations 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 informations for wavelets
  wvl_wfs <type(wvl_projector_type)>=wavefunctions informations 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)

PARENTS

      afterscfloop,energy,gstate,respfn,scfcv,vtorho

CHILDREN

      bandfft_kpt_set_ikpt,fftpac,fourwf,prep_fourwf,prtrhomxmn
      sphereboundary,symrhg,timab,wrtout,wvl_mkrho,xmpi_sum

SOURCE

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

SIDE EFFECTS

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)

PARENTS

      afterscfloop,bethe_salpeter,clnup1,mkrho,screening,sigma,vtorho

CHILDREN

      wrtout,xmpi_sum

SOURCE

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

m_mkrho/read_atomden [ Functions ]

[ Top ] [ m_mkrho ] [ Functions ]

NAME

 read_atomden

FUNCTION

COPYRIGHT

 Copyright (C) 2005-2018 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

PARENTS

      outscfcv

CHILDREN

      atomden

SOURCE

1934 subroutine read_atomden(MPI_enreg,natom,nfft,ngfft,nspden,ntypat, &
1935 &                       rhor_atm,typat,rprimd,xred,prtvol,file_prefix)
1936 
1937 
1938 !This section has been created automatically by the script Abilint (TD).
1939 !Do not modify the following lines by hand.
1940 #undef ABI_FUNC
1941 #define ABI_FUNC 'read_atomden'
1942 !End of the abilint section
1943 
1944  implicit none
1945 
1946 !Arguments ------------------------------------
1947 !scalars
1948  integer,intent(in) :: natom,nfft,nspden,ntypat,prtvol
1949 !arrays
1950  type(MPI_type),intent(in) :: MPI_enreg
1951  integer,intent(in) :: ngfft(18),typat(natom)
1952  real(dp), intent(in) :: rprimd(3,3),xred(3,natom)
1953  real(dp),intent(inout) :: rhor_atm(nfft,nspden)
1954  character(len=7), intent(in) :: file_prefix
1955 
1956 !Local variables-------------------------------
1957 !scalars
1958  character(len=500) :: message
1959  character(len=120) :: filename
1960  character(len=7) :: calctype='replace'
1961  integer :: igrid,i,i1,i2,i3,io_err,itypat,unt
1962  integer :: natomgrmax,nlines,ngrid,n1,n2,n3
1963  real(dp) :: difx,dify,difz,ucvol!,norm
1964 !arrays
1965  integer :: natomgr(ntypat)
1966  real(dp) :: a(3),b(3),c(3)
1967  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1968  real(dp),allocatable :: atomrgrid(:,:),r_vec_grid(:,:),density(:,:)
1969  real(dp),allocatable :: rho(:)
1970 
1971 ! ************************************************************************
1972 
1973 !Initialise various variables
1974  ngrid = nfft
1975  a(:) = rprimd(:,1)
1976  b(:) = rprimd(:,2)
1977  c(:) = rprimd(:,3)
1978  ABI_ALLOCATE(rho,(ngrid))
1979  if (nspden/=1) then
1980    MSG_ERROR('read_atomden: Only nspden=1 allowed.')
1981  end if
1982  rho = rhor_atm(:,1)
1983  gmet=zero;gprimd=zero;rmet=zero;ucvol=zero
1984 
1985 
1986 !Calculate the r vector (reduced coord.) of the fine gridpoints
1987  ABI_ALLOCATE(r_vec_grid,(3,ngrid))
1988  igrid = 0
1989  n1 = ngfft(1)
1990  n2 = ngfft(2)
1991  n3 = ngfft(3)
1992  do i3=0,n3-1
1993    difz=dble(i3)/dble(n3)
1994    do i2=0,n2-1
1995      dify=dble(i2)/dble(n2)
1996      do i1=0,n1-1
1997        difx=dble(i1)/dble(n1)
1998        igrid = igrid + 1
1999        r_vec_grid(1,igrid)=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
2000        r_vec_grid(2,igrid)=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
2001        r_vec_grid(3,igrid)=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
2002      end do
2003    end do
2004  end do
2005  if (igrid/=ngrid) then
2006    MSG_ERROR('read_atomden: igrid not equal to ngrid')
2007  end if
2008 
2009 !Read in atomic density data for each atom type
2010 !first check how many datapoints are in each file
2011  do itypat=1,ntypat
2012    filename='';io_err=0;
2013    if (itypat>0)  write(filename,'(a,a,i1,a)') trim(file_prefix), &
2014 &   '_density_atom_type',itypat,'.dat'
2015    if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), &
2016 &   '_density_atom_type',itypat,'.dat'
2017    if (open_file(filename, message, newunit=unt, status='old',action='read') /= 0) then
2018      write(std_out,*) 'ERROR in read_atomden: Could not open file: ',filename
2019      write(std_out,*) ' Current implementation requires this file to be present'
2020      write(std_out,*) ' for each type of atom.'
2021      write(std_out,*)trim(message)
2022      MSG_ERROR("Cannot continue")
2023    end if
2024 !  Check number of lines in file
2025    nlines = 1;io_err=0;
2026    do
2027      read(unt,*,iostat=io_err)
2028      if (io_err<0) exit
2029      nlines = nlines + 1
2030    end do
2031    close(unt)
2032    natomgr(itypat) = nlines - 2
2033  end do ! Atom type
2034 !Allocate arrays and read in data
2035  natomgrmax = maxval(natomgr)
2036  ABI_ALLOCATE(atomrgrid,(natomgrmax,ntypat))
2037  ABI_ALLOCATE(density,(natomgrmax,ntypat))
2038  atomrgrid = zero ; density = zero
2039  do itypat=1,ntypat
2040    filename='';io_err=0;
2041    if (itypat>0)  write(filename,'(a,a,i1,a)') trim(file_prefix), &
2042 &   '_density_atom_type',itypat,'.dat'
2043    if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), &
2044 &   '_density_atom_type',itypat,'.dat'
2045    if (open_file(filename,message,newunit=unt,status='old',action='read') /= 0) then
2046      MSG_ERROR(message)
2047    end if
2048    read(unt,*) ! Skip comment line
2049    do i=1,natomgr(itypat)
2050      read(unt,*) atomrgrid(i,itypat),density(i,itypat)
2051    end do
2052    close(unt)
2053    if (atomrgrid(1,itypat)/=zero) then
2054      write(std_out,*) 'ERROR in read_atomden, in file: ',filename
2055      write(std_out,*) ' First gridpoint has to be the origin.'
2056      MSG_ERROR("Cannot continue")
2057    end if
2058  end do ! Atom type
2059 
2060 !write(std_out,*) '*** --- In read_atomden before call--- ***'
2061 !write(std_out,*) '  calctype:',calctype,' natom:',natom
2062 !write(std_out,*) '    ntypat:',ntypat,' typat:',typat
2063 !write(std_out,*) '     ngrid:',ngrid
2064 !write(std_out,*) '         a:',a
2065 !write(std_out,*) '         b:',b
2066 !write(std_out,*) '         c:',c
2067 !write(std_out,*) '      xred:',xred
2068 !write(std_out,*) '   natomgr:',natomgr
2069 !write(std_out,*) 'natomgrmax:',natomgrmax
2070 !write(std_out,*) ' atomrgrid:',atomrgrid
2071 !write(std_out,*) '   density:',density
2072 !write(std_out,*) 'r_vec_grid:'
2073 !write(std_out,*) r_vec_grid
2074 
2075 !Call atomden
2076  call atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,xred, &
2077 & natomgr,natomgrmax,atomrgrid,density,prtvol,calctype)
2078 
2079 !if (prtvol>9) then ! calculate norm
2080 !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2081 !norm = SUM(rho(:))*ucvol/dble(n1*n2*n3)
2082 !write(message,'(a,F8.4)') '  In read_atomden - NORM OF DENSITY: ',norm
2083 !call wrtout(std_out,message,'COLL')
2084 !end if
2085 
2086  rhor_atm(:,1) = rho
2087 
2088  if (allocated(atomrgrid))  then
2089    ABI_DEALLOCATE(atomrgrid)
2090  end if
2091  if (allocated(density))  then
2092    ABI_DEALLOCATE(density)
2093  end if
2094  if (allocated(r_vec_grid))  then
2095    ABI_DEALLOCATE(r_vec_grid)
2096  end if
2097  if (allocated(rho))  then
2098    ABI_DEALLOCATE(rho)
2099  end if
2100 
2101  return
2102 
2103 end subroutine read_atomden