TABLE OF CONTENTS


ABINIT/m_lobpcgwf [ Functions ]

[ Top ] [ Functions ]

NAME

 m_lobpcgwf

FUNCTION

 this routine updates the whole wave functions at a given k-point,
 using the lobpcg method
 for a given spin-polarization, from a fixed hamiltonian
 but might also simply compute eigenvectors and eigenvalues at this k point.
 it will also update the matrix elements of the hamiltonian.

COPYRIGHT

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

PARENTS

      vtowfk

CHILDREN

SOURCE

 27 #if defined HAVE_CONFIG_H
 28 #include "config.h"
 29 #endif
 30 
 31 #include "abi_common.h"
 32 
 33 module m_lobpcgwf
 34 
 35  use defs_abitypes
 36  use defs_basis
 37  use m_abicore
 38  use m_lobpcg
 39  use m_xmpi
 40  use m_errors
 41  use m_time
 42  use m_xomp
 43  use m_fstrings
 44  use m_xg
 45  use m_lobpcg2
 46 
 47  use m_hamiltonian, only : gs_hamiltonian_type
 48  use m_pawcprj,     only : pawcprj_type
 49  use m_nonlop,      only : nonlop
 50  use m_prep_kgb,    only : prep_getghc, prep_nonlop
 51  use m_getghc,      only : multithreaded_getghc
 52 
 53  private
 54 
 55  integer, parameter :: l_tim_getghc=5
 56  double precision, parameter :: inv_sqrt2 = 1/sqrt2
 57  ! Use in getghc_gsc
 58  integer,save  :: l_cpopt
 59  integer,save  :: l_icplx
 60  integer,save  :: l_istwf
 61  integer,save  :: l_npw
 62  integer,save  :: l_nspinor
 63  logical,save  :: l_paw
 64  integer,save  :: l_prtvol
 65  integer,save  :: l_sij_opt
 66  real(dp), allocatable,save :: l_gvnlc(:,:)
 67  real(dp), allocatable,save ::  l_pcon(:)
 68  type(mpi_type),pointer,save :: l_mpi_enreg
 69  type(gs_hamiltonian_type),pointer,save :: l_gs_hamk
 70 
 71  public :: lobpcgwf2
 72 
 73  contains
 74 
 75 subroutine lobpcgwf2(cg,dtset,eig,enl_out,gs_hamk,kinpw,mpi_enreg,&
 76 &                   nband,npw,nspinor,prtvol,resid)
 77 
 78 
 79  use m_cgtools, only : dotprod_g
 80  use iso_c_binding
 81 
 82 !This section has been created automatically by the script Abilint (TD).
 83 !Do not modify the following lines by hand.
 84 #undef ABI_FUNC
 85 #define ABI_FUNC 'lobpcgwf2'
 86 !End of the abilint section
 87 
 88  implicit none
 89 
 90 !Arguments ------------------------------------
 91  integer,intent(in) :: nband,npw,prtvol,nspinor
 92  type(gs_hamiltonian_type),target,intent(inout) :: gs_hamk
 93  type(dataset_type)              ,intent(in   ) :: dtset
 94  type(mpi_type)           ,target,intent(inout) :: mpi_enreg
 95  real(dp)                 ,target,intent(inout) :: cg(2,nspinor*nband*npw)!,gsc(2,nspinor*nband*npw)
 96  real(dp)                        ,intent(in   ) :: kinpw(npw)
 97  real(dp)                 ,target,intent(  out) :: resid(nband)
 98  real(dp)                        ,intent(  out) :: enl_out(nband)
 99  real(dp)                 ,target,intent(  out) :: eig(nband)
100 
101 !Local variables-------------------------------
102 
103  type(xgBlock_t) :: xgx0
104  type(xgBlock_t) :: xgeigen
105  type(xgBlock_t) :: xgresidu
106  type(lobpcg_t) :: lobpcg
107 
108  integer :: space, blockdim,  nline
109  integer :: ipw
110 
111  integer :: nthreads
112 
113  double precision :: lobpcgMem(2)
114  double precision :: localmem
115 
116  integer, parameter :: tim_lobpcgwf2 = 1650
117  double precision :: cputime, walltime
118  double precision :: tsec(2)
119 
120  type(c_ptr) :: cptr
121  real(dp), pointer :: eig_ptr(:,:) => NULL()
122  real(dp), pointer :: resid_ptr(:,:) => NULL()
123 
124  ! Stupid things for NC
125  integer,parameter :: choice=1, paw_opt=0, signs=1
126  type(pawcprj_type) :: cprj_dum(gs_hamk%natom,0)
127  integer :: iband, shift
128  real(dp) :: gsc_dummy(0,0)
129 
130 ! *********************************************************************
131 
132 
133 !###########################################################################
134 !################ INITIALISATION  ##########################################
135 !###########################################################################
136 
137  call timab(tim_lobpcgwf2,1,tsec)
138  cputime = abi_cpu_time()
139  walltime = abi_wtime()
140 
141  ! Set module variables
142  l_paw = (gs_hamk%usepaw==1)
143  l_cpopt=-1;l_sij_opt=0;if (l_paw) l_sij_opt=1
144  l_istwf=gs_hamk%istwf_k
145  l_npw = npw
146  l_nspinor = nspinor
147  l_prtvol = prtvol
148  l_mpi_enreg => mpi_enreg
149  l_gs_hamk => gs_hamk
150 
151 !Variables
152  nline=dtset%nline
153  blockdim=l_mpi_enreg%nproc_band*l_mpi_enreg%bandpp
154 
155 !Depends on istwfk
156  if ( l_istwf == 2 ) then ! Real only
157    ! SPACE_CR mean that we have complex numbers but no re*im terms only re*re
158    ! and im*im so that a vector of complex is consider as a long vector of real
159    ! therefore the number of data is (2*npw*nspinor)*nband
160    ! This space is completely equivalent to SPACE_R but will correctly set and
161    ! get the array data into the xgBlock
162    space = SPACE_CR
163    l_icplx = 2
164 
165  else ! complex
166    space = SPACE_C
167    l_icplx = 1
168  end if
169 
170  ! Memory info
171  if ( prtvol >= 3 ) then
172    lobpcgMem = lobpcg_memInfo(nband,l_icplx*l_npw*l_nspinor,blockdim,space)
173    localMem = (l_npw+2*l_npw*l_nspinor*blockdim+2*nband)*kind(1.d0)
174    write(std_out,'(1x,A,F10.6,1x,A)') "Each MPI process calling lobpcg should need around ", &
175    (localMem+sum(lobpcgMem))/1e9, &
176    "GB of peak memory as follows :"
177    write(std_out,'(4x,A,F10.6,1x,A)') "Permanent memory in lobpcgwf : ", &
178    (localMem)/1e9, "GB"
179    write(std_out,'(4x,A,F10.6,1x,A)') "Permanent memory in m_lobpcg : ", &
180    (lobpcgMem(1))/1e9, "GB"
181    write(std_out,'(4x,A,F10.6,1x,A)') "Temporary memory in m_lobpcg : ", &
182    (lobpcgMem(2))/1e9, "GB"
183  end if
184 
185  !For preconditionning
186  ABI_MALLOC(l_pcon,(1:l_icplx*npw))
187  !$omp parallel do schedule(static), shared(l_pcon,kinpw)
188  do ipw=1-1,l_icplx*npw-1
189    if(kinpw(ipw/l_icplx+1)>huge(0.0_dp)*1.d-11) then
190      l_pcon(ipw+1)=0.d0
191    else
192      l_pcon(ipw+1) = (27+kinpw(ipw/l_icplx+1)*(18+kinpw(ipw/l_icplx+1)*(12+8*kinpw(ipw/l_icplx+1)))) &
193 &    / (27+kinpw(ipw/l_icplx+1)*(18+kinpw(ipw/l_icplx+1)*(12+8*kinpw(ipw/l_icplx+1))) + 16*kinpw(ipw/l_icplx+1)**4)
194    end if
195  end do
196 
197  ! Local variables for lobpcg
198  !call xg_init(xgx0,space,icplx*npw*nspinor,nband)
199  call xgBlock_map(xgx0,cg,space,l_icplx*l_npw*l_nspinor,nband,l_mpi_enreg%comm_bandspinorfft)
200  if ( l_istwf == 2 ) then ! Real only
201    ! Scale cg
202    call xgBlock_scale(xgx0,sqrt2,1)
203    ! This is possible since the memory in cg and xgx0 is the same
204    ! Don't know yet how to deal with this with xgBlock
205    if(l_mpi_enreg%me_g0 == 1) cg(:, 1:npw*nspinor*nband:npw) = cg(:, 1:npw*nspinor*nband:npw) * inv_sqrt2
206  end if
207 
208  !call xg_init(xgeigen,SPACE_R,nband,1,l_mpi_enreg%comm_bandspinorfft)
209  ! Trick the with C to change rank of arrays (:) to (:,:)
210  cptr = c_loc(eig)
211  call c_f_pointer(cptr,eig_ptr,(/ nband,1 /))
212  call xgBlock_map(xgeigen,eig_ptr,SPACE_R,nband,1)
213 
214  !call xg_init(xgresidu,SPACE_R,nband,1,l_mpi_enreg%comm_bandspinorfft)
215  ! Trick the with C to change rank of arrays (:) to (:,:)
216  cptr = c_loc(resid)
217  call c_f_pointer(cptr,resid_ptr,(/ nband,1 /))
218  call xgBlock_map(xgresidu,resid_ptr,SPACE_R,nband,1)
219 
220  ABI_MALLOC(l_gvnlc,(2,l_npw*l_nspinor*blockdim))
221 
222  call lobpcg_init(lobpcg,nband, l_icplx*l_npw*l_nspinor, blockdim,dtset%tolwfr,nline,space, l_mpi_enreg%comm_bandspinorfft)
223 
224 !###########################################################################
225 !################    RUUUUUUUN    ##########################################
226 !###########################################################################
227 
228  ! Run lobpcg
229  call lobpcg_run(lobpcg,xgx0,getghc_gsc,precond,xgeigen,xgresidu,prtvol)
230 
231  ! Free preconditionning since not needed anymore
232  ABI_FREE(l_pcon)
233 
234  ! Scale back
235  if(l_istwf == 2) then
236    call xgBlock_scale(xgx0,inv_sqrt2,1)
237    if(l_mpi_enreg%me_g0 == 1) cg(:, 1:npw*nspinor*nband:npw) = cg(:, 1:npw*nspinor*nband:npw) * sqrt2
238  end if
239 
240  ! Compute enlout (nonlocal energy for each band if necessary) This is the best
241  ! quick and dirty trick to compute this part in NC. gvnlc cannot be part of
242  ! lobpcg algorithm
243  if ( .not. l_paw ) then
244    !Check l_gvnlc size
245    !if ( size(l_gvnlc) < 2*nband*l_npw*l_nspinor ) then
246    if ( size(l_gvnlc) /= 0 ) then
247      ABI_FREE(l_gvnlc)
248      !ABI_MALLOC(l_gvnlc,(2,nband*l_npw*l_nspinor))
249      ABI_MALLOC(l_gvnlc,(0,0))
250    end if
251    !Call nonlop
252    if (mpi_enreg%paral_kgb==0) then
253 
254      call nonlop(choice,l_cpopt,cprj_dum,enl_out,l_gs_hamk,0,eig,mpi_enreg,nband,1,paw_opt,&
255 &                signs,gsc_dummy,l_tim_getghc,cg,l_gvnlc)
256 
257    else
258      do iband=1,nband/blockdim
259        shift = (iband-1)*blockdim*l_npw*l_nspinor
260       call prep_nonlop(choice,l_cpopt,cprj_dum, &
261 &       enl_out((iband-1)*blockdim+1:iband*blockdim),l_gs_hamk,0,&
262 &       eig((iband-1)*blockdim+1:iband*blockdim),blockdim,mpi_enreg,1,paw_opt,signs,&
263 &       gsc_dummy,l_tim_getghc, &
264 &       cg(:,shift+1:shift+blockdim*l_npw*l_nspinor),&
265 !&       l_gvnlc(:,shift+1:shift+blockdim*l_npw*l_nspinor),&
266 &       l_gvnlc(:,:),&
267 &       already_transposed=.false.)
268      end do
269    end if
270    !Compute enlout
271 !   do iband=1,nband
272 !     shift = npw*nspinor*(iband-1)
273 !       call dotprod_g(enl_out(iband),dprod_i,l_gs_hamk%istwf_k,npw*nspinor,1,cg(:, shift+1:shift+npw*nspinor),&
274 !  &     l_gvnlc(:, shift+1:shift+npw*nspinor),mpi_enreg%me_g0,mpi_enreg%comm_bandspinorfft)
275 !   end do
276  end if
277 
278  ABI_FREE(l_gvnlc)
279 
280  ! Free lobpcg
281  call lobpcg_free(lobpcg)
282 
283 !###########################################################################
284 !################    SORRY IT'S ALREADY FINISHED : )  ######################
285 !###########################################################################
286 
287 
288  call timab(tim_lobpcgwf2,2,tsec)
289  cputime = abi_cpu_time() - cputime
290  walltime = abi_wtime() - walltime
291  nthreads = xomp_get_num_threads(open_parallel = .true.)
292  if ( cputime/walltime/dble(nthreads) < 0.75 .and. (int(cputime/walltime)+1) /= nthreads) then
293    if ( prtvol >= 3 ) then
294      write(std_out,'(a)',advance='no') sjoin(" Lobpcg took", sec2str(cputime), "of cpu time")
295      write(std_out,*) sjoin("for a wall time of", sec2str(walltime))
296      write(std_out,'(a,f6.2)') " -> Ratio of ", cputime/walltime
297    end if
298    MSG_COMMENT(sjoin("You should set the number of threads to something close to",itoa(int(cputime/walltime)+1)))
299  end if
300 
301 
302  DBG_EXIT("COLL")
303 
304 end subroutine lobpcgwf2
305 
306 
307  subroutine getghc_gsc(X,AX,BX)
308    use m_xg, only : xg_t, xgBlock_get, xgBlock_set, xgBlock_getSize, xgBlock_t
309 #ifdef HAVE_OPENMP
310    use omp_lib
311 #endif
312 
313 !This section has been created automatically by the script Abilint (TD).
314 !Do not modify the following lines by hand.
315 #undef ABI_FUNC
316 #define ABI_FUNC 'getghc_gsc'
317 !End of the abilint section
318 
319   type(xgBlock_t), intent(inout) :: X
320   type(xgBlock_t), intent(inout) :: AX
321   type(xgBlock_t), intent(inout) :: BX
322   integer         :: blockdim
323   integer         :: spacedim
324   type(pawcprj_type) :: cprj_dum(l_gs_hamk%natom,0)
325   double precision :: dum
326   double precision, parameter :: inv_sqrt2 = 1/sqrt2
327   double precision, pointer :: cg(:,:)
328   double precision, pointer :: ghc(:,:)
329   double precision, pointer :: gsc(:,:)
330 
331   call xgBlock_getSize(X,spacedim,blockdim)
332   spacedim = spacedim/l_icplx
333 
334 
335   !call xgBlock_get(X,cg(:,1:blockdim*spacedim),0,spacedim)
336   call xgBlock_reverseMap(X,cg,l_icplx,spacedim*blockdim)
337   call xgBlock_reverseMap(AX,ghc,l_icplx,spacedim*blockdim)
338   call xgBlock_reverseMap(BX,gsc,l_icplx,spacedim*blockdim)
339 
340   ! scale back cg
341  if(l_istwf == 2) then
342    !cg(:,1:spacedim*blockdim) = cg(:,1:spacedim*blockdim) * inv_sqrt2
343    call xgBlock_scale(X,inv_sqrt2,1)
344    if(l_mpi_enreg%me_g0 == 1) cg(:, 1:spacedim*blockdim:l_npw) = cg(:, 1:spacedim*blockdim:l_npw) * sqrt2
345  end if
346 
347  if ( size(l_gvnlc) < 2*blockdim*spacedim ) then
348    ABI_FREE(l_gvnlc)
349    ABI_MALLOC(l_gvnlc,(2,blockdim*spacedim))
350  end if
351 
352   if (l_mpi_enreg%paral_kgb==0) then
353 
354     call multithreaded_getghc(l_cpopt,cg(:,1:blockdim*spacedim),cprj_dum,ghc,gsc(:,1:blockdim*spacedim),&
355       l_gs_hamk,l_gvnlc,dum, l_mpi_enreg,blockdim,l_prtvol,l_sij_opt,l_tim_getghc,0)
356 
357   else
358     call prep_getghc(cg(:,1:blockdim*spacedim),l_gs_hamk,l_gvnlc,ghc,gsc(:,1:blockdim*spacedim),dum,blockdim,l_mpi_enreg,&
359 &                     l_prtvol,l_sij_opt,l_cpopt,cprj_dum,already_transposed=.false.)
360   end if
361 
362   ! scale cg, ghc, gsc
363   if ( l_istwf == 2 ) then
364     !cg(:,1:spacedim*blockdim) = cg(:,1:spacedim*blockdim) * sqrt2
365     !ghc(:,1:spacedim*blockdim) = ghc(:,1:spacedim*blockdim) * sqrt2
366     call xgBlock_scale(X,sqrt2,1)
367     call xgBlock_scale(AX,sqrt2,1)
368     if(l_mpi_enreg%me_g0 == 1) then
369       cg(:, 1:spacedim*blockdim:l_npw) = cg(:, 1:spacedim*blockdim:l_npw) * inv_sqrt2
370       ghc(:, 1:spacedim*blockdim:l_npw) = ghc(:, 1:spacedim*blockdim:l_npw) * inv_sqrt2
371     endif
372     if(l_paw) then
373       !gsc(:,1:spacedim*blockdim) = gsc(:,1:spacedim*blockdim) * sqrt2
374       call xgBlock_scale(BX,sqrt2,1)
375       if(l_mpi_enreg%me_g0 == 1) gsc(:, 1:spacedim*blockdim:l_npw) = gsc(:, 1:spacedim*blockdim:l_npw) * inv_sqrt2
376     end if
377   end if
378 
379   if ( .not. l_paw ) call xgBlock_copy(X,BX)
380 
381   !call xgBlock_set(AX,ghc,0,spacedim)
382   !call xgBlock_set(BX,gsc(:,1:blockdim*spacedim),0,spacedim)
383  end subroutine getghc_gsc
384 
385  subroutine precond(W)
386    use m_xg, only : xg_t, xgBlock_colwiseMul
387 
388 !This section has been created automatically by the script Abilint (TD).
389 !Do not modify the following lines by hand.
390 #undef ABI_FUNC
391 #define ABI_FUNC 'precond'
392 !End of the abilint section
393 
394    type(xgBlock_t), intent(inout) :: W
395    integer :: ispinor
396    !integer :: cplx
397 
398    ! precondition resid_vec
399    do ispinor = 1,l_nspinor
400      !do cplx = 1, l_icplx
401      call xgBlock_colwiseMul(W,l_pcon,l_icplx*l_npw*(ispinor-1))
402       !end do
403    end do
404 
405  end subroutine precond
406 
407  end module m_lobpcgwf