TABLE OF CONTENTS


ABINIT/nonlop_test [ Functions ]

[ Top ] [ Functions ]

NAME

 nonlop_test

FUNCTION

 This routine is supposed to be used only for testing purpose.
 It tests the "nonlop" routine application (non-local operator) with respect to Finite Differences.
 It is not supposed to be used standalone, but via the nonlop_dfpt_test.py script to be found
 in ~abinit/scripts/post_processing/nonlop_dfpt_test directory. This Python script
 launches Abinit (several datasets) and analyse the result, in order to compare
  <Psi_i|H^(i)|Psi_j> compute with DFPT or Finite Differences.
 H^(i) is the ith derivative of the Hamiltonian with respect to one or several perturbation(s).

COPYRIGHT

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

INPUTS

  cg(2,mcg)=wavefunctions (may be read from disk file)
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  istwfk(nkpt)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis
  kpt(3,nkpt)=k points in reduced coordinates
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mgfft=maximum size of 1D FFTs
  mkmem=number of k points treated by this node.
  mpi_enreg=informations about MPI parallelization
  mpw= maximum number of plane waves
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nband(nkpt)=number of bands at each k point
  nfft=number of FFT grid points
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt=number of k points in Brillouin zone
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npwarr(nkpt)=number of planewaves in basis and boundary at each k
  nspden=Number of spin Density components
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  ntypat=number of types of atoms
  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  typat(natom)=type of each atom
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

PARENTS

      afterscfloop

CHILDREN

      destroy_hamiltonian,dotprod_g,init_hamiltonian,initylmg
      load_k_hamiltonian,load_spin_hamiltonian,mkffnl,mkkpg,nonlop

SOURCE

 65 #if defined HAVE_CONFIG_H
 66 #include "config.h"
 67 #endif
 68 
 69 #include "abi_common.h"
 70 
 71 subroutine nonlop_test(cg,eigen,istwfk,kg,kpt,mband,mcg,mgfft,mkmem,mpi_enreg,mpw,my_natom,natom,&
 72 &                      nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,ntypat,&
 73 &                       paw_ij,pawtab,ph1d,psps,rprimd,typat,xred)
 74 
 75  use defs_basis
 76  use defs_datatypes
 77  use defs_abitypes
 78  use m_profiling_abi
 79  use m_xmpi
 80  use m_errors
 81  use m_hamiltonian
 82  use m_pawtab
 83  use m_paw_ij
 84  use m_pawcprj
 85  use m_cgtools
 86 
 87 !This section has been created automatically by the script Abilint (TD).
 88 !Do not modify the following lines by hand.
 89 #undef ABI_FUNC
 90 #define ABI_FUNC 'nonlop_test'
 91  use interfaces_32_util
 92  use interfaces_56_recipspace
 93  use interfaces_66_nonlocal, except_this_one => nonlop_test
 94 !End of the abilint section
 95 
 96  implicit none
 97 
 98 !Arguments ------------------------------------
 99 !scalars
100  integer :: mband,mcg,mgfft,mkmem,mpw,my_natom,natom,nfft,nkpt,nspden,nspinor,nsppol,ntypat
101  type(MPI_type),intent(inout) :: mpi_enreg
102  type(pseudopotential_type),intent(in) :: psps
103 !arrays
104  integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol),ngfft(18),nloalg(3),npwarr(nkpt),typat(natom)
105  real(dp),intent(in) :: cg(2,mcg),eigen(mband*nkpt*nsppol),kpt(3,nkpt), ph1d(2,3*(2*mgfft+1)*natom)
106  real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
107  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
108  type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
109 !Local variables-------------------------------
110 !scalars
111  integer,parameter :: ndtset_test=6,tim_nonlop=4
112  integer,save :: idtset=0
113  integer :: bandpp,bdtot_index,blocksize,choice,cplex,cpopt,dimffnl,iatom,iatom_only
114  integer :: iband,iband_last,iband_test,iblock,icg,ider_ffnl,idir,idir_ffnl,idir_nonlop
115  integer :: ii,ikg,ikpt,ilm,inlout,isppol,istwf_k,me_distrb,my_nspinor,nband_k,nblockbd
116  integer :: nkpg,nnlout,npw_k,paw_opt,signs,spaceComm
117  logical :: ex
118  character(len=100) :: strg
119  real(dp) :: argr,argi
120  type(gs_hamiltonian_type) :: gs_hamk
121 !arrays
122  integer,allocatable :: kg_k(:,:)
123  real(dp) :: kpoint(3),rmet(3,3)
124  real(dp),allocatable :: cwavef(:,:),cwavef_out(:,:),enl(:,:,:),enlout(:),kpg_k(:,:),lambda(:)
125  real(dp),allocatable :: scwavef_out(:,:),ylm(:,:),ylmgr(:,:,:),ylm_k(:,:),ylmgr_k(:,:,:)
126  real(dp),allocatable,target :: ffnl(:,:,:,:),ph3d(:,:,:)
127  type(pawcprj_type) :: cwaveprj(1,1)
128 
129 !*************************************************************************
130 
131 !Increment dataset counter
132  idtset=idtset+1
133  if (idtset<=2) return
134 
135 !Data from parallelism
136  spaceComm=mpi_enreg%comm_kpt
137  me_distrb=mpi_enreg%me_kpt
138  my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
139 
140 !Initialize Hamiltonian datastructure
141  call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,nsppol,nspden,natom,&
142 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,&
143 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
144 & mpi_spintab=mpi_enreg%my_isppoltab,paw_ij=paw_ij,ph1d=ph1d)
145  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
146 
147 !Check for existence of files in the current directory\
148 ! and set parameters accordingly
149  choice=1 ; idir=0 ; signs=1
150  if(idtset<ndtset_test)then
151    inquire(file='config/signs1',exist=ex) ; if(ex) signs=1
152    inquire(file='config/signs2',exist=ex) ; if(ex) signs=2
153    do ii=1,100
154      if (ii< 10) write(unit=strg,fmt='(a13,i1)') "config/choice",ii
155      if (ii>=10) write(unit=strg,fmt='(a13,i2)') "config/choice",ii
156      inquire(file=trim(strg),exist=ex)  ; if(ex) choice=ii
157    end do
158    do ii=1,9
159      write(unit=strg,fmt='(a11,i1)') "config/idir",ii
160      inquire(file=trim(strg),exist=ex)  ; if(ex) idir=ii
161    end do
162  else
163    inquire(file='config/signsdfpt1',exist=ex)  ; if(ex)signs=1
164    inquire(file='config/signsdfpt2',exist=ex)  ; if(ex)signs=2
165    do ii=1,100
166      if (ii< 10) write(unit=strg,fmt='(a17,i1)') "config/choicedfpt",ii
167      if (ii>=10) write(unit=strg,fmt='(a17,i2)') "config/choicedfpt",ii
168      inquire(file=trim(strg),exist=ex)  ; if(ex) choice=ii
169    end do
170    do ii=1,36
171      if (ii< 10) write(unit=strg,fmt='(a15,i1)') "config/idirdfpt",ii
172      if (ii>=10) write(unit=strg,fmt='(a15,i2)') "config/idirdfpt",ii
173      inquire(file=trim(strg),exist=ex) ; if(ex) idir=ii
174    end do
175  end if
176  iatom=1 ; iband_test=-1
177  do ii=1,50
178    if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iatom",ii
179    if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iatom",ii
180    inquire(file=trim(strg),exist=ex)  ; if(ex) iatom=ii
181    if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iband",ii
182    if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iband",ii
183    inquire(file=trim(strg),exist=ex)  ; if(ex) iband_test=ii
184  end do
185 
186 !Set parameters for the "nonlop" routine according to users choice
187  cpopt=-1 ; paw_opt=3*psps%usepaw
188  inquire(file='config/dij',exist=ex);if(ex) paw_opt=1*psps%usepaw
189  if(signs==1)then
190    iatom_only=-1 ; idir_ffnl=0
191    idir_nonlop=0 ; cplex=1
192    if(choice==1)then
193      ider_ffnl=0
194      nnlout=1 ; inlout=1
195    end if
196    if(choice==2)then
197      ider_ffnl=0
198      nnlout=3*natom ; inlout=3*(iatom-1)+idir
199    end if
200    if(choice==3)then
201      ider_ffnl=1
202      nnlout=6 ; inlout=idir
203    end if
204    if(choice==5)then
205      ider_ffnl=1
206      nnlout=3 ; inlout=idir
207    end if
208    if(choice==51.or.choice==52)then
209      ider_ffnl=1 ; cplex=2
210      nnlout=6 ; inlout=2*idir-1
211    end if
212    if(choice==54)then
213      ider_ffnl=2 ; cplex=2
214      nnlout=18*natom ; inlout=18*(iatom-1)+2*idir-1
215    end if
216    if(choice==55)then
217      ider_ffnl=2 ; cplex=2
218      nnlout=36 ; inlout=2*idir-1
219    end if
220    if(choice==8)then
221      ider_ffnl=2
222      nnlout=6 ; inlout=idir
223    end if
224    if(choice==81)then
225      ider_ffnl=2 ; cplex=2
226      nnlout=18 ; inlout=2*idir-1
227    end if
228  else if(signs==2)then
229    nnlout=1 ; inlout =1 ; cplex=1
230    idir_nonlop=idir ; iatom_only=-1
231    if(choice==1)then
232      ider_ffnl=0 ; idir_ffnl=0
233    end if
234    if(choice==2)then
235      iatom_only=iatom
236      ider_ffnl=0 ; idir_ffnl=0
237    end if
238    if(choice==3)then
239      ider_ffnl=1 ; idir_ffnl=-7
240    end if
241    if(choice==5)then
242      ider_ffnl=1 ; idir_ffnl=4
243    end if
244    if(choice==51.or.choice==52)then
245      ider_ffnl=1 ; idir_ffnl=4 ; cplex=2
246    end if
247    if(choice==54)then
248      iatom_only=iatom
249      ider_ffnl=2 ; idir_ffnl=4 ; cplex=2
250    end if
251    if(choice==8)then
252      ider_ffnl=2 ; idir_ffnl=4
253    end if
254    if(choice==81)then
255      ider_ffnl=2 ; idir_ffnl=4 ; cplex=2
256    end if
257  end if
258 
259 !Set parameters for the "mkffnl" routine according to users choice
260  dimffnl=1+ider_ffnl
261  if (ider_ffnl==1.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=2+2*psps%useylm
262  if (ider_ffnl==2.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=3+7*psps%useylm
263  if (ider_ffnl==1.and.idir_ffnl==-7) dimffnl=2+5*psps%useylm
264  if (idir_ffnl>-7.and.idir_ffnl<0) dimffnl=2
265 
266 !Write recognizable statement in log file
267  write(std_out,'(2(a,i2),(a,i1),2(a,i2),(a,i1),(a,i2),(a,i1),2(a,i2))') &
268 & "TESTDFPT: choice=",choice,", idir(mkffnl)=",idir_ffnl,&
269 & ", ider(mkffnl)=",ider_ffnl,", dimffnl=",dimffnl,&
270 & ", idir(nonlop)=",idir_nonlop,", signs=",signs,&
271 & ", iatom=",iatom_only,", paw_opt=",paw_opt,&
272 & ", nnlout=",nnlout,", inlout=",inlout
273 
274 !Compute all spherical harmonics and gradients
275  ABI_ALLOCATE(ylm,(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm))
276  ABI_ALLOCATE(ylmgr,(mpw*mkmem,9,psps%mpsang*psps%mpsang*psps%useylm))
277  if (psps%useylm==1) then
278    call initylmg(gs_hamk%gprimd,kg,kpt,mkmem,mpi_enreg,psps%mpsang,mpw,nband,nkpt,&
279 &   npwarr,nsppol,2,rprimd,ylm,ylmgr)
280  else
281    ylm=zero ; ylmgr=zero
282  end if
283 
284 !No loop over spins; only do the first one
285  bdtot_index=0 ; icg=0 ; isppol=1
286 
287 !Continue to initialize the Hamiltonian (PAW DIJ coefficients)
288  call load_spin_hamiltonian(gs_hamk,isppol,with_nonlocal=.true.)
289 
290 !No loop over k points; only do the first one
291  ikg=0 ; ikpt=1
292 
293  nband_k=nband(ikpt+(isppol-1)*nkpt)
294  istwf_k=istwfk(ikpt)
295  npw_k=npwarr(ikpt)
296  kpoint(:)=kpt(:,ikpt)
297 
298 !My spin/kpoint or not?
299  if(.not.proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
300 
301 !  Parallelism over FFT and/or bands: define sizes and tabs
302    bandpp=mpi_enreg%bandpp
303    nblockbd=nband_k/bandpp
304    blocksize=nband_k/nblockbd
305 
306 !  Several allocations
307    ABI_ALLOCATE(lambda,(blocksize))
308    ABI_ALLOCATE(enlout,(nnlout*blocksize))
309    ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor*blocksize))
310    ABI_ALLOCATE(cwavef_out,(2,npw_k))
311    if (paw_opt>=3) then
312      ABI_ALLOCATE(scwavef_out,(2,npw_k))
313      ABI_ALLOCATE(enl,(0,0,0))
314    else
315      ABI_ALLOCATE(scwavef_out,(0,0))
316      ABI_ALLOCATE(enl,(gs_hamk%dimekb1,gs_hamk%dimekb2,gs_hamk%nspinor**2))
317      enl(:,:,:)=one
318    end if
319 
320 !  Compute (k+G) vectors and associated spherical harmonics
321    nkpg=3*nloalg(3)
322    ABI_ALLOCATE(kg_k,(3,mpw))
323    ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
324    ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
325    ABI_ALLOCATE(ylmgr_k,(npw_k,9,psps%mpsang*psps%mpsang*psps%useylm))
326    kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
327    if (nkpg>0) then
328      call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
329    end if
330    if (psps%useylm==1) then
331      do ilm=1,psps%mpsang*psps%mpsang
332        ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
333      end do
334      if (ider_ffnl>=1) then
335        do ilm=1,psps%mpsang*psps%mpsang
336          do ii=1,3+6*(ider_ffnl/2)
337            ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
338          end do
339        end do
340      end if
341    end if
342 
343 !  Compute non-local form factors
344    ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
345    call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
346 &   gs_hamk%gmet,gs_hamk%gprimd,ider_ffnl,idir_ffnl,psps%indlmn,kg_k,kpg_k,&
347 &   gs_hamk%kpt_k,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
348 &   npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
349 &   psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
350 
351 !  Load k-dependent part in the Hamiltonian datastructure
352    ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
353    call load_k_hamiltonian(gs_hamk,kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
354 &   kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.)
355 
356    do iblock=1,nblockbd
357 
358      iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
359      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
360 
361 !    Select a specific band or all
362      if (iband==iband_test.or.iband_test==-1) then
363 
364 !      Load contribution from block(n,k)
365        cwavef(:,1:npw_k*my_nspinor*blocksize)=&
366 &       cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
367        lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
368 
369 !      Call NONLOP
370        if (paw_opt<3) then
371          call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,&
372 &         mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,&
373 &         iatom_only=iatom_only,enl=enl)
374        else
375          call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,&
376 &         mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,&
377 &         iatom_only=iatom_only)
378        end if
379 
380 !      Post-processing if nonlop is called with specific options
381        if (signs==2) then
382          if (paw_opt<3) then
383            call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,cwavef_out,&
384 &           mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
385          else
386            call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,scwavef_out,&
387 &           mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
388          end if
389          enlout(inlout)=argr
390        end if
391        if (signs==1.and.choice==1) then
392          call dotprod_g(argr,argi,istwf_k,npw_k,1,cwavef,cwavef,&
393 &         mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
394          enlout(:)=enlout(:)+argr
395        end if
396 
397 !      Write recognizable statements in log file
398        if (idtset<ndtset_test) then
399          write(std_out,'(a,i3,es24.16)') "TESTDFPT_df:  ",iband,enlout(inlout)
400        else
401          write(std_out,'(a,i3,es24.16)') "TESTDFPT_dfpt:",iband,enlout(inlout)
402        end if
403 
404      end if
405 
406    end do ! End of loop on block of bands
407 
408 !  Increment indexes (not used here because only one spin/kpoint)
409    icg=icg+npw_k*my_nspinor*nband_k
410    ikg=ikg+npw_k
411 
412  end if ! Not my spin/kpoint
413 
414  bdtot_index=bdtot_index+nband_k
415 
416 !Memory deallocations
417  ABI_DEALLOCATE(enl)
418  ABI_DEALLOCATE(enlout)
419  ABI_DEALLOCATE(lambda)
420  ABI_DEALLOCATE(ph3d)
421  ABI_DEALLOCATE(ffnl)
422  ABI_DEALLOCATE(cwavef)
423  ABI_DEALLOCATE(cwavef_out)
424  ABI_DEALLOCATE(scwavef_out)
425  ABI_DEALLOCATE(kg_k)
426  ABI_DEALLOCATE(kpg_k)
427  ABI_DEALLOCATE(ylm_k)
428  ABI_DEALLOCATE(ylmgr_k)
429  ABI_DEALLOCATE(ylm)
430  ABI_DEALLOCATE(ylmgr)
431  call destroy_hamiltonian(gs_hamk)
432 
433 end subroutine nonlop_test