TABLE OF CONTENTS


ABINIT/driver [ Functions ]

[ Top ] [ Functions ]

NAME

 driver

FUNCTION

 Driver for ground state, response function, screening
 and sigma calculations. The present routine drives the following operations.
 An outer loop allows computation related to different data sets.
 For each data set, either a GS calculation, a RF calculation,
 a SUS calculation, a SCR calculation or a SIGMA calculation is made.
 In both cases, the input variables are transferred in the proper variables,
 selected big arrays are allocated, then the gstate, respfn, ...  subroutines are called.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (XG,MKV,MM,MT,FJ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

 codvsn= code version
 cpui=initial CPU time
 filnam(5)=character strings giving file names
 filstat=character strings giving name of status file
 mpi_enregs=informations about MPI parallelization
 ndtset=number of datasets
 ndtset_alloc=number of datasets, corrected for allocation of at
               least one data set.
 npsp=number of pseudopotentials
 pspheads(npsp)=<type pspheader_type>all the important information from the
   pseudopotential file header, as well as the psp file name

OUTPUT

SIDE EFFECTS

  Input/Output
 results_out(0:ndtset_alloc)=<type results_out_type>contains the results
   needed for outvars, including evolving variables
   Default values are set up in the calling routine
 dtsets(0:ndtset_alloc)=<type datasets_type>
   input: all input variables initialized from the input file.
   output: the effective set of variables used in the different datasets.
           Some variables, indeed, might have been redefined in one of the children.
           of this routine.

NOTES

 The array filnam is used for the name of input and output files,
 and roots for generic input, output or temporary files.
 Pseudopotential file names are set in pspini and pspatm,
 using another name. The name filstat will be needed beyond gstate to check
 the appearance of the "exit" flag, to make a hasty exit, as well as
 in order to output the status of the computation.

PARENTS

      abinit

CHILDREN

      abi_io_redirect,abi_linalg_finalize,abi_linalg_init,bethe_salpeter
      chkdilatmx,destroy_results_out,destroy_results_respfn,dtfil_init
      dtfil_init_img,dtset_copy,dtset_free,echo_xc_name,eph,exit_check
      f_malloc_set_status,fftw3_cleanup,fftw3_init_threads,find_getdtset
      gather_results_out,gstateimg,gwls_sternheimer,init_results_respfn
      libpaw_write_comm_set,libxc_functionals_end,libxc_functionals_init
      mkrdim,mpi_environment_set,nonlinear,nullify_wvl_data,pawang_free
      pawrad_free,pawtab_free,pawtab_nullify,psps_free,psps_init_from_dtset
      psps_init_global,respfn,screening,sigma,status,timab,wfk_analyze,wrtout
      xc_vdw_done,xc_vdw_init,xc_vdw_libxc_init,xc_vdw_memcheck,xc_vdw_read
      xc_vdw_show,xc_vdw_trigger,xc_vdw_write,xcart2xred,xg_finalize
      xgscalapack_config,xmpi_bcast,xred2xcart

SOURCE

 75 #if defined HAVE_CONFIG_H
 76 #include "config.h"
 77 #endif
 78 
 79 #include "abi_common.h"
 80 
 81 subroutine driver(codvsn,cpui,dtsets,filnam,filstat,&
 82 &                 mpi_enregs,ndtset,ndtset_alloc,npsp,pspheads,results_out)
 83 
 84  use defs_basis
 85  use defs_parameters
 86  use defs_datatypes
 87  use defs_abitypes
 88  use defs_wvltypes
 89  use m_errors
 90  use m_results_out
 91  use m_results_respfn
 92  use m_xmpi
 93  use m_xomp
 94  use m_abi_linalg
 95  use m_profiling_abi
 96  use m_exit
 97  use libxc_functionals
 98 #if defined DEV_YP_VDWXC
 99  use m_xc_vdw
100 #endif
101  use m_xgScalapack
102 
103  use m_xg,           only : xg_finalize
104  use m_libpaw_tools, only : libpaw_write_comm_set
105  use m_geometry,     only : mkrdim, xcart2xred, xred2xcart, chkdilatmx
106  use m_pawang,       only : pawang_type, pawang_free
107  use m_pawrad,       only : pawrad_type, pawrad_free
108  use m_pawtab,       only : pawtab_type, pawtab_nullify, pawtab_free
109  use m_fftw3,        only : fftw3_init_threads, fftw3_cleanup
110  use m_psps,         only : psps_init_global, psps_init_from_dtset, psps_free
111  use m_dtset,        only : dtset_copy, dtset_free, find_getdtset
112  use m_mpinfo,       only : mpi_distrib_is_ok
113 
114 #if defined HAVE_BIGDFT
115  use BigDFT_API,   only: xc_init, xc_end, XC_MIXED, XC_ABINIT,&
116 &                        mpi_environment_set,bigdft_mpi, f_malloc_set_status
117 #endif
118 
119 !This section has been created automatically by the script Abilint (TD).
120 !Do not modify the following lines by hand.
121 #undef ABI_FUNC
122 #define ABI_FUNC 'driver'
123  use interfaces_14_hidewrite
124  use interfaces_18_timing
125  use interfaces_32_util
126  use interfaces_41_xc_lowlevel
127  use interfaces_43_wvl_wrappers
128  use interfaces_95_drive, except_this_one => driver
129 !End of the abilint section
130 
131  implicit none
132 
133  !Arguments ------------------------------------
134  !scalars
135  integer,intent(in) :: ndtset,ndtset_alloc,npsp
136  real(dp),intent(in) :: cpui
137  character(len=6),intent(in) :: codvsn
138  character(len=fnlen),intent(in) :: filstat
139  type(MPI_type),intent(inout) :: mpi_enregs(0:ndtset_alloc)
140  !arrays
141  character(len=fnlen),intent(in) :: filnam(5)
142  type(dataset_type),intent(inout) :: dtsets(0:ndtset_alloc)
143  type(pspheader_type),intent(in) :: pspheads(npsp)
144  type(results_out_type),target,intent(inout) :: results_out(0:ndtset_alloc)
145 
146  !Local variables-------------------------------
147  !scalars
148 #if defined DEV_YP_VDWXC
149  character(len=fnlen) :: vdw_filnam
150 #endif
151  integer,parameter :: level=2,mdtset=9999
152  integer,save :: paw_size_old=-1
153  integer :: idtset,ierr,iexit,iget_cell,iget_occ,iget_vel,iget_xcart,iget_xred
154  integer :: ii,iimage,iimage_get,jdtset,jdtset_status,jj,kk
155  integer :: mtypalch,mu,mxnimage,nimage,openexit,paw_size,prtvol
156  real(dp) :: etotal
157  character(len=500) :: message
158  character(len=500) :: dilatmx_errmsg
159  logical :: converged,results_gathered,test_img,use_results_all
160  type(dataset_type) :: dtset
161  type(datafiles_type) :: dtfil
162  type(pawang_type) :: pawang
163  type(pseudopotential_type) :: psps
164  type(results_respfn_type) :: results_respfn
165  type(wvl_data) :: wvl
166 #if defined DEV_YP_VDWXC
167  type(xc_vdw_type) :: vdw_params
168 #endif
169  !arrays
170  integer :: mkmems(3)
171  integer,allocatable :: jdtset_(:),npwtot(:)
172  real(dp) :: acell(3),rprim(3,3),rprimd(3,3),tsec(2)
173  real(dp),allocatable :: acell_img(:,:),amu_img(:,:),rprim_img(:,:,:)
174  real(dp),allocatable :: fcart_img(:,:,:),fred_img(:,:,:)
175  real(dp),allocatable :: etotal_img(:),mixalch_img(:,:,:),strten_img(:,:),miximage(:,:)
176  real(dp),allocatable :: occ(:),xcart(:,:),xred(:,:),xredget(:,:)
177  real(dp),allocatable :: occ_img(:,:),vel_cell_img(:,:,:),vel_img(:,:,:),xred_img(:,:,:)
178  type(pawrad_type),allocatable :: pawrad(:)
179  type(pawtab_type),allocatable :: pawtab(:)
180  type(results_out_type),pointer :: results_out_all(:)
181 
182 !******************************************************************
183 
184  DBG_ENTER("COLL")
185 
186  call timab(640,1,tsec)
187  call timab(641,3,tsec)
188  call status(0,filstat,iexit,level,'enter         ')
189 
190 !Structured debugging if prtvol==-level
191  prtvol=dtsets(1)%prtvol
192  if(prtvol==-level)then
193    write(message,'(80a,a,a)')  ('=',ii=1,80),ch10,' driver : enter , debug mode '
194    call wrtout(std_out,message,'COLL')
195  end if
196 
197  if(ndtset>mdtset)then
198    write(message,'(a,i2,a,i5,a)')'  The maximal allowed ndtset is ',mdtset,' while the input value is ',ndtset,'.'
199    MSG_BUG(message)
200  end if
201 
202  mtypalch=dtsets(1)%ntypalch
203  do ii=1,ndtset_alloc
204    mtypalch=max(dtsets(ii)%ntypalch,mtypalch)
205  end do
206  call psps_init_global(mtypalch, npsp, psps, pspheads)
207 
208  ABI_ALLOCATE(jdtset_,(0:ndtset))
209  if(ndtset/=0)then
210    jdtset_(:)=dtsets(0:ndtset)%jdtset
211  else
212    jdtset_(0)=0
213  end if
214 
215  do idtset=1,ndtset_alloc
216    call init_results_respfn(dtsets,ndtset_alloc,results_respfn)
217  end do
218 
219  call timab(641,2,tsec)
220 
221 !*********************************************************************
222 !Big loop on datasets
223 
224 !Do loop on idtset (allocate statements are present)
225  do idtset=1,ndtset_alloc
226 
227    if(mpi_enregs(idtset)%me<0) cycle
228 
229    call timab(642,1,tsec)
230    call abi_io_redirect(new_io_comm=mpi_enregs(idtset)%comm_world)
231    call libpaw_write_comm_set(mpi_enregs(idtset)%comm_world)
232 
233    jdtset=dtsets(idtset)%jdtset ; if(ndtset==0)jdtset=1
234 
235    if(ndtset>=2)then
236      jdtset_status=jdtset
237    else
238      jdtset_status=0
239    end if
240 
241    call status(jdtset_status,filstat,iexit,level,'loop jdtset   ')
242 
243 !  Copy input variables into a local dtset.
244    call dtset_copy(dtset, dtsets(idtset))
245 
246 !  Set other values
247    dtset%jdtset = jdtset
248    dtset%ndtset = ndtset
249 
250 !  Print DATASET number
251    write(message,'(83a)') ch10,('=',mu=1,80),ch10,'== DATASET'
252    if (jdtset>=100) then
253      write(message,'(2a,i4,65a)') trim(message),' ',jdtset,' ',('=',mu=1,64)
254    else
255      write(message,'(2a,i2,67a)') trim(message),' ',jdtset,' ',('=',mu=1,66)
256    end if
257    write(message,'(3a,i5)') trim(message),ch10,'-   nproc =',mpi_enregs(idtset)%nproc
258    if (.not.mpi_distrib_is_ok(mpi_enregs(idtset),dtset%mband,dtset%nkpt,&
259 &   dtset%mkmem,dtset%nsppol)) then
260      write(message,'(2a)') trim(message),&
261 &     '   -> not optimal: autoparal keyword recommended in input file'
262    end if
263    write(message,'(3a)') trim(message),ch10,' '
264    call wrtout(ab_out,message,'COLL')
265    call wrtout(std_out,message,'PERS')     ! PERS is choosen to make debugging easier
266 
267    if ( dtset%np_slk == 0 ) then
268      call xgScalapack_config(SLK_DISABLED,500)
269    else if ( dtset%np_slk == 1000000 ) then
270      call xgScalapack_config(SLK_AUTO,500)
271    else if ( dtset%np_slk > 1 ) then
272      call xgScalapack_config(dtset%np_slk,500)
273    else
274      call xgScalapack_config(SLK_AUTO,500)
275    end if
276 
277 !  Copy input values
278    mkmems(1) = dtset%mkmem
279    mkmems(2) = dtset%mkqmem
280    mkmems(3) = dtset%mk1mem
281 
282 !  Initialize MPI data for the parallelism over images
283    nimage=mpi_enregs(idtset)%my_nimage
284 
285 !  Retrieve evolving arrays (most of them depend on images)
286    ABI_ALLOCATE(acell_img,(3,nimage))
287    ABI_ALLOCATE(amu_img,(dtset%ntypat,nimage))
288    ABI_ALLOCATE(mixalch_img,(dtset%npspalch,dtset%ntypalch,nimage))
289    ABI_ALLOCATE(occ_img,(dtset%mband*dtset%nkpt*dtset%nsppol,nimage))
290    ABI_ALLOCATE(rprim_img,(3,3,nimage))
291    ABI_ALLOCATE(vel_img,(3,dtset%natom,nimage))
292    ABI_ALLOCATE(vel_cell_img,(3,3,nimage))
293    ABI_ALLOCATE(xred_img,(3,dtset%natom,nimage))
294    do iimage=1,nimage
295      ii=mpi_enregs(idtset)%my_imgtab(iimage)
296      acell_img   (:  ,iimage) = dtset%acell_orig(:,ii)
297      amu_img   (:  ,iimage)   = dtset%amu_orig(1:dtset%ntypat,ii)
298      mixalch_img   (:,:,iimage) =dtset%mixalch_orig(1:dtset%npspalch,1:dtset%ntypalch,ii)
299      rprim_img   (:,:,iimage) = dtset%rprim_orig(:,:,ii)
300      vel_img     (:,:,iimage) = dtset%vel_orig(:,1:dtset%natom,ii)
301      vel_cell_img(:,:,iimage) = dtset%vel_cell_orig(:,:,ii)
302      xred_img (:,:,iimage) = dtset%xred_orig(:,1:dtset%natom,ii)
303 !    Note that occ is not supposed to depend on the image in the input file.
304 !    However, the results will depend on the image. And occ can be used
305 !    to reinitialize the next dataset, or the next timimage.
306      occ_img  (:  ,iimage) = dtset%occ_orig(1:dtset%mband*dtset%nkpt*dtset%nsppol)
307    end do
308 
309 !  ****************************************************************************
310 !  Treat the file names (get variables)
311 
312 !  In the case of multiple images, the file names will be overwritten later (for each image)
313    call dtfil_init(dtfil,dtset,filnam,filstat,idtset,jdtset_,mpi_enregs(idtset),ndtset)
314    if (dtset%optdriver==RUNL_GSTATE.and.dtset%nimage>1) then
315      call dtfil_init_img(dtfil,dtset,dtsets,idtset,jdtset_,ndtset,ndtset_alloc)
316    end if
317 
318 !  ****************************************************************************
319 !  Treat other get variables
320 
321 !  If multi dataset mode, and already the second dataset,
322 !  treatment of other get variables.
323    if( ndtset>1 .and. idtset>1 )then
324 
325 !    Check if parallelization over images is activated
326      mxnimage=maxval(dtsets(1:ndtset_alloc)%nimage)
327      ABI_ALLOCATE(miximage,(mxnimage,mxnimage))
328      test_img=(mxnimage/=1.and.maxval(dtsets(:)%npimage)>1)
329      use_results_all=.false.
330      if (test_img.and.mpi_enregs(idtset)%me_cell==0) then
331        use_results_all=.true.
332        ABI_DATATYPE_ALLOCATE(results_out_all,(0:ndtset_alloc))
333      else
334        results_out_all => results_out
335      end if
336      results_gathered=.false.
337 
338      call find_getdtset(dtsets,dtset%getocc,'getocc',idtset,iget_occ,miximage,mxnimage,ndtset_alloc)
339      if(iget_occ/=0)then
340 !      Gather contributions to results_out from images, if needed
341        if (test_img.and.mpi_enregs(idtset)%me_cell==0.and.(.not.results_gathered)) then
342          call gather_results_out(dtsets,mpi_enregs,results_out,results_out_all,use_results_all, &
343 &         allgather=.true.,only_one_per_img=.true.)
344          results_gathered=.true.
345        end if
346        if ((.not.test_img).or.mpi_enregs(idtset)%me_cell==0) then
347          do iimage=1,nimage
348            ii=mpi_enregs(idtset)%my_imgtab(iimage)
349            occ_img(:,iimage)=zero
350            do iimage_get=1,dtsets(iget_occ)%nimage
351              do jj=1,dtset%mband*dtset%nkpt*dtset%nsppol
352                occ_img(jj,iimage)=occ_img(jj,iimage)+miximage(ii,iimage_get) &
353 &               *results_out_all(iget_occ)%occ(jj,iimage_get)
354              end do
355            end do
356          end do
357        end if
358      end if
359 
360 !    Getcell has to be treated BEFORE getxcart since acell and rprim will be used
361      call find_getdtset(dtsets,dtset%getcell,'getcell',idtset,iget_cell,miximage,mxnimage,ndtset_alloc)
362      if(iget_cell/=0)then
363 !      Gather contributions to results_out from images, if needed
364        if (test_img.and.mpi_enregs(idtset)%me_cell==0.and.(.not.results_gathered)) then
365          call gather_results_out(dtsets,mpi_enregs,results_out,results_out_all,use_results_all, &
366 &         allgather=.true.,only_one_per_img=.true.)
367          results_gathered=.true.
368        end if
369        if ((.not.test_img).or.mpi_enregs(idtset)%me_cell==0) then
370          do iimage=1,nimage
371            ii=mpi_enregs(idtset)%my_imgtab(iimage)
372            acell_img(:,iimage)=zero
373            rprim_img(:,:,iimage)=zero
374            do iimage_get=1,dtsets(iget_cell)%nimage
375              do jj=1,3
376                acell_img(jj,iimage)=acell_img(jj,iimage)+&
377 &               miximage(ii,iimage_get)*results_out_all(iget_cell)%acell(jj,iimage_get)
378                do kk=1,3
379                  rprim_img(kk,jj,iimage)=rprim_img(kk,jj,iimage)+&
380 &                 miximage(ii,iimage_get)*results_out_all(iget_cell)%rprim(kk,jj,iimage_get)
381                end do
382              end do
383              dtset%rprimd_orig(:,:,iimage)=dtsets(iget_cell)%rprimd_orig(:,:,iimage)
384 !            Check that the new acell and rprim are consistent with the input dilatmx
385              call mkrdim(acell_img(:,iimage),rprim_img(:,:,iimage),rprimd)
386              call chkdilatmx(dtset%chkdilatmx,dtset%dilatmx,rprimd,&
387 &             dtset%rprimd_orig(1:3,1:3,iimage), dilatmx_errmsg)
388              _IBM6("dilatxm_errmsg: "//TRIM(dilatmx_errmsg))
389              if (LEN_TRIM(dilatmx_errmsg) /= 0) then
390                acell_img(1,iimage) = sqrt(sum(rprimd(:,1)**2))
391                acell_img(2,iimage) = sqrt(sum(rprimd(:,2)**2))
392                acell_img(3,iimage) = sqrt(sum(rprimd(:,3)**2))
393                rprim_img(:,1,iimage) = rprimd(:,1) / acell_img(1,iimage)
394                rprim_img(:,2,iimage) = rprimd(:,2) / acell_img(1,iimage)
395                rprim_img(:,3,iimage) = rprimd(:,3) / acell_img(1,iimage)
396                MSG_WARNING(dilatmx_errmsg)
397              end if
398            end do
399          end do
400        end if
401      end if
402 
403      call find_getdtset(dtsets,dtset%getxred,'getxred',idtset,iget_xred,miximage,mxnimage,ndtset_alloc)
404      if(iget_xred/=0)then
405 !      Gather contributions to results_out from images, if needed
406        if (test_img.and.mpi_enregs(idtset)%me_cell==0.and.(.not.results_gathered)) then
407          call gather_results_out(dtsets,mpi_enregs,results_out,results_out_all,use_results_all, &
408 &         allgather=.true.,only_one_per_img=.true.)
409          results_gathered=.true.
410        end if
411        if ((.not.test_img).or.mpi_enregs(idtset)%me_cell==0) then
412          do iimage=1,nimage
413            ii=mpi_enregs(idtset)%my_imgtab(iimage)
414            xred_img(:,:,iimage)=zero
415            do iimage_get=1,dtsets(iget_xred)%nimage
416              do jj=1,dtset%natom
417                do kk=1,3
418                  xred_img(kk,jj,iimage)=xred_img(kk,jj,iimage)+&
419 &                 miximage(ii,iimage_get)*results_out_all(iget_xred)%xred(kk,jj,iimage_get)
420                end do
421              end do
422            end do
423          end do
424        end if
425      end if
426 
427      call find_getdtset(dtsets,dtset%getxcart,'getxcart',idtset,iget_xcart,miximage,mxnimage,ndtset_alloc)
428      if(iget_xcart/=0)then
429 !      Gather contributions to results_out from images, if needed
430        if (test_img.and.mpi_enregs(idtset)%me_cell==0.and.(.not.results_gathered)) then
431          call gather_results_out(dtsets,mpi_enregs,results_out,results_out_all,use_results_all, &
432 &         allgather=.true.,only_one_per_img=.true.)
433          results_gathered=.true.
434        end if
435        if ((.not.test_img).or.mpi_enregs(idtset)%me_cell==0) then
436          ABI_ALLOCATE(xcart,(3,dtset%natom))
437          ABI_ALLOCATE(xredget,(3,dtset%natom))
438          do iimage=1,nimage
439            ii=mpi_enregs(idtset)%my_imgtab(iimage)
440            xred_img(:,:,iimage)=zero
441            do iimage_get=1,dtsets(iget_xcart)%nimage
442 !            Compute xcart of the previous dataset
443              call  mkrdim(results_out_all(iget_xcart)%acell(:,iimage_get),&
444 &             results_out_all(iget_xcart)%rprim(:,:,iimage_get),rprimd)
445              do jj=1,dtset%natom
446                do kk=1,3
447                  xredget (kk,jj)=results_out_all(iget_xcart)%xred(kk,jj,iimage_get)
448                end do
449              end do
450              call xred2xcart(dtset%natom,rprimd,xcart,xredget)
451 !            xcart from previous dataset is computed. Now, produce xred for the new dataset,
452 !            with the new acell and rprim ...
453              call mkrdim(acell_img(:,iimage),rprim_img(:,:,iimage),rprimd)
454              call xcart2xred(dtset%natom,rprimd,xcart,xredget(:,:))
455              do jj=1,dtset%natom
456                do kk=1,3
457                  xred_img(kk,jj,iimage)=xred_img(kk,jj,iimage)+miximage(ii,iimage_get)*xredget(kk,jj)
458                end do
459              end do
460            end do
461          end do
462          ABI_DEALLOCATE(xcart)
463          ABI_DEALLOCATE(xredget)
464        end if
465      end if
466 
467      call find_getdtset(dtsets,dtset%getvel,'getvel',idtset,iget_vel,miximage,mxnimage,ndtset_alloc)
468      if(iget_vel/=0)then
469 !      Gather contributions to results_out from images, if needed
470        if (test_img.and.mpi_enregs(idtset)%me_cell==0.and.(.not.results_gathered)) then
471          call gather_results_out(dtsets,mpi_enregs,results_out,results_out_all,use_results_all, &
472 &         allgather=.true.,only_one_per_img=.true.)
473          results_gathered=.true.
474        end if
475        if ((.not.test_img).or.mpi_enregs(idtset)%me_cell==0) then
476          do iimage=1,nimage
477            ii=mpi_enregs(idtset)%my_imgtab(iimage)
478            vel_img(:,:,iimage)=zero
479            vel_cell_img(:,:,iimage)=zero
480            do iimage_get=1,dtsets(iget_vel)%nimage
481              do jj=1,dtset%natom
482                do kk=1,3
483                  vel_img(kk,jj,iimage)=vel_img(kk,jj,iimage)+&
484 &                 miximage(ii,iimage_get)*results_out_all(iget_vel)%vel(kk,jj,iimage_get)
485                  vel_cell_img(kk,jj,iimage)=vel_cell_img(kk,jj,iimage)+&
486 &                 miximage(ii,iimage_get)*results_out_all(iget_vel)%vel_cell(kk,jj,iimage_get)
487                end do
488              end do
489            end do
490          end do
491        end if
492      end if
493 
494 !    In the case of parallelization over images, has to distribute data
495      if (test_img) then
496        if (iget_occ/=0) then
497          call xmpi_bcast(occ_img,0,mpi_enregs(idtset)%comm_cell,ierr)
498        end if
499        if (iget_cell/=0) then
500          call xmpi_bcast(acell_img,0,mpi_enregs(idtset)%comm_cell,ierr)
501          call xmpi_bcast(rprim_img,0,mpi_enregs(idtset)%comm_cell,ierr)
502        end if
503        if (iget_vel/=0) then
504          call xmpi_bcast(vel_img,0,mpi_enregs(idtset)%comm_cell,ierr)
505          call xmpi_bcast(vel_cell_img,0,mpi_enregs(idtset)%comm_cell,ierr)
506        end if
507        if (iget_xred/=0.or.iget_xcart/=0) then
508          call xmpi_bcast(xred_img,0,mpi_enregs(idtset)%comm_cell,ierr)
509        end if
510      end if
511 
512 !    Clean memory
513      ABI_DEALLOCATE(miximage)
514      if (test_img.and.mpi_enregs(idtset)%me_cell==0) then
515        if (results_gathered) then
516          call destroy_results_out(results_out_all)
517        end if
518        ABI_DATATYPE_DEALLOCATE(results_out_all)
519      else
520        nullify(results_out_all)
521      end if
522 
523    end if
524 
525 !  ****************************************************************************
526 !  Treat the pseudopotentials : initialize the psps/PAW variable
527 
528    call status(jdtset_status,filstat,iexit,level,'init psps     ')
529    call psps_init_from_dtset(dtset, idtset, psps, pspheads)
530 
531 !  The correct dimension of pawrad/tab is ntypat. In case of alchemical psps
532 !  pawrad/tab(ipsp) is invoked with ipsp<=npsp. So, in order to avoid any problem,
533 !  declare pawrad/tab at paw_size=max(ntypat,npsp).
534    paw_size=0;if (psps%usepaw==1) paw_size=max(dtset%ntypat,dtset%npsp)
535    if (paw_size/=paw_size_old) then
536      if (paw_size_old/=-1) then
537        call pawrad_free(pawrad)
538        call pawtab_free(pawtab)
539        ABI_DATATYPE_DEALLOCATE(pawrad)
540        ABI_DATATYPE_DEALLOCATE(pawtab)
541      end if
542      ABI_DATATYPE_ALLOCATE(pawrad,(paw_size))
543      ABI_DATATYPE_ALLOCATE(pawtab,(paw_size))
544      call pawtab_nullify(pawtab)
545      paw_size_old=paw_size
546    end if
547 
548 !  ****************************************************************************
549 !  WVL allocations.
550 
551 !  Nullify wvl_data. It is important to do so irregardless of the value of usewvl:
552    call nullify_wvl_data(wvl)
553 
554 !  Set up mpi informations from the dataset
555    if (dtset%usewvl == 1) then
556 #if defined HAVE_BIGDFT
557      call f_malloc_set_status(iproc=mpi_enregs(idtset)%me_wvl)
558      call mpi_environment_set(bigdft_mpi,mpi_enregs(idtset)%me_wvl,&
559 &     mpi_enregs(idtset)%nproc_wvl,mpi_enregs(idtset)%comm_wvl,&
560 &     mpi_enregs(idtset)%nproc_wvl)
561 #endif
562    end if
563 
564 !  ****************************************************************************
565 !  At this stage, all the data needed for the treatment of one dataset
566 !  have been transferred from multi-dataset arrays.
567 
568    iexit=0
569 
570 !  Smaller integer arrays
571    etotal=zero
572    ABI_ALLOCATE(npwtot,(dtset%nkpt))
573    npwtot = 0
574    ABI_ALLOCATE(xred,(3,dtset%natom))
575    ABI_ALLOCATE(occ,(dtset%mband*dtset%nkpt*dtset%nsppol))
576    if(dtset%optdriver/=RUNL_GSTATE)then
577      occ(:)=occ_img(:,1)
578      acell(:)=acell_img(:,1)
579      rprim(:,:)=rprim_img(:,:,1)
580      xred(:,:)=xred_img(:,:,1)
581    end if
582 
583 !  ****************************************************************************
584 !  Exchange-correlation
585 
586    call echo_xc_name (dtset%ixc)
587 
588    if (dtset%ixc<0) then
589      call libxc_functionals_init(dtset%ixc,dtset%nspden)
590 
591 #if defined DEV_YP_VDWXC
592      if ( (dtset%vdw_xc > 0) .and. (dtset%vdw_xc < 3) ) then
593        vdw_params%functional = dtset%vdw_xc
594        vdw_params%acutmin = dtset%vdw_df_acutmin
595        vdw_params%aratio = dtset%vdw_df_aratio
596        vdw_params%damax = dtset%vdw_df_damax
597        vdw_params%damin = dtset%vdw_df_damin
598        vdw_params%dcut = dtset%vdw_df_dcut
599        vdw_params%dratio = dtset%vdw_df_dratio
600        vdw_params%dsoft = dtset%vdw_df_dsoft
601        vdw_params%gcut = dtset%vdw_df_gcut
602        vdw_params%ndpts = dtset%vdw_df_ndpts
603        vdw_params%ngpts = dtset%vdw_df_ngpts
604        vdw_params%nqpts = dtset%vdw_df_nqpts
605        vdw_params%nrpts = dtset%vdw_df_nrpts
606        vdw_params%nsmooth = dtset%vdw_df_nsmooth
607        vdw_params%phisoft = dtset%vdw_df_phisoft
608        vdw_params%qcut = dtset%vdw_df_qcut
609        vdw_params%qratio = dtset%vdw_df_qratio
610        vdw_params%rcut = dtset%vdw_df_rcut
611        vdw_params%rsoft = dtset%vdw_df_rsoft
612        vdw_params%tolerance = dtset%vdw_df_tolerance
613        vdw_params%tweaks = dtset%vdw_df_tweaks
614        vdw_params%zab = dtset%vdw_df_zab
615        write(message,'(a,1x,a)') ch10,'[vdW-DF] *** Before init ***'
616        call wrtout(std_out,message,'COLL')
617        call xc_vdw_show(std_out,vdw_params)
618        if ( dtset%irdvdw == 1 ) then
619          write(vdw_filnam,'(a,a)') trim(filnam(3)),'_VDW.nc'
620          call xc_vdw_read(vdw_filnam)
621        else
622          call xc_vdw_init(vdw_params)
623        end if
624        call xc_vdw_libxc_init(vdw_params%functional)
625        write(message,'(a,1x,a)') ch10,'[vdW-DF] *** After init ***'
626        call wrtout(std_out,message,'COLL')
627 !      call xc_vdw_get_params(vdw_params)
628        call xc_vdw_memcheck(std_out)
629        call xc_vdw_show(std_out)
630        call xc_vdw_show(ab_out)
631 
632        write (message,'(a,1x,a,e10.3,a)')ch10,&
633 &       '[vdW-DF] activation threshold: vdw_df_threshold=',dtset%vdw_df_threshold,ch10
634        call xc_vdw_trigger(.false.)
635        call wrtout(std_out,message,'COLL')
636      end if
637 #else
638      if ( (dtset%vdw_xc > 0) .and. (dtset%vdw_xc < 3) ) then
639        write(message,'(3a)')&
640 &       'vdW-DF functionals are not fully operational yet.',ch10,&
641 &       'Action : modify vdw_xc'
642        MSG_ERROR(message)
643      end if
644 #endif
645    end if
646 
647 !  FFTW3 threads initialization
648    if (dtset%ngfft(7)/100==FFT_FFTW3)then
649      call fftw3_init_threads()
650    end if
651 
652 !  linalg initialisation:
653    call abi_linalg_init(mpi_enregs(idtset)%comm_bandspinorfft,dtset%np_slk,3*maxval(dtset%nband(:)), mpi_enregs(idtset)%me)
654 
655    call timab(642,2,tsec)
656 
657    select case(dtset%optdriver)
658 
659    case(RUNL_GSTATE)
660 
661      ABI_ALLOCATE(fcart_img,(3,dtset%natom,nimage))
662      ABI_ALLOCATE(fred_img,(3,dtset%natom,nimage))
663      ABI_ALLOCATE(etotal_img,(nimage))
664      ABI_ALLOCATE(strten_img,(6,nimage))
665 
666      call status(jdtset_status,filstat,iexit,level,'call gstateimg')
667      call gstateimg(acell_img,amu_img,codvsn,cpui,dtfil,dtset,etotal_img,fcart_img,&
668 &     fred_img,iexit,mixalch_img,mpi_enregs(idtset),nimage,npwtot,occ_img,&
669 &     pawang,pawrad,pawtab,psps,rprim_img,strten_img,vel_cell_img,vel_img,wvl,xred_img,&
670 &     filnam,filstat,idtset,jdtset_,ndtset)
671 
672    case(RUNL_RESPFN)
673 
674      call status(jdtset_status,filstat,iexit,level,'call respfn')
675      call respfn(codvsn,cpui,dtfil,dtset,etotal,iexit,mkmems,mpi_enregs(idtset),&
676 &     npwtot,occ,pawang,pawrad,pawtab,psps,results_respfn,xred)
677 
678    case(RUNL_SCREENING)
679 
680      call status(jdtset_status,filstat,iexit,level,'call screening')
681      call screening(acell,codvsn,dtfil,dtset,pawang,pawrad,pawtab,psps,rprim)
682 
683    case(RUNL_SIGMA)
684 
685      call status(jdtset_status,filstat,iexit,level,'call sigma')
686      call sigma(acell,codvsn,dtfil,dtset,pawang,pawrad,pawtab,psps,rprim,converged)
687 
688    case(RUNL_NONLINEAR)
689 
690      call status(jdtset_status,filstat,iexit,level,'call nonlinear')
691      call nonlinear(codvsn,dtfil,dtset,etotal,iexit,&
692 &     dtset%mband,dtset%mgfft,dtset%mkmem,mpi_enregs(idtset),dtset%mpw,dtset%natom,dtset%nfft,dtset%nkpt,npwtot,dtset%nspden,&
693 &     dtset%nspinor,dtset%nsppol,dtset%nsym,occ,pawrad,pawtab,psps,xred)
694 
695    case(6)
696 
697      write(message,'(3a)')&
698 &     'The optdriver value 6 has been disabled since ABINITv6.0.',ch10,&
699 &     'Action : modify optdriver in the input file.'
700      MSG_ERROR(message)
701 
702    case (RUNL_BSE)
703 
704      call status(jdtset_status,filstat,iexit,level,'call bethe_salpeter')
705      call bethe_salpeter(acell,codvsn,dtfil,dtset,pawang,pawrad,pawtab,psps,rprim,xred)
706 
707    case(RUNL_GWLS) ! For running G0W0 calculations with Lanczos basis for dielectric operator and Sternheimer equation for avoiding the use of conduction states (MC+JJL)
708 
709      ABI_ALLOCATE(fcart_img,(3,dtset%natom,nimage))
710      ABI_ALLOCATE(fred_img,(3,dtset%natom,nimage))
711      ABI_ALLOCATE(etotal_img,(nimage))
712      ABI_ALLOCATE(strten_img,(6,nimage))
713 
714      call status(jdtset_status,filstat,iexit,level,'call gwls_sternheimer')
715      call gwls_sternheimer(acell_img,amu_img,codvsn,cpui,dtfil,dtset,etotal_img,fcart_img,&
716 &     fred_img,iexit,mixalch_img,mpi_enregs(idtset),nimage,npwtot,occ_img,&
717 &     pawang,pawrad,pawtab,psps,rprim_img,strten_img,vel_cell_img,vel_img,xred_img,&
718 &     filnam,filstat,idtset,jdtset_,ndtset)
719 
720      ABI_DEALLOCATE(etotal_img)
721      ABI_DEALLOCATE(fcart_img)
722      ABI_DEALLOCATE(fred_img)
723      ABI_DEALLOCATE(strten_img)
724 
725    case (RUNL_WFK)
726      call status(jdtset_status,filstat,iexit,level,'call wfk_analyze')
727      call wfk_analyze(acell,codvsn,dtfil,dtset,pawang,pawrad,pawtab,psps,rprim,xred)
728 
729    case (RUNL_EPH)
730      call status(jdtset_status,filstat,iexit,level,'call eph      ')
731      call eph(acell,codvsn,dtfil,dtset,pawang,pawrad,pawtab,psps,rprim,xred)
732 
733    case default ! Bad value for optdriver
734      write(message,'(a,i0,4a)')&
735 &     'Unknown value for the variable optdriver: ',dtset%optdriver,ch10,&
736 &     'This is not allowed. ',ch10,&
737 &     'Action: modify optdriver in the input file.'
738      MSG_ERROR(message)
739    end select
740 
741    call timab(643,1,tsec)
742 !  ****************************************************************************
743 
744 !  Transfer of multi dataset outputs from temporaries :
745 !  acell, xred, occ rprim, and vel might be modified from their
746 !  input values
747 !  etotal, fcart, fred, and strten have been computed
748 !  npwtot was already computed before, but is stored only now
749 
750    if(dtset%optdriver==RUNL_GSTATE)then
751      do iimage=1,nimage
752        results_out(idtset)%etotal(iimage)                 =etotal_img(iimage)
753        results_out(idtset)%acell(:,iimage)                =acell_img(:,iimage)
754        results_out(idtset)%amu(1:dtset%ntypat,iimage)     =amu_img(:,iimage)
755        results_out(idtset)%rprim(:,:,iimage)              =rprim_img(:,:,iimage)
756        results_out(idtset)%strten(:,iimage)                =strten_img(:,iimage)
757        results_out(idtset)%fcart(1:3,1:dtset%natom,iimage)=fcart_img(:,:,iimage)
758        results_out(idtset)%fred(1:3,1:dtset%natom,iimage) =fred_img(:,:,iimage)
759        results_out(idtset)%mixalch(1:dtset%npspalch,1:dtset%ntypalch,iimage) &
760 &       =mixalch_img(1:dtset%npspalch,1:dtset%ntypalch,iimage)
761        results_out(idtset)%npwtot(1:dtset%nkpt,iimage)    =npwtot(1:dtset%nkpt)
762        results_out(idtset)%occ(1:dtset%mband*dtset%nkpt*dtset%nsppol,iimage)=&
763 &       occ_img(1:dtset%mband*dtset%nkpt*dtset%nsppol,iimage)
764        results_out(idtset)%vel(:,1:dtset%natom,iimage)    =vel_img(:,:,iimage)
765        results_out(idtset)%vel_cell(:,:,iimage)           =vel_cell_img(:,:,iimage)
766        results_out(idtset)%xred(:,1:dtset%natom,iimage)   =xred_img(:,:,iimage)
767      end do
768      ABI_DEALLOCATE(etotal_img)
769      ABI_DEALLOCATE(fcart_img)
770      ABI_DEALLOCATE(fred_img)
771      ABI_DEALLOCATE(strten_img)
772    else
773      results_out(idtset)%acell(:,1)                =acell(:)
774      results_out(idtset)%etotal(1)                 =etotal
775      results_out(idtset)%rprim(:,:,1)              =rprim(:,:)
776      results_out(idtset)%npwtot(1:dtset%nkpt,1)    =npwtot(1:dtset%nkpt)
777      results_out(idtset)%occ(1:dtset%mband*dtset%nkpt*dtset%nsppol,1)=&
778 &     occ(1:dtset%mband*dtset%nkpt*dtset%nsppol)
779      results_out(idtset)%xred(:,1:dtset%natom,1)   =xred(:,:)
780    end if
781    ABI_DEALLOCATE(acell_img)
782    ABI_DEALLOCATE(amu_img)
783    ABI_DEALLOCATE(mixalch_img)
784    ABI_DEALLOCATE(occ_img)
785    ABI_DEALLOCATE(rprim_img)
786    ABI_DEALLOCATE(vel_img)
787    ABI_DEALLOCATE(vel_cell_img)
788    ABI_DEALLOCATE(xred_img)
789 
790    if (dtset%ngfft(7)/100==FFT_FFTW3) then
791      call fftw3_cleanup()
792    end if
793 
794    if (dtset%ixc<0) then
795      call libxc_functionals_end()
796 
797 #if defined DEV_YP_VDWXC
798      if ( (dtset%vdw_xc > 0) .and. (dtset%vdw_xc < 10) ) then
799        if ( dtset%prtvdw /= 0 ) then
800          write(vdw_filnam,'(a,a)') trim(filnam(4)),'_VDW.nc'
801          call xc_vdw_write(vdw_filnam)
802        end if
803        call xc_vdw_show(std_out,vdw_params)
804        call xc_vdw_done(vdw_params)
805      end if
806 #endif
807    end if
808 
809 !  MG: There are routines such as GW and Berry phase that can change/compute
810 !  entries in the datatype at run-time. These changes won't be visibile
811 !  in the main output file since we are passing a copy of dtsets.
812 !  I tried to update the results with the call below but this creates
813 !  several problems in outvars since one should take into account
814 !  the new dimensions (e.g. nkptgw) and their maximum value.
815 !  For the time being, we continue to pass a copy of dtsets(idtset).
816 #if 0
817    call dtset_free(dtsets(idtset))
818    call dtset_copy(dtsets(idtset),dtset)
819 #endif
820 
821    call dtset_free(dtset)
822 
823    ABI_DEALLOCATE(occ)
824    ABI_DEALLOCATE(xred)
825    ABI_DEALLOCATE(npwtot)
826 
827    call abi_linalg_finalize()
828    call xg_finalize()
829 
830    ! Check whether exiting was required by the user.
831    ! If found then beat a hasty exit from time steps
832    openexit=1; if(dtset%chkexit==0) openexit=0
833    call exit_check(zero,dtfil%filnam_ds(1),iexit,ab_out,mpi_enregs(idtset)%comm_cell,openexit)
834 
835    call timab(643,2,tsec)
836 
837    if (iexit/=0) exit
838 
839  end do ! idtset (allocate statements are present - an exit statement is present)
840 
841 !*********************************************************************
842 
843  call timab(644,1,tsec)
844 
845 !PSP deallocation
846  call psps_free(psps)
847 
848 !XG 121126 : One should not use dtset or idtset in this section, as these might not be defined for all processors.
849 
850 !PAW deallocation
851  if (allocated(pawrad)) then
852    call pawrad_free(pawrad)
853    ABI_DATATYPE_DEALLOCATE(pawrad)
854  end if
855  if (allocated(pawtab)) then
856    call pawtab_free(pawtab)
857    ABI_DATATYPE_DEALLOCATE(pawtab)
858  end if
859  call pawang_free(pawang)
860 
861  ABI_DEALLOCATE(jdtset_)
862 
863 !Results_respfn deallocation
864  call destroy_results_respfn(results_respfn)
865 
866 #if defined HAVE_BIGDFT
867 !XG 121126 : NOTE that the next debugging section was quite problematic : indeed we are
868 !outside the loop over datasets, so the value of dtset%usewvl, that is referring to the last
869 !dtset, might not be initialized, if the dataset is NOT treated by the processor.
870 !See line 253 for the initialisation of dtset, while this might not happen if
871 !if(mpi_enregs(idtset)%me<0) cycle  , as defined by line 216
872 !if (dtset%usewvl == 1) then
873 !!  WVL - debugging
874 !call memocc_abi(0,mpi_enregs(1)%me_wvl,'count','stop')
875 !call wvl_timing(mpi_enregs(1)%me_wvl,'WFN_OPT','PR')
876 !call wvl_timing(mpi_enregs(1)%me_wvl,'              ','RE')
877 !end if
878 #endif
879 
880  call status(0,filstat,iexit,level,'exit')
881 
882  call timab(644,2,tsec)
883  call timab(640,2,tsec)
884 
885  DBG_EXIT("COLL")
886 
887 end subroutine driver