Line data Source code
1 : !!****m* ABINIT/m_nonlop_test
2 : !! NAME
3 : !! m_nonlop_test
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2017-2026 ABINIT group (MT)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_nonlop_test
22 :
23 : implicit none
24 :
25 : private
26 : !!***
27 :
28 : public :: nonlop_test
29 : !!***
30 :
31 : contains
32 : !!***
33 :
34 : !!****f* ABINIT/nonlop_test
35 : !! NAME
36 : !! nonlop_test
37 : !!
38 : !! FUNCTION
39 : !! This routine is supposed to be used only for testing purpose.
40 : !! It tests the "nonlop" routine application (non-local operator) with respect to Finite Differences.
41 : !! It is not supposed to be used standalone, but via the nonlop_dfpt_test.py script to be found
42 : !! in ~abinit/scripts/post_processing/nonlop_dfpt_test directory. This Python script
43 : !! launches Abinit (several datasets) and analyse the result, in order to compare
44 : !! <Psi_i|H^(i)|Psi_j> compute with DFPT or Finite Differences.
45 : !! H^(i) is the ith derivative of the Hamiltonian with respect to one or several perturbation(s).
46 : !!
47 : !! INPUTS
48 : !! cg(2,mcg)=wavefunctions (may be read from disk file)
49 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
50 : !! istwfk(nkpt)=input option parameter that describes the storage of wfs
51 : !! kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis
52 : !! kpt(3,nkpt)=k points in reduced coordinates
53 : !! mband=maximum number of bands
54 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
55 : !! mgfft=maximum size of 1D FFTs
56 : !! mkmem=number of k points treated by this node.
57 : !! mpi_enreg=information about MPI parallelization
58 : !! mpw= maximum number of plane waves
59 : !! my_natom=number of atoms treated by current processor
60 : !! natom=number of atoms in cell.
61 : !! nband(nkpt)=number of bands at each k point
62 : !! nfft=number of FFT grid points
63 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
64 : !! nkpt=number of k points in Brillouin zone
65 : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
66 : !! npwarr(nkpt)=number of planewaves in basis and boundary at each k
67 : !! nspden=Number of spin Density components
68 : !! nspinor=number of spinorial components of the wavefunctions
69 : !! nsppol=1 for unpolarized, 2 for spin-polarized
70 : !! ntypat=number of types of atoms
71 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
72 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
73 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
74 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
75 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
76 : !! typat(natom)=type of each atom
77 : !! xred(3,natom)=reduced dimensionless atomic coordinates
78 : !!
79 : !! OUTPUT
80 : !!
81 : !! SOURCE
82 :
83 0 : subroutine nonlop_test(cg,eigen,istwfk,kg,kpt,mband,mcg,mgfft,mkmem,mpi_enreg,mpw,my_natom,natom,&
84 0 : & nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,ntypat,&
85 0 : & paw_ij,pawtab,ph1d,psps,rprimd,typat,xred)
86 :
87 : use defs_basis
88 : use m_abicore
89 : use m_xmpi
90 : use m_errors
91 : use m_hamiltonian
92 : use m_pawtab
93 : use m_paw_ij
94 : use m_pawcprj
95 : use m_cgtools
96 :
97 :
98 : use defs_datatypes, only : pseudopotential_type
99 : use defs_abitypes, only : MPI_type
100 : use m_kg, only : mkkpg
101 : use m_initylmg, only : initylmg
102 : use m_mkffnl, only : mkffnl
103 : use m_mpinfo, only : proc_distrb_cycle
104 : use m_nonlop, only : nonlop
105 :
106 : !Arguments ------------------------------------
107 : !scalars
108 : integer :: mband,mcg,mgfft,mkmem,mpw,my_natom,natom,nfft,nkpt,nspden,nspinor,nsppol,ntypat
109 : type(MPI_type),intent(inout) :: mpi_enreg
110 : type(pseudopotential_type),intent(in) :: psps
111 : !arrays
112 : integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol),ngfft(18),nloalg(3),npwarr(nkpt),typat(natom)
113 : real(dp),intent(in) :: cg(2,mcg),eigen(mband*nkpt*nsppol),kpt(3,nkpt), ph1d(2,3*(2*mgfft+1)*natom)
114 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
115 : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
116 : type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
117 : !Local variables-------------------------------
118 : !scalars
119 : integer,parameter :: ndtset_test=6,tim_nonlop=4
120 : integer,save :: idtset=0
121 : integer :: bandpp,bdtot_index,blocksize,choice,cplex,cpopt,dimffnl,iatm,iatom,iatom_only
122 : integer :: iband,iband_last,iband_test,iblock,icg,ider_ffnl,idir,idir_ffnl,idir_nonlop
123 : integer :: ii,ikg,ikpt,ilm,inlout,isppol,istwf_k,me_distrb,my_nspinor,nband_k,nblockbd
124 : integer :: nkpg,nnlout,npw_k,paw_opt,signs,spaceComm
125 : logical :: ex
126 : character(len=100) :: strg
127 : real(dp) :: argr,argi
128 0 : type(gs_hamiltonian_type) :: gs_hamk
129 : !arrays
130 0 : integer,allocatable :: kg_k(:,:)
131 : real(dp) :: kpoint(3),rmet(3,3)
132 0 : real(dp),allocatable :: cwavef(:,:),cwavef_out(:,:),enl(:,:,:,:),enlout(:),kpg_k(:,:),lambda(:)
133 0 : real(dp),allocatable :: scwavef_out(:,:),ylm(:,:),ylmgr(:,:,:),ylm_k(:,:),ylmgr_k(:,:,:)
134 0 : real(dp),allocatable,target :: ffnl(:,:,:,:),ph3d(:,:,:)
135 0 : type(pawcprj_type) :: cwaveprj(1,1)
136 :
137 : !*************************************************************************
138 :
139 : !Increment dataset counter
140 0 : idtset=idtset+1
141 0 : if (idtset<=2) return
142 :
143 : !Data from parallelism
144 0 : spaceComm=mpi_enreg%comm_kpt
145 0 : me_distrb=mpi_enreg%me_kpt
146 0 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
147 :
148 : !Initialize Hamiltonian datastructure
149 : call gs_hamk%init(psps,pawtab,nspinor,nsppol,nspden,natom,&
150 : & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,&
151 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
152 0 : & mpi_spintab=mpi_enreg%my_isppoltab,paw_ij=paw_ij,ph1d=ph1d)
153 0 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
154 :
155 : !Check for existence of files in the current directory\
156 : ! and set parameters accordingly
157 0 : choice=1 ; idir=0 ; signs=1
158 0 : if(idtset<ndtset_test)then
159 0 : inquire(file='config/signs1',exist=ex) ; if(ex) signs=1
160 0 : inquire(file='config/signs2',exist=ex) ; if(ex) signs=2
161 0 : do ii=1,100
162 0 : if (ii< 10) write(unit=strg,fmt='(a13,i1)') "config/choice",ii
163 0 : if (ii>=10) write(unit=strg,fmt='(a13,i2)') "config/choice",ii
164 0 : inquire(file=trim(strg),exist=ex) ; if(ex) choice=ii
165 : end do
166 0 : do ii=1,9
167 0 : write(unit=strg,fmt='(a11,i1)') "config/idir",ii
168 0 : inquire(file=trim(strg),exist=ex) ; if(ex) idir=ii
169 : end do
170 : else
171 0 : inquire(file='config/signsdfpt1',exist=ex) ; if(ex)signs=1
172 0 : inquire(file='config/signsdfpt2',exist=ex) ; if(ex)signs=2
173 0 : do ii=1,100
174 0 : if (ii< 10) write(unit=strg,fmt='(a17,i1)') "config/choicedfpt",ii
175 0 : if (ii>=10) write(unit=strg,fmt='(a17,i2)') "config/choicedfpt",ii
176 0 : inquire(file=trim(strg),exist=ex) ; if(ex) choice=ii
177 : end do
178 0 : do ii=1,36
179 0 : if (ii< 10) write(unit=strg,fmt='(a15,i1)') "config/idirdfpt",ii
180 0 : if (ii>=10) write(unit=strg,fmt='(a15,i2)') "config/idirdfpt",ii
181 0 : inquire(file=trim(strg),exist=ex) ; if(ex) idir=ii
182 : end do
183 : end if
184 0 : iatom=1 ; iband_test=-1
185 0 : do ii=1,50
186 0 : if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iatom",ii
187 0 : if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iatom",ii
188 0 : inquire(file=trim(strg),exist=ex) ; if(ex) iatom=ii
189 0 : if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iband",ii
190 0 : if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iband",ii
191 0 : inquire(file=trim(strg),exist=ex) ; if(ex) iband_test=ii
192 : end do
193 :
194 : !Set parameters for the "nonlop" routine according to users choice
195 0 : cpopt=-1 ; paw_opt=3*psps%usepaw ; iatm=gs_hamk%atindx(iatom)
196 0 : inquire(file='config/dij',exist=ex);if(ex) paw_opt=1*psps%usepaw
197 0 : if(signs==1)then
198 0 : iatom_only=-1 ; idir_ffnl=0
199 0 : idir_nonlop=0 ; cplex=1
200 0 : if(choice==1)then
201 0 : ider_ffnl=0
202 0 : nnlout=1 ; inlout=1
203 : end if
204 0 : if(choice==2)then
205 0 : ider_ffnl=0
206 0 : nnlout=3*natom ; inlout=3*(iatm-1)+idir ! Atoms are type-sorted in enlout()
207 : end if
208 0 : if(choice==3)then
209 0 : ider_ffnl=1
210 0 : nnlout=6 ; inlout=idir
211 : end if
212 0 : if(choice==5)then
213 0 : ider_ffnl=1
214 0 : nnlout=3 ; inlout=idir
215 : end if
216 0 : if(choice==51.or.choice==52)then
217 0 : ider_ffnl=1 ; cplex=2
218 0 : nnlout=6 ; inlout=2*idir-1
219 : end if
220 0 : if(choice==54)then
221 0 : ider_ffnl=2 ; cplex=2
222 0 : nnlout=18*natom ; inlout=18*(iatm-1)+2*idir-1 ! Atoms are type-sorted in enlout()
223 : end if
224 0 : if(choice==55)then
225 0 : ider_ffnl=2 ; cplex=2
226 0 : nnlout=36 ; inlout=2*idir-1
227 : end if
228 0 : if(choice==8)then
229 0 : ider_ffnl=2
230 0 : nnlout=6 ; inlout=idir
231 : end if
232 0 : if(choice==81)then
233 0 : ider_ffnl=2 ; cplex=2
234 0 : nnlout=18 ; inlout=2*idir-1
235 : end if
236 0 : else if(signs==2)then
237 0 : nnlout=1 ; inlout =1 ; cplex=1
238 0 : idir_nonlop=idir ; iatom_only=-1
239 0 : if(choice==1)then
240 0 : ider_ffnl=0 ; idir_ffnl=0
241 : end if
242 0 : if(choice==2)then
243 0 : iatom_only=iatom
244 0 : ider_ffnl=0 ; idir_ffnl=0
245 : end if
246 0 : if(choice==3)then
247 0 : ider_ffnl=1 ; idir_ffnl=-7
248 : end if
249 0 : if(choice==5)then
250 0 : ider_ffnl=1 ; idir_ffnl=4
251 : end if
252 0 : if(choice==51.or.choice==52)then
253 0 : ider_ffnl=1 ; idir_ffnl=4 ; cplex=2
254 : end if
255 0 : if(choice==54)then
256 0 : iatom_only=iatom
257 0 : ider_ffnl=2 ; idir_ffnl=4 ; cplex=2
258 : end if
259 0 : if(choice==8)then
260 0 : ider_ffnl=2 ; idir_ffnl=4
261 : end if
262 0 : if(choice==81)then
263 0 : ider_ffnl=2 ; idir_ffnl=4 ; cplex=2
264 : end if
265 : end if
266 :
267 : !Set parameters for the "mkffnl" routine according to users choice
268 0 : dimffnl=1+ider_ffnl
269 0 : if (ider_ffnl==1.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=2+2*psps%useylm
270 0 : if (ider_ffnl==2.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=3+7*psps%useylm
271 0 : if (ider_ffnl==1.and.idir_ffnl==-7) dimffnl=2+5*psps%useylm
272 0 : if (idir_ffnl>-7.and.idir_ffnl<0) dimffnl=2
273 :
274 : !Write recognizable statement in log file
275 : write(std_out,'(2(a,i2),(a,i1),2(a,i2),(a,i1),(a,i2),(a,i1),2(a,i2))') &
276 0 : & "TESTDFPT: choice=",choice,", idir(mkffnl)=",idir_ffnl,&
277 0 : & ", ider(mkffnl)=",ider_ffnl,", dimffnl=",dimffnl,&
278 0 : & ", idir(nonlop)=",idir_nonlop,", signs=",signs,&
279 0 : & ", iatom=",iatom_only,", paw_opt=",paw_opt,&
280 0 : & ", nnlout=",nnlout,", inlout=",inlout
281 :
282 : !Compute all spherical harmonics and gradients
283 0 : ABI_MALLOC(ylm,(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm))
284 0 : ABI_MALLOC(ylmgr,(mpw*mkmem,9,psps%mpsang*psps%mpsang*psps%useylm))
285 0 : if (psps%useylm==1) then
286 : call initylmg(gs_hamk%gprimd,kg,kpt,mkmem,mpi_enreg,psps%mpsang,mpw,nband,nkpt,&
287 0 : & npwarr,nsppol,2,rprimd,ylm,ylmgr)
288 : else
289 0 : ylm=zero ; ylmgr=zero
290 : end if
291 :
292 : !No loop over spins; only do the first one
293 0 : bdtot_index=0 ; icg=0 ; isppol=1
294 :
295 : !Continue to initialize the Hamiltonian (PAW DIJ coefficients)
296 0 : call gs_hamk%load_spin(isppol,with_nonlocal=.true.)
297 :
298 : !No loop over k points; only do the first one
299 0 : ikg=0 ; ikpt=1
300 :
301 0 : nband_k=nband(ikpt+(isppol-1)*nkpt)
302 0 : istwf_k=istwfk(ikpt)
303 0 : npw_k=npwarr(ikpt)
304 0 : kpoint(:)=kpt(:,ikpt)
305 :
306 : !My spin/kpoint or not?
307 0 : if(.not.proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
308 :
309 : ! Parallelism over FFT and/or bands: define sizes and tabs
310 0 : bandpp=mpi_enreg%bandpp
311 0 : nblockbd=nband_k/bandpp
312 0 : blocksize=nband_k/nblockbd
313 :
314 : ! Several allocations
315 0 : ABI_MALLOC(lambda,(blocksize))
316 0 : ABI_MALLOC(enlout,(nnlout*blocksize))
317 0 : ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize))
318 0 : ABI_MALLOC(cwavef_out,(2,npw_k))
319 0 : if (paw_opt>=3) then
320 0 : ABI_MALLOC(scwavef_out,(2,npw_k))
321 0 : ABI_MALLOC(enl,(0,0,0,0))
322 : else
323 0 : ABI_MALLOC(scwavef_out,(0,0))
324 0 : ABI_MALLOC(enl,(gs_hamk%dimekb1,gs_hamk%dimekb2,gs_hamk%nspinor**2,1))
325 0 : enl(:,:,:,:)=one
326 : end if
327 :
328 : ! Compute (k+G) vectors and associated spherical harmonics
329 0 : nkpg=3*nloalg(3)
330 0 : ABI_MALLOC(kg_k,(3,mpw))
331 0 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
332 0 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
333 0 : ABI_MALLOC(ylmgr_k,(npw_k,9,psps%mpsang*psps%mpsang*psps%useylm))
334 0 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
335 0 : if (nkpg>0) then
336 0 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
337 : end if
338 0 : if (psps%useylm==1) then
339 0 : do ilm=1,psps%mpsang*psps%mpsang
340 0 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
341 : end do
342 0 : if (ider_ffnl>=1) then
343 0 : do ilm=1,psps%mpsang*psps%mpsang
344 0 : do ii=1,3+6*(ider_ffnl/2)
345 0 : ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
346 : end do
347 : end do
348 : end if
349 : end if
350 :
351 : ! Compute non-local form factors
352 0 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
353 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
354 : & gs_hamk%gmet,gs_hamk%gprimd,ider_ffnl,idir_ffnl,psps%indlmn,kg_k,kpg_k,&
355 : & gs_hamk%kpt_k,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
356 : & npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
357 0 : & psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
358 :
359 : ! Load k-dependent part in the Hamiltonian datastructure
360 0 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
361 : call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
362 0 : & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.)
363 :
364 0 : do iblock=1,nblockbd
365 :
366 0 : iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
367 0 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
368 :
369 : ! Select a specific band or all
370 0 : if (iband==iband_test.or.iband_test==-1) then
371 :
372 : ! Load contribution from block(n,k)
373 : cwavef(:,1:npw_k*my_nspinor*blocksize)=&
374 0 : & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
375 0 : lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
376 :
377 : ! Call NONLOP
378 0 : if (paw_opt<3) then
379 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,&
380 : & mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,&
381 0 : & iatom_only=iatom_only,enl=enl)
382 : else
383 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,&
384 : & mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,&
385 0 : & iatom_only=iatom_only)
386 : end if
387 :
388 : ! Post-processing if nonlop is called with specific options
389 0 : if (signs==2) then
390 0 : if (paw_opt<3) then
391 : call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,cwavef_out,&
392 0 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
393 : else
394 : call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,scwavef_out,&
395 0 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
396 : end if
397 0 : enlout(inlout)=argr
398 : end if
399 0 : if (signs==1.and.choice==1) then
400 : call dotprod_g(argr,argi,istwf_k,npw_k,1,cwavef,cwavef,&
401 0 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
402 0 : enlout(:)=enlout(:)+argr
403 : end if
404 :
405 : ! Write recognizable statements in log file
406 0 : if (idtset<ndtset_test) then
407 0 : write(std_out,'(a,i3,es24.16)') "TESTDFPT_df: ",iband,enlout(inlout)
408 : else
409 0 : write(std_out,'(a,i3,es24.16)') "TESTDFPT_dfpt:",iband,enlout(inlout)
410 : end if
411 :
412 : end if
413 :
414 : end do ! End of loop on block of bands
415 :
416 : ! Increment indexes (not used here because only one spin/kpoint)
417 0 : icg=icg+npw_k*my_nspinor*nband_k
418 0 : ikg=ikg+npw_k
419 :
420 : end if ! Not my spin/kpoint
421 :
422 0 : bdtot_index=bdtot_index+nband_k
423 :
424 : !Memory deallocations
425 0 : ABI_FREE(enl)
426 0 : ABI_FREE(enlout)
427 0 : ABI_FREE(lambda)
428 0 : ABI_FREE(ph3d)
429 0 : ABI_FREE(ffnl)
430 0 : ABI_FREE(cwavef)
431 0 : ABI_FREE(cwavef_out)
432 0 : ABI_FREE(scwavef_out)
433 0 : ABI_FREE(kg_k)
434 0 : ABI_FREE(kpg_k)
435 0 : ABI_FREE(ylm_k)
436 0 : ABI_FREE(ylmgr_k)
437 0 : ABI_FREE(ylm)
438 0 : ABI_FREE(ylmgr)
439 0 : call gs_hamk%free()
440 :
441 0 : end subroutine nonlop_test
442 : !!***
443 :
444 : end module m_nonlop_test
445 : !!***
|
