LCOV - code coverage report
Current view: top level - src/57_iovars - m_memeval.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 90.3 % 1042 941
Test Date: 2026-09-20 18:56:22 Functions: 100.0 % 6 6

            Line data    Source code
       1              : !!****m* ABINIT/m_memeval
       2              : !! NAME
       3              : !! m_memeval
       4              : !!
       5              : !! FUNCTION
       6              : !!  Functions to estimate memory requirements from the calculation parameters.
       7              : !!
       8              : !! COPYRIGHT
       9              : !!  Copyright (C) 2008-2026 ABINIT group (XG, DC, DW)
      10              : !!  This file is distributed under the terms of the
      11              : !!  GNU General Public License, see ~abinit/COPYING
      12              : !!  or http://www.gnu.org/copyleft/gpl.txt .
      13              : !!
      14              : !! SOURCE
      15              : 
      16              : #if defined HAVE_CONFIG_H
      17              : #include "config.h"
      18              : #endif
      19              : 
      20              : #include "abi_common.h"
      21              : 
      22              : MODULE m_memeval
      23              : 
      24              :  use defs_basis
      25              :  use m_abicore
      26              :  use m_xmpi
      27              :  use m_errors
      28              :  use m_dtset
      29              : 
      30              :  use defs_datatypes, only : pspheader_type
      31              :  use defs_abitypes,   only : MPI_type
      32              :  use m_geometry,      only : mkradim, mkrdim, xred2xcart, metric
      33              :  use m_matrix,        only : mati3inv
      34              :  use m_symtk,         only : littlegroup_q
      35              :  use m_spgdata,       only : prtspgroup
      36              :  use m_fftcore,       only : getng
      37              :  use m_kg,            only : getmpw
      38              :  use m_libpaw_tools,  only : libpaw_write_comm_set
      39              : 
      40              :  implicit none
      41              : 
      42              :  private
      43              : !!***
      44              : 
      45              :  public :: memory_eval   ! Main entry point
      46              :  public :: getdim_nloc   ! Determine the dimensions of arrays with non-local projectors: ekb, ffspl, indlmn
      47              :  public :: setmqgrid     ! Sets the number of points needed to represent the pseudopotentials in q-space
      48              : 
      49              : contains
      50              : 
      51              : !!****f* m_memeval/memory_eval
      52              : !! NAME
      53              : !! memory_eval
      54              : !!
      55              : !! FUNCTION
      56              : !! Big loop on the datasets:
      57              : !! - for each of the datasets, write one line about the crystallographic data
      58              : !! - compute the memory needs for this data set.
      59              : !!
      60              : !! INPUTS
      61              : !!  dtsets(0:ndtset_alloc)=<type datafiles_type>contains all input variables
      62              : !!  iout=unit number of output file
      63              : !!  mpi_enregs=information about MPI parallelization
      64              : !!  ndtset= number of datasets to be read; if 0, no multi-dataset mode
      65              : !!  ndtset_alloc=number of datasets, corrected for allocation of at least
      66              : !!      one data set.
      67              : !!  npsp=number of pseudopotentials
      68              : !!  pspheads(npsp)=<type pspheader_type>all the important information from the
      69              : !!   pseudopotential file header, as well as the psp file name
      70              : !!
      71              : !! OUTPUT
      72              : !!   printing only
      73              : !!
      74              : !! SOURCE
      75              : 
      76         1427 : subroutine memory_eval(dtsets,iout,mpi_enregs,ndtset,ndtset_alloc,npsp,pspheads)
      77              : 
      78              : !Arguments ------------------------------------
      79              : !scalars
      80              :  integer,intent(in) :: iout,ndtset,ndtset_alloc,npsp
      81              :  type(MPI_type),intent(inout) :: mpi_enregs(0:ndtset_alloc)
      82              : !arrays
      83              :  type(dataset_type),intent(inout) :: dtsets(0:ndtset_alloc)
      84              :  type(pspheader_type),intent(in) :: pspheads(npsp)
      85              : 
      86              : !Local variables -------------------------------
      87              : !scalars
      88              :  integer :: cplex,exchn2n3d,extrapwf,getcell,idtset,ii,intxc,densfor_pred,iprcel
      89              :  integer :: iscf,isym,jdtset,lmnmax,mem_test
      90              :  integer :: lmnmax_eff,lmnmaxso,lnmax,lnmax_eff,lnmaxso,mband
      91              :  integer :: me_fft,mffmem,mgfftdiel,mgfftf,mkmem,mpsang,mpspso
      92              :  integer :: mpssoang,mpw,mqgrid,mqgriddg,mqgrid_ff,mqgrid_vl,n1xccc,natom
      93              :  integer :: nfftdiel,nfftf,nkpt,nproc_fft,nptsgvec,npulayit,npwdiel,nspden,nspinor
      94              :  integer :: nsppol,nsym,ntypat,occopt,optddk,optforces,optphon,optstress
      95              :  integer :: optstrs,paral_fft,pawcpxocc,pawmixdg,pawnhatxc,pawspnorb,pawstgylm,prtvol,ptgroupma,response
      96              :  integer :: spgroup,timrev,usepaw,useylm,gpu_option,xclevel
      97              :  real(dp) :: diecut,dilatmx,ecut,ecut_eff,ecutdg_eff,ecutsus,ucvol
      98              : !arrays
      99              :  integer :: bravais(11),mkmems(3),ngfftdiel(18)
     100              :  integer :: ngfftf(18),nloalg(3)
     101         1427 :  integer,allocatable :: nband(:),symq(:,:,:),symrec(:,:,:),symrel(:,:,:)
     102              :  real(dp),parameter :: k0(3)=(/zero,zero,zero/)
     103              :  real(dp) :: genafm(3),gmet(3,3),gprimd(3,3),kpt_diel(3),qphon(3),rmet(3,3),rprimd(3,3)
     104              : 
     105              : !*************************************************************************
     106              : 
     107         8385 :  do idtset=1,ndtset_alloc
     108         6958 :    if(mpi_enregs(idtset)%me<0) cycle
     109         6958 :    call abi_io_redirect(new_io_comm=mpi_enregs(idtset)%comm_world)
     110         6958 :    call libpaw_write_comm_set(mpi_enregs(idtset)%comm_world)
     111              : 
     112              : !  Initialisations
     113        83496 :    bravais(:)=dtsets(idtset)%bravais(:)
     114         6958 :    exchn2n3d=dtsets(idtset)%exchn2n3d
     115         6958 :    extrapwf=dtsets(idtset)%extrapwf
     116        27832 :    genafm(:) =dtsets(idtset)%genafm(:)
     117         6958 :    getcell=dtsets(idtset)%getcell
     118         6958 :    intxc=dtsets(idtset)%intxc
     119         6958 :    densfor_pred=dtsets(idtset)%densfor_pred
     120         6958 :    iprcel=dtsets(idtset)%iprcel
     121         6958 :    iscf=dtsets(idtset)%iscf
     122         6958 :    jdtset=dtsets(idtset)%jdtset ; if(ndtset==0)jdtset=0
     123         6958 :    me_fft=mpi_enregs(idtset)%me_fft
     124         6958 :    mffmem=dtsets(idtset)%mffmem
     125         6958 :    mpw=dtsets(idtset)%mpw
     126         6958 :    mqgrid=dtsets(idtset)%mqgrid
     127         6958 :    mqgriddg=dtsets(idtset)%mqgriddg
     128         6958 :    natom=dtsets(idtset)%natom
     129         6958 :    nkpt  =dtsets(idtset)%nkpt
     130        27832 :    nloalg(:)=dtsets(idtset)%nloalg(:)
     131         6958 :    nproc_fft=mpi_enregs(idtset)%nproc_fft
     132         6958 :    npulayit=dtsets(idtset)%npulayit
     133         6958 :    nspden=dtsets(idtset)%nspden
     134         6958 :    nspinor=dtsets(idtset)%nspinor
     135         6958 :    nsppol=dtsets(idtset)%nsppol
     136         6958 :    nsym     =dtsets(idtset)%nsym
     137         6958 :    ntypat=dtsets(idtset)%ntypat
     138         6958 :    occopt=dtsets(idtset)%occopt
     139         6958 :    optforces=dtsets(idtset)%optforces
     140         6958 :    paral_fft=mpi_enregs(idtset)%paral_kgb
     141         6958 :    pawcpxocc=dtsets(idtset)%pawcpxocc
     142         6958 :    pawmixdg=dtsets(idtset)%pawmixdg
     143         6958 :    pawnhatxc=dtsets(idtset)%pawnhatxc
     144         6958 :    pawspnorb=dtsets(idtset)%pawspnorb
     145         6958 :    pawstgylm=dtsets(idtset)%pawstgylm
     146         6958 :    prtvol=dtsets(idtset)%prtvol
     147         6958 :    ptgroupma =dtsets(idtset)%ptgroupma
     148        27832 :    qphon(:)=dtsets(idtset)%qptn(:)
     149         6958 :    spgroup   =dtsets(idtset)%spgroup
     150         6958 :    usepaw=dtsets(idtset)%usepaw
     151         6958 :    useylm=dtsets(idtset)%useylm
     152         6958 :    gpu_option=dtsets(idtset)%gpu_option
     153         6958 :    xclevel=dtsets(idtset)%xclevel
     154              : 
     155        20874 :    ABI_MALLOC(symrel,(3,3,nsym))
     156      2199265 :    symrel(:,:,1:nsym)=dtsets(idtset)%symrel(:,:,1:nsym)
     157              : 
     158              : !  Space group output
     159         6958 :    call prtspgroup(bravais,genafm,std_out,jdtset,ptgroupma,spgroup)
     160         6958 :    call prtspgroup(bravais,genafm,iout,jdtset,ptgroupma,spgroup)
     161              : 
     162         6958 :    if (dtsets(idtset)%toldff>tol16.and.optforces==0) optforces=1
     163         6958 :    if (dtsets(idtset)%tolrff>tol16.and.optforces==0) optforces=1
     164         6958 :    if (dtsets(idtset)%ionmov>tol16.and.optforces==0) optforces=1
     165         6958 :    if (dtsets(idtset)%imgmov>tol16.and.optforces==0) optforces=1
     166         6958 :    optstress=dtsets(idtset)%optstress
     167         6958 :    optddk=0;optphon=0;optstrs=0
     168         6958 :    if (dtsets(idtset)%rfddk>0.or.dtsets(idtset)%rf2_dkdk>0.or.dtsets(idtset)%rf2_dkde>0) optddk=1
     169              :    if (dtsets(idtset)%rfelfd>0.or.dtsets(idtset)%d3e_pert1_elfd>0.or.&
     170         6958 : &   dtsets(idtset)%d3e_pert2_elfd>0.or.dtsets(idtset)%d3e_pert3_elfd>0) optddk=1
     171              :    if (dtsets(idtset)%rfphon>0.or.dtsets(idtset)%d3e_pert1_phon>0.or.&
     172         6958 : &   dtsets(idtset)%d3e_pert2_phon>0.or.dtsets(idtset)%d3e_pert3_phon>0) optphon=1
     173         6958 :    if (dtsets(idtset)%rfstrs>0) optstrs=1
     174              : 
     175        20874 :    ABI_MALLOC(nband,(nkpt*nsppol))
     176       124357 :    nband(1:nkpt*nsppol)=dtsets(idtset)%nband(1:nkpt*nsppol)
     177       124357 :    mband=maxval(nband(1:nkpt*nsppol))
     178         6958 :    dtsets(idtset)%mband=mband
     179              : 
     180              : !  mpsang=max(maxval(pspheads(1:npsp)%lmax)+1,1) ! Likely problems with the HP compiler
     181              : !  n1xccc=maxval(pspheads(1:npsp)%xccc)
     182         6958 :    mpsang=1
     183         6958 :    n1xccc=pspheads(1)%xccc
     184        15563 :    do ii=1,npsp
     185         8605 :      mpsang=max(pspheads(ii)%lmax+1,mpsang)
     186        15563 :      n1xccc=max(pspheads(ii)%xccc,n1xccc)
     187              :    end do
     188              : 
     189              : !  Determine the maximum number of projectors, for the set of pseudo atom
     190              :    call getdim_nloc(lmnmax,lmnmaxso,lnmax,lnmaxso,dtsets(idtset)%mixalch_orig,&
     191         6958 : &   dtsets(idtset)%nimage,npsp,dtsets(idtset)%npspalch,ntypat,dtsets(idtset)%ntypalch,pspheads)
     192              : 
     193              : !  Treatment of the effect of using a spin-orbit part
     194              : !  Warning: mpspso is different for each dataset; not relevant for PAW
     195         6958 :    mpspso=1
     196         6958 :    if (dtsets(idtset)%usepaw==0) then
     197        12607 :      do ii=1,npsp
     198        12607 :        if(nspinor/=1)then
     199          273 :          if(pspheads(ii)%pspso/=0)then
     200          204 :            if(dtsets(idtset)%so_psp(ii)/=0)then
     201         6829 :              mpspso=2
     202              :            end if
     203              :          end if
     204              :        end if
     205              :      end do
     206              :    end if
     207              : !  In case of no spin-orbit
     208         5778 :    if(mpspso==1)then
     209         6861 :      mpssoang=mpsang ; lmnmax_eff =lmnmax; lnmax_eff =lnmax
     210              :    else ! spin-orbit will be used
     211           97 :      mpssoang=2*mpsang-1 ; lmnmax_eff =lmnmaxso ; lnmax_eff =lnmaxso
     212              :    end if
     213              : !  lmnmax is not used if the Ylm are not used
     214         6958 :    if (useylm==0) lmnmax_eff =lnmax_eff
     215              : 
     216         6958 :    ecut     =dtsets(idtset)%ecut
     217         6958 :    dilatmx  =dtsets(idtset)%dilatmx
     218         6958 :    ecut_eff=ecut*dilatmx**2
     219         6958 :    ecutdg_eff=dtsets(idtset)%pawecutdg*dtsets(idtset)%dilatmx**2
     220              : 
     221              : !  Compute mgfft,mpw,nfft for this data set
     222         6958 :    call mkrdim(dtsets(idtset)%acell_orig(1:3,1),dtsets(idtset)%rprim_orig(1:3,1:3,1),rprimd)
     223         6958 :    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
     224              : 
     225         6958 :    if (usepaw==0) then
     226       109782 :      mgfftf=dtsets(idtset)%mgfft;nfftf=dtsets(idtset)%nfft;ngfftf(:)=dtsets(idtset)%ngfft(:)
     227              :    else
     228        22420 :      mgfftf=dtsets(idtset)%mgfftdg;nfftf=dtsets(idtset)%nfftdg;ngfftf(:)=dtsets(idtset)%ngfftdg(:)
     229              :    end if
     230         6958 :    response=0
     231              :    if(dtsets(idtset)%rfddk/=0  .or. dtsets(idtset)%rf2_dkdk/=0 .or. dtsets(idtset)%rf2_dkde/=0 .or. &
     232              : &   dtsets(idtset)%rfphon/=0 .or. dtsets(idtset)%rfelfd/=0 .or. &
     233         6958 : &   dtsets(idtset)%rfstrs/=0 .or. dtsets(idtset)%rfmagn/=0    ) response=1
     234              : 
     235              : !  Compute mgfftdiel,npwdiel,nfftdiel for this data set
     236         6958 :    if((modulo(iprcel,100)>=20 .and.modulo(iprcel,100) < 71).or. iscf==-1)then
     237              : !    Get diecut, and the fft grid to be used for the susceptibility computation
     238           64 :      diecut=abs(dtsets(idtset)%diecut)
     239           64 :      if( dtsets(idtset)%diecut < zero )then
     240              :        ecutsus=ecut
     241              :      else
     242           64 :        ecutsus= ( sqrt(ecut) *0.5_dp + sqrt(diecut) *0.25_dp )**2
     243              :      end if
     244              : !    Beware, for the dielectric matrix fftalg=ngfftdiel(7) is default here
     245          960 :      ngfftdiel(1:3)=0 ; ngfftdiel(7)=101 ; ngfftdiel(8:18)=dtsets(idtset)%ngfft(8:18)
     246           64 :      if(iscf==-1)ngfftdiel(7)=102
     247           64 :      ecut_eff=ecutsus*dilatmx**2
     248              :      call getng(dtsets(idtset)%boxcutmin,dtsets(idtset)%chksymtnons,&
     249              : &     ecut_eff,gmet,k0,me_fft,mgfftdiel,nfftdiel,&
     250              : &     ngfftdiel,nproc_fft,nsym,paral_fft,symrel,dtsets(idtset)%tnons,&
     251           64 : &     gpu_option=dtsets(idtset)%gpu_option)
     252              : !    Compute the size of the dielectric matrix : npwdiel
     253           64 :      kpt_diel(1:3)=(/ 0.0_dp, 0.0_dp, 0.0_dp /)
     254           64 :      ecut_eff=diecut*dilatmx**2
     255           64 :      call getmpw(ecut_eff,exchn2n3d,gmet,(/1/),kpt_diel,mpi_enregs(idtset),npwdiel,1)
     256              :    else
     257        62046 :      npwdiel=1 ; mgfftdiel=1 ; nfftdiel=1 ; ngfftdiel(1:8)=1
     258              :    end if
     259              : 
     260              : !  Special treatment for the value of mqgrid to be fed in memory.F90
     261              : 
     262         6958 :    nptsgvec         = 200 ! At present, this has to be chosen once and for all ...
     263         6958 :    if ( dtsets(idtset)%usewvl == 0) then
     264         6958 :      call setmqgrid(mqgrid,mqgriddg,ecut_eff,ecutdg_eff,gprimd,nptsgvec,usepaw)
     265              :    else
     266            0 :      call setmqgrid(mqgrid,mqgriddg,one,one,gprimd,nptsgvec,usepaw)
     267              :    end if
     268         6958 :    mqgrid_ff=mqgrid
     269         6958 :    if (usepaw==0) mqgrid_vl=mqgrid
     270         6958 :    if (usepaw==1) mqgrid_vl=mqgriddg
     271              : 
     272              : !  Compute the memory needs for this data set.
     273         6958 :    if(response==0)then
     274              : 
     275         6240 :      if (dtsets(idtset)%usewvl == 0) then
     276         6240 :        mkmem=dtsets(idtset)%mkmem
     277        67273 :        mband=maxval(dtsets(idtset)%nband(1:nkpt*nsppol))
     278              : 
     279              :        ! Don't perform memory tests if MBPT.
     280         6240 :        mem_test = dtsets(idtset)%mem_test
     281        35538 :        if (any(dtsets(idtset)%optdriver == [RUNL_SIGMA, RUNL_SCREENING, RUNL_BSE, RUNL_EPH, RUNL_GWR])) mem_test = 0
     282              : 
     283              :        call memory(n1xccc,extrapwf,getcell,idtset,dtsets(idtset)%icoulomb,&
     284              : &       intxc,dtsets(idtset)%ionmov,iout,densfor_pred,&
     285              : &       iprcel,iscf,jdtset,lmnmax_eff,lnmax_eff,mband,mffmem,dtsets(idtset)%mgfft,mgfftdiel,mgfftf,mkmem,&
     286              : &       mpi_enregs(idtset),mpsang,mpssoang,mpw,mqgrid_ff,mqgrid_vl,natom,nband,dtsets(idtset)%nfft,nfftdiel,nfftf,&
     287              : &       dtsets(idtset)%ngfft,ngfftdiel,ngfftf,dtsets(idtset)%nimage,nkpt,nloalg,npsp,npulayit,npwdiel,nspden,nspinor,&
     288              : &       nsppol,nsym,ntypat,occopt,optforces,mem_test,optstress,pawcpxocc,pawmixdg,&
     289              : &       pawnhatxc,pawspnorb,pawstgylm,prtvol,pspheads,qphon,dtsets(idtset)%tfkinfunc,&
     290         6240 : &       dtsets(idtset)%typat,ucvol,usepaw,useylm,gpu_option,xclevel)
     291            0 :      else if( dtsets(idtset)%usepaw==0) then
     292            0 :        if (mpi_enregs(idtset)%me == 0) then
     293            0 :          call wvl_memory(dtsets(idtset), idtset, mpi_enregs(idtset), npsp, 1, pspheads)
     294              :        end if
     295              :      end if
     296              : 
     297              :    else
     298              : !    Compute the value of cplex, for which one needs symrec
     299         2154 :      ABI_MALLOC(symq,(4,2,nsym))
     300         1436 :      ABI_MALLOC(symrec,(3,3,nsym))
     301        17871 :      do isym=1,nsym
     302        17871 :        call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
     303              :      end do
     304          718 :      call littlegroup_q(nsym,qphon,symq,symrec,dtsets(idtset)%symafm,timrev)
     305          718 :      cplex=2-timrev
     306          718 :      ABI_FREE(symq)
     307          718 :      ABI_FREE(symrec)
     308          718 :      mkmems(1)=dtsets(idtset)%mkmem
     309          718 :      mkmems(2)=dtsets(idtset)%mkqmem
     310          718 :      mkmems(3)=dtsets(idtset)%mk1mem
     311              : 
     312          718 :      mem_test = dtsets(idtset)%mem_test
     313              : 
     314              :      call memorf(cplex,n1xccc,getcell,idtset,intxc,iout,iprcel,&
     315              : &     iscf,jdtset,lmnmax_eff,lnmax_eff,mband,mffmem,dtsets(idtset)%mgfft,&
     316              : &     mkmems,mpi_enregs(idtset),mpsang,mpssoang,mpw,mqgrid_ff,natom,nband,dtsets(idtset)%nfft,&
     317              : &     dtsets(idtset)%ngfft,nkpt,nloalg,nspden,nspinor,nsppol,nsym,&
     318         1436 : &     ntypat,occopt,optddk,optphon,mem_test,optstrs,prtvol,useylm,gpu_option,xclevel)
     319              :    end if
     320              : 
     321              : !  Deallocate temporary arrays (when they will really be temporary !)
     322         6958 :    ABI_FREE(nband)
     323         8385 :    ABI_FREE(symrel)
     324              : 
     325              :  end do ! idtset
     326              : 
     327         1427 : end subroutine memory_eval
     328              : !!***
     329              : 
     330              : !!****f* m_memeval/memory
     331              : !! NAME
     332              : !! memory
     333              : !!
     334              : !! FUNCTION
     335              : !! Estimation of the memory needed for a ground-state job.
     336              : !! According to the value of the option variable,
     337              : !! might also try to allocate this amount of memory, and if it fails,
     338              : !! might estimate the available memory.
     339              : !!
     340              : !! INPUTS
     341              : !!  extrapwf=flag controlling the extrapolation of wave functions during MD or relaxation
     342              : !!  getcell=if non-zero, the values of acell and rprim are taken from
     343              : !!   the output of another dataset
     344              : !!  idtset=number of the current dataset
     345              : !!  icoulomb=0 for periodic Fourier calculation of Hartree potential; 1 for isolated system using Poisson solver.
     346              : !!  intxc=control xc quadrature
     347              : !!  ionmov=control force calculations
     348              : !!  iout=unit number for output of formatted data.
     349              : !!  densfor_pred=govern the choice of density prediction and/or forces correction
     350              : !!  iprcel=govern the choice of preconditioner for the SCF cycle
     351              : !!  iscf=governs the choice of SCF algorithm, or non-SCF calculation.
     352              : !!  jdtset=index of the current dataset
     353              : !!  lmnmax=max. number of (l,m,n) components over all type of psps
     354              : !!  lnmax =max. number of (l,n)   components over all type of psps
     355              : !!  mband =maximum number of bands
     356              : !!  mffmem =governs the number of FFT arrays which are fit in core memory
     357              : !!  mgfftf =maximum single fft dimension (fine grid, if PAW)
     358              : !!  mgfft  =maximum single fft dimension (coarse grid, if PAW)
     359              : !!  mgfftdiel =maximum single fft dimension for susceptibility and dielectric
     360              : !!   matrices.
     361              : !!  mkmem =maximum number of k points which can fit in core memory
     362              : !!  mpi_enreg=information about MPI parallelization
     363              : !!  mpssang is 1+maximum angular momentum for nonlocal pseudopotential
     364              : !!  mpssoang is 1+maximum (spin*angular momentum) for nonlocal pseudopotential
     365              : !!  mpw   =maximum number of planewaves in basis sphere (large number)
     366              : !!  mqgrid_ff=dimension of q (or G) grid for nl form factors (array ffspl)
     367              : !!  mqgrid_vl=dimension of q (or G) grid for Vloc (array vlspl)
     368              : !!  natom =number of atoms in unit cell
     369              : !!  nband(nkpt*nsppol)=number of bands at each k point, for each polarization
     370              : !!  nfftf =number of fft grid points for density        (fine grid, if PAW)
     371              : !!  nfft  =number of fft grid points for wavefunctions  (coarse grid, if PAW)
     372              : !!  nfftdiel  =maximum number of fft grid points for susceptibility
     373              : !!    and dielectric matrices
     374              : !!  ngfftf(18)=contain all needed information about 3D FFT (fine grid, if PAW)
     375              : !!  ngfft(18) =contain all needed information about 3D FFT (coarse grid, if PAW)
     376              : !!  ngfftdiel(18)=contain all needed information about 3D FFT, dielectric case,
     377              : !!                 see ~abinit/doc/variables/vargs.htm#ngfft
     378              : !!    for susceptibility and dielectric matrices
     379              : !!  nimage=number of images (replicas) of the cell
     380              : !!  nkpt  =number of k points
     381              : !!  npsp=number of different pseudopotentials
     382              : !!  npwdiel=number of plane wave for susceptibility and dielectric matrix
     383              : !!  npulayit=number of iterations used in Pulay SCF mixing
     384              : !!  nloalg(3)=governs the choice of the algorithm for non-local operator.
     385              : !!  nspden=number of spin-density components
     386              : !!  nspinor=number of spinorial components of the wavefunctions
     387              : !!  nsppol=number of channels for spin-polarization (1 or 2)
     388              : !!  nsym  =number of symmetry elements in space group
     389              : !!  ntypat =number of types of atoms
     390              : !!  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
     391              : !!  occopt=option for occupation numbers. If 3<=occopt<=8, varying occupation
     392              : !!  optforces=1 if forces are computed during run
     393              : !!  option : if 0 , no test of available memory
     394              : !!           if 1 , the routine tries to allocate the estimated memory, for testing
     395              : !!                    purposes, and if a failure occurs, the routine stops.
     396              : !!           if 2 , like 1, but before stopping, the routine will provide
     397              : !!                    an estimation of the available memory.
     398              : !!  optstress=1 if stresses are computed during run
     399              : !!  pawcpxocc=2 if PAW occupancies (rhoij) are complex
     400              : !!  pawmixdg=1 if mixing (in PAW) is done on the fine grid
     401              : !!  pawnhatxc=1 if nhat PAW density has to be analytically included in XC
     402              : !!  pawspnorb=1 when spin-orbit is activated within PAW
     403              : !!  pawstgylm=1 if g_l(r).Y_lm(r) factors are stored in memory (PAW)
     404              : !!  prtvol=control print volume
     405              : !!  pspheads(npsp)=<type pspheader_type>all the important information from the header
     406              : !!  tfkinfun=flag controling the use of Thomas-Fermi algorithme (without WF)
     407              : !!  typat(natom)=type of each atom
     408              : !!  ucvol= unit cell volume
     409              : !!  usepaw= 0 for non paw calculation; =1 for paw calculation
     410              : !!  useylm=governs the way the nonlocal operator is to be applied:
     411              : !!         1=using Ylm, 0=using Legendre polynomials
     412              : !!  gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
     413              : !!  xclevel=XC functional level
     414              : !!
     415              : !! OUTPUT
     416              : !!  (only writing)
     417              : !!
     418              : !! NOTES
     419              : !! for the estimation, it is only taken into account those
     420              : !! arrays that have some probability of being larger than 1000*8 bytes :
     421              : !! - All the arrays that have large numbers as one of their dimensions
     422              : !! (mqgrid, mpw, nfft, ngfft(4)*ngfft(5)*ngfft(6),
     423              : !!                     ngfftdiel(4)*ngfftdiel(5)*ngfftdiel(6), n1xccc
     424              : !!                                      or a constant larger than 1000)
     425              : !! - All the arrays that have a product of two moderately large numbers
     426              : !! (potential size above 30  : mband, mgfft, mkmem, natom, nkpt, nsym,
     427              : !!  or a constant larger than 30)
     428              : !! After this estimation, an amount of (176 + 55 + 6*natom) Kbytes is added
     429              : !! to take into account the static arrays declared
     430              : !! in rhotoxc and daughter routines (at maximum 22*1000 dp numbers),
     431              : !! as well as other arrays like
     432              : !! character(len=500) :: message (present in about 100 routines), or the different
     433              : !! arrays allocated in move.f, brdmin.f, gstate.f (xf array) or pspini.f
     434              : !! In the case 3<=occopt<=8 this amount is increased by 760 Kbytes
     435              : !! to take into account the arrays smdfun, occfun, entfun, workfun and xgrid,
     436              : !! declared in getnel.
     437              : !!
     438              : !! The current version takes into account
     439              : !! 1) and 2) the "main chain" in its two slightly different versions :
     440              : !! driver - gstate - (move or brdmin) - scfcv - vtorho - vtowfk -
     441              : ! !!     cgwf - getghc - fourwf or (nonlop+opernl)
     442              : !! 3) the xc chain :
     443              : !! driver - gstate - (move or brdmin) - scfcv - (vresfo) - rhotoxc - xcden
     444              : !! 4) the mkrho chain :
     445              : !! driver - gstate - (move or brdmin) - scfcv - vtorho - mkrho
     446              : !! 5) the fourdp chain :
     447              : !! driver - gstate - (move or brdmin) - scfcv - vtorho
     448              : !!         ( + ftofr - fourdp - symrhg )
     449              : !! 6) the parallel k-point chain :
     450              : !! driver - gstate - (move or brdmin) - scfcv - vtorho - MPI_ALLREDUCE
     451              : !! 7) the newvtr chain :
     452              : !! driver - gstate - (move or brdmin) - scfcv - newvtr
     453              : !! 8) the susceptibility chain :
     454              : !! driver - gstate - (move or brdmin) - scfcv - vtorho - suscep - suskmm
     455              : !! 9) the dielectric chain :
     456              : !! driver - gstate - (move or brdmin) - scfcv - vtorho - dielmt
     457              : !! 10) the tddft chain :
     458              : !! driver - gstate - (move or brdmin) - scfcv - vtorho - tddft
     459              : !!
     460              : !! It is valid for all values of iscf, but not for nstep=0 (when the chain
     461              : !! goes through energy instead of vtorho).
     462              : !!
     463              : !! Also, it is assumed that the potentials are non-local, even if there
     464              : !! are local ! It would be necessary to update this routine
     465              : !! now that the beginning of psp files is read before the present call (XG 980502)
     466              : !!
     467              : !! One might also estimate if there must be a chain arriving at:
     468              : !!  strnps , mkffnl, mkcore, mklocl, mkrho, prcpot, irrzg, initro, clnup1.
     469              : !! This is because there are allocated arrays in these routines.
     470              : !!
     471              : !! SOURCE
     472              : 
     473         6240 : subroutine memory(n1xccc,extrapwf,getcell,idtset,icoulomb,intxc,ionmov,iout,densfor_pred,iprcel,&
     474              : & iscf,jdtset,lmnmax,lnmax,&
     475              : & mband,mffmem,mgfft,mgfftdiel,mgfftf,mkmem,mpi_enreg,mpsang,mpssoang,mpw,mqgrid_ff,mqgrid_vl,&
     476         6240 : & natom,nband,nfft,nfftdiel,nfftf,ngfft,ngfftdiel,ngfftf,nimage,&
     477              : & nkpt,nloalg,npsp,npulayit,npwdiel,nspden,nspinor,nsppol,nsym,ntypat,&
     478              : & occopt,optforces,option,optstress,pawcpxocc,pawmixdg,pawnhatxc,pawspnorb,pawstgylm,&
     479         6240 : & prtvol,pspheads,qphon,tfkinfunc,typat,ucvol,usepaw,useylm,gpu_option,xclevel)
     480              : 
     481              : !Arguments ------------------------------------
     482              : !scalars
     483              :  integer,intent(in) :: extrapwf,getcell,icoulomb,idtset,intxc,ionmov,iout,densfor_pred
     484              :  integer,intent(in) :: iprcel,iscf,jdtset,lmnmax,lnmax,mband,mffmem,mgfft
     485              :  integer,intent(in) :: mgfftdiel,mgfftf,mkmem,mpsang,mpssoang,mpw,mqgrid_ff
     486              :  integer,intent(in) :: mqgrid_vl,n1xccc,natom,nfft,nfftdiel,nfftf,nimage,nkpt,npsp
     487              :  integer,intent(in) :: npulayit,npwdiel,nspden,nspinor,nsppol,nsym,ntypat
     488              :  integer,intent(in) :: occopt,optforces,option,optstress
     489              :  integer,intent(in) :: pawcpxocc,pawmixdg,pawnhatxc,pawspnorb,pawstgylm
     490              :  integer,intent(in) :: prtvol,tfkinfunc,usepaw,useylm,gpu_option,xclevel
     491              :  real(dp) :: ucvol
     492              :  type(MPI_type),intent(in) :: mpi_enreg
     493              : !arrays
     494              :  integer,intent(in) :: nband(nkpt*nsppol),ngfft(18),ngfftdiel(18),ngfftf(18)
     495              :  integer,intent(in) :: nloalg(3),typat(natom)
     496              :  real(dp),intent(in) :: qphon(3)
     497              :  type(pspheader_type) :: pspheads(npsp)
     498              : 
     499              : !Local variables-------------------------------
     500              : !marrays=maximal number of arrays to be monitored (or group of arrays)
     501              : !cmpw(marrays)=count of blocks of size mpw bytes
     502              : !cfft(marrays) =count of blocks of size nfft bytes (coarse grid, if PAW)
     503              : !cfftf(marrays)=count of blocks of size nfft bytes (fine grid, if PAW)
     504              : !cadd(marrays)=count of additional storage needed (in bytes)
     505              : !dttyp(marrays)=datatype of the array : 4 for integers, 8 for real(dp)
     506              : !nchain=number of different chains of routines
     507              : !chain(marrays,nchain)=different chains of routines
     508              :  ! The cfoo arrays are used to store the allocated memory in the different
     509              :  ! routines of the program. Each stack of the program can allocate some
     510              :  ! memory and the amount is estimated and stored in cfoo(i). The lower i,
     511              :  ! the higher routine. cfft is memory used by FFT handling, cmpw for
     512              :  ! plane waves storage and cadd for miscellaneous memory occupation.
     513              :  ! The unit is the multiplier of the size of nfft for cfft, the multiplier
     514              :  ! of mpw for cmpw and the actually allocated memory for cadd.
     515              :  ! This array stores the size of each chunk of memory (8 for double
     516              :  ! floating point precision, 4 for integers and so on).
     517              :  ! This array defines if the chain defined above allocate or not the
     518              :  ! memory (depending on options).
     519              : !scalars
     520              :  integer,parameter :: marrays=150,nchain=10
     521              :  integer :: fftalgb,histsz,ii,iscf10,jj,l_max,l_size_max,matblk,mblk,mincat,mu
     522              :  integer :: my_natom,n_fftgr,narr_fourdp,nbnd_in_blk,ndiel4,ndiel456,ndiel5,ndiel6
     523              :  integer :: ngrad,nprocwf,nspgrad,qphase_rhoij,rhoij_nspden
     524              :  real(dp) :: mbcg,mbdiskpd,mbdiskwf,mbf_fftgr,mbgylm
     525              :  character(len=500) :: msg
     526              : ! character(len=1) :: firstchar
     527              : !arrays
     528        12480 :  integer :: dttyp(marrays),nattyp(ntypat)
     529         6240 :  integer,allocatable :: basis_size(:),l_size(:),lmn2_size(:),lmn_size(:)
     530         6240 :  integer,allocatable :: mesh_size(:),my_nattyp(:),pawver(:),shape_type(:)
     531              :  real(dp) :: cadd(marrays),cfft(marrays),cfftf(marrays),cmpw(marrays)
     532         6240 :  real(dp),allocatable :: rshp(:)
     533              :  logical :: chain(marrays,nchain)
     534              : 
     535              : ! **************************************************************************
     536              : 
     537         6240 :  if(option<0 .or. option>2)then
     538            0 :    write(msg,'(A,I0,A)')'option=',option,' while the only allowed values are 0, 1, or 2.'
     539            0 :    ABI_BUG(msg)
     540              :  end if
     541              : 
     542              : !firstchar=' ';if (gpu_option/=0) firstchar='_'
     543         6240 :  cmpw(:)=zero ; cfft(:)=zero ; cfftf(:)=zero ; cadd(:)=zero
     544         6240 :  dttyp(:)=0
     545              : 
     546         6240 :  my_natom=natom;if (mpi_enreg%nproc_atom>1) my_natom=mpi_enreg%my_natom
     547              : 
     548         6240 :  call wrtout(std_out,'memory: analysis of memory needs ')
     549              : 
     550         6240 :  if(jdtset>=100)then
     551        20493 :    write(msg,'(80a,a,a,i5,a)')('=',mu=1,80),ch10,&
     552        20746 :     ' Values of the parameters that define the memory need for DATASET',jdtset,'.'
     553         5987 :  else if(jdtset/=0)then
     554       454815 :    write(msg,'(80a,a,a,i3,a)')('=',mu=1,80),ch10,&
     555       460430 :     ' Values of the parameters that define the memory need for DATASET',jdtset,'.'
     556              :  else
     557        30132 :    write(msg,'(80a,a,a)')('=',mu=1,80),ch10,&
     558        30504 :     ' Values of the parameters that define the memory need of the present run '
     559              :  end if
     560         6240 :  call wrtout(iout,msg)
     561         6240 :  call wrtout(std_out,msg)
     562              : 
     563              :  write(msg,'( 4(a,i8),a,4(a,i8) )' ) &
     564         6240 : & '     intxc =',intxc   ,'    ionmov =',ionmov,&
     565         6240 : & '      iscf =',iscf    ,'    lmnmax =',lmnmax,ch10,&
     566         6240 : & '     lnmax =',lnmax   ,'     mgfft =',mgfft,&
     567        12480 : & '  mpssoang =',mpssoang,'    mqgrid =',mqgrid_vl
     568         6240 :  call wrtout(iout,msg)
     569         6240 :  call wrtout(std_out,msg)
     570              : 
     571              :  write(msg,'( 4(a,i8),a,4(a,i8),a,4(a,i8) )' ) &
     572         6240 : & '     natom =',natom  ,'  nloc_mem =',nloalg(2)*(nloalg(3)+1),&
     573         6240 : & '    nspden =',nspden ,'   nspinor =',nspinor,ch10,&
     574         6240 : & '    nsppol =',nsppol ,'      nsym =',nsym,&
     575         6240 : & '    n1xccc =',n1xccc ,'    ntypat =',ntypat,ch10,&
     576        12480 : & '    occopt =',occopt ,'   xclevel =',xclevel
     577         6240 :  call wrtout(iout,msg)
     578         6240 :  call wrtout(std_out,msg)
     579              : 
     580              :  write(msg,'(4(3(a,i12),a))') &
     581         6240 : & '-    mband =',mband  ,'        mffmem =',mffmem,&
     582         6240 : & '         mkmem =',mkmem  ,ch10,&
     583         6240 : & '       mpw =',mpw    ,'          nfft =',nfft ,&
     584        12480 : & '          nkpt =',nkpt
     585         6240 :  call wrtout(iout,msg)
     586         6240 :  call wrtout(std_out,msg)
     587              : 
     588         6240 :  if (my_natom/=natom)then
     589          426 :    write(msg,'(a,i10)') 'Pmy_natom=',my_natom
     590          426 :    call wrtout(iout,msg)
     591          426 :    call wrtout(std_out,msg)
     592              :  end if
     593              : 
     594              : !Additional information if imgmov is activated (use of replicas of the cell)
     595         6240 :  if (nimage>1) then
     596           50 :    write(msg,'(1(a,i10))' ) '  nimage =',nimage
     597           50 :    call wrtout(iout,msg)
     598           50 :    call wrtout(std_out,msg)
     599              :  end if
     600              : 
     601              : !Additional information on FFT grids if PAW
     602         6240 :  if (usepaw==1) then
     603              :    write(msg, '(a,a,a,i10,a,i10)' )&
     604         1068 : &   ' PAW method is used; the additional fine FFT grid is defined by:',ch10,&
     605         2136 : &   '   mgfftf=',mgfftf,'    nfftf =',nfftf
     606         1068 :    call wrtout(iout,msg)
     607         1068 :    call wrtout(std_out,msg)
     608              :  end if
     609              : 
     610              : !Additional information if GPU
     611              :  if (gpu_option/=ABI_GPU_DISABLED) then
     612              : !  write(msg, '(a)' )' GPU method is used'
     613              : !  call wrtout(iout,msg)
     614              : !  call wrtout(std_out,msg)
     615              :  end if
     616              : 
     617              : !Additional information needed for the susceptibility and dielectric matrices
     618         6240 :  if((modulo(iprcel,100)>=20.and.modulo(iprcel,100)<70) .or. iscf==-1)then
     619              : 
     620              : !  Compute the number of bands in blocks (nbnd_in_blk) from mband (see suskmm.f)
     621              : !  Consider that if the number of bands is large, there are at most 8 blocks
     622           61 :    if(mband>=48)then
     623            0 :      mblk=8
     624            0 :      nbnd_in_blk=(mband-1)/mblk+1
     625              : !    If the number of bands is medium, place 6 bands per block
     626           61 :    else if(mband>=12)then
     627           47 :      nbnd_in_blk=6
     628              : !    Otherwise, must have at least 2 blocks
     629              :    else
     630           14 :      mblk=2
     631           14 :      nbnd_in_blk=(mband-1)/mblk+1
     632              :    end if
     633              : 
     634              :    write(msg, '(a,a,a,i10,a,i6,a,i10,a,i10)' )&
     635           61 : &   ' For the susceptibility and dielectric matrices, or tddft :',ch10,&
     636           61 : &   '   mgfft =',mgfftdiel,'  nbnd_in_blk=',nbnd_in_blk,'    nfft =',nfftdiel,&
     637          122 : &   '     npw =',npwdiel
     638           61 :    call wrtout(iout,msg)
     639           61 :    call wrtout(std_out,msg)
     640           61 :    ndiel4=ngfftdiel(4) ; ndiel5=ngfftdiel(5) ; ndiel6=ngfftdiel(6)
     641           61 :    ndiel456=ndiel4*ndiel5*ndiel6
     642              :  else
     643              : !  To be sure of initialisation.
     644              :    ndiel456 = 1
     645              :  end if
     646              : 
     647       505440 :  write(msg,'(80a)') ('=',mu=1,80)
     648         6240 :  call wrtout(iout,msg)
     649         6240 :  call wrtout(std_out,msg)
     650              : 
     651         6240 :  if(getcell>0 .or. (getcell<0 .and. idtset+getcell>0) )then
     652           11 :    write(msg,'(a,a,a,a,a,a,i3,a,i3,a,a,a,a,a,a)' )ch10,&
     653           11 : &   ' memory : COMMENT -',ch10,&
     654           11 : &   '  The determination of memory needs at this stage is meaningless,',ch10,&
     655           11 : &   '  since getcell = ',getcell,' is non-zero, while idtset=',idtset,'.',ch10,&
     656           11 : &   '  The following numbers are obtained by supposing that acell and rprim',ch10,&
     657           22 : &   '  are NOT taken from a previous dataset. You cannot rely on them.',ch10
     658           11 :    call wrtout(iout,msg)
     659           11 :    call wrtout(std_out,msg)
     660              :  end if
     661              : 
     662              : !Compute number of atoms per type for current proc
     663        13749 :  nattyp(:)=0
     664        42667 :  do ii=1,natom
     665        42667 :    nattyp(typat(ii))=nattyp(typat(ii))+1
     666              :  end do
     667              : 
     668              : !PAW: store useful dims
     669         6240 :  if (usepaw==1) then
     670         3204 :    ABI_MALLOC(basis_size,(npsp))
     671         2136 :    ABI_MALLOC(l_size,(npsp))
     672         2136 :    ABI_MALLOC(lmn_size,(npsp))
     673         2136 :    ABI_MALLOC(lmn2_size,(npsp))
     674         2136 :    ABI_MALLOC(mesh_size,(npsp))
     675         2136 :    ABI_MALLOC(shape_type,(npsp))
     676         2136 :    ABI_MALLOC(pawver,(npsp))
     677         3204 :    ABI_MALLOC(rshp,(npsp))
     678         2646 :    do ii=1,npsp
     679         1578 :      basis_size(ii)=pspheads(ii)%pawheader%basis_size
     680         1578 :      mesh_size(ii)=pspheads(ii)%pawheader%mesh_size
     681         1578 :      l_size(ii)=pspheads(ii)%pawheader%l_size
     682         1578 :      lmn_size(ii)=pspheads(ii)%pawheader%lmn_size
     683         1578 :      lmn2_size(ii)=lmn_size(ii)*(lmn_size(ii)+1)/2
     684         1578 :      pawver(ii)=pspheads(ii)%pawheader%pawver
     685         1578 :      rshp(ii)=pspheads(ii)%pawheader%rshp
     686         2646 :      shape_type(ii)=pspheads(ii)%pawheader%shape_type
     687              :    end do
     688         2646 :    l_max=maxval(pspheads(:)%lmax)
     689         2646 :    l_size_max=maxval(pspheads(:)%pawheader%l_size)
     690         1068 :    rhoij_nspden=nspden;if (pawspnorb>0) rhoij_nspden=4
     691         3204 :    ABI_MALLOC(my_nattyp,(ntypat))
     692         1068 :    if ((mpi_enreg%nproc_atom<=1).or.(.not.associated(mpi_enreg%my_atmtab))) then
     693         2234 :      my_nattyp=nattyp
     694              :    else
     695         1054 :      my_nattyp=0
     696          806 :      do ii=1,my_natom
     697          380 :        jj=typat(mpi_enreg%my_atmtab(ii))
     698          806 :        my_nattyp(jj)=my_nattyp(jj)+1
     699              :      end do
     700              :    end if
     701         4272 :    qphase_rhoij=merge(2,1,any(qphon(:)>tol8))
     702              :  else
     703              : !  Do the allocation to avoid uninitialised variables.
     704         5172 :    ABI_MALLOC(my_nattyp,(1))
     705         5172 :    ABI_MALLOC(basis_size,(1))
     706         5172 :    ABI_MALLOC(l_size,(1))
     707         5172 :    ABI_MALLOC(lmn_size,(1))
     708         5172 :    ABI_MALLOC(lmn2_size,(1))
     709         5172 :    ABI_MALLOC(mesh_size,(1))
     710         5172 :    ABI_MALLOC(shape_type,(1))
     711         5172 :    ABI_MALLOC(pawver,(1))
     712         5172 :    ABI_MALLOC(rshp,(1))
     713         5172 :    rhoij_nspden=nspden
     714         5172 :    l_size_max=1
     715         5172 :    l_max=1
     716         5172 :    qphase_rhoij=1
     717              :  end if
     718              : 
     719         6240 :  n_fftgr=1;iscf10=mod(iscf,10)
     720         6240 :  if(iscf10==1)              n_fftgr=5
     721         6237 :  if(iscf10==2)              n_fftgr=3
     722         6227 :  if(iscf10==3)              n_fftgr=4
     723         6224 :  if(iscf10==4)              n_fftgr=6
     724         6238 :  if(iscf10==5.or.iscf10==6) n_fftgr=10
     725         6216 :  if(iscf10==7)              n_fftgr=2+2*npulayit
     726              : 
     727              : !work1 and work2 in fourdp : take into account approximately fftalgb
     728         6240 :  fftalgb=mod(ngfft(7),100)/10
     729              :  if(fftalgb==0)narr_fourdp=2*2
     730         6240 :  if(fftalgb==1)narr_fourdp=2
     731              : 
     732         6240 :  ngrad=1;if(xclevel==2.or.tfkinfunc>10)ngrad=2
     733              : 
     734              : !(1)                     in main, driver, gstate and brdmin ----------------
     735              : !in move, nothing interesting is allocated.
     736              : !kg (gstate)
     737         6240 :  cmpw(1)=3*mkmem               ; dttyp(1)=4
     738              : !indsym (gstate)
     739         6240 :  cadd(3)=4*nsym*natom          ; dttyp(3)=4
     740              : !irrzon  (gstate)
     741         6240 :  if(nsym/=1)then
     742         5385 :    cfft(4)=2*((nspden/nsppol)-3*(nspden/4))    ; dttyp(4)=4
     743              :  end if
     744              : !ylm (gstate)
     745         6240 :  cmpw(5)=mkmem*mpsang*mpsang*useylm ; dttyp(5)=8
     746              : !
     747              : !rhor,rhog (gstate)
     748         6240 :  cfftf(5)=nspden+2              ; dttyp(5)=8
     749              : !cg (gstate)
     750         6240 :  cmpw(6)=2*nspinor*mband*mkmem*nsppol  ; dttyp(6)=8
     751              : !eigen,resid,occ (occ is initialized in abinit, and not in driver)
     752         6240 :  cadd(7)=3*mband*nkpt*nsppol   ; dttyp(7)=8
     753              : !qgrid_vl,qgrid_ff,vlspl,ffspl
     754              :  cadd(8)=mqgrid_vl*(1+2*ntypat)   &
     755              : & +mqgrid_ff*(1+2*ntypat*lnmax)   &
     756         6240 : & ; dttyp(8)=8
     757              : !ph1d (actually allocated in scfcv !!)
     758         6240 :  cadd(9)=2*3*(2*mgfft+1)*natom ; dttyp(9)=8
     759         6240 :  cadd(9)=cadd(9)+2*3*(2*mgfftf+1)*natom*usepaw  !Additional ph1df for PAW
     760              : !phnons (in gstate)
     761         6240 :  if(nsym/=1)then
     762         5385 :    cfft(10)=2*((nspden/nsppol)-3*(nspden/4))    ; dttyp(10)=8
     763              :  end if
     764              : !xccc1d (in driver)
     765         6240 :  cadd(11)=n1xccc*6*ntypat      ; dttyp(11)=8
     766              : 
     767              : !hessin in brdmin
     768         6240 :  if(ionmov==2)then
     769          162 :    cadd(15)=3*natom*3*natom      ; dttyp(15)=8
     770              :  end if
     771              : 
     772              : !Additional PAW arrays
     773              : !PAW datasets (pawtab)
     774         6240 :  if (usepaw==1) then
     775         1068 :    dttyp(16)=8 ; dttyp(17)=4
     776         2646 :    do ii=1,npsp
     777         1578 :      cadd(16)=cadd(16)+2*mesh_size(ii)*basis_size(ii)   !phi,tphi
     778              :      cadd(16)=cadd(16)+2*mesh_size(ii)*basis_size(ii)&  !phiphj,tphiphj
     779         1578 : &    *(basis_size(ii)+1)/2
     780         1578 :      cadd(16)=cadd(16)+mesh_size(ii)*l_size(ii)         !shapefunc
     781         1578 :      cadd(16)=cadd(16)+lmn2_size(ii)*l_size(ii)**2      !qijl
     782         1578 :      cadd(16)=cadd(16)+l_size(ii)*5                     !gnorm,shape_a,shape_q
     783         1578 :      cadd(16)=cadd(16)+lmn2_size(ii)*(4+lmn2_size(ii))  !eijkl,dltij,dij0,rhoij0,sij
     784         1578 :      cadd(17)=cadd(17)+lmn2_size(ii)*8                  !indklmn
     785         1578 :      cadd(16)=cadd(16)+mesh_size(ii)*5                  !coreden,tcoreden,rad,radfact,simfact
     786         1578 :      if (shape_type(ii)==-1) cadd(16)=cadd(16)+4*mesh_size(ii)*l_size(ii)  !dshpfunc
     787         1578 :      cadd(16)=cadd(16)+mqgrid_vl*2                      !tncorespl
     788         2646 :      if (pawver(ii)>=4) cadd(16)=cadd(16)+mqgrid_vl*2   !tnvalespl
     789              :    end do
     790              : !  additional arrays
     791         1068 :    cadd(16)=cadd(16)+l_size_max*2*l_max*nsym                 !zarot
     792         1068 :    cadd(16)=cadd(16)+(2*l_max-1)**2*l_max**2*(l_max**2+1)/2  !realgnt
     793         1068 :    cadd(17)=cadd(17)+nfft+nfftf                              ! fintocoa,coatofin
     794         2646 :    do ii=1,ntypat
     795         1578 :      cadd(16)=cadd(16)+my_nattyp(ii)*lmn2_size(ii)*rhoij_nspden*pawcpxocc ! Rhoij and related data
     796         2646 :      cadd(17)=cadd(17)+my_nattyp(ii)*(2+lmn2_size(ii))    ! (rhoijselect, ...)
     797              :    end do
     798              :    !PAW:cprj
     799         2646 :    do ii=1,ntypat
     800         2646 :      cadd(16)=cadd(16)+2*nattyp(ii)*nkpt*nspinor*mband*nsppol*lmn_size(ii)/max(mpi_enreg%nproc_band,1)
     801              :    end do
     802              :  end if
     803              : 
     804              : !SCF history (if selected)
     805         6240 :  if (abs(densfor_pred)==5.or.abs(densfor_pred)==6) then          ! scf_history...
     806          624 :    histsz=2
     807          624 :    cfftf(18)=nspden*(histsz+1)+1      ; dttyp(18)=8  ! %deltarhor, %atmrho_last, %rhor_last
     808          624 :    cadd(19)=3*natom*2*histsz          ; dttyp(19)=8  ! %xreddiff,xred_last
     809          624 :    dttyp(20)=4
     810          624 :    if (usepaw==1) then
     811         1044 :      do ii=1,ntypat
     812          619 :        cadd(19)=cadd(19)+histsz*2*my_nattyp(ii)*lmn2_size(ii)*rhoij_nspden*qphase_rhoij*pawcpxocc ! %pawrhoij()%rhoijp
     813         1044 :        cadd(20)=cadd(20)+histsz*2*my_nattyp(ii)*(2+lmn2_size(ii))*nspden ! %pawrhoij()%rhoijselect
     814              :      end do
     815              :    end if
     816          624 :    if (extrapwf>0) then
     817            1 :      cadd(19)=cadd(19)+histsz*2*nspinor*mband*mkmem*nsppol  ; dttyp(19)=8  ! %cg
     818              :    end if
     819              :  end if
     820              : 
     821              : !(2)                     in scfcv----------------------------------------
     822              : 
     823              : !vhartr,vpsp,vtrial,vxc
     824         6240 :  cfftf(21)=2+2*nspden           ; dttyp(21)=8
     825              : !kxc
     826         6240 :  if (abs(densfor_pred)>0.and.iscf>=10) then
     827          999 :    cfftf(21)=cfftf(21)+3*nspden
     828          999 :    if (densfor_pred<0.and.xclevel==2) cfftf(21)=cfftf(21)+20*nspden
     829              :  end if
     830         6240 :  if(iscf>0)then
     831              : !  f_fftgr
     832         5726 :    if (pawmixdg==1) then
     833          188 :      cfftf(22)=nspden*n_fftgr*mffmem; dttyp(22)=8
     834              :    else
     835         5538 :      cfft(22)=nspden*n_fftgr*mffmem; dttyp(22)=8
     836              :    end if
     837              :  end if
     838         6240 :  if( iscf>0 .and. (modulo(iprcel,100)>=20.and.modulo(iprcel,100)<70))then
     839              : !  dielinv, susmat
     840           39 :    cadd(23)=4*(npwdiel*min(nspden,2))**2; dttyp(23)=8
     841              :  end if
     842              : !Kernel of Poisson's solver
     843         6240 :  if (icoulomb == 1) then
     844            0 :    cadd(24) = ngfft(4)*ngfft(5)*ngfft(6) ; dttyp(24) = 8
     845              :  end if
     846         6240 :  if( (iscf>0 .and. modulo(iprcel,100)>=20 .and. modulo(iprcel,100)<70) .or. iscf==-1 )then
     847              : !  kg_diel
     848           61 :    cadd(27)=3*npwdiel             ; dttyp(27)=4
     849           61 :    if(nsym/=1)then
     850              : !    irrzondiel
     851           47 :      cadd(27)=cadd(27)+2*nfftdiel*(nspden/nsppol)
     852              : !    phnonsdiel
     853           47 :      cadd(28)=2*nfftdiel*(nspden/nsppol)   ; dttyp(28)=8
     854              :    end if
     855              :  end if
     856         6240 :  if(n1xccc/=0)then
     857              : !  xccc3d
     858         3153 :    cfftf(29)=1                    ; dttyp(29)=8
     859              :  end if
     860              : 
     861              : !Additional PAW arrays
     862         6240 :  dttyp(25)=8 ; dttyp(26)=4
     863         6240 :  if (usepaw==1) then
     864         2646 :    do ii=1,ntypat
     865         1578 :      jj=(1+int(nfftf*four_pi*rshp(ii)**3/(three*ucvol)))        ! pawfgrtab
     866         1578 :      cadd(26)=cadd(26)+my_nattyp(ii)*jj                         !   %ifftsph
     867         1578 :      cadd(25)=cadd(25)+my_nattyp(ii)*jj*(1-pawstgylm)*3         !   %rfgd (if pawstgylm=0)
     868         1578 :      cadd(25)=cadd(25)+my_nattyp(ii)*jj*pawstgylm*l_size(ii)**2 !   %gylm (if pawstgylm=1)
     869         1578 :      if (optforces==1) cadd(25)=cadd(25)+my_nattyp(ii)*jj&      !   %gylmgr,%rfgd (if pawstgylm=1)
     870          401 : &    *pawstgylm*(3*l_size(ii)**2+3*optstress)
     871         1578 :      cadd(26)=cadd(26)+my_nattyp(ii)*l_size(ii)**2/32           ! lmselect  !now a boolean
     872         1578 :      cadd(25)=cadd(25)+my_nattyp(ii)*lmn2_size(ii)*nspinor**3   ! dij
     873         2646 :      if (iscf>0) then
     874         1549 :        cadd(25)=cadd(25)+my_nattyp(ii)*lmn2_size(ii)*rhoij_nspden*pawcpxocc                ! rhoijres
     875         1549 :        cadd(25)=cadd(25)+my_nattyp(ii)*lmn2_size(ii)*rhoij_nspden*pawcpxocc*n_fftgr*mffmem ! f_paw
     876              :      end if
     877              :    end do
     878              : !   cadd(25)=cadd(25)+(1+3*pawnhatxc*(ngrad/2))*nspden*nfftf       !nhat,nhatgr
     879         1068 :    cfftf(29)=cfftf(29)+(1+3*pawnhatxc*(ngrad/2))*nspden       !nhat,nhatgr
     880              :  end if
     881              : 
     882              : !(3)                     in rhotoxc, xcden -------------------------------
     883              : 
     884         6240 :  if(xclevel/=0)then
     885         6194 :    if(n1xccc/=0)then
     886              : !    rhocorval
     887         3136 :      cfftf(31)=nspden               ; dttyp(31)=8
     888              :    end if
     889              : !  dnexcdn, rhonow
     890         6194 :    nspgrad=nspden*ngrad
     891         6194 :    if(nspden==2 .and. ngrad==2)nspgrad=5
     892         6194 :    cfftf(32)=nspden*ngrad*ngrad+nspgrad  ; dttyp(32)=8
     893         6194 :    if(intxc==1 .or. ngrad==2)then
     894              : !    wkcmpx,work in xcden +work1,work2 in fourdp
     895         1851 :      cfftf(33)=3+narr_fourdp        ; dttyp(33)=8
     896         1851 :      cadd(33)=narr_fourdp*(ngfftf(4)*ngfftf(5)*ngfftf(6)-nfftf)
     897              :    end if
     898         1851 :    if(ngrad==2)then
     899              : !    workgr in xcden
     900         1313 :      cfftf(34)=2                    ; dttyp(34)=8
     901              :    end if
     902              :  end if
     903         6240 :  if(iscf>0)then
     904              : !  In this case, rhotoxc is called from rhotov also,
     905              : !  for which vresid was allocated in scfcv
     906              : !  vresid
     907         5726 :    cfftf(35)=nspden               ; dttyp(35)=8
     908              :  end if
     909              : !Poisson's solver with zero padding
     910         6240 :  if (icoulomb == 1) then
     911            0 :    cfft(36) = 8                   ; dttyp(36) = 8
     912            0 :    cadd(36) = ngfft(4) * ngfft(5) * ngfft(6) - nfft
     913              :  end if
     914              : 
     915              : !Note : in hartre, called by rhotoxc, one uses
     916              : !2 dp arrays of total size 3*nfft,
     917              : !and 2 arrays of total size 4*n4*n5*n6 for fourdp
     918              : !This will be smaller than the total use for symrhg
     919              : 
     920              : !(4)                     in newvtr/newrho --------------------------------------
     921              : 
     922         6240 :  if(iscf>0)then
     923              : !  vresid (allocated in scfcv) and vrespc
     924         5726 :    if (pawmixdg==1) then
     925          188 :      cfftf(41)=2*nspden             ; dttyp(41)=8
     926              :    else
     927         5538 :      cfft(41)=2*nspden             ; dttyp(41)=8
     928              :    end if
     929         5726 :    if(mffmem==0)then
     930              : !    f_fftgr_disk
     931            0 :      if (pawmixdg==1) then
     932            0 :        cfftf(42)=nspden*n_fftgr       ; dttyp(42)=8
     933              :      else
     934            0 :        cfft(42)=nspden*n_fftgr       ; dttyp(42)=8
     935              :      end if
     936              : !    f_paw_disk
     937            0 :      if (usepaw==1) then
     938            0 :        dttyp(43)=8
     939            0 :        do ii=1,ntypat
     940            0 :          cadd(43)=cadd(43)+my_nattyp(ii)*lmn2_size(ii)*nspden*n_fftgr
     941              :        end do
     942              :      end if
     943              :    end if
     944              : !  rhoupdn, n(v)resid0, vtrialg, rhog2, magng
     945         5726 :    if (pawmixdg==1) then
     946          188 :      cfftf(43)=2*nspden       ; dttyp(43)=8
     947              :    else
     948         5538 :      cfft(43)=2*nspden       ; dttyp(43)=8
     949         5538 :      if (nspden>1) cfftf(43)=2*(nspden-1)
     950              :    end if
     951              :  end if
     952              : 
     953              : !(5-6)                   in vtorho-----------------------------------------
     954              : 
     955              : !Note : (5) is for the arrays inside the spin and k-point loop
     956              : !they belong to the main chain
     957              : !(6) is for the arrays after the spin and k-point loop
     958              : !(6a) is for the arrays after that loop, for the parallel k-point chain
     959              : !(6b) is for the arrays in mkrho, for the mkrho chain
     960              : !(6c) is for the arrays in symrhg, for the fourdp chain
     961              : !(6d) is for the arrays in suscep, for the suscep chain, see (10)
     962              : !(6e) is for the arrays in dielmt, for the dielmt chain, see (11)
     963              : !(6f) is for the arrays in pawmkrhoij
     964              : 
     965              : !eknlk, enlxnk, grnlnk
     966              :  cadd(51)=(11+3*natom)*mband*nkpt*nsppol &
     967         6240 : & ; dttyp(51)=8
     968              : !kg_k
     969         6240 :  cmpw(52)=3                    ; dttyp(52)=4
     970              : !rhoaug,vlocal
     971         6240 :  cfft(53)=2                    ; dttyp(53)=8
     972         6240 :  cadd(53)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
     973              : !rhowfr,rhowfg
     974         6240 :  cfft(53)=cfft(53)+2+nspden
     975         6240 :  if(mkmem==0)then
     976              : !  cg_disk
     977            0 :    cmpw(54)=2*nspinor*mband      ; dttyp(54)=8
     978              :  end if
     979              : !eig_k, ek_k, enlx_k, grnl_k, occ_k, resid_k
     980         6240 :  cadd(56)=(14+3*natom)*mband   ; dttyp(56)=8
     981              : !ylm_k
     982         6240 :  cmpw(57)=mpsang*mpsang*useylm ; dttyp(57)=8
     983              : !!PAW:cprj
     984              : ! if (usepaw==1) then
     985              : !   dttyp(58)=8
     986              : !   do ii=1,ntypat
     987              : !     cadd(58)=cadd(58)+2*nattyp(ii)*nkpt*nspinor*mband*nsppol*lmn_size(ii)/max(mpi_enreg%nproc_band,1)
     988              : !   end do
     989              : ! end if
     990              : 
     991              : !(6)                     in vtorho----------------------------------------
     992              : 
     993              : !doccde
     994         6240 :  cadd(60)=mband*nkpt*nsppol    ; dttyp(60)=8
     995              : 
     996              : !(6a)                    in vtorho----------------------------------------
     997              :  if(xmpi_paral==1)then
     998              : !  Parallel case
     999              : !  buffer1
    1000              : !  buffer2
    1001         6240 :    if(occopt>=3 .and. occopt <=8) then
    1002         2364 :      dttyp(61)=8
    1003         2364 :      if(nsppol*nfft >= (13+3*natom)*mband*nkpt*nspden)then
    1004         1554 :        cfft(61)=2*nspden
    1005              :      else
    1006          810 :        cadd(61)=(13+3*natom)*mband*nkpt*nspden
    1007              :      end if
    1008              :    else
    1009         3876 :      cfft(61)=2*nspden             ; dttyp(61)=8
    1010         3876 :      cadd(61)=9+3*natom+2+2*mband*nkpt*nspden
    1011              :    end if
    1012              :  end if
    1013              : 
    1014              : 
    1015              : !(6b)                    in mkrho, called by vtorho--------------------------
    1016         6240 :  if(occopt>=3 .and. occopt <=8)then
    1017         2364 :    if(mkmem==0)then
    1018              : !    cg_disk
    1019            0 :      cmpw(62)=2*nspinor*mband      ; dttyp(62)=8
    1020              :    end if
    1021              : !  cwavef
    1022         2364 :    cmpw(65)=2*nspinor            ; dttyp(65)=8
    1023              : 
    1024              : !  rhoaug, wfraug, work1 in fourwf
    1025         2364 :    cfft(66)=5                    ; dttyp(66)=8
    1026         2364 :    cadd(66)=5*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    1027              :  end if
    1028              : 
    1029              : !(6c)                    in symrhg, called by vtorho--------------------------
    1030         6240 :  if(iscf>0)then
    1031         5726 :    cfft(67)=narr_fourdp          ; dttyp(67)=8
    1032         5726 :    cadd(67)=narr_fourdp*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    1033         5726 :    if(nsym>1)then
    1034              : !    work1  in symrhg
    1035         4959 :      cfft(68)=2                    ; dttyp(68)=8
    1036         4959 :      cadd(68)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    1037              :    end if
    1038              :  end if
    1039              : 
    1040              : 
    1041              : !(6d) and (6e)           in suscep and dielmt, called by vtorho,
    1042              : !see (10) and (11) -------------------------------
    1043              : 
    1044              : !(6f)  in pawmkrhoij or pawrhoij_symrhoij called by pawmkrho, called by vtorho--------
    1045              : !only when paralellim over atoms is activated
    1046         6240 :  dttyp(63)=8
    1047         6240 :  if((usepaw==1) .and. ((iscf>0) .or. (iscf == -3) .and. mpi_enreg%nproc_atom>1 ))then
    1048         2594 :    do ii=1,ntypat
    1049         1549 :      cadd(63)=cadd(63)+nattyp(ii)*lmn2_size(ii)*rhoij_nspden*pawcpxocc*qphase_rhoij ! Rhoij_gather and related data
    1050         7789 :      cadd(63)=cadd(63)+nattyp(ii)*(2+lmn2_size(ii)) ! Rhoij_gather (rhoijselect, ...)
    1051              :    end do
    1052              :  end if
    1053              : 
    1054              : !(7)                     in vtowfk----------------------------------------
    1055              : 
    1056              : !evec
    1057         6240 :  cadd(71)=2*mband*mband        ; dttyp(71)=8
    1058              : !subham, subvnlx(if not PAW or if usefock_ACE)
    1059         6240 :  cadd(72)=(1+usepaw)*mband*(mband+1)    ; dttyp(72)=8
    1060              : !gkpsq
    1061         6240 :  cmpw(73)=1                    ; dttyp(73)=8
    1062              : !ffnl
    1063         6240 :  cmpw(74)=2*ntypat*lmnmax      ; dttyp(74)=8
    1064              : !ph3d
    1065        13749 :  matblk=min(NLO_MINCAT,maxval(nattyp))
    1066         6240 :  if(nloalg(2)<=0)matblk=natom
    1067         6240 :  cmpw(75)=2*matblk             ; dttyp(75)=8
    1068              : !gsc(if PAW)
    1069              : ! cmpw(76)=2*mband*nspinor*usepaw          ; dttyp(76)=8
    1070              : !Note : matvnl and mat1 do not belong to a chain defined until now
    1071              : !
    1072         6240 :  if(occopt<3 .and. iscf>0)then
    1073              : !  cwavef
    1074         3430 :    cmpw(77)=2*nspinor            ; dttyp(77)=8
    1075              : !  wfraug
    1076         3430 :    cfft(78)=2                    ; dttyp(78)=8
    1077         3430 :    cadd(78)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    1078              : !  work1 in fourwf
    1079         3430 :    cfft(79)=2                    ; dttyp(79)=8
    1080         3430 :    cadd(79)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    1081              :  end if
    1082              : 
    1083              : 
    1084              : !(8)                     in cgwf_cprj-------------------------------------
    1085              : 
    1086              : !conjgr, direc, direc_tmp, gvnlx
    1087         6240 :  cmpw(81)=2*4*nspinor          ; dttyp(81)=8
    1088              : ! cwavef_r,direc_r
    1089         6240 :  cfft(82)=2*2*nspinor          ; dttyp(82)=8
    1090              : 
    1091              : !(8)                     in cgwf------------------------------------------
    1092              : 
    1093              : !!conjgr, cwavef, direc, gh_direc, gvnlx_direc
    1094              : ! cmpw(81)=2*5*nspinor          ; dttyp(81)=8
    1095              : !!ghc,gvnlxc
    1096              : ! cmpw(82)=2*2*nspinor          ; dttyp(82)=8
    1097              : !!PAW: scwavef,direc_tmp,ghc_all
    1098              : ! cmpw(83)=2*(2+mband)*nspinor*usepaw  ; dttyp(83)=8
    1099              : 
    1100              : 
    1101              : !(9a)                    in getghc and fourwf----------------------------
    1102              : 
    1103              : !work (in getghc)
    1104         6240 :  cfft(91)=2                    ; dttyp(91)=8
    1105         6240 :  cadd(92)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    1106              : !work1 (in fourwf)
    1107         6240 :  cfft(92)=2                    ; dttyp(92)=8
    1108              :  cadd(92)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    1109              : 
    1110              : !(9b)                    in getghc, nonlop and opernl--------------------
    1111         6240 :  mincat=min(NLO_MINCAT,natom-ntypat+1)
    1112         6240 :  if (useylm==0) then                          ! ===== nonlop_pl
    1113              : !  gxa  (in nonlop)
    1114         5098 :    cadd(94)=2*20*mincat*2       ; dttyp(94)=8
    1115              : !  dgxdt  (in nonlop)            !MT20072002: not allocated in getghc !!
    1116         5098 :    if (optforces==1) then
    1117          388 :      cadd(95)=2*3*20*mincat*2    ; dttyp(95)=8
    1118              :    end if
    1119              : !  teffv (in opernl4 - no distinction is made for opernl, opernl2 or opernl3)
    1120              : !  kpgx, ffkg
    1121              : !  here, evaluate an upper value, with nproj=2, p,d and f orbitals, but not
    1122              : !  considering the stress, since it will be called outside of the main chain
    1123         5098 :    cadd(97)=NLO_MBLKPW*40        ; dttyp(97)=8
    1124              : !  kpg if nloalg(3)=1
    1125         5098 :    cadd(98)=3*mpw*nloalg(3)      ; dttyp(98)=8
    1126              :  else                                        ! ===== nonlop_ylm
    1127              : !  gx + gxfac + gxfac_sij
    1128              : !   cadd(94)=2*lmnmax*mincat*(mpw+1+usepaw)    ; dttyp(94)=8
    1129         1142 :    cmpw(94)=2*lmnmax*mincat                   ; dttyp(94)=8
    1130         1142 :    cadd(99)=2*lmnmax*mincat*(1+usepaw)        ; dttyp(99)=8
    1131              : !  kpg
    1132         1142 :    cmpw(95)=3             ; dttyp(95)=8
    1133              : !  indlmn_typ, ffnl_typ
    1134         1142 :    cadd(96)=lmnmax*6; dttyp(96)=4
    1135              : !  ffnl_typ
    1136         1142 :    cmpw(97)=lmnmax; dttyp(97)=8
    1137              : !  opernla_ylm: scalar,scali
    1138         1142 :    cmpw(98)=2; dttyp(98)=8
    1139              :  end if
    1140              : 
    1141              : !(10)                    in suscep and suskmm ----------------------------
    1142              : 
    1143         6240 :  if(modulo(iprcel,100)>=20.and.modulo(iprcel,100)<70)then
    1144              : !  Variables allocated in suscep
    1145           39 :    if(mkmem==0)then
    1146              : !    cg_disk
    1147            0 :      cmpw(101)=2*mband             ; dttyp(101)=8
    1148              :    end if
    1149           39 :    if(occopt>=3)then
    1150              : !    drhode
    1151           24 :      cadd(103)=2*npwdiel*nsppol    ; dttyp(103)=8
    1152              :    end if
    1153              : !  rhoextrap (always included, although it appears only when extrap==1)
    1154           39 :    cadd(104)=ndiel456            ; dttyp(104)=8
    1155              : 
    1156              : !  Variables allocated in suskmm
    1157              : !  cwavef
    1158           39 :    cmpw(106)=2                   ; dttyp(106)=8
    1159              : !  rhoaug, wfraug
    1160           39 :    cadd(107)=3*ndiel456          ; dttyp(107)=8
    1161              : !  wfprod
    1162           39 :    cadd(108)=2*npwdiel           ; dttyp(108)=8
    1163              : !  wfrspa1, wfrspa2
    1164           39 :    cadd(109)=4*ndiel456*nbnd_in_blk ; dttyp(109)=8
    1165              : 
    1166              :  end if
    1167              : 
    1168              : !(11)                    in dielmt ---------------------------------------
    1169              : 
    1170         6240 :  if(modulo(iprcel,100)>=20.and.modulo(iprcel,100)<70)then
    1171              : !  dielh,dielvec,eig_diel,zhpev1,zhpev2
    1172              :    cadd(111)=3*npwdiel*npwdiel                   &
    1173           39 : &   +9*npwdiel           ; dttyp(111)=8
    1174              :  end if
    1175              : 
    1176              : !(12)                    in tddft  ---------------------------------------
    1177              : 
    1178         6240 :  if(iscf==-1)then
    1179           22 :    if(mkmem/=0)then
    1180              : !    cg_disk
    1181           22 :      cmpw(121)=2*mband            ; dttyp(121)=8
    1182              :    end if
    1183              : !  cwavef
    1184           22 :    cmpw(124)=2*mband             ; dttyp(124)=8
    1185              : !  rhoaug,wfraug,wfrspa
    1186           22 :    cadd(125)=(2+mband)*ndiel456  ; dttyp(125)=8
    1187              :  end if
    1188              : 
    1189              : !--------------------------------------------------------------------------
    1190              : 
    1191      9428640 :  chain(:,:)=.true.
    1192              : 
    1193              : !Define the main chain version a (fourwf)
    1194       131040 :  chain(31:50,1)=.false.
    1195        74880 :  chain(60:70,1)=.false.
    1196        31200 :  chain(77:80,1)=.false.
    1197        56160 :  chain(93:100,1)=.false.
    1198       318240 :  chain(101:marrays,1)=.false.
    1199              : 
    1200              : !Define the main chain version b (nonlop+opernl)
    1201       131040 :  chain(31:50,2)=.false.
    1202        74880 :  chain(60:70,2)=.false.
    1203        31200 :  chain(77:80,2)=.false.
    1204        18720 :  chain(91:92,2)=.false.
    1205       318240 :  chain(101:marrays,2)=.false.
    1206              : 
    1207              : !Define the XC chain ( 31:40 belong only to this chain)
    1208       692640 :  chain(41:marrays,3)=.false.
    1209              : 
    1210              : !Define the mkrho chain ( 62:66 and 76:77 belong only to this chain)
    1211              : !is it sure that they have to be summed ?)
    1212       131040 :  chain(31:50,4)=.false.
    1213        62400 :  chain(51:59,4)=.false.
    1214         6240 :  chain(61   ,4)=.false.
    1215        31200 :  chain(67:70,4)=.false.
    1216       505440 :  chain(71:marrays,4)=.false.
    1217        31200 :  chain(77:80,4)=.true.
    1218              : 
    1219              : !Define the fourdp chain ( 67:70 belong only to this chain)
    1220       131040 :  chain(31:50,5)=.false.
    1221       106080 :  chain(51:66,5)=.false.
    1222         6240 :  chain(60   ,5)=.true.
    1223       505440 :  chain(71:marrays,5)=.false.
    1224              : 
    1225              : !Define the parallel k-point chain ( 61 belong only to this chain )
    1226       131040 :  chain(31:50,6)=.false.
    1227        62400 :  chain(51:59,6)=.false.
    1228        62400 :  chain(62:70,6)=.false.
    1229       505440 :  chain(71:marrays,6)=.false.
    1230              : 
    1231              : !Define the newvtr chain ( 41:50 belong only to this chain)
    1232        68640 :  chain(31:40,7)=.false.
    1233       630240 :  chain(51:marrays,7)=.false.
    1234              : 
    1235              : !Define the suscep chain ( 101:110 belong only to this chain)
    1236       755040 :  chain(31:marrays,8)=.false.
    1237         6240 :  chain(60    ,8)=.true.
    1238        68640 :  chain(101:110,8)=.true.
    1239              : 
    1240              : !Define the dielmt chain ( 111:120 belong only to this chain)
    1241       755040 :  chain(31:marrays,9)=.false.
    1242         6240 :  chain(60    ,9)=.true.
    1243        68640 :  chain(111:120,9)=.true.
    1244              : 
    1245              : !Define the tddft chain ( 121:130 belong only to this chain)
    1246       755040 :  chain(31:marrays,10)=.false.
    1247         6240 :  chain(60    ,10)=.true.
    1248        68640 :  chain(121:130,10)=.true.
    1249              : 
    1250              : !The memory needed for each chain has been computed
    1251              : !-------------------------------------------------------------------------
    1252              : !Still need some auxiliary data : estimate the disk space
    1253              : !or the maximum segment size.
    1254              : 
    1255              : !XG030513 : MPIWF need to multiply mbdiskwf by the number of processors
    1256              : !in the WF group. For the time being, nprocwf=1
    1257         6240 :  nprocwf=mpi_enreg%nproc_fft
    1258              : 
    1259        67273 :  mbdiskwf=(8*two*mpw*nprocwf*sum(nband(1:nkpt*nsppol)))/1024._dp**2 + 0.002_dp
    1260         6240 :  mbdiskpd=(8*nfftf*nsppol)/1024._dp**2 + 0.002_dp
    1261              : 
    1262              : !Determine the largest array out of cg (cg_disk), f_fftgr (f_fftgr_disk), or pawfgrtab%gylm
    1263         6240 :  if(mkmem==0)then
    1264            0 :    mbcg=(8*2*mpw*mband)/1024._dp**2 + 0.002_dp
    1265              :  else
    1266         6240 :    mbcg=(8*2*mpw*mband*mkmem*nsppol)/1024._dp**2 + 0.002_dp
    1267              :  end if
    1268         6240 :  if(mffmem==0)then
    1269            0 :    if (pawmixdg==1) then
    1270            0 :      mbf_fftgr=(8*nfftf*n_fftgr)/1024._dp**2 + 0.002_dp
    1271              :    else
    1272            0 :      mbf_fftgr=(8*nfft*n_fftgr)/1024._dp**2 + 0.002_dp
    1273              :    end if
    1274              :  else
    1275         6240 :    if (pawmixdg==1) then
    1276          189 :      mbf_fftgr=(8*nfftf*n_fftgr*nsppol*mffmem)/1024._dp**2 + 0.002_dp
    1277              :    else
    1278         6051 :      mbf_fftgr=(8*nfft*n_fftgr*nsppol*mffmem)/1024._dp**2 + 0.002_dp
    1279              :    end if
    1280              :  end if
    1281         6240 :  if(usepaw==1)then
    1282         1068 :    mbgylm=0
    1283         2646 :    do ii=1,ntypat                                        ! pawfgrtab
    1284         1578 :      jj=(1+int(nfftf*four_pi/(three*ucvol)*rshp(ii)**3))
    1285              :      mbgylm=mbgylm+my_nattyp(ii)*jj &
    1286              : &     *( l_size(ii)**2*pawstgylm &                              !   %gylm   (if pawstgylm=1)
    1287              : &    +3*max((optforces+1)/2,optstress)*l_size(ii)**2*pawstgylm& !   %gylmgr (if pawstgylm=1)
    1288              : &    +3*optstress*pawstgylm&                                    !   %rfgd   (if pawstgylm=1)
    1289         2646 : &    +3*(1-pawstgylm) )                                         !   %rfgd   (if pawstgylm=0)
    1290              :    end do
    1291         1068 :    mbgylm=8*mbgylm/1024._dp**2 + 0.002_dp
    1292              :  else
    1293         5172 :    mbgylm=0
    1294              :  end if
    1295              : 
    1296              : !-------------------------------------------------------------------------
    1297         6240 :  ABI_FREE(my_nattyp)
    1298         6240 :  ABI_FREE(basis_size)
    1299         6240 :  ABI_FREE(l_size)
    1300         6240 :  ABI_FREE(lmn_size)
    1301         6240 :  ABI_FREE(lmn2_size)
    1302         6240 :  ABI_FREE(mesh_size)
    1303         6240 :  ABI_FREE(pawver)
    1304         6240 :  ABI_FREE(shape_type)
    1305         6240 :  ABI_FREE(rshp)
    1306              : 
    1307              : !---------------------------------------------------------------------
    1308              : !Now, analyze the data
    1309              : 
    1310              :  call memana(cadd,cfft,cfftf,chain,cmpw,dttyp,iout,iprcel,iscf,&
    1311              : & marrays,mbcg,mbdiskpd,mbdiskwf,mbf_fftgr,mbgylm,mffmem,&
    1312         6240 : & mpw,natom,nchain,nfft,nfftf,occopt,option,prtvol)
    1313              : 
    1314         6240 : end subroutine memory
    1315              : !!***
    1316              : 
    1317              : !!****f* m_memeval/memana
    1318              : !! NAME
    1319              : !! memana
    1320              : !!
    1321              : !! FUNCTION
    1322              : !! Analysis of the memory and disk space needed for the job,
    1323              : !! thanks to the data computed in the calling routine: for each
    1324              : !! array, the number of blocks of size mpw or nfft bytes, and the
    1325              : !! additional memory occupation;
    1326              : !! the list of arrays that are used for each chain.
    1327              : !!
    1328              : !! According to the value of the option variable,
    1329              : !! the routine will eventually try to allocate this amount of memory,
    1330              : !! and if it fails, estimate the maximum value nfft compatible with
    1331              : !! the available memory.
    1332              : !!
    1333              : !! INPUTS
    1334              : !!  cadd(marrays)= count of bytes needed in addition of cmpw, cfftc and cfft.
    1335              : !!  cfft(marrays) =for each array, count of blocks of size nfft bytes (coarse grid, if PAW)
    1336              : !!  cfftf(marrays)=for each array, count of blocks of size nfft bytes (fine grid, if PAW)
    1337              : !!  chain(marrays,nchain)=logical variable, that informs whether an array
    1338              : !!    belongs to a given chain.
    1339              : !!  cmpw(marrays)=for each array, count of blocks of size mpw bytes.
    1340              : !!  dttyp(marrays)=datatype of the array : 4 for integers, 8 for real(dp)
    1341              : !!  iout=unit number for output of formatted data.
    1342              : !!  iprcel=govern the choice of preconditioner for the SCF cycle
    1343              : !!  iscf=governs the choice of SCF algorithm, or non-SCF calculation.
    1344              : !!  marrays=maximal number of arrays (or group of arrays) to be monitored.
    1345              : !!  mbcg=number of MB needed for the cg array.
    1346              : !!  mbdiskpd=number of MB needed to store a density or potential file on disk
    1347              : !!  mbdiskwf=number of MB needed to store a wavefunction file on disk
    1348              : !!  mbf_fftgr=number of MB needed for the f_fftgr array.
    1349              : !!  mbgylm=number of MB needed for the pawfgrtab%gylm array (paw only)
    1350              : !!  mffmem =governs the number of FFT arrays which are fit in core memory
    1351              : !!  mpw   =maximum number of planewaves in basis sphere (large number)
    1352              : !!  natom =number of atoms in unit cell
    1353              : !!  nchain=number of chains to be used in the estimation of memory.
    1354              : !!  nfft =(effective) number of FFT grid points (for one processor) (coarse grid, if PAW)
    1355              : !!  nfftf=(effective) number of FFT grid points (for one processor) (fine grid, if PAW)
    1356              : !!  occopt=option for occupation numbers. If 3<=occopt<=8, varying occupation
    1357              : !!  option : if 0 , no test of available memory
    1358              : !!           if 1 , the routine tries to allocate the estimated memory, for testing
    1359              : !!                    purposes, and if a failure occurs, the routine stops.
    1360              : !!           if 2 , like 1, but before stopping, the routine will provide
    1361              : !!                    an estimation of the available memory.
    1362              : !!  prtvol=control print volume
    1363              : !!
    1364              : !! OUTPUT
    1365              : !!  (only writing)
    1366              : !!
    1367              : !! SOURCE
    1368              : 
    1369         6958 : subroutine memana(cadd,cfft,cfftf,chain,cmpw,dttyp,iout,iprcel,iscf,&
    1370              : & marrays,mbcg,mbdiskpd,mbdiskwf,mbf_fftgr,mbgylm,mffmem,&
    1371              : & mpw,natom,nchain,nfft,nfftf,occopt,option,prtvol)
    1372              : 
    1373              : !Arguments ------------------------------------
    1374              : !scalars
    1375              :  integer,intent(in) :: iout,iprcel,iscf,marrays,mffmem,mpw,natom,nchain
    1376              :  integer,intent(in) :: nfft,nfftf,occopt,option,prtvol
    1377              :  real(dp),intent(in) :: mbcg,mbdiskpd,mbdiskwf,mbf_fftgr,mbgylm
    1378              : !arrays
    1379              :  integer,intent(in) :: dttyp(marrays)
    1380              :  logical,intent(in) :: chain(marrays,nchain)
    1381              :  real(dp),intent(in) :: cadd(marrays),cfft(marrays),cfftf(marrays),cmpw(marrays)
    1382              : 
    1383              : !Local variables-------------------------------
    1384              : !scalars
    1385              :  integer :: biggest,ichain,ier,ier1,ier2,ier3,ier4,ier5,ier6,ier7,ier8,ii
    1386              : !integer :: jj,kk
    1387              :  integer :: mu,nmbytes,nquarter_mbytes,quit
    1388              :  real(dp) :: mbbigarr,mbbiggest
    1389              :  character(len=500) :: msg
    1390              : !arrays
    1391         6958 :  real(dp),allocatable :: bigarray(:,:),bigarray1(:,:),bigarray2(:,:)
    1392         6958 :  real(dp),allocatable :: bigarray3(:,:),bigarray4(:,:),bigarray5(:,:)
    1393         6958 :  real(dp),allocatable :: bigarray6(:,:),bigarray7(:,:),bigarray8(:,:)
    1394         6958 :  real(dp),allocatable :: cdpadd(:),cdpfft(:),cdpfftf(:),cdpmpw(:)
    1395         6958 :  real(dp),allocatable :: cintfft(:),cintfftf(:),cintmpw(:),cintadd(:)
    1396         6958 :  real(dp),allocatable :: mbdpadd(:),mbdpfft(:),mbdpfftf(:)
    1397         6958 :  real(dp),allocatable :: mbdpmpw(:),mbintadd(:),mbintfft(:),mbintfftf(:)
    1398         6958 :  real(dp),allocatable :: mbintmpw(:),mbother(:),mbtot(:)
    1399              : 
    1400              : ! **************************************************************************
    1401              : 
    1402              : !write(std_out,*)' memana : nchain=',nchain
    1403              : 
    1404        20874 :  ABI_MALLOC(cdpfftf,(nchain))
    1405        13916 :  ABI_MALLOC(cdpfft,(nchain))
    1406        13916 :  ABI_MALLOC(cdpmpw,(nchain))
    1407        13916 :  ABI_MALLOC(cintfftf,(nchain))
    1408        13916 :  ABI_MALLOC(cintfft,(nchain))
    1409        13916 :  ABI_MALLOC(cintmpw,(nchain))
    1410        13916 :  ABI_MALLOC(cdpadd,(nchain))
    1411        13916 :  ABI_MALLOC(cintadd,(nchain))
    1412        13916 :  ABI_MALLOC(mbdpadd,(nchain))
    1413        13916 :  ABI_MALLOC(mbdpfftf,(nchain))
    1414        13916 :  ABI_MALLOC(mbdpfft,(nchain))
    1415        13916 :  ABI_MALLOC(mbdpmpw,(nchain))
    1416        13916 :  ABI_MALLOC(mbintadd,(nchain))
    1417        13916 :  ABI_MALLOC(mbintfftf,(nchain))
    1418        13916 :  ABI_MALLOC(mbintfft,(nchain))
    1419        13916 :  ABI_MALLOC(mbintmpw,(nchain))
    1420        13916 :  ABI_MALLOC(mbother,(nchain))
    1421        13916 :  ABI_MALLOC(mbtot,(nchain))
    1422              : 
    1423         6958 :  biggest=0
    1424         6958 :  mbbiggest=0.0_dp
    1425              : 
    1426              : !For each chain, compute the number of bytes
    1427        70794 :  do ichain=1,nchain
    1428              : 
    1429              : !  First, the number of integer or real(dp), fft, mpw or add blocks
    1430      9639236 :    cdpmpw(ichain) =sum(cmpw(:),MASK=(dttyp(:)==8).and.chain(:,ichain))
    1431      9639236 :    cintmpw(ichain)=sum(cmpw(:),MASK=(dttyp(:)==4).and.chain(:,ichain))
    1432      9639236 :    cdpfftf(ichain) =sum(cfftf(:),MASK=(dttyp(:)==8).and.chain(:,ichain))
    1433      9639236 :    cintfftf(ichain)=sum(cfftf(:),MASK=(dttyp(:)==4).and.chain(:,ichain))
    1434      9639236 :    cdpfft(ichain) =sum(cfft(:),MASK=(dttyp(:)==8).and.chain(:,ichain))
    1435      9639236 :    cintfft(ichain)=sum(cfft(:),MASK=(dttyp(:)==4).and.chain(:,ichain))
    1436      9639236 :    cdpadd(ichain) =sum(cadd(:),MASK=(dttyp(:)==8).and.chain(:,ichain))
    1437      9639236 :    cintadd(ichain)=sum(cadd(:),MASK=(dttyp(:)==4).and.chain(:,ichain))
    1438              : 
    1439              : !  Compute the corresponding number of Mbytes
    1440        63836 :    mbdpmpw(ichain) =8*cdpmpw(ichain) *dble(mpw) /1024._dp**2
    1441        63836 :    mbintmpw(ichain)=4*cintmpw(ichain)*dble(mpw) /1024._dp**2
    1442        63836 :    mbdpfftf(ichain) =8*cdpfftf(ichain) *dble(nfftf)/1024._dp**2
    1443        63836 :    mbintfftf(ichain)=4*cintfftf(ichain)*dble(nfftf)/1024._dp**2
    1444        63836 :    mbdpfft(ichain) =8*cdpfft(ichain) *dble(nfft)/1024._dp**2
    1445        63836 :    mbintfft(ichain)=4*cintfft(ichain)*dble(nfft)/1024._dp**2
    1446        63836 :    mbdpadd(ichain) =8*cdpadd(ichain)              /1024._dp**2
    1447        63836 :    mbintadd(ichain)=4*cintadd(ichain)             /1024._dp**2
    1448        63836 :    mbother(ichain) =dble(231+6*natom)/1024._dp
    1449        63836 :    if(3<=occopt .and. occopt<=8)mbother(ichain)=dble(991+natom)/1024._dp
    1450              : 
    1451              : !  Compute the total number of Mbytes
    1452              :    mbtot(ichain)=mbdpmpw(ichain)+mbintmpw(ichain)&
    1453              : &   +mbdpfftf(ichain)+mbintfftf(ichain)&
    1454              : &   +mbdpfft(ichain)+mbintfft(ichain)&
    1455        63836 : &   +mbdpadd(ichain)+mbintadd(ichain)+mbother(ichain)
    1456              : 
    1457              : !  Select the biggest chain
    1458        70794 :    if(mbtot(ichain)>mbbiggest)then
    1459        11035 :      mbbiggest=mbtot(ichain)
    1460        11035 :      biggest=ichain
    1461              :    end if
    1462              :  end do
    1463              : !When iprcel<20, the biggest chains cannot be number 8 or 9 ...
    1464         6958 :  if(modulo(iprcel,100)<20 .and. (biggest==8 .or. biggest==9))then
    1465            0 :    write(msg,'(a,a,a,a,i3,a,a,a)') ch10,&
    1466            0 : &   ' memana: BUG -',ch10,&
    1467            0 : &   '  The biggest chain is number',biggest,' while iprcel==20.',ch10,&
    1468            0 : &   '  This is not allowed.'
    1469            0 :    call wrtout(std_out,msg)
    1470              :  end if
    1471              : 
    1472              :  write(msg, '(a,f11.3,a)' ) &
    1473         6958 : & 'P This job should need less than                 ',&
    1474        13916 : & mbbiggest+tol10,' Mbytes of memory. '
    1475         6958 :  call wrtout(std_out,msg)
    1476         6958 :  call wrtout(iout,msg)
    1477              : 
    1478         6958 :  if(prtvol>=10)then
    1479          245 :    if(biggest==1)write(msg,'(a)')'P Max. in main chain + fourwf.f '
    1480          245 :    if(biggest==2)write(msg,'(a)')'P Max. in main chain + nonlop.f + opernl.f '
    1481          245 :    if(biggest==3)write(msg,'(a)')'P Max. in XC chain '
    1482          245 :    if(biggest==4)write(msg,'(a)')'P Max. in mkrho chain '
    1483          245 :    if(biggest==5)write(msg,'(a)')'P Max. in fourdp chain '
    1484          245 :    if(biggest==6)write(msg,'(a)')'P Max. in parallel k-point chain '
    1485          245 :    if(biggest==7)write(msg,'(a)')'P Max. in newvtr chain '
    1486          245 :    if(biggest==8)write(msg,'(a)')'P Max. in suscep chain '
    1487          245 :    if(biggest==9)write(msg,'(a)')'P Max. in dielmt chain '
    1488          245 :    if(biggest==10)write(msg,'(a)')'P Max. in tddft chain '
    1489          245 :    call wrtout(iout,msg)
    1490              : 
    1491              :    write(msg, '(a,i13,a,f11.3,a)' )&
    1492          245 : &   'P',nint(cintmpw(biggest)),' blocks of mpw  integer numbers, for',&
    1493          490 : &   mbintmpw(biggest)+tol10,' Mbytes. '
    1494          245 :    call wrtout(iout,msg)
    1495              :    write(msg, '(a,i13,a,f11.3,a)' )&
    1496          245 : &   'P',nint(cdpmpw(biggest)),' blocks of mpw  real(dp)  numbers, for',&
    1497          490 : &   mbdpmpw(biggest)+tol10,' Mbytes. '
    1498          245 :    call wrtout(iout,msg)
    1499          245 :    if (nfft==nfftf) then
    1500          225 :      if(mbintfft(biggest)+mbintfftf(biggest)>0.001)then
    1501              :        write(msg, '(a,i13,a,f11.3,a)' )&
    1502          149 : &       'P',nint(cintfft(biggest)+cintfftf(biggest)),' blocks of nfft integer numbers, for',&
    1503          298 : &       mbintfft(biggest)+mbintfftf(biggest)+tol10,' Mbytes. '
    1504          149 :        call wrtout(iout,msg)
    1505              :      end if
    1506              :      write(msg, '(a,i13,a,f11.3,a)' )&
    1507          225 : &     'P',nint(cdpfft(biggest)+cdpfftf(biggest)),' blocks of nfft real(dp)  numbers, for',&
    1508          450 : &     mbdpfft(biggest)+mbdpfftf(biggest)+tol10,' Mbytes. '
    1509          225 :      call wrtout(iout,msg)
    1510              :    else
    1511           20 :      if(mbintfftf(biggest)>0.001)then
    1512              :        write(msg, '(a,i13,a,f11.3,a)' )&
    1513            0 : &       'P',nint(cintfftf(biggest)),' blocks of nfft (fine grid) integer numbers, for',&
    1514            0 : &       mbintfftf(biggest)+tol10,' Mbytes. '
    1515            0 :        call wrtout(iout,msg)
    1516              :      end if
    1517              :      write(msg, '(a,i13,a,f11.3,a)' )&
    1518           20 : &     'P',nint(cdpfftf(biggest)),' blocks of nfft (fine grid) real(dp)  numbers, for',&
    1519           40 : &     mbdpfftf(biggest)+tol10,' Mbytes. '
    1520           20 :      call wrtout(iout,msg)
    1521           20 :      if(mbintfft(biggest)>0.001)then
    1522              :        write(msg, '(a,i13,a,f11.3,a)' )&
    1523           20 : &       'P',nint(cintfft(biggest)),' blocks of nfft (coarse grid) integer numbers, for',&
    1524           40 : &       mbintfft(biggest)+tol10,' Mbytes. '
    1525           20 :        call wrtout(iout,msg)
    1526              :      end if
    1527              :      write(msg, '(a,i13,a,f11.3,a)' )&
    1528           20 : &     'P',nint(cdpfft(biggest)),' blocks of nfft (coarse grid) real(dp)  numbers, for',&
    1529           40 : &     mbdpfft(biggest)+tol10,' Mbytes. '
    1530           20 :      call wrtout(iout,msg)
    1531              :    end if
    1532          245 :    if(mbintadd(biggest)>0.001)then
    1533          107 :      write(msg, '(a,13x,a,f11.3,a)' )'P',' Additional     integer numbers, for',mbintadd(biggest)+tol10,' Mbytes. '
    1534          107 :      call wrtout(iout,msg)
    1535              :    end if
    1536          245 :    write(msg, '(a,13x,a,f11.3,a)' )'P',' Additional     real(dp)  numbers, for',mbdpadd(biggest)+tol10,' Mbytes. '
    1537          245 :    call wrtout(iout,msg)
    1538          245 :    write(msg, '(a,13x,a,f11.3,a)' )'P',' With residue estimated to be       ',mbother(biggest)+tol10,' Mbytes. '
    1539          245 :    call wrtout(iout,msg)
    1540          245 :    write(msg, '(a)' )'P'
    1541          245 :    call wrtout(iout,msg)
    1542          245 :    write(msg, '(a)' )'P Comparison of the memory needs of different chains'
    1543          245 :    call wrtout(iout,msg)
    1544              : 
    1545          245 :    write(msg, '(a,f11.3,a)' )'P Main chain + fourwf.f           ',mbtot(1)+tol10,' Mbytes. '
    1546          245 :    call wrtout(iout,msg)
    1547          245 :    write(msg, '(a,f11.3,a)' )'P Main chain + nonlop.f + opernl.f',mbtot(2)+tol10,' Mbytes. '
    1548          245 :    call wrtout(iout,msg)
    1549              : 
    1550              : !  The next chains are not defined in the RF case.
    1551          245 :    if(nchain>2)then
    1552          194 :      write(msg, '(a,f11.3,a)' )'P XC chain                        ',mbtot(3)+tol10,' Mbytes. '
    1553          194 :      call wrtout(iout,msg)
    1554              :      write(msg, '(a,f11.3,a)' )&
    1555          194 : &     'P mkrho chain                     ',mbtot(4)+tol10,' Mbytes. '
    1556          194 :      call wrtout(iout,msg)
    1557              :      write(msg, '(a,f11.3,a)' )&
    1558          194 : &     'P fourdp chain                    ',mbtot(5)+tol10,' Mbytes. '
    1559          194 :      call wrtout(iout,msg)
    1560              :      if(xmpi_paral==1)then
    1561              :        write(msg, '(a,f11.3,a)' )&
    1562          194 : &       '- parallel k-point chain          ',mbtot(6)+tol10,' Mbytes. '
    1563          194 :        call wrtout(iout,msg)
    1564              :      end if
    1565              :      write(msg, '(a,f11.3,a)' )&
    1566          194 : &     'P newvtr chain                    ',mbtot(7)+tol10,' Mbytes. '
    1567          194 :      call wrtout(iout,msg)
    1568          194 :      if(modulo(iprcel,100)>=20.and.modulo(iprcel,100)<70)then
    1569              :        write(msg, '(a,f11.3,a)' )&
    1570           16 : &       'P suscep chain                    ',mbtot(8)+tol10,' Mbytes. '
    1571           16 :        call wrtout(iout,msg)
    1572              :        write(msg, '(a,f11.3,a)' )&
    1573           16 : &       'P dielmt chain                    ',mbtot(9)+tol10,' Mbytes. '
    1574           16 :        call wrtout(iout,msg)
    1575              :      end if
    1576          194 :      if(iscf==-1)then
    1577              :        write(msg, '(a,f11.3,a)' )&
    1578            0 : &       'P tddft  chain                    ',mbtot(10)+tol10,' Mbytes. '
    1579              :      end if
    1580              :    end if ! nchain>2
    1581              : 
    1582              :  end if
    1583              : 
    1584              : !--------------------------------------------------------------------
    1585              : 
    1586         6958 :  write(msg, '(a)' ) '  Rough estimation (10% accuracy) of disk space for files :'
    1587         6958 :  call wrtout(iout,msg)
    1588         6958 :  call wrtout(std_out,msg)
    1589              : 
    1590              :  write(msg, '(a,f11.3,a,a,f11.3,a)' ) &
    1591         6958 : & '_ WF disk file :',mbdiskwf+tol10,' Mbytes ;',&
    1592        13916 : & ' DEN or POT disk file :',mbdiskpd+tol10,' Mbytes.'
    1593         6958 :  call wrtout(iout,msg)
    1594         6958 :  call wrtout(std_out,msg)
    1595              : 
    1596         6958 :  if(mffmem==0 .and. iscf>0)then
    1597              :    if(iscf==1)then
    1598              :      write(msg, '(a,a,a)' )&
    1599            0 : &     '  mffmem==0, iscf==1 => use of 1 FFT temporary disk file,',ch10,&
    1600            0 : &     '                       5 times bigger than a DEN file.'
    1601              :    else if(iscf==2.or.iscf==12)then
    1602              :      write(msg, '(a,a,a)' )&
    1603            0 : &     '  mffmem==0, iscf==2 => use of 1 FFT temporary disk file,',ch10,&
    1604            0 : &     '                       3 times bigger than a DEN file.'
    1605              :    else if(iscf==3.or.iscf==13)then
    1606              :      write(msg, '(a,a,a)' )&
    1607            0 : &     '  mffmem==0, iscf==3 => use of 1 FFT temporary disk file,',ch10,&
    1608            0 : &     '                       4 times bigger than a DEN file.'
    1609              :    else if(iscf==4.or.iscf==14)then
    1610              :      write(msg, '(a,a,a)' )&
    1611            0 : &     '  mffmem==0, iscf==4 => use of 1 FFT temporary disk file,',ch10,&
    1612            0 : &     '                       6 times bigger than a DEN file.'
    1613              :    else if(iscf==5)then
    1614              :      write(msg, '(a,a,a)' )&
    1615            0 : &     '  mffmem==0, iscf==5 => use of 1 FFT temporary disk file,',ch10,&
    1616            0 : &     '                       10 times bigger than a DEN file.'
    1617              :    else if(iscf==6)then
    1618              :      write(msg, '(a,a,a)' )&
    1619            0 : &     '  mffmem==0, iscf==6 => use of 1 FFT temporary disk file,',ch10,&
    1620            0 : &     '                       10 times bigger than a DEN file.'
    1621              :    else if(iscf==7.or.iscf==17)then
    1622              :      write(msg, '(a,a,a)' )&
    1623            0 : &     '  mffmem==0, iscf==7 => use of 1 FFT temporary disk file,',ch10,&
    1624            0 : &     '                       (2+2*npulayit) times bigger than a DEN file.'
    1625              :    end if
    1626            0 :    call wrtout(iout,msg)
    1627            0 :    call wrtout(std_out,msg)
    1628              :  end if
    1629              : 
    1630              : !Temporary msg - estimation of PAW specific data has to be done...
    1631              : !Have to add the usepaw argument to use this.
    1632              : !if (usepaw==1) then
    1633              : !write(msg,'(5a)') '  WARNING: You are using PAW formalism;',ch10,&
    1634              : !&       '           Above estimations do not take PAW',ch10,&
    1635              : !&       '           specific data into account !'
    1636              : !call wrtout(iout,msg)
    1637              : !call wrtout(std_out,msg)
    1638              : !end if
    1639              : 
    1640       563598 :  write(msg,'(80a,a)') ('=',mu=1,80),ch10
    1641         6958 :  call wrtout(iout,msg)
    1642         6958 :  call wrtout(std_out,msg)
    1643              : 
    1644              : !--------------------------------------------------------------------
    1645              : !Here, each processor must test its memory, so use
    1646              : !the PERS mode for error msgs, followed by synchronisation
    1647              : 
    1648         6958 :  mbbigarr=max(mbf_fftgr,mbcg,mbgylm)
    1649         6958 :  if(mbbigarr==mbcg) then
    1650         1317 :    write(msg, '(a,f12.4,a)' ) ' Biggest array : cg(disk), with',mbcg+tol10,' MBytes.'
    1651         5641 :  else if (mbbigarr==mbf_fftgr) then
    1652         5164 :    write(msg, '(a,f12.4,a)' ) ' Biggest array : f_fftgr(disk), with',mbf_fftgr+tol10,' MBytes.'
    1653          477 :  else if (mbbigarr==mbgylm)then
    1654          477 :    write(msg, '(a,f12.4,a)' ) ' Biggest array : pawfgrtab%gylm(gr), with',mbgylm+tol10,' MBytes.'
    1655              :  end if
    1656         6958 :  call wrtout(std_out,msg)
    1657              : 
    1658              : !if (mpi_enreg%my_nimage>1) then
    1659              : !write(msg, '(a,f12.4,a)' ) &
    1660              : !&   ' These estimations take the distribution over replicas (images) of the cell into account.'
    1661              : !call wrtout(std_out,msg)
    1662              : !end if
    1663              : 
    1664         6958 :  quit=0
    1665              : 
    1666         6958 :  if(option>=1)then
    1667              : 
    1668              : !  Test the ability to allocate the biggest array
    1669         6467 :    nquarter_mbytes=4.0_dp*mbbigarr+1.0_dp
    1670        19401 :    ABI_STAT_MALLOC(bigarray,(32*1024,nquarter_mbytes), ier)
    1671         6467 :    if(ier/=0)then
    1672              :      write(msg,'(a,f11.3,a,a,a,a,a,a,a)')&
    1673            0 : &     'Test failed to allocate an array of',mbbigarr,' Mbytes',ch10,&
    1674            0 : &     'It is not worth to continue ',ch10,&
    1675            0 : &     'Action: modify input variable to fit the available memory,',ch10,&
    1676            0 : &     'increase limit on maximal array size or set mem_test to 0 to disable this test.'
    1677            0 :      call wrtout(std_out,msg,'PERS')
    1678            0 :      if(option==1)then
    1679            0 :        ABI_ERROR_CLASS(msg, "MemanaError")
    1680              :      else
    1681            0 :        ABI_WARNING(msg)
    1682              :        quit=1
    1683              :      end if
    1684              :    end if
    1685              :    if(quit==0)then
    1686         6467 :      write(msg,'(a,f11.3,a)')' memana : allocated an array of',mbbigarr+tol10,' Mbytes, for testing purposes. '
    1687         6467 :      call wrtout(std_out,msg)
    1688              :    end if
    1689              : 
    1690         6467 :    ABI_SFREE(bigarray)
    1691              : 
    1692              : !  Test the ability to allocate the needed total memory : use 8 segments,
    1693              : !  hoping that the maximal segment size is not so much smaller than the
    1694              : !  total memory
    1695         6467 :    nquarter_mbytes=0.5_dp*mbbiggest+1.0_dp
    1696        19401 :    ABI_STAT_MALLOC(bigarray1,(32*1024,nquarter_mbytes), ier1)
    1697        12934 :    ABI_STAT_MALLOC(bigarray2,(32*1024,nquarter_mbytes), ier2)
    1698        12934 :    ABI_STAT_MALLOC(bigarray3,(32*1024,nquarter_mbytes), ier3)
    1699        12934 :    ABI_STAT_MALLOC(bigarray4,(32*1024,nquarter_mbytes), ier4)
    1700        12934 :    ABI_STAT_MALLOC(bigarray5,(32*1024,nquarter_mbytes), ier5)
    1701        12934 :    ABI_STAT_MALLOC(bigarray6,(32*1024,nquarter_mbytes), ier6)
    1702        12934 :    ABI_STAT_MALLOC(bigarray7,(32*1024,nquarter_mbytes), ier7)
    1703        12934 :    ABI_STAT_MALLOC(bigarray8,(32*1024,nquarter_mbytes), ier8)
    1704              : 
    1705         6467 :    if(ier1/=0 .or. ier2/=0 .or. ier3/=0 .or. ier4/=0 .or. ier5/=0 .or. ier6/=0 .or. ier7/=0 .or. ier8/=0) then
    1706              :      write(msg,'(a,f11.3,a,a,a,a,a,a,a)')&
    1707            0 : &     'Test failed to allocate ',mbbiggest,' Mbytes',ch10,&
    1708            0 : &     'It is not worth to continue ',ch10,&
    1709            0 : &     'Action: modify input variables or submission parameters to fit the available memory,',ch10,&
    1710            0 : &     'increase limit on available memory or set mem_test to 0 to disable this test.'
    1711            0 :      if(option==1)then
    1712            0 :        ABI_ERROR_CLASS(msg, "MemanaError")
    1713              :      else
    1714            0 :        ABI_WARNING(msg)
    1715              :        quit=1
    1716              :      end if
    1717              :    end if
    1718              : 
    1719         6467 :    if(quit==0)then
    1720              :      write(msg,'(a,f11.3,a,a,a)')&
    1721         6467 : &     ' memana: allocated ',mbbiggest,'Mbytes, for testing purposes. ',ch10,&
    1722        12934 : &     ' The job will continue.'
    1723         6467 :      call wrtout(std_out,msg)
    1724              :    end if
    1725         6467 :    ABI_SFREE(bigarray1)
    1726         6467 :    ABI_SFREE(bigarray2)
    1727         6467 :    ABI_SFREE(bigarray3)
    1728         6467 :    ABI_SFREE(bigarray4)
    1729         6467 :    ABI_SFREE(bigarray5)
    1730         6467 :    ABI_SFREE(bigarray6)
    1731         6467 :    ABI_SFREE(bigarray7)
    1732         6467 :    ABI_SFREE(bigarray8)
    1733              : 
    1734              :  end if
    1735              : 
    1736              : !--------------------------------------------------------------------
    1737              : 
    1738         6958 :  if(option==2 .and. quit==1 )then
    1739              : 
    1740              : !  Estimation of the available memory
    1741              : !
    1742              : !  A quarter of Mbyte is 256*1024/8 real(dp) numbers,
    1743              : !  that is 32*1024 dp numbers.
    1744              : !  One begins with the allocation of 4 Mbytes. If successful,
    1745              : !  one increases that number, until the allocation is not successfull
    1746              : !  any more. Unfortunately, on a P6 with the pghpf compiler, the
    1747              : !  allocate instruction generate a core dump, instead of returning
    1748              : !  an error code, so that this part of code has been made optional.
    1749              : 
    1750            0 :    nquarter_mbytes=16
    1751            0 :    nmbytes=nquarter_mbytes/4.0_dp
    1752              : 
    1753              : !  With an increase ratio of 1.25_dp (see below), ii=5 leads to 9 MB,
    1754              : !  ii=10 leads to 28 MB, ii=15 leads to 85 MB, ii=18 leads to 165 MB,
    1755              : !  ii=30 is over 2 GB
    1756            0 :    do ii=1,30
    1757            0 :      ABI_STAT_MALLOC(bigarray,(32*1024,nquarter_mbytes), ier)
    1758            0 :      if(ier/=0)then
    1759            0 :        write(msg,'(a,i0,a)')' memana : failed to allocate ',nmbytes,' Mbytes'
    1760            0 :        call wrtout(std_out,msg,'PERS')
    1761            0 :        exit
    1762              :      end if
    1763            0 :      write(msg,'(a,i0,a)')' memana : succeeded to allocate ',nmbytes,' Mbytes'
    1764            0 :      call wrtout(std_out,msg,'PERS')
    1765              : !    Here really test the space
    1766              : !    do kk=1,nquarter_mbytes
    1767              : !    do jj=1,32*1024,37
    1768              : !    bigarray(jj,kk)=0.0_dp
    1769              : !    end do
    1770              : !    write(std_out,*)' memana : wrote ',kk,' quarter of mbytes'
    1771              : !    end do
    1772            0 :      ABI_FREE(bigarray)
    1773            0 :      nquarter_mbytes=dble(nquarter_mbytes)*1.25_dp
    1774            0 :      nmbytes=nquarter_mbytes/4.0_dp
    1775              :    end do
    1776            0 :    ABI_SFREE(bigarray)
    1777              : 
    1778            0 :    ABI_ERROR_CLASS("in memana with option==2 .and. quit==1", "MemanaError")
    1779              :  end if !  End the test of the available memory
    1780              : 
    1781              : !--------------------------------------------------------------------
    1782              : 
    1783         6958 :  ABI_FREE(cdpfftf)
    1784         6958 :  ABI_FREE(cdpfft)
    1785         6958 :  ABI_FREE(cdpmpw)
    1786         6958 :  ABI_FREE(cintfftf)
    1787         6958 :  ABI_FREE(cintfft)
    1788         6958 :  ABI_FREE(cintmpw)
    1789         6958 :  ABI_FREE(cdpadd)
    1790         6958 :  ABI_FREE(cintadd)
    1791         6958 :  ABI_FREE(mbdpadd)
    1792         6958 :  ABI_FREE(mbdpfftf)
    1793         6958 :  ABI_FREE(mbdpfft)
    1794         6958 :  ABI_FREE(mbdpmpw)
    1795         6958 :  ABI_FREE(mbintadd)
    1796         6958 :  ABI_FREE(mbintfftf)
    1797         6958 :  ABI_FREE(mbintfft)
    1798         6958 :  ABI_FREE(mbintmpw)
    1799         6958 :  ABI_FREE(mbother)
    1800         6958 :  ABI_FREE(mbtot)
    1801              : 
    1802         6958 : end subroutine memana
    1803              : !!***
    1804              : 
    1805              : !!****f* m_memeval/memorf
    1806              : !! NAME
    1807              : !! memorf
    1808              : !!
    1809              : !! FUNCTION
    1810              : !! Estimation of the memory needed for a response-function job.
    1811              : !! According to the value of the option variable,
    1812              : !! might also try to allocate this amount of memory, and if it fails,
    1813              : !! might estimate the available memory.
    1814              : !!
    1815              : !! INPUTS
    1816              : !!  cplex=1 or 2, indicate whether the den and pot functions are real or complex
    1817              : !!  getcell=if non-zero, the values of acell and rprim are taken from
    1818              : !!   the output of another dataset
    1819              : !!  idtset=number of the current dataset
    1820              : !!  intxc=control xc quadrature
    1821              : !!  iout=unit number for output of formatted data.
    1822              : !!  iprcel=govern the choice of preconditioner for the SCF cycle
    1823              : !!  iscf=governs the choice of SCF algorithm, or non-SCF calculation.
    1824              : !!  jdtset=index of the current dataset
    1825              : !!  lmnmax=max. number of (l,m,n) components over all type of psps
    1826              : !!  lnmax =max. number of (l,n)   components over all type of psps
    1827              : !!  mband =maximum number of bands
    1828              : !!  mffmem =governs the number of FFT arrays which are fit in core memory
    1829              : !!  mgfft =maximum single fft dimension
    1830              : !!  mkmems=number of k points which can fit in memory; set to 0 if use disk
    1831              : !!    the three values correspond to mkmem, mkqmem and mk1mem
    1832              : !!  mpi_enreg=information about MPI parallelization
    1833              : !!  mpssang is 1+maximum angular momentum for nonlocal pseudopotential
    1834              : !!  mpssoang is 1+maximum (spin*angular momentum) for nonlocal pseudopotential
    1835              : !!  mpw   =maximum number of planewaves in basis sphere (large number)
    1836              : !!  mqgrid=maximum dimension of grid of q values for psp representations
    1837              : !!  natom =number of atoms in unit cell
    1838              : !!  nband(nkpt*nsppol)=number of bands at each k point, for each polarization
    1839              : !!  nfft  =(effective) number of FFT grid points (for one processor)
    1840              : !!  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
    1841              : !!  nkpt  =number of k points
    1842              : !!  nloalg(3)=governs the choice of the algorithm for non-local operator.
    1843              : !!  nspden=number of spin-density components
    1844              : !!  nspinor=number of spinorial components of the wavefunctions
    1845              : !!  nsppol=number of channels for spin-polarization (1 or 2)
    1846              : !!  nsym  =number of symmetry elements in space group
    1847              : !!  ntypat=number of types of atoms
    1848              : !!  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
    1849              : !!  occopt=option for occupation numbers. If 3<=occopt<=7, varying occupation
    1850              : !!  optddk=1 if ddk is computed during run
    1851              : !!  optphon=1 if phonons are computed during run
    1852              : !!  option : if 0 , no test of available memory
    1853              : !!           if 1 , the routine tries to allocate the estimated memory, for testing
    1854              : !!                    purposes, and if a failure occurs, the routine stops.
    1855              : !!           if 2 , like 1, but before stopping, the routine will provide
    1856              : !!                    an estimation of the available memory.
    1857              : !!  optstrs=1 if strain perturbation is computing during run
    1858              : !!  prtvol=control print volume
    1859              : !!  useylm=governs the way the nonlocal operator is to be applied:
    1860              : !!         1=using Ylm, 0=using Legendre polynomials
    1861              : !!  gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
    1862              : !!  xclevel= level of the XC functional
    1863              : !!
    1864              : !! OUTPUT
    1865              : !!  (only writing)
    1866              : !!
    1867              : !! NOTES
    1868              : !! for the estimation, it is only taken into account those
    1869              : !! arrays that have some probability of being larger than 1000*8 bytes :
    1870              : !! - All the arrays that have large numbers as one of their dimensions
    1871              : !! (mqgrid, mpw, nfft, ngfft(4)*ngfft(5)*ngfft(6),n1xccc
    1872              : !!                                      or a constant larger than 1000)
    1873              : !! - All the arrays that have a product of two moderately large numbers
    1874              : !! (potential size above 30  : mband, mgfft, mkmems, natom, nkpt, nsym,
    1875              : !!  or a constant larger than 30)
    1876              : !! After this estimation, an amount of (176 + 55 + 6*natom) Kbytes is added
    1877              : !! to take into account the static arrays declared
    1878              : !! in rhotoxc and daughter routines (at maximum 22*1000 dp numbers),
    1879              : !! as well as other arrays like
    1880              : !! character(len=500) :: msg (present in about 100 routines), or the different
    1881              : !! arrays allocated in move.f, brdmin.f, gstate.f (xf array) or pspini.f
    1882              : !! In the case 3<=occopt<=7 this amount is increased by 760 Kbytes
    1883              : !! to take into account the arrays smdfun, occfun, entfun, workfun and xgrid,
    1884              : !! declared in getnel
    1885              : !!
    1886              : !! The current version takes into account only :
    1887              : !! 1) and 2) the "main chain" in its two slightly different versions :
    1888              : !! driver - respfn - dfpt_looppert - dfpt_scfcv - dfpt_vtorho - dfpt_vtowfk -
    1889              : !!     dfpt_cgwf - getghc - fourwf or (nonlop+opernl)
    1890              : !!
    1891              : !! Also, it is assumed that the potentials are non-local, even if they
    1892              : !!     are local ! It would be necessary to update this routine
    1893              : !!     now that the beginning of psp files is read before
    1894              : !!     the present call (XG 980502)
    1895              : !!
    1896              : !! Some BIG approximations, not present in the GS corresponding routine
    1897              : !!  have been done : nsym=nsym1, nkpt=nkpt_rbz, mpw=mpw1 ...
    1898              : !!
    1899              : !! SOURCE
    1900              : 
    1901          718 : subroutine memorf(cplex,n1xccc,getcell,idtset,intxc,iout,iprcel,&
    1902              : & iscf,jdtset,lmnmax,lnmax,mband,mffmem,mgfft,&
    1903              : & mkmems,mpi_enreg,mpsang,mpssoang,mpw,mqgrid,&
    1904          718 : & natom,nband,nfft,ngfft,&
    1905              : & nkpt,nloalg,nspden,nspinor,nsppol,nsym,ntypat,&
    1906              : & occopt,optddk,optphon,option,optstrs,prtvol,useylm,gpu_option,xclevel)
    1907              : 
    1908              : !Arguments ------------------------------------
    1909              : !scalars
    1910              :  integer,intent(in) :: cplex,getcell,idtset,intxc,iout,iprcel,iscf
    1911              :  integer,intent(in) :: jdtset,lmnmax,lnmax,mband,mffmem,mgfft,mpsang
    1912              :  integer,intent(in) :: mpssoang,mpw,mqgrid,n1xccc,natom,nfft,nkpt
    1913              :  integer,intent(in) :: nspden,nspinor,nsppol,nsym,ntypat,occopt
    1914              :  integer,intent(in) :: optddk,option,optphon,optstrs,prtvol,useylm
    1915              :  integer,intent(in) :: gpu_option,xclevel
    1916              :  type(MPI_type),intent(in) :: mpi_enreg
    1917              : !arrays
    1918              :  integer,intent(in) :: mkmems(3),nband(nkpt*nsppol),ngfft(18)
    1919              :  integer,intent(in) :: nloalg(3)
    1920              : 
    1921              : !Local variables-------------------------------
    1922              : !marrays= maximal number of arrays to be monitored (or group of arrays)
    1923              : !cmpw(marrays)=count of blocks of size mpw bytes
    1924              : !cfft(marrays)=number of blocks of size nfft bytes
    1925              : !cadd(marrays)=additional storage needed (in bytes)
    1926              : !dttyp(marrays)=datatype of the array : 4 for integers, 8 for real(dp)
    1927              : !nchain= number of different chains of routines
    1928              : !chain(marrays,nchain)=different chains of routines
    1929              : !scalars
    1930              :  integer,parameter :: marrays=150,nchain=2
    1931              :  integer :: fftalgb,matblk,maxmkmem,mincat,mk1mem,mkmem,mkqmem,mu,n_fftgr
    1932              :  integer :: narr_fourdp,ngrad,nprocwf
    1933              :  integer :: my_natom
    1934              :  integer :: my_mband
    1935              :  real(dp) :: mbcg,mbdiskpd,mbdiskwf,mbf_fftgr,mbgylm
    1936              :  character(len=500) :: msg
    1937              :  character(len=1) :: firstchar
    1938              : !arrays
    1939              :  integer :: dttyp(marrays)
    1940              :  real(dp) :: cadd(marrays),cfft(marrays),cmpw(marrays)
    1941          718 :  real(dp),allocatable :: cfft_dum(:)
    1942              :  logical :: chain(marrays,nchain)
    1943              : 
    1944              : ! **************************************************************************
    1945              : 
    1946          718 :  if(option<0 .or. option>2)then
    1947            0 :    write(msg, '(a,i0,a)')'option= ',option,' while the only allowed values are 0, 1, or 2.'
    1948            0 :    ABI_BUG(msg)
    1949              :  end if
    1950              : 
    1951          718 :  firstchar=' ';if (gpu_option/=ABI_GPU_DISABLED) firstchar='_'
    1952          718 :  cmpw(:)=zero ; cfft(:)=zero ; cadd(:)=zero
    1953          718 :  dttyp(:)=0
    1954              : 
    1955          718 :  call wrtout(std_out,' memorf : analysis of memory needs ')
    1956              : 
    1957          718 :  if(jdtset>=100)then
    1958            0 :    write(msg,'(80a,a,a,i5,a)')('=',mu=1,80),ch10,&
    1959            0 :    ' Values of the parameters that define the memory need for DATASET',jdtset,' (RF).'
    1960          718 :  else if(jdtset/=0)then
    1961        55971 :    write(msg,'(80a,a,a,i3,a)')('=',mu=1,80),ch10,&
    1962        56662 :    ' Values of the parameters that define the memory need for DATASET',jdtset,' (RF).'
    1963              :  else
    1964         2187 :    write(msg,'(80a,a,a,a)')('=',mu=1,80),ch10,&
    1965         2214 :    ' Values of the parameters that define the memory need of the present run',' (RF).'
    1966              :  end if
    1967          718 :  call wrtout(iout,msg)
    1968          718 :  call wrtout(std_out,msg)
    1969              : 
    1970          718 :  mkmem=mkmems(1)
    1971          718 :  mkqmem=mkmems(2)
    1972          718 :  mk1mem=mkmems(3)
    1973          718 :  my_natom=natom;if (mpi_enreg%nproc_atom>1) my_natom=mpi_enreg%my_natom
    1974              :  ! TODO: this does not work yet, as nproc_band is not initialized in the DFPT case (paralkgb is 0)
    1975              :  ! solution:  around src/57_iovars/m_mpi_setup.F90 line 1341, update to give nproc_band a value
    1976          718 :  my_mband = CEILING(dble(mband)/ max(mpi_enreg%nproc_band,1))
    1977              : 
    1978              :  write(msg,'( 4(a,i8),a,4(a,i8) )' ) &
    1979          718 : & '     intxc =',intxc   ,'      iscf =',iscf,&
    1980          718 : & '    lmnmax =',lmnmax  ,'     lnmax =',lnmax,ch10,&
    1981          718 : & '     mgfft =',mgfft,'  mpssoang =',mpssoang,&
    1982         1436 : & '    mqgrid =',mqgrid,'     natom =',natom
    1983          718 :  call wrtout(iout,msg)
    1984          718 :  call wrtout(std_out,msg)
    1985              : 
    1986              :  write(msg,'( 4(a,i8),a,4(a,i8),a,4(a,i8) )' ) &
    1987          718 : & '  nloc_mem =',nloalg(2)*(nloalg(3)+1),'    nspden =',nspden ,&
    1988          718 : & '   nspinor =',nspinor,'    nsppol =',nsppol ,ch10,&
    1989          718 : & '      nsym =',nsym,'    n1xccc =',n1xccc ,&
    1990          718 : & '    ntypat =',ntypat,'    occopt =',occopt ,ch10,&
    1991         1436 : & '   xclevel =',xclevel
    1992          718 :  call wrtout(iout,msg)
    1993          718 :  call wrtout(std_out,msg)
    1994              : 
    1995              :  write(msg,'(4(3(a,i12),a))') &
    1996          718 : & '-    mband =',mband  ,'        mffmem =',mffmem,&
    1997          718 : & '         mkmem =',mkmem  ,ch10,&
    1998          718 : & '-   mkqmem =',mkqmem ,'        mk1mem =',mk1mem,&
    1999          718 : & '           mpw =',mpw  ,ch10,&
    2000         1436 : & '      nfft =',nfft ,'          nkpt =',nkpt
    2001          718 :  call wrtout(iout,msg)
    2002          718 :  call wrtout(std_out,msg)
    2003              : 
    2004              :  write(msg,'(4(3(a,i12),a))') &
    2005          718 : & '- my_mband =',my_mband
    2006          718 :  call wrtout(iout,msg)
    2007          718 :  call wrtout(std_out,msg)
    2008              : 
    2009          718 :  if (my_natom/=natom)then
    2010           36 :    write(msg,'(a,i10)') 'Pmy_natom=',my_natom
    2011           36 :    call wrtout(iout,msg)
    2012           36 :    call wrtout(std_out,msg)
    2013              :  end if
    2014              : 
    2015        58158 :  write(msg,'(80a)') ('=',mu=1,80)
    2016          718 :  call wrtout(iout,msg)
    2017          718 :  call wrtout(std_out,msg)
    2018              : 
    2019          718 :  if(getcell>0 .or. (getcell<0 .and. idtset+getcell>0) )then
    2020            3 :    write(msg,'(a,a,a,a,a,a,i3,a,i3,a,a,a,a,a,a)' )ch10,&
    2021            3 : &   ' memorf : COMMENT -',ch10,&
    2022            3 : &   '  The determination of memory needs at this stage is meaningless,',ch10,&
    2023            3 : &   '  since getcell = ',getcell,' is non-zero, while idtset=',idtset,'.',ch10,&
    2024            3 : &   '  The following numbers are obtained by supposing that acell and rprim',ch10,&
    2025            6 : &   '  are NOT taken from a previous dataset. You cannot rely on them.',ch10
    2026            3 :    call wrtout(iout,msg)
    2027            3 :    call wrtout(std_out,msg)
    2028              :  end if
    2029              : 
    2030          718 :  n_fftgr=1
    2031          718 :  if(iscf==1)            n_fftgr=5
    2032          718 :  if(iscf==2.or.iscf==3) n_fftgr=4
    2033          718 :  if(iscf==5.or.iscf==6) n_fftgr=10
    2034              : 
    2035              : !work1 and work2 in fourdp : take into account approximately fftalgb
    2036          718 :  fftalgb=mod(ngfft(7),100)/10
    2037              :  if(fftalgb==0)narr_fourdp=2*2
    2038          718 :  if(fftalgb==1)narr_fourdp=2
    2039              : 
    2040              :  ngrad=1
    2041          718 :  if(xclevel==2)ngrad=2
    2042              : 
    2043              : !(0)                     in main, driver, and respfn -------------------
    2044              : !indsym (respfn)
    2045          718 :  cadd(1)=4*nsym*natom          ; dttyp(1)=4
    2046              : !rhor,rhog (respfn)
    2047          718 :  cfft(2)=nspden+2              ; dttyp(2)=8
    2048              : !occ (driver), doccde (respfn)
    2049          718 :  cadd(3)=2*mband*nkpt*nsppol   ; dttyp(3)=8
    2050              : !qgrid,vlspl,ffspl (driver)
    2051              :  cadd(4)=mqgrid*(1+2*ntypat*(1+lnmax))   &
    2052          718 : & ; dttyp(4)=8
    2053              : !xccc1d (driver)
    2054          718 :  cadd(5)=n1xccc*6*ntypat       ; dttyp(5)=8
    2055              : !vtrial (respfn)
    2056          718 :  cfft(6)=nspden                ; dttyp(6)=8
    2057              : !kxc (respfn)
    2058          718 :  cfft(7)=2*nspden-1            ; dttyp(7)=8
    2059              : 
    2060              : !(1-2)                   in dfpt_looppert --------------------------------------
    2061              : !ph1d
    2062          718 :  cadd(11)=2*3*(2*mgfft+1)*natom ; dttyp(11)=8
    2063              : !vpsp1
    2064          718 :  cfft(12)=cplex                ; dttyp(12)=8
    2065              : !indsy1  assume that nsym=nsym1
    2066          718 :  cadd(13)=4*nsym*natom         ; dttyp(13)=4
    2067              : !irrzonr1 and phnons1  assume that nsym=nsym1
    2068          718 :  if(nsym/=1)then
    2069          555 :    cfft(14)=(2+(nspden/4))*((nspden/nsppol)-3*nspden/3)     ; dttyp(14)=4
    2070          555 :    cfft(15)=2*((nspden/nsppol)-3*nspden/3)                  ; dttyp(15)=8
    2071              :  end if
    2072              : !doccde_rbz, eigen0, eigenq, occ_rbz, docckqde, occkq, resid
    2073              : !assume than nkpt=nkpt_rbz
    2074          718 :  cadd(16)=7*mband*nkpt*nsppol  ; dttyp(16)=8
    2075              : !kg
    2076          718 :  cmpw(18)=3*mkmem              ; dttyp(18)=4
    2077              : !cg
    2078          718 :  cmpw(19)=2*nspinor*my_mband*mkmem*nsppol  ; dttyp(19)=8
    2079              : !kg1
    2080          718 :  cmpw(21)=3*mk1mem             ; dttyp(21)=4
    2081              : !cgq
    2082          718 :  cmpw(22)=2*nspinor*my_mband*mkqmem*nsppol  ; dttyp(22)=8
    2083              : !cg1
    2084          718 :  cmpw(23)=2*nspinor*my_mband*mk1mem*nsppol  ; dttyp(23)=8
    2085              : !rhor1,rhog1
    2086          718 :  cfft(24)=cplex*nspden+2       ; dttyp(24)=8
    2087              : !eigen1
    2088              : !assume than nkpt=nkpt_rbz
    2089          718 :  cadd(25)=2*mband*mband*nkpt*nsppol      ; dttyp(25)=8
    2090              : !ylm
    2091          718 :  cmpw(26)=mkmem*mpsang*mpsang*useylm     ; dttyp(26)=8
    2092              : 
    2093              : !(3)                     in dfpt_scfcv --------------------------------------
    2094              : 
    2095              : !vhartr1,vtrial1,vxc
    2096          718 :  cfft(31)=cplex+cplex*nspden+nspden      ; dttyp(31)=8
    2097          718 :  if(iscf>0)then
    2098              : !  f_fftgr
    2099          521 :    cfft(32)=cplex*nspden*n_fftgr*mffmem    ; dttyp(32)=8
    2100              :  end if
    2101              : 
    2102              : !(4)                   in dfpt_vtorho----------------------------------------
    2103              : 
    2104              : !proc_distrb
    2105          718 :  cadd(41)=nkpt*mband*nsppol    ; dttyp(41)=4
    2106              : !kg_k,kg1_k
    2107          718 :  cmpw(42)=6                    ; dttyp(42)=4
    2108              : !rhoaug1, vlocal, vlocal1
    2109          718 :  cfft(43)=2*cplex+1            ; dttyp(43)=8
    2110          718 :  cadd(43)=(2*cplex+1)*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    2111              : 
    2112          718 :  if(mkqmem==0)then
    2113              : !  cgq_disk
    2114            0 :    cmpw(45)=2*nspinor*my_mband      ; dttyp(45)=8
    2115              :  end if
    2116              : !doccde_k,doccde_kq,eig0_k, ..., eig1_k, rocceig
    2117          718 :  cadd(47)=(14+3*mband)*mband   ; dttyp(47)=8
    2118              : !ylm_k,ylm1_k
    2119          718 :  cmpw(49)=2*mpsang*mpsang*useylm  ; dttyp(49)=8
    2120              : 
    2121              : !(5)                     in dfpt_vtowfk --------------------------------------
    2122              : 
    2123              : !dkinpw,kinpw1
    2124          718 :  cmpw(51)=2                    ; dttyp(51)=8
    2125              : !ffnlk,ffnl1,ffnlkq
    2126          718 :  cmpw(52)=2*(ntypat+2)*lmnmax  ; dttyp(52)=8
    2127              : !ghc,gvnlxc,gvnlx1
    2128          718 :  cmpw(53)=6*nspinor            ; dttyp(53)=8
    2129              : !ph3d
    2130          718 :  matblk=NLO_MINCAT
    2131          718 :  if(nloalg(2)<=0)matblk=natom
    2132          718 :  cmpw(54)=2*matblk             ; dttyp(54)=8
    2133              : !wfraug,wfraug1,rhoaug
    2134          718 :  cfft(55)=5                    ; dttyp(55)=8
    2135          718 :  cadd(55)=5*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    2136              : !cwavef,cwave0,cwave1
    2137          718 :  cmpw(56)=6*nspinor            ; dttyp(56)=8
    2138              : 
    2139              : !(6)                     in dfpt_cgwf ----------------------------------------
    2140              : 
    2141              : !gh1, gh_direc, gvnlx_direc, conjgr, direc, vresid, cwaveq
    2142          718 :  cmpw(61)=14*nspinor            ; dttyp(61)=8
    2143              : 
    2144              : !(9a)                    in getghc and fourwf----------------------------
    2145              : 
    2146              : !work (in getghc)
    2147          718 :  cfft(91)=2                    ; dttyp(91)=8
    2148          718 :  cadd(92)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    2149              : !work1 (in fourwf)
    2150          718 :  cfft(92)=2                    ; dttyp(92)=8
    2151              :  cadd(92)=2*(ngfft(4)*ngfft(5)*ngfft(6)-nfft)
    2152              : 
    2153              : !(9b)                    in getghc, nonlop and opernl--------------------
    2154          718 :  mincat=min(NLO_MINCAT,natom-ntypat+1)
    2155          718 :  if (useylm==0) then                          ! ===== nonlop_pl
    2156              : !  gxa  (in nonlop)
    2157          544 :    cadd(94)=2*20*mincat*2       ; dttyp(94)=8
    2158              : !  dgxdt  (in nonlop)
    2159          544 :    cadd(95)=2*3*20*mincat*2    ; dttyp(95)=8
    2160              : !  dgxds  (in nonlop)
    2161          544 :    cadd(96)=2*56*mincat*2      ; dttyp(96)=8
    2162              : !  teffv (in opernl4 - no distinction is made for opernl, opernl2 or opernl3)
    2163              : !  kpgx, ffkg
    2164              : !  here, evaluate an upper value, with nproj=2, p,d and f orbitals, but not
    2165              : !  considering the stress, since it will be called outside of the main chain
    2166          544 :    cadd(97)=NLO_MBLKPW*40        ; dttyp(97)=8
    2167              : !  kpg if nloalg(3)=1
    2168          544 :    cadd(98)=3*mpw*nloalg(3)     ; dttyp(98)=8
    2169              :  else                                        ! ===== nonlop_ylm
    2170              : !  gx + gxfac
    2171          174 :    cadd(94)=2* 2*mpw*lmnmax*mincat    ; dttyp(94)=8
    2172              : !  dgxdt + dgxdtfac + d2gxdt
    2173          174 :    if (optddk>0.and.optphon==0.and.optstrs==0) cadd(95)=2*2*mpw*lmnmax*mincat
    2174          174 :    if (optphon>0) cadd(95)=12*2*mpw*lmnmax*mincat
    2175          174 :    if (optstrs>0) cadd(95)=72*2*mpw*lmnmax*mincat
    2176          174 :    dttyp(95)=8
    2177              : !  kpg
    2178          174 :    cadd(96)=2*3*mpw       ; dttyp(96)=8
    2179          174 :    if (optphon>0) cadd(96)=cadd(96)+2*6*mpw
    2180              : !  miscelaneous: indlmn_typ, ffnl_typ
    2181          174 :    cadd(97)=lmnmax*(6+mpw*(2+optstrs)); dttyp(97)=8
    2182              : !  opernla_ylm: scalar,scali,scalarr,scalari
    2183          174 :    cadd(98)=2*mpw+2*mpw
    2184          174 :    if (optddk>0.and.optstrs==0) cadd(98)=cadd(98)+2*mpw
    2185          174 :    if (optstrs>0) cadd(98)=cadd(98)+9*2*mpw
    2186          174 :    dttyp(98)=8
    2187              :  end if
    2188              : 
    2189              : !--------------------------------------------------------------------------
    2190              : 
    2191       217554 :  chain(:,:)=.true.
    2192              : 
    2193              : !Define the main chain version a (fourwf)
    2194         6462 :  chain(93:100,1)=.false.
    2195              : 
    2196              : !Define the main chain version b (nonlop+opernl)
    2197         2154 :  chain(91:92,2)=.false.
    2198              : 
    2199              : !The memory needed for each chain has been computed
    2200              : !-------------------------------------------------------------------------
    2201              : !Still need some auxiliary data : estimate the disk space
    2202              : !or the maximum segment size.
    2203              : 
    2204              : !XG030513 : MPIWF need to multiply mbdiskwf by the number of processors
    2205              : !in the WF group. For the time being, nprocwf=1
    2206          718 :  nprocwf=mpi_enreg%nproc_fft
    2207              : 
    2208        57084 :  mbdiskwf=(8*2*mpw*nprocwf*sum(nband(1:nkpt*nsppol)))/1024._dp**2 + 0.002_dp
    2209          718 :  mbdiskpd=(8*nfft*nsppol)/1024._dp**2 + 0.002_dp
    2210              : 
    2211              : !Determine the largest array out of cg,cg1,cgq, cg_disk or f_fftgr (f_fftgr_disk)
    2212          718 :  if(mkmem==0 .and. mk1mem==0 .and. mkqmem==0)then
    2213            0 :    mbcg=(8*2*mpw*nspinor*my_mband)/1024._dp**2 + 0.002_dp
    2214              :  else
    2215         2872 :    maxmkmem=maxval(mkmems(:))
    2216          718 :    mbcg=(8*2*mpw*nspinor*my_mband*maxmkmem*nsppol)/1024._dp**2 + 0.002_dp
    2217              :  end if
    2218          718 :  if(mffmem==0)then
    2219            0 :    mbf_fftgr=(8*cplex*nfft*n_fftgr)/1024._dp**2 + 0.002_dp
    2220              :  else
    2221          718 :    mbf_fftgr=(8*cplex*nfft*n_fftgr*nspden*mffmem)/1024._dp**2 + 0.002_dp
    2222              :  end if
    2223              : 
    2224              : !---------------------------------------------------------------------
    2225              : !Now, analyze the data
    2226              : 
    2227              : !DEBUG
    2228              : !write(std_out,*)' memorf : nchain=',nchain
    2229              : !ENDDEBUG
    2230              : 
    2231          718 :  ABI_MALLOC(cfft_dum,(marrays))
    2232       108418 :  cfft_dum=zero
    2233          718 :  mbgylm=zero
    2234              :  call memana(cadd,cfft,cfft_dum,chain,cmpw,dttyp,iout,iprcel,iscf,&
    2235              : & marrays,mbcg,mbdiskpd,mbdiskwf,mbf_fftgr,mbgylm,mffmem,&
    2236          718 : & mpw,natom,nchain,nfft,nfft,occopt,option,prtvol)
    2237          718 :  ABI_FREE(cfft_dum)
    2238              : 
    2239          718 : end subroutine memorf
    2240              : !!***
    2241              : 
    2242              : !!****f* m_memeval/getdim_nloc
    2243              : !! NAME
    2244              : !! getdim_nloc
    2245              : !!
    2246              : !! FUNCTION
    2247              : !! Determine the dimensions of arrays that contain
    2248              : !! the definition of non-local projectors : ekb, ffspl, indlmn
    2249              : !!
    2250              : !! INPUTS
    2251              : !!  mixalch(npspalch,ntypalch,nimage)=alchemical mixing coefficients
    2252              : !!  nimage=number of images
    2253              : !!  npsp=number of pseudopotentials
    2254              : !!  npspalch=number of pseudopotentials for alchemical purposes
    2255              : !!  ntypat=number of types of pseudo atoms
    2256              : !!  ntypalch=number of types of alchemical pseudo atoms
    2257              : !!  pspheads(npsp)=<type pspheader_type>all the important information from the
    2258              : !!   pseudopotential file headers, as well as the psp file names
    2259              : !!
    2260              : !! OUTPUT
    2261              : !!  lmnmax=maximum number of l,m,n projectors, not taking into account the spin-orbit
    2262              : !!  lmnmaxso=maximum number of l,m,n projectors, taking into account the spin-orbit
    2263              : !!  lnmax=maximum number of l,n projectors, not taking into account the spin-orbit
    2264              : !!  lnmaxso=maximum number of l,n projectors, taking into account the spin-orbit
    2265              : !!
    2266              : !! SOURCE
    2267              : 
    2268        12243 : subroutine getdim_nloc(lmnmax,lmnmaxso,lnmax,lnmaxso,mixalch,nimage,npsp,npspalch,&
    2269        12243 : & ntypat,ntypalch,pspheads)
    2270              : 
    2271              : !Arguments ------------------------------------
    2272              : !scalars
    2273              :  integer,intent(in) :: nimage,npsp,npspalch,ntypalch,ntypat
    2274              :  integer,intent(out) :: lmnmax,lmnmaxso,lnmax,lnmaxso
    2275              : !arrays
    2276              :  real(dp),intent(in) :: mixalch(npspalch,ntypalch,nimage)
    2277              :  type(pspheader_type),intent(in) :: pspheads(npsp)
    2278              : 
    2279              : !Local variables-------------------------------
    2280              : !scalars
    2281              :  integer :: ilang,ipsp,ipspalch,itypalch,itypat,ntyppure,max_l
    2282              : !integer :: llmax
    2283              :  character(len=500) :: msg
    2284              : !arrays
    2285        12243 :  integer,allocatable :: lmnproj_typat(:),lmnprojso_typat(:),lnproj_typat(:)
    2286        12243 :  integer,allocatable :: lnprojso_typat(:),nproj_typat(:,:),nprojso_typat(:,:)
    2287              : 
    2288              : ! *************************************************************************
    2289              : 
    2290              : !write(std_out,*)' getdim_nloc: 'pspheads(1)%nproj(0:3)=',pspheads(1)%nproj(0:3)
    2291        12243 :  max_l=0
    2292        27766 :  do ipsp=1,npsp
    2293        27766 :    max_l=max(max_l,pspheads(ipsp)%lmax)
    2294              :  enddo
    2295              : 
    2296        36729 :  ABI_MALLOC(lmnproj_typat,(ntypat))
    2297        24486 :  ABI_MALLOC(lmnprojso_typat,(ntypat))
    2298        24486 :  ABI_MALLOC(lnproj_typat,(ntypat))
    2299        24486 :  ABI_MALLOC(lnprojso_typat,(ntypat))
    2300        48972 :  ABI_MALLOC(nproj_typat,(0:max_l,ntypat))
    2301        48972 :  ABI_MALLOC(nprojso_typat,(max_l,ntypat))
    2302        55452 :  lmnproj_typat(:)=0 ; lmnprojso_typat(:)=0
    2303        55452 :  lnproj_typat(:)=0 ; lnprojso_typat(:)=0
    2304       112457 :  nproj_typat(:,:)=0 ; nprojso_typat(:,:)=0
    2305              : 
    2306        12243 :  ntyppure=ntypat-ntypalch
    2307              : 
    2308              : !For each type of pseudo atom, compute the number of projectors
    2309              : !First, pure pseudo atoms
    2310        12243 :  if(ntyppure>0)then
    2311        27666 :    do itypat=1,ntyppure
    2312        50746 :      nproj_typat(0:pspheads(itypat)%lmax,itypat)=pspheads(itypat)%nproj(0:pspheads(itypat)%lmax)
    2313        47530 :      nprojso_typat(1:pspheads(itypat)%lmax,itypat)=pspheads(itypat)%nprojso(1:pspheads(itypat)%lmax)
    2314              :    end do
    2315              :  end if
    2316              : 
    2317              : !Then, alchemical pseudo atoms
    2318        12243 :  if(ntypalch>0)then
    2319           80 :    do itypat=ntyppure+1,ntypat
    2320           40 :      itypalch=itypat-ntyppure
    2321          160 :      do ipsp=ntyppure+1,npsp
    2322           80 :        ipspalch=ipsp-ntyppure
    2323              : !      If there is some mixing, must accumulate the projectors
    2324          224 :        if(sum(abs(mixalch(ipspalch,itypalch,:)))>tol10)then
    2325              :          nproj_typat(0:pspheads(ipsp)%lmax,itypat)=nproj_typat(0:pspheads(ipsp)%lmax,itypat)+&
    2326          242 : &                                                  pspheads(ipsp)%nproj(0:pspheads(ipsp)%lmax)
    2327              :          nprojso_typat(1:pspheads(ipsp)%lmax,itypat)=nprojso_typat(1:pspheads(ipsp)%lmax,itypat)+&
    2328          176 : &                                                    pspheads(ipsp)%nprojso(:)
    2329              :        end if
    2330              :      end do
    2331              :    end do
    2332              :  end if
    2333              : 
    2334              : !Now that the number of projectors is known, accumulate the dimensions
    2335        27726 :  do itypat=1,ntypat
    2336        51727 :    do ilang=0,max_l
    2337        36244 :      lnproj_typat(itypat)=lnproj_typat(itypat)+nproj_typat(ilang,itypat)
    2338        51727 :      lmnproj_typat(itypat)=lmnproj_typat(itypat)+nproj_typat(ilang,itypat)*(2*ilang+1)
    2339              :    end do
    2340        15483 :    lnprojso_typat(itypat)=lnproj_typat(itypat)
    2341        15483 :    lmnprojso_typat(itypat)=lmnproj_typat(itypat)
    2342        48487 :    do ilang=1,max_l
    2343        20761 :      lnprojso_typat(itypat)=lnprojso_typat(itypat)+nprojso_typat(ilang,itypat)
    2344        36244 :      lmnprojso_typat(itypat)=lmnprojso_typat(itypat)+nprojso_typat(ilang,itypat)*(2*ilang+1)
    2345              :    end do
    2346              :  end do
    2347              : 
    2348              : !Compute the maximal bounds, at least equal to 1, even for local psps
    2349        12243 :  lmnmax=1;lmnmaxso=1;lnmax=1;lnmaxso=1
    2350        27726 :  do itypat=1,ntypat
    2351        15483 :    lmnmax  =max(lmnmax  ,lmnproj_typat  (itypat))
    2352        15483 :    lmnmaxso=max(lmnmaxso,lmnprojso_typat(itypat))
    2353        15483 :    lnmax   =max(lnmax   ,lnproj_typat   (itypat))
    2354        27726 :    lnmaxso =max(lnmaxso ,lnprojso_typat (itypat))
    2355              :  end do
    2356              : !The initial coding (below) was not totally portable (MT 110215)
    2357              : !lmnmax=max(maxval(lmnproj_typat(1:ntypat)),1)
    2358              : !lmnmaxso=max(maxval(lmnprojso_typat(1:ntypat)),1)
    2359              : !lnmax=max(maxval(lnproj_typat(1:ntypat)),1)
    2360              : !lnmaxso=max(maxval(lnprojso_typat(1:ntypat)),1)
    2361              : 
    2362        27726 :  if(maxval(lmnproj_typat(1:ntypat))==0)then
    2363              :    write(msg, '(3a)' )&
    2364         2724 :     'Despite there is only a local part to pseudopotential(s),',ch10,&
    2365         5448 :     'lmnmax and lnmax are set to 1.'
    2366         2724 :    ABI_COMMENT(msg)
    2367              :  end if
    2368              : 
    2369              : !XG040806 : These lines make modifications of lnmax and lmnmax
    2370              : !that are unjustified in many cases, according to the many tests cases
    2371              : !where they produce a changes, while the test case was working properly.
    2372              : !One should understand better the needs, and code more appropriate changes ...
    2373              : !lnmax/lmnmax has to be bigger than 1+lmax (for compatibility reasons)
    2374              : !llmax=maxval(pspheads(1:ntypat)%lmax)+1 ! And this line might have trouble with HP compiler
    2375              : !if (lnmax   <llmax) lnmax=llmax
    2376              : !if (lnmaxso <llmax) lnmaxso=llmax
    2377              : !if (lmnmax  <llmax) lmnmax=llmax
    2378              : !if (lmnmaxso<llmax) lmnmaxso=llmax
    2379              : 
    2380        12243 :  write(msg, '(a,a,i4,a,i4,3a,i4,a,i4,a)' ) ch10,&
    2381        12243 :  ' getdim_nloc: deduce lmnmax  =',lmnmax,', lnmax  =',lnmax,',',ch10,&
    2382        24486 :  '                      lmnmaxso=',lmnmaxso,', lnmaxso=',lnmaxso,'.'
    2383        12243 :  call wrtout(std_out,msg)
    2384              : 
    2385        12243 :  ABI_FREE(lmnproj_typat)
    2386        12243 :  ABI_FREE(lmnprojso_typat)
    2387        12243 :  ABI_FREE(lnproj_typat)
    2388        12243 :  ABI_FREE(lnprojso_typat)
    2389        12243 :  ABI_FREE(nproj_typat)
    2390        12243 :  ABI_FREE(nprojso_typat)
    2391              : 
    2392        12243 : end subroutine getdim_nloc
    2393              : !!***
    2394              : 
    2395              : !!****f* m_memeval/setmqgrid
    2396              : !! NAME
    2397              : !!  setmqgrid
    2398              : !!
    2399              : !! FUNCTION
    2400              : !!  Sets the number of points needed to represent the pseudopotentials in
    2401              : !!  reciprocal space for a specified resolution.
    2402              : !!
    2403              : !! INPUTS
    2404              : !!  ecut=cutoff energy for the wavefunctions
    2405              : !!  ecutdg=cutoff energy for the fine grid in case usepaw==1
    2406              : !!  gprimd=primitive translation vectors for reciprocal space
    2407              : !!  nptsgvec=number of points along the smallest primitive translation vector
    2408              : !!    of the reciprocal space
    2409              : !!  usepaw=1 if PAW is used, 0 otherwise
    2410              : !!
    2411              : !! OUTPUT
    2412              : !!
    2413              : !! SOURCE
    2414              : 
    2415        12243 : subroutine setmqgrid(mqgrid,mqgriddg,ecut,ecutdg,gprimd,nptsgvec,usepaw)
    2416              : 
    2417              : !Arguments ------------------------------------
    2418              :  integer , intent(inout)  :: mqgrid,mqgriddg
    2419              :  integer , intent(in)  :: nptsgvec,usepaw
    2420              :  real(dp), intent(in) :: ecut,ecutdg
    2421              :  real(dp), intent(in) :: gprimd(3,3)
    2422              : 
    2423              : !Local variables-------------------------------
    2424              :  integer :: mqgrid2,mqgriddg2
    2425              :  real(dp) :: gmax,gmaxdg,gvecnorm
    2426              :  character(len=500) :: msg
    2427              : 
    2428              : ! *************************************************************************
    2429              : 
    2430              :  gvecnorm=sqrt(min(dot_product(gprimd(:,1),gprimd(:,1)), &
    2431              : & dot_product(gprimd(:,2),gprimd(:,2)), &
    2432       122430 : & dot_product(gprimd(:,3),gprimd(:,3))))
    2433        12243 :  gmax=one/(sqrt2*pi)*sqrt(ecut)
    2434              : 
    2435        12243 :  if (mqgrid == 0) then
    2436        12239 :    mqgrid2=ceiling(gmax/gvecnorm*nptsgvec)
    2437        12239 :    mqgrid=max(mqgrid2,3001)
    2438              :    write(msg, '(5a,i0,a)' )&
    2439        12239 : &   'The number of points "mqgrid" in reciprocal space used for the',ch10,&
    2440        12239 : &   'description of the pseudopotentials has been set automatically',ch10,&
    2441        24478 : &   'by abinit to: ',mqgrid,'.'
    2442              :    !ABI_COMMENT(msg)
    2443              :  else
    2444            4 :    mqgrid2=ceiling(gmax/gvecnorm*nptsgvec)
    2445            4 :    if (mqgrid2>mqgrid) then
    2446              :      write(msg, '(3a,i8,3a,i8,3a)' )&
    2447            0 : &     'The number of points "mqgrid" in reciprocal space used for the',ch10,&
    2448            0 : &     'description of the pseudopotentials is : ',mqgrid,'.',ch10,&
    2449            0 : &     'It would be better to increase it to at least ',mqgrid2,', or',ch10,&
    2450            0 : &     'let abinit choose it automatically by setting mqgrid = 0.'
    2451            0 :      ABI_WARNING(msg)
    2452              :    end if
    2453              :  end if
    2454              : 
    2455        12243 :  if (usepaw==1) then
    2456         2356 :    if(ecutdg<tol6)then
    2457            0 :      write(msg,'(a)')'The value of (paw)ecutdg is zero or negative, which is forbidden.'
    2458            0 :      ABI_ERROR(msg)
    2459              :    end if
    2460         2356 :    gmaxdg=one/(sqrt2*pi)*sqrt(ecutdg)
    2461         2356 :    if (mqgriddg == 0) then
    2462         2354 :      mqgriddg2=ceiling(gmaxdg/gvecnorm*nptsgvec)
    2463         2354 :      mqgriddg=max(mqgriddg2,3001)
    2464              :      write(msg, '(5a,i0,a)' )&
    2465         2354 : &     'The number of points "mqgriddg" in reciprocal space used for the',ch10,&
    2466         2354 : &     'description of the pseudopotentials has been set automatically',ch10,&
    2467         4708 : &     'by abinit to: ',mqgriddg,'.'
    2468              :      !ABI_COMMENT(msg)
    2469              :    else
    2470            2 :      mqgriddg2=ceiling(gmax/gvecnorm*nptsgvec)
    2471            2 :      if (mqgriddg2>mqgriddg) then
    2472              :        write(msg, '(3a,i8,3a,i8,3a)' )&
    2473            0 : &       'The number of points "mqgriddg" in reciprocal space used for the',ch10,&
    2474            0 : &       'description of the pseudopotentials (fine grid) is :',mqgriddg,'.',ch10,&
    2475            0 : &       'It would be better to increase it to at least ',mqgriddg2,', or',ch10,&
    2476            0 : &       'let abinit choose it automatically by setting mqgrid = 0.'
    2477            0 :        ABI_WARNING(msg)
    2478              :      end if
    2479              :    end if
    2480              :  end if
    2481              : 
    2482        12243 : end subroutine setmqgrid
    2483              : !!***
    2484              : 
    2485              : !!****f* m_memeval/wvl_memory
    2486              : !! NAME
    2487              : !! wvl_memory
    2488              : !!
    2489              : !! FUNCTION
    2490              : !! Estimation of the memory needed for waelet based computation job.
    2491              : !! According to the value of the option variable,
    2492              : !! might also try to allocate this amount of memory, and if it fails,
    2493              : !! might estimate the available memory.
    2494              : !!
    2495              : !! INPUTS
    2496              : !!  dtset=<type datafiles_type>contains all input variables.
    2497              : !!  idtset=number of the current dataset
    2498              : !!  mpi_enreg=information about MPI parallelization
    2499              : !!  npsp=number of pseudopotentials
    2500              : !!  option: if 0, no test of available memory
    2501              : !!          if 1, the routine tries to allocate the estimated memory, for testing
    2502              : !!                purposes, and if a failure occurs, the routine stops.
    2503              : !!          if 2, like 1, but before stopping, the routine will provide
    2504              : !!                an estimation of the available memory.
    2505              : !!  pspheads(npsp)=<type pspheader_type>all the important information from the
    2506              : !!   pseudopotential file header, as well as the psp file name
    2507              : !!
    2508              : !! OUTPUT
    2509              : !!  (only writing)
    2510              : !!
    2511              : !! NOTES
    2512              : !! The estimator is the one provided by BigDFT.
    2513              : !!
    2514              : !! SOURCE
    2515              : 
    2516            0 : subroutine wvl_memory(dtset, idtset, mpi_enreg, npsp, option, pspheads)
    2517              : 
    2518              :  use defs_wvltypes
    2519              :  use m_abi2big, only : wvl_setBoxGeometry
    2520              :  use m_wvl_descr_psp,    only : wvl_descr_free, wvl_descr_atoms_set
    2521              : 
    2522              : #if defined HAVE_BIGDFT
    2523              :  use BigDFT_API, only: MemoryEstimator, createWavefunctionsDescriptors, deallocate_lr, &
    2524              :       & atomic_info, memory_estimation
    2525              : #endif
    2526              : 
    2527              : !Arguments ------------------------------------
    2528              :   !scalars
    2529              :   integer,intent(in) :: idtset, npsp, option
    2530              :   type(dataset_type),intent(in) :: dtset
    2531              :   type(MPI_type),intent(in) :: mpi_enreg
    2532              :   !arrays
    2533              :   type(pspheader_type),intent(in) :: pspheads(npsp)
    2534              : 
    2535              : !Local variables-------------------------------
    2536              : #if defined HAVE_BIGDFT
    2537              :   !scalars
    2538              :   integer :: ityp, i, mu, nstates, me, nproc, comm
    2539              :   character(len=500) :: msg
    2540              :   real(dp) :: ehomo, radfine
    2541              :   type(wvl_internal_type) :: wvl
    2542              :   type(memory_estimation) :: peakmem
    2543              :   !arrays
    2544              :   real(dp) :: acell(3), rprimd(3,3), rprim(3,3)
    2545              :   real(dp), allocatable :: radii_cf(:,:)
    2546              :   real(dp), allocatable :: xred(:,:), xcart(:,:)
    2547              : #endif
    2548              : 
    2549              : ! **************************************************************************
    2550              : 
    2551              : #if defined HAVE_BIGDFT
    2552              : 
    2553              :  comm=mpi_enreg%comm_wvl
    2554              :  me=xmpi_comm_rank(comm)
    2555              :  nproc=xmpi_comm_size(comm)
    2556              : 
    2557              :  if(option<0 .or. option>2)then
    2558              :    write(msg, '(A,A,A,A,I0,A)') ch10,&
    2559              : &   ' wvl_memory : BUG -',ch10,&
    2560              : &   '  option=',option,' while the only allowed values are 0, 1, or 2.'
    2561              :    call wrtout(std_out,msg)
    2562              :  end if
    2563              : 
    2564              :  wvl%paw%usepaw=0 !no PAW here
    2565              :  nullify(wvl%rholoc%d)
    2566              :  nullify(wvl%rholoc%msz)
    2567              :  nullify(wvl%rholoc%rad)
    2568              :  nullify(wvl%rholoc%radius)
    2569              :  nullify(wvl%paw%spsi)
    2570              :  nullify(wvl%paw%indlmn)
    2571              :  nullify(wvl%paw%spsi)
    2572              :  nullify(wvl%paw%indlmn)
    2573              : 
    2574              :  write(msg,*)' wvl_memory : analysis of memory needs '
    2575              :  call wrtout(std_out,msg)
    2576              : 
    2577              :  if(idtset>=100)then
    2578              :    write(msg,'(80a,a,a,i5,a)')('=',mu=1,80),ch10,&
    2579              : &   ' Values of the parameters that define the memory need for DATASET', idtset,&
    2580              : &   ' (WVL).'
    2581              :  else if(idtset/=0)then
    2582              :    write(msg,'(80a,a,a,i3,a)')('=',mu=1,80),ch10,&
    2583              : &   ' Values of the parameters that define the memory need for DATASET', idtset,&
    2584              : &   ' (WVL).'
    2585              :  else
    2586              :    write(msg,'(80a,a,a,a)')('=',mu=1,80),ch10,&
    2587              : &   ' Values of the parameters that define the memory need of the present run',&
    2588              : &   ' (WVL).'
    2589              :  end if
    2590              :  call wrtout(ab_out,msg)
    2591              :  call wrtout(std_out,msg)
    2592              : 
    2593              :  write(msg,'( a,f7.3,a,i7,2(a,F7.3),a,a,f7.3,a,i7 )' ) &
    2594              : & '  wvl_hgrid =', dtset%wvl_hgrid , '   nwfshist =', dtset%nwfshist, &
    2595              : & ' wvl_crmult =', dtset%wvl_crmult, ' wvl_frmult =', dtset%wvl_frmult, ch10,&
    2596              : & '  tl_radius =', dtset%tl_radius , '  tl_nprccg =', dtset%tl_nprccg
    2597              :  call wrtout(ab_out,msg)
    2598              :  call wrtout(std_out,msg)
    2599              : 
    2600              :  if (dtset%nsppol == 2) then
    2601              :    nstates = dtset%nelect
    2602              :  else
    2603              :    nstates = dtset%mband
    2604              :  end if
    2605              :  write(msg,'(4(a,i7))')&
    2606              : & '      natom =', dtset%natom, '     ntypat =', dtset%ntypat, &
    2607              : & '    nstates =', nstates,     '     nsppol =', dtset%nsppol
    2608              :  call wrtout(ab_out,msg)
    2609              :  call wrtout(std_out,msg)
    2610              : 
    2611              :  write(msg,'(80a)') ('=',mu=1,80)
    2612              :  call wrtout(ab_out,msg)
    2613              :  call wrtout(std_out,msg)
    2614              : 
    2615              : !First, use eleconf to get radii_cf().
    2616              :  ABI_MALLOC(radii_cf,(npsp, 3))
    2617              :  do ityp = 1, npsp, 1
    2618              :    call atomic_info(int(pspheads(ityp)%znuclpsp), int(pspheads(ityp)%zionpsp), ehomo = ehomo)
    2619              : 
    2620              : !  new method for assigning the radii
    2621              :    radii_cf(ityp, 1) = one / sqrt(abs(two * ehomo))
    2622              :    radfine = 100.d0
    2623              :    do i = 0, 4, 1
    2624              :      if (pspheads(ityp)%GTHradii(i) /= zero) then
    2625              :        radfine = min(radfine, pspheads(ityp)%GTHradii(i))
    2626              :      end if
    2627              :    end do
    2628              :    radii_cf(ityp,2) = radfine
    2629              :  end do
    2630              : 
    2631              : !Compute the shifted positions and acell
    2632              :  acell = dtset%acell_orig(1:3,1)
    2633              :  call wvl_descr_atoms_set(acell, dtset%icoulomb, dtset%natom, dtset%ntypat, dtset%typat, wvl)
    2634              :  ABI_MALLOC(xred,(3, dtset%natom))
    2635              :  xred = dtset%xred_orig(:,:,1)
    2636              :  rprimd = dtset%rprimd_orig(1:3,1:3,1)
    2637              :  wvl%h(:) = dtset%wvl_hgrid
    2638              :  call wvl_setBoxGeometry(1, radii_cf, rprimd, xred, &
    2639              : & wvl, dtset%wvl_crmult, dtset%wvl_frmult)
    2640              : !Compute acell and rprim from rprimd
    2641              :  call mkradim(acell,rprim,rprimd)
    2642              :  ABI_MALLOC(xcart,(3, dtset%natom))
    2643              :  call xred2xcart(dtset%natom, rprimd, xcart, xred)
    2644              :  call createWavefunctionsDescriptors(me, wvl%h(1), wvl%h(2), wvl%h(3), &
    2645              : & wvl%atoms, xcart, radii_cf, dtset%wvl_crmult, dtset%wvl_frmult, wvl%Glr)
    2646              :  call MemoryEstimator(nproc, dtset%nwfshist, wvl%Glr, &
    2647              : & dtset%mband, dtset%nspinor, dtset%nkpt, 0, dtset%nsppol, &
    2648              : & 0, dtset%iscf, peakmem)
    2649              : 
    2650              :  call deallocate_lr(wvl%Glr)
    2651              :  call wvl_descr_free(wvl)
    2652              :  ABI_FREE(radii_cf)
    2653              :  ABI_FREE(xred)
    2654              :  ABI_FREE(xcart)
    2655              : 
    2656              :  write(msg,'(80a,a)') ('=',mu=1,80), ch10
    2657              :  call wrtout(ab_out,msg)
    2658              :  call wrtout(std_out,msg)
    2659              : 
    2660              : #else
    2661            0 :  BIGDFT_NOTENABLED_ERROR()
    2662              :  if (.false.) write(std_out,*) idtset,npsp,option,dtset%nstep,mpi_enreg%nproc,pspheads(1)%zionpsp
    2663              : #endif
    2664              : 
    2665              : end subroutine wvl_memory
    2666              : !!***
    2667              : 
    2668              : end module m_memeval
    2669              : !!***
        

Generated by: LCOV version 2.3-1