Branch data Line data Source code
1 : : !!****m* ABINIT/m_paw_optics
2 : : !! NAME
3 : : !! m_paw_optics
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains several routines related to conductivity:
7 : : !! optical conductivity, X spectroscopy, linear susceptibility, ...
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2018-2026 ABINIT group (SM,VR,FJ,MT,NB,PGhosh)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : ! nvtx related macro definition
24 : : #include "nvtx_macros.h"
25 : :
26 : : MODULE m_paw_optics
27 : :
28 : : use defs_basis
29 : : use m_xmpi
30 : : use m_errors
31 : : use m_wffile
32 : : use m_abicore
33 : : use m_hdr
34 : : use m_dtset
35 : : use m_dtfil
36 : : use m_nctk
37 : : use netcdf
38 : :
39 : : use defs_datatypes, only : pseudopotential_type
40 : : use defs_abitypes, only : MPI_type
41 : : use m_time, only : timab
42 : : use m_io_tools, only : open_file,get_unit,close_unit
43 : : use m_pawpsp, only : pawpsp_init_core
44 : : use m_paw_atomorb, only : atomorb_type,destroy_atomorb,copy_atomorb
45 : : use m_pawrad, only : pawrad_type,pawrad_deducer0,simp_gen,nderiv_gen,poisson
46 : : use m_pawtab, only : pawtab_type
47 : : use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_get, &
48 : : & pawcprj_free,pawcprj_mpi_allgather,pawcprj_pack
49 : : use m_pawang, only : pawang_type
50 : : use m_paw_denpot, only : pawdensities,pawkindensities,pawdenpot
51 : : use m_paw_an, only : paw_an_type,paw_an_init,paw_an_free,paw_an_copy
52 : : use m_pawrhoij, only : pawrhoij_type
53 : : use m_paw_ij, only : paw_ij_type
54 : : use m_paw_onsite, only : pawnabla_init,pawnabla_core_init
55 : : use m_paw_sphharm, only : setnabla_ylm
56 : : use m_pawxc, only : pawxc,pawxcm,pawxc_get_xclevel,pawxc_get_usekden
57 : : use m_rcpaw, only : rcpaw_type
58 : : use m_mpinfo, only : destroy_mpi_enreg,nullify_mpi_enreg,initmpi_seq,proc_distrb_cycle
59 : : use m_numeric_tools,only : kramerskronig
60 : : use m_geometry, only : metric
61 : : use m_hide_lapack, only : matrginv
62 : : use m_abi_linalg, only : abi_xgemm
63 : : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
64 : : #if defined(HAVE_GPU_MARKERS)
65 : : use m_nvtx
66 : : use m_nvtx_data
67 : : #endif
68 : :
69 : : implicit none
70 : :
71 : : private
72 : :
73 : : !public procedures.
74 : : public :: optics_paw
75 : : public :: optics_paw_core
76 : : public :: linear_optics_paw
77 : :
78 : : !I/O parameters
79 : : !Set to true to force the use of netCDF when available
80 : : ! overriding the value of dtset%iomode
81 : : logical,parameter :: use_netcdf_forced=.true.
82 : : !Set to true to use netcdf-MPIIO when available
83 : : logical,parameter :: use_netcdf_mpiio=.true.
84 : : !Set to true to compute/write only half of the (n,m) dipoles
85 : : logical,parameter :: compute_half_dipoles=.true.
86 : : !Set to true to use unlimited dimensions in netCDF file (experimental)
87 : : ! This is not mandatory because we know exactly the amount of data to write
88 : : ! and seems to impact performances negatively...
89 : : ! Not compatible with compute_half_dipoles=.true.
90 : : logical,parameter :: use_netcdf_unlimited=.false.
91 : :
92 : : !Computation parameters
93 : : !Target memory (in bytes) of the work arrays for one block of bands in optics_paw_dipoles_k
94 : : real(dp),parameter :: dipoles_blk_mem=512._dp*1024._dp**2
95 : :
96 : : CONTAINS !========================================================================================
97 : : !!***
98 : :
99 : : !----------------------------------------------------------------------
100 : :
101 : : !!****f* m_paw_optics/optics_paw
102 : : !! NAME
103 : : !! optics_paw
104 : : !!
105 : : !! FUNCTION
106 : : !! Compute matrix elements need for optical conductivity (in the PAW context) and store them in a file
107 : : !! Matrix elements = <Phi_i|Nabla|Phi_j>
108 : : !!
109 : : !! INPUTS
110 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
111 : : !! cg(2,mcg)=planewave coefficients of wavefunctions.
112 : : !! cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk>
113 : : !! and each |p_lmn> non-local projector
114 : : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
115 : : !! dtfil <type(datafiles_type)>=variables related to files
116 : : !! dtset <type(dataset_type)>=all input variables for this dataset
117 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
118 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
119 : : !! mband=maximum number of bands
120 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
121 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
122 : : !! mkmem =number of k points treated by this node.
123 : : !! mpi_enreg=information about MPI parallelization
124 : : !! mpsang =1+maximum angular momentum for nonlocal pseudopotentials
125 : : !! mpw=maximum dimensioned size of npw.
126 : : !! natom=number of atoms in cell.
127 : : !! nkpt=number of k points.
128 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
129 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
130 : : !! pawang <type(pawang_type)>= PAW ANGular mesh discretization and related data
131 : : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
132 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
133 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
134 : : !! znucl(ntypat)=atomic number of atom type
135 : : !!
136 : : !! OUTPUT
137 : : !! psinablapsi_out=contains the matrix elements
138 : : !! (The size of the psinablapsi_out decide whether we compute the full matrix or
139 : : !! only the diagonal part)
140 : : !! (if not present only writing in a file)
141 : : !!
142 : : !! SIDE EFFECTS
143 : : !!
144 : : !! NOTES
145 : : !!
146 : : !! SOURCE
147 : :
148 : 2316 : subroutine optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen0,gprimd,hdr,kg,&
149 : 1158 : & mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,&
150 [ + + + - : 1158 : & pawang,pawrad,pawrhoij,pawtab,znucl,psinablapsi_out)
- + ]
151 : :
152 : : !Arguments ------------------------------------
153 : : !scalars
154 : : integer,intent(in) :: mband,mcg,mcprj,mkmem,mpsang,mpw,natom,nkpt,nsppol
155 : : type(MPI_type),intent(in) :: mpi_enreg
156 : : type(datafiles_type),intent(in) :: dtfil
157 : : type(dataset_type),intent(in) :: dtset
158 : : type(hdr_type),intent(inout) :: hdr
159 : : type(pawang_type),intent(in) :: pawang
160 : : real(dp),optional,target,intent(out) :: psinablapsi_out(:,:,:,:)
161 : : !arrays
162 : : integer,intent(in) :: atindx1(natom),dimcprj(natom),npwarr(nkpt)
163 : : integer,intent(in),target :: kg(3,mpw*mkmem)
164 : : real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
165 : : real(dp),intent(in) :: gprimd(3,3),znucl(dtset%ntypat)
166 : : real(dp),intent(inout) :: cg(2,mcg)
167 : : type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
168 : : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
169 : : type(pawrhoij_type),intent(inout) :: pawrhoij(mpi_enreg%my_natom)
170 : : type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
171 : :
172 : : !Local variables-------------------------------
173 : : !scalars
174 : : integer,parameter :: master=0
175 : : integer :: bsize,iomode,bdtot_index,cplex,etiq,fformopt,iatom,ib,ibg
176 : : integer :: icg,ierr,ikg,ikpt,ount,ncid,varid
177 : : integer :: iorder_cprj,ipw,isppol,istwf_k,itypat
178 : : integer :: gpu_option,iblk,iproj,j0,j0s,j1,jb,jpaw1,jpaw2,jpw1,jpw2,lmn_size,lmnmax,mband_cprj
179 : : integer :: me_band_io,my_j0,my_j1,nbj_io,nent,nproc_band_io,option_core
180 : : integer :: my_nspinor,nband_k,nband_cprj_k,nproj,npw_k,sender,me,me_fft,master_spfftband,pnp_size
181 : : integer :: spaceComm_band,spaceComm_bandspinorfft,spaceComm_fft,spaceComm_kpt
182 : : integer :: spaceComm_spinorfft,spaceComm_w
183 : : integer, parameter :: NO_FILE_OUT=-1
184 : : logical :: already_has_nabla,cprj_paral_band,mykpt,iomode_etsf_mpiio
185 : : logical :: i_am_master,i_am_master_kpt,i_am_master_band,i_am_master_spfft,nc_unlimited,store_half_dipoles
186 : : logical :: diag_only
187 : : real(dp) :: el_temp
188 : : character(len=500) :: msg
189 : : type(nctkdim_t) :: nctkdim
190 : : !arrays
191 : : integer :: nc_count_5(5),nc_count_6(6),nc_start_5(5),nc_start_6(6),nc_stride_5(5),nc_stride_6(6),tmp_shape(3)
192 : 1158 : integer, contiguous, pointer :: kg_k(:,:)
193 : 1158 : integer,allocatable :: atom_off(:),lmn_at(:),nlmn(:),row_atom(:)
194 : : real(dp) :: kpoint(3),tsec(2)
195 : 1158 : real(dp),allocatable :: cprj_pk(:,:),kpg_k(:,:),nabla_at(:,:,:,:),soc_at(:,:,:,:,:,:)
196 : 1158 : real(dp),pointer :: psinablapsi(:,:,:),psinablapsi_w(:,:,:),psi_blk(:,:,:)
197 : 1158 : type(coeff5_type),allocatable :: phisocphj(:)
198 : 1158 : type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
199 : : type(nctkarr_t) :: nctk_arrays(1)
200 : :
201 : : ! ************************************************************************
202 : :
203 : : DBG_ENTER("COLL")
204 : 1158 : call timab(1184,1,tsec)
205 : : ABI_NVTX_START_RANGE(NVTX_OPTICS_PAW)
206 : :
207 : : !Compatibility tests
208 [ - + ]: 1158 : ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
209 [ - + - - ]: 1158 : ABI_CHECK(mpi_enreg%paral_spinor==0.or.dtset%pawspnorb==0,"spinor parallelization not supported with SOC!")
210 : : ! MJV 6/12/2008: looks like mpi_enreg may not be completely initialized here
211 : : if (xmpi_paral==1) then
212 [ + + ]: 4632 : tmp_shape = shape(mpi_enreg%proc_distrb)
213 [ - + ]: 1158 : if (nkpt > tmp_shape(1)) then
214 : 0 : ABI_BUG('problem with proc_distrb!')
215 : : end if
216 : : end if
217 : :
218 : : !Init parallelism
219 : 1158 : spaceComm_w=mpi_enreg%comm_cell
220 [ + + ]: 1158 : if (mpi_enreg%paral_kgb==1) then
221 : 844 : spaceComm_kpt=mpi_enreg%comm_kpt
222 : 844 : spaceComm_fft=mpi_enreg%comm_fft
223 : 844 : spaceComm_band=mpi_enreg%comm_band
224 : 844 : spaceComm_spinorfft=mpi_enreg%comm_spinorfft
225 : 844 : spaceComm_bandspinorfft=mpi_enreg%comm_bandspinorfft
226 : : else
227 : 314 : spaceComm_kpt=mpi_enreg%comm_kpt
228 : 314 : spaceComm_fft=xmpi_comm_self
229 : 314 : spaceComm_band=mpi_enreg%comm_band
230 : 314 : spaceComm_spinorfft=xmpi_comm_self
231 : : !Plane-wave and on-site contributions are reduced together over this communicator
232 : : !(comm_band is xmpi_comm_self here, except for band parallelism (paralbd=1))
233 : 314 : spaceComm_bandspinorfft=spaceComm_band
234 : : end if
235 : 1158 : me=xmpi_comm_rank(spaceComm_w)
236 : 1158 : i_am_master=(me==master)
237 : 1158 : i_am_master_kpt=(xmpi_comm_rank(spaceComm_kpt)==master)
238 : 1158 : i_am_master_band=(xmpi_comm_rank(spaceComm_band)==master)
239 : 1158 : i_am_master_spfft=(xmpi_comm_rank(spaceComm_spinorfft)==master)
240 : 1158 : me_fft=xmpi_comm_rank(spaceComm_fft)
241 : 1158 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
242 : :
243 : : !GPU: only OpenMP offload variant supported
244 : : gpu_option=dtset%gpu_option
245 : : if (gpu_option/=ABI_GPU_OPENMP) gpu_option=ABI_GPU_DISABLED
246 : : #ifndef HAVE_OPENMP_OFFLOAD
247 : 1158 : gpu_option=ABI_GPU_DISABLED
248 : : #endif
249 : :
250 : : !Check wether we write in file or save the matrix elements in psinablapsi_out
251 : : !and if we need to compute the full matrix or only the diagonal part
252 : 1158 : iomode=dtset%iomode
253 : 1158 : diag_only = .false.
254 : 1158 : store_half_dipoles = .false.
255 [ + + ]: 1158 : if (present(psinablapsi_out)) then
256 : 1121 : iomode = NO_FILE_OUT
257 : 1121 : iomode_etsf_mpiio = .false.
258 : 1121 : nc_unlimited = .false.
259 : 1121 : bsize = size(psinablapsi_out(1,1,:,1))
260 [ + - ]: 1121 : if (bsize == mband) then
261 : 1121 : diag_only = .true.
262 : 1121 : store_half_dipoles = .true. !FB ??
263 [ # # ]: 0 : else if (bsize == mband*(mband+1)/2) then
264 : 0 : diag_only = .false.
265 : 0 : store_half_dipoles = .true.
266 [ # # ]: 0 : else if (bsize == mband**2) then
267 : 0 : diag_only = .false.
268 : 0 : store_half_dipoles = .false.
269 : : else
270 : 0 : msg = "Wrong dimensions of psinablapsi_out!"
271 : 0 : ABI_ERROR(msg)
272 : : end if
273 [ + + + + : 2906753 : psinablapsi_out = zero
+ + + + ]
274 : : end if
275 : :
276 : : !----------------------------------------------------------------------------------
277 : : !1- Opening of OPT file and header writing
278 : : !----------------------------------------------------------------------------------
279 : :
280 : : #if defined(HAVE_GPU_MARKERS)
281 : : call nvtxStartRange("header_write")
282 : : #endif
283 [ + + ]: 1158 : if (iomode /= NO_FILE_OUT) then
284 : : ! I/O mode is netCDF or Fortran
285 : : iomode=merge(IO_MODE_ETSF,IO_MODE_FORTRAN_MASTER,dtset%iomode==IO_MODE_ETSF)
286 : 37 : if (use_netcdf_forced) iomode=IO_MODE_ETSF
287 : :
288 : : ! (master proc only)
289 [ + + ]: 37 : if (i_am_master) then
290 : 13 : fformopt=610 ; if (compute_half_dipoles) fformopt=620
291 : : ! ====> NETCDF format
292 : : if (iomode==IO_MODE_ETSF) then
293 : : ! Open/create nc file
294 [ - + ]: 13 : NCF_CHECK(nctk_open_create(ncid,nctk_ncify(dtfil%fnameabo_app_opt),xmpi_comm_self))
295 : : ! Write header data
296 [ - + ]: 13 : NCF_CHECK(hdr%ncwrite(ncid,fformopt,nc_define=.true.))
297 : : ! Define dims and array for dipole variables
298 : 13 : nctk_arrays(1)%name="dipole_valence_valence"
299 : 13 : nctk_arrays(1)%dtype="dp"
300 : 13 : nc_unlimited=(use_netcdf_unlimited.and.(.not.(nctk_has_mpiio.and.use_netcdf_mpiio)))
301 : : if (nc_unlimited) then
302 : : nctkdim%name="unlimited_bands"
303 : : nctkdim%value=NF90_UNLIMITED
304 : : NCF_CHECK(nctk_def_dims(ncid,nctkdim))
305 : : nctk_arrays(1)%shape_str=&
306 : : & "complex,number_of_cartesian_directions,max_number_of_states,number_of_kpoints,number_of_spins,unlimited_bands"
307 : : else if (compute_half_dipoles) then
308 : 13 : nctkdim%name="max_number_of_state_pairs"
309 : 13 : nctkdim%value=(mband*(mband+1))/2
310 [ - + ]: 13 : NCF_CHECK(nctk_def_dims(ncid,nctkdim))
311 : : nctk_arrays(1)%shape_str=&
312 : 13 : & "complex,number_of_cartesian_directions,max_number_of_state_pairs,number_of_kpoints,number_of_spins"
313 : : else
314 : : nctk_arrays(1)%shape_str=&
315 : : & "complex,number_of_cartesian_directions,max_number_of_states,max_number_of_states,number_of_kpoints,number_of_spins"
316 : : end if
317 [ - + ]: 13 : NCF_CHECK(nctk_def_arrays(ncid, nctk_arrays))
318 [ - + ]: 13 : NCF_CHECK(nctk_set_atomic_units(ncid, "dipole_valence_valence"))
319 : : ! Write eigenvalues
320 [ - + ]: 13 : NCF_CHECK(nctk_set_datamode(ncid))
321 : 13 : varid=nctk_idname(ncid,"eigenvalues")
322 [ + + - + : 52 : NCF_CHECK(nf90_put_var(ncid,varid,reshape(eigen0,[mband,nkpt,nsppol])))
- - ]
323 : : !Close file here because the rest has possibly to be written with collective I/O
324 [ - + ]: 13 : NCF_CHECK(nf90_close(ncid))
325 : : ! ====> Standard FORTRAN binary format
326 : : else if (iomode==IO_MODE_FORTRAN_MASTER) then
327 : : if (open_file(dtfil%fnameabo_app_opt,msg,newunit=ount,form="unformatted",status="unknown")/= 0) then
328 : : ABI_ERROR(msg)
329 : : end if
330 : : call hdr%fort_write(ount,fformopt,ierr,rewind=.true.)
331 : : write(ount)(eigen0(ib),ib=1,mband*nkpt*nsppol)
332 : : else
333 : : msg = "Wrong OPT file format!"
334 : : ABI_BUG(msg)
335 : : end if ! File format
336 : : end if ! master node
337 : 37 : call xmpi_bcast(iomode,master,spaceComm_w,ierr) ! Seems mandatory; why ?
338 [ + - - + ]: 37 : iomode_etsf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
339 : 37 : nc_unlimited=(iomode==IO_MODE_ETSF.and.use_netcdf_unlimited.and.(.not.iomode_etsf_mpiio)) ! UNLIMITED not compatible with mpi-io
340 : 37 : store_half_dipoles=(compute_half_dipoles.and.(.not.nc_unlimited))
341 : : end if
342 : : #if defined(HAVE_GPU_MARKERS)
343 : : call nvtxEndRange()
344 : : #endif
345 : :
346 : : !----------------------------------------------------------------------------------
347 : : !2- Computation of on-site contribution: <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>
348 : : !----------------------------------------------------------------------------------
349 : :
350 : : #if defined(HAVE_GPU_MARKERS)
351 : : call nvtxStartRange("pawnabla_init")
352 : : #endif
353 [ + - - + ]: 1158 : already_has_nabla=all(pawtab(:)%has_nabla==2)
354 : 1158 : call pawnabla_init(mpsang,dtset%ntypat,pawrad,pawtab)
355 : :
356 : : !Get electronic temperature from dtset
357 [ - + - - : 1158 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
- - ]
358 : : #if defined(HAVE_GPU_MARKERS)
359 : : call nvtxEndRange()
360 : : #endif
361 : :
362 : : !Compute spin-orbit contributions if necessary
363 [ + + ]: 1158 : if (dtset%pawspnorb==1) then
364 : 5 : option_core=0
365 : : call pawnabla_soc_init(el_temp,phisocphj,option_core,dtset%ixc,mpi_enreg%my_natom,natom,&
366 : : & dtset%nspden,dtset%ntypat,pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,&
367 : : & dtset%spnorbscl,dtset%typat,dtset%xc_denpos,dtset%xc_taupos,znucl,&
368 : 5 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
369 : : else
370 [ - + ]: 1153 : ABI_MALLOC(phisocphj,(0))
371 : : end if
372 : :
373 : : !----------------------------------------------------------------------------------
374 : : !3- Computation of <psi_n|-i.nabla|psi_m> for each k
375 : : !----------------------------------------------------------------------------------
376 : :
377 : : #if defined(HAVE_GPU_MARKERS)
378 : : call nvtxStartRange("prepare_netcdf")
379 : : #endif
380 : : !Prepare valence-valence dipoles writing
381 : : !In case of netCDF access to OPT file, prepare collective I/O
382 [ + + ]: 1158 : if (iomode == IO_MODE_ETSF) then
383 [ + - ]: 37 : if (iomode_etsf_mpiio) then
384 [ + + ]: 37 : if (i_am_master_spfft) then
385 [ - + ]: 23 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt),spaceComm_band))
386 : 23 : varid=nctk_idname(ncid,"dipole_valence_valence")
387 [ + + ]: 23 : if (xmpi_comm_size(spaceComm_w)>1) then
388 [ - + ]: 18 : NCF_CHECK(nctk_set_collective(ncid,varid))
389 : : end if
390 [ - + ]: 23 : NCF_CHECK(nctk_set_datamode(ncid))
391 : : end if
392 [ # # ]: 0 : else if (i_am_master) then
393 [ # # ]: 0 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt),xmpi_comm_self))
394 : 0 : varid=nctk_idname(ncid,"dipole_valence_valence")
395 : : if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
396 : : NCF_CHECK(nctk_set_collective(ncid,varid))
397 : : end if
398 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
399 : : end if
400 : : end if
401 [ + + ]: 1158 : if (iomode /= NO_FILE_OUT) then
402 [ + - ]: 37 : if (iomode_etsf_mpiio) then
403 : : !If MPI-IO, bands are treated by blocks of nbj_io columns: store only these columns
404 : : !psinablapsi holds the block written by current proc, psinablapsi_w is a work buffer
405 : 37 : nproc_band_io=xmpi_comm_size(spaceComm_band)
406 : 37 : me_band_io=xmpi_comm_rank(spaceComm_band)
407 : 37 : nbj_io=int(dipoles_blk_mem/(96._dp*real(mband,dp)))
408 : 37 : nbj_io=max(1,min(nbj_io,mband))
409 [ - + + - : 111 : ABI_MALLOC(psinablapsi,(2,3,nbj_io*mband))
+ - + - -
+ - + ]
410 [ - + + - : 74 : ABI_MALLOC(psinablapsi_w,(2,3,nbj_io*mband))
- + - + ]
411 : : else
412 : : !If not, store all (ib,jb) pairs (or half)
413 [ # # ]: 0 : bsize=mband**2 ; if (store_half_dipoles) bsize=(mband*(mband+1))/2
414 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,bsize))
# # # # #
# # # ]
415 [ # # # # : 0 : psinablapsi=zero
# # ]
416 : : end if
417 [ + + ]: 148 : pnp_size=size(psinablapsi)
418 : : end if
419 : : #if defined(HAVE_GPU_MARKERS)
420 : : call nvtxEndRange()
421 : : #endif
422 : : !Note: psinablapsi accumulates both plane-wave and on-site (PAW+SOC) contributions
423 : :
424 : : #if defined(HAVE_GPU_MARKERS)
425 : : call nvtxStartRange("alloc")
426 : : #endif
427 : : !Determine if cprj datastructure is distributed over bands
428 : 1158 : mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
429 : 1158 : cprj_paral_band=(mband_cprj<mband)
430 : :
431 : : !Store on-site matrix elements for each atom in plain arrays (for optics_paw_dipoles_k)
432 [ - + + - : 3474 : ABI_MALLOC(lmn_at,(natom))
+ - + - -
+ - + ]
433 [ + + ]: 3496 : do iatom=1,natom
434 : 3496 : lmn_at(iatom)=pawtab(dtset%typat(iatom))%lmn_size
435 : : end do
436 [ + + ]: 3496 : lmnmax=maxval(lmn_at)
437 [ - + + - : 5790 : ABI_MALLOC(nabla_at,(3,lmnmax,lmnmax,natom))
- + + - -
+ + - + -
+ - - + -
+ ]
438 [ - + + - : 5795 : ABI_MALLOC(soc_at,(2,2,3,lmnmax,lmnmax,merge(natom,1,dtset%pawspnorb==1)))
- + + - +
+ - + + -
+ - + - -
+ - + ]
439 [ + + + + : 2287693 : nabla_at=zero ; soc_at=zero
+ + + + +
+ + + + +
+ + + + +
+ ]
440 [ + + ]: 3496 : do iatom=1,natom
441 : 2338 : lmn_size=lmn_at(iatom)
442 [ + + + + : 620218 : nabla_at(:,1:lmn_size,1:lmn_size,iatom)=pawtab(dtset%typat(iatom))%nabla_ij(:,1:lmn_size,1:lmn_size)
+ + ]
443 [ + + ]: 3496 : if (dtset%pawspnorb==1) then
444 [ + + + + : 28340 : soc_at(:,:,:,1:lmn_size,1:lmn_size,iatom)=phisocphj(iatom)%value(:,:,:,1:lmn_size,1:lmn_size)
+ + + + +
+ ]
445 : : end if
446 : : end do
447 : : #ifdef HAVE_OPENMP_OFFLOAD
448 : : !$OMP TARGET ENTER DATA MAP(to:lmn_at,nabla_at,soc_at) IF(gpu_option==ABI_GPU_OPENMP)
449 : : #endif
450 : : #if defined(HAVE_GPU_MARKERS)
451 : : call nvtxEndRange()
452 : : #endif
453 : :
454 : : !LOOP OVER SPINS
455 : 1158 : ibg=0;icg=0
456 : 1158 : bdtot_index=0
457 [ + + ]: 2316 : do isppol=1,nsppol
458 : :
459 : : ! LOOP OVER k POINTS
460 : 1158 : ikg=0
461 [ + + ]: 38443 : do ikpt=1,nkpt
462 : : #if defined(HAVE_GPU_MARKERS)
463 : : call nvtxStartRange("k-point step")
464 : : #endif
465 : :
466 [ + + + - ]: 36127 : if (iomode == NO_FILE_OUT) psinablapsi => psinablapsi_out(:,:,:,ikpt)
467 : :
468 : 36127 : etiq=ikpt+(isppol-1)*nkpt
469 : 36127 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
470 [ + + ]: 326207 : master_spfftband=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
471 : :
472 : : ! Select k-points for current proc
473 : 36127 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt))
474 [ + + ]: 36127 : if (mykpt) then
475 : :
476 : : ! Data depending on k-point
477 : 16493 : npw_k=npwarr(ikpt)
478 : 16493 : istwf_k=dtset%istwfk(ikpt)
479 [ + + ]: 16493 : cplex=2;if (istwf_k>1) cplex=1
480 [ + + ]: 65972 : kpoint(:)=dtset%kptns(:,ikpt)
481 : :
482 : : #if defined(HAVE_GPU_MARKERS)
483 : : call nvtxStartRange("pawcprj_alloc")
484 : : #endif
485 : : ! Extract cprj for this k-point
486 [ + + ]: 16493 : nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
487 [ + + ]: 16493 : if (mkmem*nsppol/=1) then
488 : 16481 : iorder_cprj=0
489 [ - + + - : 383696 : ABI_MALLOC(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
- + + - +
- + - - +
- + + + +
+ ]
490 : 16481 : call pawcprj_alloc(cprj_k_loc,0,dimcprj)
491 : : call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
492 : : & mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
493 [ - + ]: 16481 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
494 : : else
495 : 12 : cprj_k_loc => cprj
496 : : end if
497 : :
498 : : ! if cprj are distributed over bands, gather them (because we need to mix bands)
499 [ + + ]: 16493 : if (cprj_paral_band) then
500 [ - + + - : 187016 : ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k))
- + + - +
- + - - +
- + + + +
+ ]
501 : 6616 : call pawcprj_alloc(cprj_k,0,dimcprj)
502 : : call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom, &
503 : : & my_nspinor*nband_cprj_k,my_nspinor*mpi_enreg%bandpp,&
504 : 6616 : & dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
505 : : else
506 : 9877 : cprj_k => cprj_k_loc
507 : : end if
508 : :
509 : : ! Pack cprj into a contiguous array: cprj_pk(2,nproj,my_nspinor*nband_k)
510 : : ! (projectors of all atoms concatenated, in cprj atom order)
511 [ - + + - : 49479 : ABI_MALLOC(nlmn,(natom))
+ - + - -
+ - + -
+ ]
512 [ + + ]: 49501 : nlmn(:)=cprj_k(:,1)%nlmn
513 [ + + ]: 49501 : nproj=sum(nlmn)
514 [ - + + - : 49479 : ABI_MALLOC(cprj_pk,(2,nproj*my_nspinor*nband_k))
+ - + - -
+ - + -
+ ]
515 : 16493 : call pawcprj_pack(nlmn,cprj_k(:,1:my_nspinor*nband_k),cprj_pk)
516 : : ! Offset of each atom in cprj_pk and atom of each projector
517 [ - + + - : 32986 : ABI_MALLOC(atom_off,(natom))
- + - + ]
518 [ - + + - : 49479 : ABI_MALLOC(row_atom,(nproj))
+ - + - -
+ - + ]
519 : 49501 : iproj=0
520 [ + + ]: 49501 : do iatom=1,natom
521 : 33008 : atom_off(iatom)=iproj
522 [ + + ]: 297088 : row_atom(iproj+1:iproj+nlmn(iatom))=iatom
523 : 49501 : iproj=iproj+nlmn(iatom)
524 : : end do
525 : :
526 : : ! Compute k+G in cartesian coordinates
527 [ - + + - : 49479 : ABI_MALLOC(kpg_k,(3,npw_k*dtset%nspinor))
+ - + - -
+ - + ]
528 : 16493 : kg_k => kg(:,1+ikg:npw_k+ikg)
529 [ + + ]: 2330285 : do ipw=1,npw_k
530 : : kpg_k(1,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1) &
531 : : & +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2) &
532 : 2313792 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3)
533 : : kpg_k(2,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1) &
534 : : & +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2) &
535 : 2313792 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3)
536 : : kpg_k(3,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1) &
537 : : & +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2) &
538 : 2330285 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3)
539 : : end do
540 [ + + + + ]: 9271661 : kpg_k(:,1:npw_k)=two_pi*kpg_k(:,1:npw_k)
541 [ + + + + : 21445 : if (dtset%nspinor==2) kpg_k(1:3,npw_k+1:2*npw_k)=kpg_k(1:3,1:npw_k)
+ + ]
542 : : #if defined(HAVE_GPU_MARKERS)
543 : : call nvtxEndRange()
544 : : #endif
545 : :
546 : : ! Transfer data of this k-point on GPU
547 : : #ifdef HAVE_OPENMP_OFFLOAD
548 : : !$OMP TARGET ENTER DATA MAP(to:cg(:,icg+1:icg+npw_k*my_nspinor*nband_k)) IF(gpu_option==ABI_GPU_OPENMP)
549 : : !$OMP TARGET ENTER DATA MAP(to:cprj_pk,kpg_k,atom_off,row_atom) IF(gpu_option==ABI_GPU_OPENMP)
550 : : #endif
551 : :
552 : : ! Select bands for current proc
553 : : ! jpw1:jpw2 = bands jb for which the plane-wave contribution is computed
554 : : ! jpaw1:jpaw2 = bands jb for which the on-site (PAW+SOC) contribution is computed
555 : 16493 : jpw1=1 ; jpw2=nband_k ; jpaw1=1 ; jpaw2=nband_k
556 [ + + ]: 16493 : if (mpi_enreg%paral_kgb==1) then
557 : : ! cprj are known for all bands: distribute contiguous blocks of bands
558 : 6716 : jpaw1=(mpi_enreg%me_band*nband_k)/mpi_enreg%nproc_band+1
559 : 6716 : jpaw2=((mpi_enreg%me_band+1)*nband_k)/mpi_enreg%nproc_band
560 : : else if (xmpi_paral==1) then
561 : : ! Bands of this k-point treated by current proc (contiguous, see distrb2)
562 : 9777 : jpaw1=nband_k+1 ; jpaw2=0
563 [ + + ]: 88081 : do jb=1,nband_k
564 [ + - ]: 88081 : if (mpi_enreg%proc_distrb(ikpt,jb,isppol)==mpi_enreg%me_kpt) then
565 : 78304 : jpaw1=min(jpaw1,jb) ; jpaw2=max(jpaw2,jb)
566 : : end if
567 : : end do
568 [ + + + - : 88081 : if (count(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)==mpi_enreg%me_kpt)/=max(0,jpaw2-jpaw1+1)) then
- + ]
569 : 0 : ABI_BUG('Bands of a k-point are not distributed in contiguous blocks!')
570 : : end if
571 [ + + ]: 9777 : if (.not.iomode_etsf_mpiio) then
572 : 9760 : jpw1=jpaw1 ; jpw2=jpaw2
573 : : end if
574 : : end if
575 : : ! On-site contribution is the same on all FFT procs: only one of them adds it
576 [ + + ]: 16493 : if (me_fft/=0) then
577 : 94 : jpaw1=1 ; jpaw2=0
578 : : end if
579 : :
580 [ + + ]: 16493 : if (iomode_etsf_mpiio) then
581 : :
582 : : #if defined(HAVE_GPU_MARKERS)
583 : : call nvtxStartRange("jband_loop_mpiio",5)
584 : : #endif
585 : : ! MPI-IO: bands are treated by super-blocks made of one block of nbj_io columns
586 : : ! per band proc. Each block is computed by all procs and reduced to its owner;
587 : : ! then all owners write their block with a single collective call.
588 : : ! Columns j0:j1 are contiguous in the OPT file, for both full and half storage.
589 [ - + - - : 410 : do j0s=1,nband_k,nbj_io*nproc_band_io
- + - + ]
590 : 205 : my_j0=1 ; my_j1=0
591 : : #if defined(HAVE_GPU_MARKERS)
592 : : call nvtxStartRange("mpiio_blk",2)
593 : : #endif
594 [ + + ]: 410 : do iblk=0,nproc_band_io-1
595 : 293 : j0=j0s+iblk*nbj_io
596 [ + + ]: 293 : if (j0>nband_k) exit
597 : 205 : j1=min(nband_k,j0+nbj_io-1)
598 [ + + ]: 205 : if (iblk==me_band_io) then
599 : 151 : psi_blk => psinablapsi ; my_j0=j0 ; my_j1=j1
600 : : else
601 : 54 : psi_blk => psinablapsi_w
602 : : end if
603 : :
604 : : #if defined(HAVE_GPU_MARKERS)
605 : : call nvtxStartRange("optics_paw_dipoles",6)
606 : : #endif
607 : : call optics_paw_dipoles_k(cg(:,icg+1:icg+npw_k*my_nspinor*nband_k),cprj_pk,kpg_k,&
608 : : & nabla_at,soc_at,lmn_at,atom_off,row_atom,cplex,istwf_k,j0,j1,jpw1,jpw2,jpaw1,jpaw2,&
609 : : & lmnmax,mband,my_nspinor,natom,nband_k,nproj,npw_k,dtset%pawspnorb,&
610 : 205 : & (.not.store_half_dipoles),diag_only,store_half_dipoles,gpu_option,psi_blk)
611 : :
612 : : #if defined(HAVE_GPU_MARKERS)
613 : : call nvtxEndRange()
614 : : call nvtxStartRange("xmpi_sum",17)
615 : : #endif
616 : : ! Reduction of the block on its owner
617 [ - + ]: 205 : nent=merge((j1*(j1+1))/2-(j0*(j0-1))/2,(j1-j0+1)*mband,store_half_dipoles)
618 : 205 : call timab(48,1,tsec)
619 [ + + ]: 205 : if (mpi_enreg%paral_kgb==1) then
620 : 188 : call xmpi_sum_master(psi_blk(:,:,1:nent),master,spaceComm_spinorfft,ierr)
621 : : end if
622 [ + + ]: 205 : if (i_am_master_spfft) then
623 : 111 : call xmpi_sum_master(psi_blk(:,:,1:nent),iblk,spaceComm_band,ierr)
624 : : end if
625 : 410 : call timab(48,2,tsec)
626 : : #if defined(HAVE_GPU_MARKERS)
627 : : call nvtxEndRange()
628 : : #endif
629 : : end do ! iblk
630 : : #if defined(HAVE_GPU_MARKERS)
631 : : call nvtxEndRange()
632 : : #endif
633 : :
634 : : #if defined(HAVE_GPU_MARKERS)
635 : : call nvtxStartRange("mpiio_write")
636 : : #endif
637 : : ! Write to OPT file: each owner writes its block
638 [ + + ]: 410 : if (i_am_master_spfft) then
639 [ + - ]: 111 : if (store_half_dipoles) then
640 [ + + ]: 666 : nc_start_5=[1,1,(my_j0*(my_j0-1))/2+1,ikpt,isppol];nc_stride_5=[1,1,1,1,1]
641 : 111 : nc_count_5=[0,0,0,0,0]
642 [ + + + + ]: 456 : if (my_j1>=my_j0) nc_count_5=[2,3,(my_j1*(my_j1+1))/2-(my_j0*(my_j0-1))/2,1,1]
643 [ - + ]: 111 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
644 : : else
645 [ # # ]: 0 : nc_start_6=[1,1,1,my_j0,ikpt,isppol];nc_stride_6=[1,1,1,1,1,1]
646 : 0 : nc_count_6=[0,0,0,0,0,0]
647 [ # # # # ]: 0 : if (my_j1>=my_j0) nc_count_6=[2,3,mband,my_j1-my_j0+1,1,1]
648 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
649 : : end if
650 : : end if
651 : : #if defined(HAVE_GPU_MARKERS)
652 : : call nvtxEndRange()
653 : : #endif
654 : :
655 : : end do ! j0s
656 : : #if defined(HAVE_GPU_MARKERS)
657 : : call nvtxEndRange()
658 : : #endif
659 : :
660 : : else
661 : :
662 : : #if defined(HAVE_GPU_MARKERS)
663 : : call nvtxStartRange("optics_paw_dipoles_no_mpiio",2)
664 : : #endif
665 : : ! No MPI-IO: compute all (ib,jb) pairs at once
666 : : call optics_paw_dipoles_k(cg(:,icg+1:icg+npw_k*my_nspinor*nband_k),cprj_pk,kpg_k,&
667 : : & nabla_at,soc_at,lmn_at,atom_off,row_atom,cplex,istwf_k,1,nband_k,jpw1,jpw2,jpaw1,jpaw2,&
668 : : & lmnmax,mband,my_nspinor,natom,nband_k,nproj,npw_k,dtset%pawspnorb,.false.,&
669 : 16288 : & diag_only,store_half_dipoles,gpu_option,psinablapsi)
670 : : #if defined(HAVE_GPU_MARKERS)
671 : : call nvtxEndRange()
672 : : #endif
673 : :
674 : : end if
675 : :
676 : : ! Release data of this k-point on GPU
677 : : #ifdef HAVE_OPENMP_OFFLOAD
678 : : !$OMP TARGET EXIT DATA MAP(delete:cg(:,icg+1:icg+npw_k*my_nspinor*nband_k)) IF(gpu_option==ABI_GPU_OPENMP)
679 : : !$OMP TARGET EXIT DATA MAP(delete:cprj_pk,kpg_k,atom_off,row_atom) IF(gpu_option==ABI_GPU_OPENMP)
680 : : #endif
681 : :
682 [ + - ]: 16493 : if (mkmem/=0) then
683 : 16493 : ibg = ibg + my_nspinor*nband_cprj_k
684 : 16493 : icg = icg + npw_k*my_nspinor*nband_k
685 : 16493 : ikg = ikg + npw_k
686 : : end if
687 : :
688 : : #if defined(HAVE_GPU_MARKERS)
689 : : call nvtxStartRange("free")
690 : : #endif
691 [ - + ]: 16493 : ABI_FREE(cprj_pk)
692 [ - + ]: 16493 : ABI_FREE(nlmn)
693 [ - + ]: 16493 : ABI_FREE(atom_off)
694 [ - + ]: 16493 : ABI_FREE(row_atom)
695 [ + + ]: 16493 : if (cprj_paral_band) then
696 : 6616 : call pawcprj_free(cprj_k)
697 [ - + ]: 6616 : ABI_FREE(cprj_k)
698 : : end if
699 [ + + ]: 16493 : if (mkmem*nsppol/=1) then
700 : 16481 : call pawcprj_free(cprj_k_loc)
701 [ - + ]: 16481 : ABI_FREE(cprj_k_loc)
702 : : end if
703 [ - + ]: 16493 : ABI_FREE(kpg_k)
704 : : #if defined(HAVE_GPU_MARKERS)
705 : : call nvtxEndRange()
706 : : #endif
707 : :
708 : : ! Write to OPT file if not MPI-IO
709 [ + + ]: 16493 : if (iomode /= NO_FILE_OUT) then
710 : : #if defined(HAVE_GPU_MARKERS)
711 : : call nvtxStartRange("write")
712 : : #endif
713 : : ! >>> Reduction in case of parallelism
714 [ - + ]: 205 : if (.not.iomode_etsf_mpiio) then
715 : 0 : call timab(48,1,tsec)
716 : 0 : call xmpi_sum_master(psinablapsi,master,spaceComm_bandspinorfft,ierr)
717 : 0 : call timab(48,2,tsec)
718 : : end if
719 : :
720 : : ! >>> This my kpt and I am the master node: I write the data
721 : : if (.not.iomode_etsf_mpiio) then
722 [ # # ]: 0 : if (i_am_master) then
723 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
724 : : if (nc_unlimited) then
725 : : nc_start_6=[1,1,1,ikpt,isppol,1] ; nc_count_6=[2,3,mband,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1]
726 : : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
727 [ # # ]: 0 : else if (.not.store_half_dipoles) then
728 [ # # # # ]: 0 : nc_start_6=[1,1,1,1,ikpt,isppol] ; nc_count_6=[2,3,mband,mband,1,1] ; nc_stride_6=[1,1,1,1,1,1]
729 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
730 : : else
731 [ # # # # ]: 0 : nc_start_5=[1,1,1,ikpt,isppol] ; nc_count_5=[2,3,(mband*(mband+1))/2,1,1] ; nc_stride_5=[1,1,1,1,1]
732 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
733 : : end if
734 : : else
735 [ # # ]: 0 : bsize=nband_k**2;if (store_half_dipoles) bsize=(nband_k*(nband_k+1))/2
736 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,1,ib),ib=1,bsize)
737 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,2,ib),ib=1,bsize)
738 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,3,ib),ib=1,bsize)
739 : : end if
740 : :
741 : : ! >>> This my kpt and I am not the master node: I send the data
742 [ # # ]: 0 : else if (i_am_master_band.and.i_am_master_spfft) then
743 [ # # ]: 0 : if (mpi_enreg%me_kpt/=master_spfftband) then
744 : 0 : ABI_BUG('Problem with band communicator!')
745 : : end if
746 : 0 : call xmpi_exch(psinablapsi,pnp_size,mpi_enreg%me_kpt,psinablapsi,master,spaceComm_kpt,etiq,ierr)
747 : : end if
748 : : end if
749 : : #if defined(HAVE_GPU_MARKERS)
750 : : call nvtxEndRange()
751 : : #endif
752 : : else
753 : : ! >>> Reduction in case of parallelism
754 : 16288 : call timab(48,1,tsec)
755 : 16288 : call xmpi_sum(psinablapsi,spaceComm_bandspinorfft,ierr)
756 : 16288 : call timab(48,2,tsec)
757 : : end if ! no_file_out
758 : :
759 : : ! >>> This is not my kpt and I am the master node: I receive the data and I write
760 [ + + + - ]: 19634 : elseif ((.not.iomode_etsf_mpiio).and.i_am_master.and.(iomode/=NO_FILE_OUT)) then ! mykpt
761 : 0 : sender=master_spfftband
762 : 0 : call xmpi_exch(psinablapsi,pnp_size,sender,psinablapsi,master,spaceComm_kpt,etiq,ierr)
763 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
764 : : if (nc_unlimited) then
765 : : nc_start_6=[1,1,1,ikpt,isppol,1] ; nc_count_6=[2,3,mband,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1]
766 : : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
767 [ # # ]: 0 : else if (.not.store_half_dipoles) then
768 [ # # # # ]: 0 : nc_start_6=[1,1,1,1,ikpt,isppol] ; nc_count_6=[2,3,mband,mband,1,1] ; nc_stride_6=[1,1,1,1,1,1]
769 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
770 : : else
771 [ # # # # ]: 0 : nc_start_5=[1,1,1,ikpt,isppol] ; nc_count_5=[2,3,(mband*(mband+1))/2,1,1] ; nc_stride_5=[1,1,1,1,1]
772 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
773 : : end if
774 : : else
775 [ # # ]: 0 : bsize=nband_k**2;if (store_half_dipoles) bsize=(nband_k*(nband_k+1))/2
776 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,1,ib),ib=1,bsize)
777 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,2,ib),ib=1,bsize)
778 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,3,ib),ib=1,bsize)
779 : : end if
780 : : end if ! mykpt
781 : :
782 : 37285 : bdtot_index=bdtot_index+nband_k
783 : :
784 : : ! End loop on spin,kpt
785 : : #if defined(HAVE_GPU_MARKERS)
786 : : call nvtxEndRange()
787 : : #endif
788 : : end do ! ikpt
789 : : end do !isppol
790 : :
791 : : ! >>> Last reduction over k-points in case of parallelism
792 : : ! >>> if no output to file
793 [ + + ]: 1158 : if (iomode == NO_FILE_OUT) then
794 : 1121 : call xmpi_sum(psinablapsi_out,spaceComm_kpt,ierr)
795 : : end if
796 : :
797 : : !Close file
798 [ + + + + : 1158 : if ((i_am_master.or.(iomode_etsf_mpiio.and.i_am_master_spfft)).and.iomode /= NO_FILE_OUT) then
+ + ]
799 [ + - ]: 23 : if (iomode==IO_MODE_ETSF) then
800 [ - + ]: 23 : NCF_CHECK(nf90_close(ncid))
801 : : else
802 : 0 : ierr=close_unit(ount,msg)
803 [ # # ]: 0 : ABI_CHECK(ierr==0,"Error while closing OPT file")
804 : : end if
805 : : end if
806 : :
807 : : !Datastructures deallocations
808 [ + + ]: 1158 : if (iomode /= NO_FILE_OUT) then
809 [ - + ]: 37 : ABI_FREE(psinablapsi)
810 : : end if
811 [ + + ]: 1158 : if (iomode_etsf_mpiio) then
812 [ - + ]: 37 : ABI_FREE(psinablapsi_w)
813 : : end if
814 [ + + ]: 1158 : if (dtset%pawspnorb==1) then
815 [ + + ]: 25 : do iatom=1,natom
816 [ - + ]: 25 : ABI_FREE(phisocphj(iatom)%value)
817 : : end do
818 : : end if
819 [ + - + + : 1178 : ABI_FREE(phisocphj)
- + - + ]
820 : : #ifdef HAVE_OPENMP_OFFLOAD
821 : : !$OMP TARGET EXIT DATA MAP(delete:lmn_at,nabla_at,soc_at) IF(gpu_option==ABI_GPU_OPENMP)
822 : : #endif
823 [ - + ]: 1158 : ABI_FREE(lmn_at)
824 [ - + ]: 1158 : ABI_FREE(nabla_at)
825 [ - + ]: 1158 : ABI_FREE(soc_at)
826 [ + - ]: 1158 : if (.not.already_has_nabla) then
827 [ + + ]: 2318 : do itypat=1,dtset%ntypat
828 [ + - ]: 2318 : if (allocated(pawtab(itypat)%nabla_ij)) then
829 : 1160 : ABI_FREE(pawtab(itypat)%nabla_ij)
830 : 1160 : pawtab(itypat)%has_nabla=0
831 : : end if
832 : : end do
833 : : end if
834 : :
835 : : ABI_NVTX_END_RANGE()
836 : 1158 : call timab(1184,2,tsec)
837 : : DBG_EXIT("COLL")
838 : :
839 [ - + - + ]: 2316 : end subroutine optics_paw
840 : : !!***
841 : :
842 : : !----------------------------------------------------------------------
843 : :
844 : : !!****f* m_paw_optics/optics_paw_dipoles_k
845 : : !! NAME
846 : : !! optics_paw_dipoles_k
847 : : !!
848 : : !! FUNCTION
849 : : !! Compute, for one k-point and a range of bands jb1:jb2, the matrix elements
850 : : !! <Psi_ib|-i.nabla|Psi_jb> as the sum of:
851 : : !! - plane-wave part: <Psi_tild_ib|-i.nabla|Psi_tild_jb>
852 : : !! - on-site part: <Psi_ib|p_i><p_j|Psi_jb> (<phi_i|-i.nabla|phi_j>-<tphi_i|-i.nabla|tphi_j>)
853 : : !! + spin-orbit part (if pawspnorb==1)
854 : : !! No MPI communication and no I/O are done here.
855 : : !!
856 : : !! Both parts are sesquilinear forms X^H.D.X computed with GEMMs over blocks of columns jb:
857 : : !! - plane-wave: X=C (plane-wave coefficients), D=diag(k+G)_idir
858 : : !! - on-site: X=P (packed cprj), D=K_idir, block-diagonal over atoms, with
859 : : !! K_idir^{ss'}_ij = -i.nabla_ij.delta_ss' + (Sigma^ss' X g_ij)_idir
860 : : !!
861 : : !! INPUTS
862 : : !! cg_k(2,npw_k*my_nspinor,nband_k)=planewave coefficients of wavefunctions at this k-point
863 : : !! cprj_pk(2,nproj*my_nspinor,nband_k)= <p_lmn|Cnk> coefficients at this k-point (all bands),
864 : : !! packed with pawcprj_pack (projectors of all atoms concatenated, in cprj atom order)
865 : : !! kpg_k(3,npw_k*my_nspinor)=k+G vectors in cartesian coordinates (times 2pi)
866 : : !! nabla_at(3,lmnmax,lmnmax,natom)=<phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for each atom
867 : : !! soc_at(2,2,3,lmnmax,lmnmax,:)=spin-orbit on-site matrix elements for each atom (if pawspnorb==1)
868 : : !! lmn_at(natom)=number of (l,m,n) channels of each atom
869 : : !! atom_off(natom)=offset of each atom in the projector dimension of cprj_pk
870 : : !! row_atom(nproj)=atom corresponding to each projector of cprj_pk
871 : : !! cplex=1 if cprj are real (istwfk>1), 2 otherwise
872 : : !! istwf_k=storage mode of wavefunctions at this k-point
873 : : !! jb1,jb2=range of bands jb to compute
874 : : !! jpw1,jpw2=range of bands jb for which the plane-wave part is computed by this proc
875 : : !! jpaw1,jpaw2=range of bands jb for which the on-site part is computed by this proc
876 : : !! lmnmax=max. number of (l,m,n) channels over atoms
877 : : !! mband=maximum number of bands
878 : : !! my_nspinor=number of spinorial components treated by this proc
879 : : !! natom=number of atoms in cell
880 : : !! nband_k=number of bands at this k-point
881 : : !! nproj=total number of projectors (per spinor)
882 : : !! npw_k=number of planewaves at this k-point
883 : : !! pawspnorb=1 if spin-orbit contribution is included
884 : : !! all_rows=if true, compute all ib for each jb, without filling the second half of the matrix
885 : : !! (full storage only, used for MPI-IO); if false, compute only ib<=jb
886 : : !! diag_only=if true, compute only ib=jb elements (jb1 must be 1)
887 : : !! store_half=if true, use packed storage for the (ib<=jb) half matrix
888 : : !! Note: psinablapsi only stores columns jb1:jb2: element (ib,jb) is stored at index
889 : : !! jb(jb-1)/2+ib-ishift (store_half) or (jb-1)*mband+ib-ishift (full storage),
890 : : !! where ishift is the index of the last element before column jb1.
891 : : !! The second half of the full matrix is filled only if jb1=1.
892 : : !! gpu_option=ABI_GPU_OPENMP to run on GPU (all input arrays must then be mapped on device)
893 : : !!
894 : : !! OUTPUT
895 : : !! psinablapsi(2,3,:)=sum of plane-wave and on-site (PAW+SOC) contributions
896 : : !!
897 : : !! SOURCE
898 : :
899 : 49479 : subroutine optics_paw_dipoles_k(cg_k,cprj_pk,kpg_k,nabla_at,soc_at,lmn_at,atom_off,row_atom,&
900 : : & cplex,istwf_k,jb1,jb2,jpw1,jpw2,jpaw1,jpaw2,lmnmax,mband,my_nspinor,natom,nband_k,&
901 [ - + - + ]: 16493 : & nproj,npw_k,pawspnorb,all_rows,diag_only,store_half,gpu_option,psinablapsi)
902 : :
903 : : !Arguments ------------------------------------
904 : : !scalars
905 : : integer,intent(in) :: cplex,gpu_option,istwf_k,jb1,jb2,jpw1,jpw2,jpaw1,jpaw2,lmnmax,mband
906 : : integer,intent(in) :: my_nspinor,natom,nband_k,nproj,npw_k,pawspnorb
907 : : logical,intent(in) :: all_rows,diag_only,store_half
908 : : !arrays
909 : : integer,intent(in) :: atom_off(natom),lmn_at(natom),row_atom(nproj)
910 : : real(dp),intent(in) :: cg_k(2,npw_k*my_nspinor,nband_k),cprj_pk(2,nproj*my_nspinor,nband_k)
911 : : real(dp),intent(in) :: kpg_k(3,npw_k*my_nspinor),nabla_at(3,lmnmax,lmnmax,natom)
912 : : real(dp),intent(in) :: soc_at(:,:,:,:,:,:)
913 : : real(dp),intent(inout) :: psinablapsi(:,:,:)
914 : :
915 : : !Local variables-------------------------------
916 : : !scalars
917 : : integer :: iatom,ib,icol,idir,ii,ilmn,ioff,iproj0,ipw,irow,ispinor,j,j0,jb,jbhi,jblo,jlmn
918 : : #ifdef HAVE_OPENMP_OFFLOAD
919 : : integer :: ib,ii
920 : : #endif
921 : : integer :: ishift,lmn_size,nbj,nbj_blk,nbjmax,nprojsp,npwsp,nrows,nzero
922 : : logical :: gpu
923 : : real(dp) :: aa,bb,cc,nab,re1,im1,x1,y1,x2,y2
924 : : !arrays
925 : 16493 : real(dp),allocatable :: cprj_r(:,:),tpaw(:,:,:),wpw(:,:,:),zblk(:,:,:),zr(:,:,:)
926 : :
927 : : ! ************************************************************************
928 : :
929 [ - + - - : 16493 : ABI_CHECK(.not.(all_rows.and.(store_half.or.diag_only)),"all_rows requires full storage!")
- - ]
930 [ - + - - : 16493 : ABI_CHECK(jb1==1.or.all_rows.or.store_half,"Full matrix storage requires jb1=1!")
- - ]
931 [ - + - - ]: 16493 : ABI_CHECK(jb1==1.or.(.not.diag_only),"diag_only requires jb1=1!")
932 [ + + + - : 16493 : ABI_CHECK(pawspnorb==0.or.(cplex==2.and.my_nspinor==2),"SOC requires complex cprj and both spinors!")
- + ]
933 : :
934 : 16493 : gpu=(gpu_option==ABI_GPU_OPENMP)
935 : 16493 : npwsp=npw_k*my_nspinor
936 : 16493 : nprojsp=nproj*my_nspinor
937 : :
938 : : !Index of the last element before column jb1 in the storage, and number of stored elements
939 [ + + ]: 16493 : if (diag_only) then
940 : 16288 : ishift=0 ; nzero=nband_k
941 [ + - ]: 205 : else if (store_half) then
942 : 205 : ishift=(jb1*(jb1-1))/2 ; nzero=(jb2*(jb2+1))/2-ishift
943 : : else
944 : 0 : ishift=(jb1-1)*mband ; nzero=jb2*mband-ishift
945 : : end if
946 : :
947 : : !Fill output array with zeros
948 [ + - ]: 16493 : if (.not.gpu) then
949 [ + + + + : 1486973 : psinablapsi(:,:,1:nzero)=zero
+ + ]
950 : : else
951 : : #ifdef HAVE_OPENMP_OFFLOAD
952 : : !$OMP TARGET ENTER DATA MAP(alloc:psinablapsi)
953 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) MAP(to:psinablapsi) PRIVATE(ii,idir)
954 : : do ii=1,nzero
955 : : do idir=1,3
956 : : psinablapsi(1,idir,ii)=zero ; psinablapsi(2,idir,ii)=zero
957 : : end do
958 : : end do
959 : : #endif
960 : : end if
961 : :
962 : : !Number of columns jb treated at once (limited by the memory of work arrays)
963 : : !(each column jb needs 2x3 reals per plane wave, per band ib and per projector)
964 : 16493 : nbj_blk=int(dipoles_blk_mem/(48._dp*real(npwsp+nband_k+nprojsp,dp)))
965 : 16493 : nbj_blk=max(1,min(nbj_blk,jb2-jb1+1))
966 : :
967 : : !2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m>
968 : : !----------------------------------------------------------------------------------
969 : : !Note: <psi_nk|-i.nabla|psi_mk> => Sum_g[ <G|Psi_nk>^* <G|Psi_mk> G ]
970 : : ! computed as C^H.W with W(g,jb,idir)=(k+G)_idir*C(g,jb)
971 : : ! For istwfk>1: Sum_g>=g0[ 2i.Im(<G|Psi_nk>^* <G|Psi_mk> G) ]
972 : : ! G=k+g=0 term is included but do not contribute
973 : :
974 : 16493 : jblo=max(jb1,jpw1) ; jbhi=min(jb2,jpw2)
975 [ + - ]: 16493 : if (jbhi>=jblo) then
976 : 16493 : nbjmax=min(nbj_blk,jbhi-jblo+1)
977 [ - + + - : 65972 : ABI_MALLOC(wpw,(2,npwsp,3*nbjmax))
- + + - +
- + - - +
- + - + ]
978 [ - + + - : 65972 : ABI_MALLOC(zblk,(2,nband_k,3*nbjmax))
- + + - +
- + - - +
- + - + ]
979 [ - + + - : 65972 : ABI_MALLOC(zr,(1,nband_k,3*nbjmax))
- + + - +
- + - - +
- + - + ]
980 : : #ifdef HAVE_OPENMP_OFFLOAD
981 : : !$OMP TARGET ENTER DATA MAP(alloc:wpw,zblk,zr) IF(gpu)
982 : : #endif
983 : :
984 [ + - + - ]: 32986 : do j0=jblo,jbhi,nbjmax
985 : 16493 : nbj=min(nbjmax,jbhi-j0+1)
986 : : !If all_rows, compute all ib ; if not, compute only ib<=jb (trapezoidal block)
987 [ + - ]: 16493 : nrows=merge(nband_k,j0+nbj-1,all_rows)
988 : :
989 : : ! Build W=(k+G)_idir*C_jb (or (k+G)_idir*(Im C_jb,-Re C_jb) for istwfk>1)
990 [ + - ]: 16493 : if (.not.gpu) then
991 : : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(SHARED) PRIVATE(idir,j,jb,icol,ipw)
992 [ + + ]: 65972 : do idir=1,3
993 [ + + ]: 464396 : do j=1,nbj
994 : 398424 : jb=j0+j-1 ; icol=(idir-1)*nbj+j
995 [ + + ]: 447903 : if (istwf_k>1) then
996 [ + + ]: 23058 : do ipw=1,npwsp
997 : 22932 : wpw(1,ipw,icol)= kpg_k(idir,ipw)*cg_k(2,ipw,jb)
998 : 23058 : wpw(2,ipw,icol)=-kpg_k(idir,ipw)*cg_k(1,ipw,jb)
999 : : end do
1000 : : else
1001 [ + + ]: 56227110 : do ipw=1,npwsp
1002 : 55828812 : wpw(1,ipw,icol)=kpg_k(idir,ipw)*cg_k(1,ipw,jb)
1003 : 56227110 : wpw(2,ipw,icol)=kpg_k(idir,ipw)*cg_k(2,ipw,jb)
1004 : : end do
1005 : : end if
1006 : : end do
1007 : : end do
1008 : : !$OMP END PARALLEL DO
1009 : : else
1010 : : #ifdef HAVE_OPENMP_OFFLOAD
1011 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) MAP(to:wpw,kpg_k,cg_k) &
1012 : : !$OMP& PRIVATE(idir,j,jb,icol,ipw)
1013 : : do idir=1,3
1014 : : do j=1,nbj
1015 : : do ipw=1,npwsp
1016 : : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1017 : : if (istwf_k>1) then
1018 : : wpw(1,ipw,icol)= kpg_k(idir,ipw)*cg_k(2,ipw,jb)
1019 : : wpw(2,ipw,icol)=-kpg_k(idir,ipw)*cg_k(1,ipw,jb)
1020 : : else
1021 : : wpw(1,ipw,icol)=kpg_k(idir,ipw)*cg_k(1,ipw,jb)
1022 : : wpw(2,ipw,icol)=kpg_k(idir,ipw)*cg_k(2,ipw,jb)
1023 : : end if
1024 : : end do
1025 : : end do
1026 : : end do
1027 : : #endif
1028 : : end if
1029 : :
1030 : : ! Z=C^H.W (only diagonal elements if diag_only)
1031 [ + + ]: 16493 : if (diag_only) then
1032 [ + - ]: 16288 : if (.not.gpu) then
1033 : : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(SHARED) PRIVATE(idir,j,jb,icol,ipw,re1,im1)
1034 [ + + ]: 65152 : do idir=1,3
1035 [ + + ]: 456064 : do j=1,nbj
1036 : 390912 : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1037 : 390912 : re1=zero ; im1=zero
1038 [ + - ]: 390912 : if (istwf_k>1) then
1039 [ # # ]: 0 : do ipw=1,npwsp
1040 : 0 : im1=im1+cg_k(1,ipw,jb)*wpw(1,ipw,icol)+cg_k(2,ipw,jb)*wpw(2,ipw,icol)
1041 : : end do
1042 : 0 : im1=two*im1
1043 : : else
1044 [ + + ]: 55490880 : do ipw=1,npwsp
1045 : 55099968 : re1=re1+cg_k(1,ipw,jb)*wpw(1,ipw,icol)+cg_k(2,ipw,jb)*wpw(2,ipw,icol)
1046 : 55490880 : im1=im1+cg_k(1,ipw,jb)*wpw(2,ipw,icol)-cg_k(2,ipw,jb)*wpw(1,ipw,icol)
1047 : : end do
1048 : : end if
1049 : 439776 : zblk(1,jb,icol)=re1 ; zblk(2,jb,icol)=im1
1050 : : end do
1051 : : end do
1052 : : !$OMP END PARALLEL DO
1053 : : else
1054 : : #ifdef HAVE_OPENMP_OFFLOAD
1055 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:cg_k,wpw,zblk) PRIVATE(idir,j,jb,icol,re1,im1)
1056 : : do idir=1,3
1057 : : do j=1,nbj
1058 : : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1059 : : re1=zero ; im1=zero
1060 : : if (istwf_k>1) then
1061 : : !$OMP PARALLEL DO REDUCTION(+:im1) PRIVATE(ipw)
1062 : : do ipw=1,npwsp
1063 : : im1=im1+cg_k(1,ipw,jb)*wpw(1,ipw,icol)+cg_k(2,ipw,jb)*wpw(2,ipw,icol)
1064 : : end do
1065 : : im1=two*im1
1066 : : else
1067 : : !$OMP PARALLEL DO REDUCTION(+:re1,im1) PRIVATE(ipw)
1068 : : do ipw=1,npwsp
1069 : : re1=re1+cg_k(1,ipw,jb)*wpw(1,ipw,icol)+cg_k(2,ipw,jb)*wpw(2,ipw,icol)
1070 : : im1=im1+cg_k(1,ipw,jb)*wpw(2,ipw,icol)-cg_k(2,ipw,jb)*wpw(1,ipw,icol)
1071 : : end do
1072 : : end if
1073 : : zblk(1,jb,icol)=re1 ; zblk(2,jb,icol)=im1
1074 : : end do
1075 : : end do
1076 : : #endif
1077 : : end if
1078 [ + + ]: 205 : else if (istwf_k>1) then
1079 : : call abi_xgemm('T','N',nrows,3*nbj,2*npwsp,cmplx(two,zero,kind=dp),cg_k,2*npwsp,&
1080 : 3 : & wpw,2*npwsp,czero,zr,nband_k,x_cplx=1,gpu_option=gpu_option)
1081 : 3 : call optics_paw_zr2zblk(zr,zblk,nrows,3*nbj,gpu)
1082 : : else
1083 : : call abi_xgemm('C','N',nrows,3*nbj,npwsp,cone,cg_k,npwsp,&
1084 : 202 : & wpw,npwsp,czero,zblk,nband_k,x_cplx=2,gpu_option=gpu_option)
1085 : : end if
1086 : :
1087 : : call optics_paw_scatter(zblk,j0,nbj,nrows,ishift,mband,all_rows,diag_only,store_half,&
1088 : 32986 : & gpu,psinablapsi)
1089 : :
1090 : : end do ! j0
1091 : :
1092 : : #ifdef HAVE_OPENMP_OFFLOAD
1093 : : !$OMP TARGET EXIT DATA MAP(delete:wpw,zblk,zr) IF(gpu)
1094 : : #endif
1095 [ - + ]: 16493 : ABI_FREE(wpw)
1096 [ - + ]: 16493 : ABI_FREE(zblk)
1097 [ - + ]: 16493 : ABI_FREE(zr)
1098 : : end if
1099 : :
1100 : : !2-B Computation of <psi_n|p_i><p_j|psi_m>(<phi_i|-i.nabla|phi_j>-<tphi_i|-i.nabla|tphi_j>)
1101 : : !2-C + Spin-orbit coupling contribution:
1102 : : ! Sum_ij,ss'[<psi_n,s|p_i><p_j|psi_m,s'>(<phi_i|1/4 Alpha^2 (Sigma^ss' X dV/dr)|phi_j>]
1103 : : !----------------------------------------------------------------------------------
1104 : : !Note: computed as P^H.T with T=K.P, K block-diagonal over atoms:
1105 : : ! Non relativistic part: K^ss'_ij = -i.(<Phi_i|Nabla|Phi_j>-<tPhi_i|Nabla|tPhi_j>) delta_ss'
1106 : : ! Spin-orbit part (not compatible with parallelization over spinors):
1107 : : ! K^ss'_ij = (Sigma^ss' X g_ij), where g_ij = <Phi_i| 1/4 Alpha^2 dV/dr vec(r)/r) |Phi_j>
1108 : : ! With A=Re(Sigma^up-up X g_ij), B=Re(Sigma^dn-up X g_ij), C=Im(Sigma^dn-up X g_ij):
1109 : : ! K^up-up=A, K^dn-dn=-A, K^dn-up=B+iC, K^up-dn=B-iC
1110 : : ! soc_at(re:im,1,idir,ilmn,jlmn,iatom) is (Sigma^up-up X g_ij)
1111 : : ! soc_at(re:im,2,idir,ilmn,jlmn,iatom) is (Sigma^dn-up X g_ij)
1112 : : !Note: nabla_ij is antisymmetric (enforced in pawnabla_init), so ib>jb needs no special care
1113 : : !Note: if cplex=1 (istwfk>1, hence nspinor=1), only the real part of cprj is used
1114 : :
1115 : 16493 : jblo=max(jb1,jpaw1) ; jbhi=min(jb2,jpaw2)
1116 [ + + ]: 16493 : if (jbhi>=jblo) then
1117 : 16399 : nbjmax=min(nbj_blk,jbhi-jblo+1)
1118 [ - + + - : 65596 : ABI_MALLOC(zblk,(2,nband_k,3*nbjmax))
- + + - +
- + - - +
- + - + ]
1119 [ - + + - : 81995 : ABI_MALLOC(tpaw,(cplex,nprojsp,3*nbjmax))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1120 [ - + + - : 65596 : ABI_MALLOC(zr,(1,nband_k,3*nbjmax))
- + + - +
- + - - +
- + - + ]
1121 [ - + + - : 65599 : ABI_MALLOC(cprj_r,(nprojsp,merge(nband_k,1,cplex==1))) ! Only used if cplex=1
+ + - + +
- + - + -
- + - + -
+ ]
1122 : : #ifdef HAVE_OPENMP_OFFLOAD
1123 : : !$OMP TARGET ENTER DATA MAP(alloc:zblk,tpaw,zr,cprj_r) IF(gpu)
1124 : : #endif
1125 [ + + ]: 16399 : if (cplex==1) then
1126 [ + - ]: 3 : if (.not.gpu) then
1127 [ + + + + ]: 1005 : cprj_r(:,:)=cprj_pk(1,:,:)
1128 : : else
1129 : : #ifdef HAVE_OPENMP_OFFLOAD
1130 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) MAP(to:cprj_r,cprj_pk) PRIVATE(ib,irow)
1131 : : do ib=1,nband_k
1132 : : do irow=1,nprojsp
1133 : : cprj_r(irow,ib)=cprj_pk(1,irow,ib)
1134 : : end do
1135 : : end do
1136 : : #endif
1137 : : end if
1138 : : end if
1139 : :
1140 [ + - + - ]: 32798 : do j0=jblo,jbhi,nbjmax
1141 : 16399 : nbj=min(nbjmax,jbhi-j0+1)
1142 : : !If all_rows, compute all ib ; if not, compute only ib<=jb (trapezoidal block)
1143 [ + - ]: 16399 : nrows=merge(nband_k,j0+nbj-1,all_rows)
1144 : :
1145 : : ! Build T=K.P_jb (one element of T per projector row, spinor, column and direction)
1146 [ + - ]: 16399 : if (.not.gpu) then
1147 : : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(SHARED) &
1148 : : !$OMP PRIVATE(idir,j,jb,icol,ispinor,ioff,irow,iatom,iproj0,ilmn,lmn_size,jlmn, &
1149 : : !$OMP nab,aa,bb,cc,re1,im1,x1,y1,x2,y2)
1150 [ + + ]: 65596 : do idir=1,3
1151 [ + + ]: 380734 : do j=1,nbj
1152 [ + + ]: 679593 : do ispinor=1,my_nspinor
1153 [ + + ]: 5680572 : do irow=1,nproj
1154 : 5050176 : jb=j0+j-1 ; icol=(idir-1)*nbj+j ; ioff=(ispinor-1)*nproj
1155 : 5050176 : iatom=row_atom(irow) ; iproj0=atom_off(iatom) ; ilmn=irow-iproj0
1156 [ - + ]: 5050176 : lmn_size=merge(lmn_at(iatom),0,ilmn<=lmn_at(iatom))
1157 : 5050176 : re1=zero ; im1=zero
1158 [ + + ]: 5365434 : if (cplex==1) then
1159 [ + + ]: 30528 : do jlmn=1,lmn_size
1160 : 30528 : re1=re1+nabla_at(idir,ilmn,jlmn,iatom)*cprj_r(ioff+iproj0+jlmn,jb)
1161 : : end do
1162 : 2880 : tpaw(1,ioff+irow,icol)=re1
1163 : : else
1164 [ + + ]: 45430272 : do jlmn=1,lmn_size
1165 : 40382976 : nab=nabla_at(idir,ilmn,jlmn,iatom)
1166 : 40382976 : re1=re1+nab*cprj_pk(2,ioff+iproj0+jlmn,jb)
1167 : 45430272 : im1=im1-nab*cprj_pk(1,ioff+iproj0+jlmn,jb)
1168 : : end do
1169 [ + + ]: 5047296 : if (pawspnorb==1) then
1170 [ + + ]: 69120 : do jlmn=1,lmn_size
1171 : 61440 : aa=soc_at(1,1,idir,ilmn,jlmn,iatom)
1172 : 61440 : bb=soc_at(1,2,idir,ilmn,jlmn,iatom)
1173 : 61440 : cc=soc_at(2,2,idir,ilmn,jlmn,iatom)
1174 : 61440 : x1=cprj_pk(1,iproj0+jlmn,jb) ; y1=cprj_pk(2,iproj0+jlmn,jb)
1175 : 61440 : x2=cprj_pk(1,nproj+iproj0+jlmn,jb) ; y2=cprj_pk(2,nproj+iproj0+jlmn,jb)
1176 [ + + ]: 69120 : if (ispinor==1) then
1177 : : !Spin up: A.P^up + (B-iC).P^dn
1178 : 30720 : re1=re1+aa*x1+bb*x2+cc*y2
1179 : 30720 : im1=im1+aa*y1+bb*y2-cc*x2
1180 : : else
1181 : : !Spin down: (B+iC).P^up - A.P^dn
1182 : 30720 : re1=re1+bb*x1-cc*y1-aa*x2
1183 : 30720 : im1=im1+bb*y1+cc*x1-aa*y2
1184 : : end if
1185 : : end do
1186 : : end if
1187 : 5047296 : tpaw(1,ioff+irow,icol)=re1 ; tpaw(2,ioff+irow,icol)=im1
1188 : : end if
1189 : : end do ! irow
1190 : : end do ! ispinor
1191 : : end do ! j
1192 : : end do ! idir
1193 : : !$OMP END PARALLEL DO
1194 : : else
1195 : : #ifdef HAVE_OPENMP_OFFLOAD
1196 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(4) &
1197 : : !$OMP& MAP(to:tpaw,cprj_pk,cprj_r,nabla_at,soc_at,lmn_at,atom_off,row_atom) &
1198 : : !$OMP& PRIVATE(idir,j,jb,icol,ispinor,ioff,irow,iatom,iproj0,ilmn,lmn_size,jlmn, &
1199 : : !$OMP& nab,aa,bb,cc,re1,im1,x1,y1,x2,y2)
1200 : : do idir=1,3
1201 : : do j=1,nbj
1202 : : do ispinor=1,my_nspinor
1203 : : do irow=1,nproj
1204 : : jb=j0+j-1 ; icol=(idir-1)*nbj+j ; ioff=(ispinor-1)*nproj
1205 : : iatom=row_atom(irow) ; iproj0=atom_off(iatom) ; ilmn=irow-iproj0
1206 : : lmn_size=merge(lmn_at(iatom),0,ilmn<=lmn_at(iatom))
1207 : : re1=zero ; im1=zero
1208 : : if (cplex==1) then
1209 : : do jlmn=1,lmn_size
1210 : : re1=re1+nabla_at(idir,ilmn,jlmn,iatom)*cprj_r(ioff+iproj0+jlmn,jb)
1211 : : end do
1212 : : tpaw(1,ioff+irow,icol)=re1
1213 : : else
1214 : : do jlmn=1,lmn_size
1215 : : nab=nabla_at(idir,ilmn,jlmn,iatom)
1216 : : re1=re1+nab*cprj_pk(2,ioff+iproj0+jlmn,jb)
1217 : : im1=im1-nab*cprj_pk(1,ioff+iproj0+jlmn,jb)
1218 : : end do
1219 : : if (pawspnorb==1) then
1220 : : do jlmn=1,lmn_size
1221 : : aa=soc_at(1,1,idir,ilmn,jlmn,iatom)
1222 : : bb=soc_at(1,2,idir,ilmn,jlmn,iatom)
1223 : : cc=soc_at(2,2,idir,ilmn,jlmn,iatom)
1224 : : x1=cprj_pk(1,iproj0+jlmn,jb) ; y1=cprj_pk(2,iproj0+jlmn,jb)
1225 : : x2=cprj_pk(1,nproj+iproj0+jlmn,jb) ; y2=cprj_pk(2,nproj+iproj0+jlmn,jb)
1226 : : if (ispinor==1) then
1227 : : !Spin up: A.P^up + (B-iC).P^dn
1228 : : re1=re1+aa*x1+bb*x2+cc*y2
1229 : : im1=im1+aa*y1+bb*y2-cc*x2
1230 : : else
1231 : : !Spin down: (B+iC).P^up - A.P^dn
1232 : : re1=re1+bb*x1-cc*y1-aa*x2
1233 : : im1=im1+bb*y1+cc*x1-aa*y2
1234 : : end if
1235 : : end do
1236 : : end if
1237 : : tpaw(1,ioff+irow,icol)=re1 ; tpaw(2,ioff+irow,icol)=im1
1238 : : end if
1239 : : end do ! irow
1240 : : end do ! ispinor
1241 : : end do ! j
1242 : : end do ! idir
1243 : : #endif
1244 : : end if
1245 : :
1246 : : ! Z=P^H.T (only diagonal elements if diag_only)
1247 [ + + ]: 16399 : if (diag_only) then
1248 [ + - ]: 16288 : if (.not.gpu) then
1249 : : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(SHARED) PRIVATE(idir,j,jb,icol,irow,re1,im1)
1250 [ + + ]: 65152 : do idir=1,3
1251 [ + + ]: 377728 : do j=1,nbj
1252 : 312576 : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1253 : 312576 : re1=zero ; im1=zero
1254 [ + - ]: 312576 : if (cplex==1) then
1255 [ # # ]: 0 : do irow=1,nprojsp
1256 : 0 : im1=im1-cprj_r(irow,jb)*tpaw(1,irow,icol)
1257 : : end do
1258 : : else
1259 [ + + ]: 5313792 : do irow=1,nprojsp
1260 : 5001216 : re1=re1+cprj_pk(1,irow,jb)*tpaw(1,irow,icol)+cprj_pk(2,irow,jb)*tpaw(2,irow,icol)
1261 : 5313792 : im1=im1+cprj_pk(1,irow,jb)*tpaw(2,irow,icol)-cprj_pk(2,irow,jb)*tpaw(1,irow,icol)
1262 : : end do
1263 : : end if
1264 : 361440 : zblk(1,jb,icol)=re1 ; zblk(2,jb,icol)=im1
1265 : : end do
1266 : : end do
1267 : : !$OMP END PARALLEL DO
1268 : : else
1269 : : #ifdef HAVE_OPENMP_OFFLOAD
1270 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:cprj_pk,cprj_r,tpaw,zblk) PRIVATE(idir,j,jb,icol,re1,im1)
1271 : : do idir=1,3
1272 : : do j=1,nbj
1273 : : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1274 : : re1=zero ; im1=zero
1275 : : if (cplex==1) then
1276 : : !$OMP PARALLEL DO REDUCTION(+:im1) PRIVATE(irow)
1277 : : do irow=1,nprojsp
1278 : : im1=im1-cprj_r(irow,jb)*tpaw(1,irow,icol)
1279 : : end do
1280 : : else
1281 : : !$OMP PARALLEL DO REDUCTION(+:re1,im1) PRIVATE(irow)
1282 : : do irow=1,nprojsp
1283 : : re1=re1+cprj_pk(1,irow,jb)*tpaw(1,irow,icol)+cprj_pk(2,irow,jb)*tpaw(2,irow,icol)
1284 : : im1=im1+cprj_pk(1,irow,jb)*tpaw(2,irow,icol)-cprj_pk(2,irow,jb)*tpaw(1,irow,icol)
1285 : : end do
1286 : : end if
1287 : : zblk(1,jb,icol)=re1 ; zblk(2,jb,icol)=im1
1288 : : end do
1289 : : end do
1290 : : #endif
1291 : : end if
1292 [ + + ]: 111 : else if (cplex==1) then
1293 : : call abi_xgemm('T','N',nrows,3*nbj,nprojsp,-cone,cprj_r,nprojsp,&
1294 : 3 : & tpaw,nprojsp,czero,zr,nband_k,x_cplx=1,gpu_option=gpu_option)
1295 : 3 : call optics_paw_zr2zblk(zr,zblk,nrows,3*nbj,gpu)
1296 : : else
1297 : : call abi_xgemm('C','N',nrows,3*nbj,nprojsp,cone,cprj_pk,nprojsp,&
1298 : 108 : & tpaw,nprojsp,czero,zblk,nband_k,x_cplx=2,gpu_option=gpu_option)
1299 : : end if
1300 : :
1301 : : call optics_paw_scatter(zblk,j0,nbj,nrows,ishift,mband,all_rows,diag_only,store_half,&
1302 : 32798 : & gpu,psinablapsi)
1303 : :
1304 : : end do ! j0
1305 : :
1306 : : #ifdef HAVE_OPENMP_OFFLOAD
1307 : : !$OMP TARGET EXIT DATA MAP(delete:zblk,tpaw,zr,cprj_r) IF(gpu)
1308 : : #endif
1309 [ - + ]: 16399 : ABI_FREE(zblk)
1310 [ - + ]: 16399 : ABI_FREE(tpaw)
1311 [ - + ]: 16399 : ABI_FREE(zr)
1312 [ - + ]: 16399 : ABI_FREE(cprj_r)
1313 : : end if
1314 : :
1315 : : !Get result back from device
1316 : : if (gpu) then
1317 : : #ifdef HAVE_OPENMP_OFFLOAD
1318 : : !$OMP TARGET UPDATE FROM(psinablapsi(:,:,1:nzero))
1319 : : !$OMP TARGET EXIT DATA MAP(delete:psinablapsi)
1320 : : #endif
1321 : : end if
1322 : :
1323 [ - + - + : 16493 : end subroutine optics_paw_dipoles_k
- + - + -
+ ]
1324 : : !!***
1325 : :
1326 : : !----------------------------------------------------------------------
1327 : :
1328 : : !!****f* m_paw_optics/optics_paw_zr2zblk
1329 : : !! NAME
1330 : : !! optics_paw_zr2zblk
1331 : : !!
1332 : : !! FUNCTION
1333 : : !! Store a real block (imaginary part of the matrix elements, computed with real GEMMs
1334 : : !! when istwfk>1) into the complex work array used in optics_paw_dipoles_k.
1335 : : !!
1336 : : !! INPUTS
1337 : : !! zr(1,:,ncol)=real block (imaginary part)
1338 : : !! nrows,ncol=size of the block
1339 : : !! gpu=true if arrays are on GPU
1340 : : !!
1341 : : !! OUTPUT
1342 : : !! zblk(2,:,ncol)=complex block (real part set to zero)
1343 : : !!
1344 : : !! SOURCE
1345 : :
1346 [ - + - + ]: 6 : subroutine optics_paw_zr2zblk(zr,zblk,nrows,ncol,gpu)
1347 : :
1348 : : !Arguments ------------------------------------
1349 : : !scalars
1350 : : integer,intent(in) :: nrows,ncol
1351 : : logical,intent(in) :: gpu
1352 : : !arrays
1353 : : real(dp),intent(in) :: zr(:,:,:)
1354 : : real(dp),intent(inout) :: zblk(:,:,:)
1355 : :
1356 : : !Local variables-------------------------------
1357 : : !scalars
1358 : : integer :: ib,icol
1359 : :
1360 : : ! ************************************************************************
1361 : :
1362 [ + - ]: 6 : if (.not.gpu) then
1363 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(icol,ib)
1364 [ + + ]: 258 : do icol=1,ncol
1365 [ + + ]: 3930 : do ib=1,nrows
1366 : 3924 : zblk(1,ib,icol)=zero ; zblk(2,ib,icol)=zr(1,ib,icol)
1367 : : end do
1368 : : end do
1369 : : !$OMP END PARALLEL DO
1370 : : else
1371 : : #ifdef HAVE_OPENMP_OFFLOAD
1372 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) MAP(to:zr,zblk) PRIVATE(icol,ib)
1373 : : do icol=1,ncol
1374 : : do ib=1,nrows
1375 : : zblk(1,ib,icol)=zero ; zblk(2,ib,icol)=zr(1,ib,icol)
1376 : : end do
1377 : : end do
1378 : : #endif
1379 : : end if
1380 : :
1381 : 6 : end subroutine optics_paw_zr2zblk
1382 : : !!***
1383 : :
1384 : : !----------------------------------------------------------------------
1385 : :
1386 : : !!****f* m_paw_optics/optics_paw_scatter
1387 : : !! NAME
1388 : : !! optics_paw_scatter
1389 : : !!
1390 : : !! FUNCTION
1391 : : !! Add a block of matrix elements Z(ib,jb,idir), for jb=j0:j0+nbj-1, from the dense
1392 : : !! work array used in optics_paw_dipoles_k to the (2,3,:) storage of the output array.
1393 : : !!
1394 : : !! INPUTS
1395 : : !! zblk(2,:,3*nbj)=dense block; column (idir-1)*nbj+j holds jb=j0+j-1
1396 : : !! j0=first band jb of the block
1397 : : !! nbj=number of bands jb in the block
1398 : : !! nrows=number of rows ib computed in the block
1399 : : !! ishift=index of the last element before the first column stored in psinablapsi
1400 : : !! mband=maximum number of bands
1401 : : !! all_rows=if true, store all ib computed for each jb (full storage, no second half)
1402 : : !! diag_only=if true, store only ib=jb elements
1403 : : !! store_half=if true, use packed storage for the (ib<=jb) half matrix;
1404 : : !! if false, also fill the second half of the matrix
1405 : : !! gpu=true if arrays are on GPU
1406 : : !!
1407 : : !! SIDE EFFECTS
1408 : : !! psinablapsi(2,3,:)=output array, Z block is added to it
1409 : : !!
1410 : : !! SOURCE
1411 : :
1412 : 98676 : subroutine optics_paw_scatter(zblk,j0,nbj,nrows,ishift,mband,all_rows,diag_only,store_half,&
1413 [ - + - + ]: 32892 : & gpu,psinablapsi)
1414 : :
1415 : : !Arguments ------------------------------------
1416 : : !scalars
1417 : : integer,intent(in) :: j0,nbj,nrows,ishift,mband
1418 : : logical,intent(in) :: all_rows,diag_only,store_half,gpu
1419 : : !arrays
1420 : : real(dp),intent(in) :: zblk(:,:,:)
1421 : : real(dp),intent(inout) :: psinablapsi(:,:,:)
1422 : :
1423 : : !Local variables-------------------------------
1424 : : !scalars
1425 : : integer :: ib,icol,idir,ii,j,jb
1426 : :
1427 : : ! ************************************************************************
1428 : :
1429 : : !Note: in the full matrix case (if not all_rows), element (ib,jb) with ib<jb is also stored
1430 : : ! at (jb,ib); these two locations never overlap between different (ib,jb) pairs.
1431 : : ! This is only done if ishift=0 (all columns are stored).
1432 : :
1433 [ + + ]: 32892 : if (diag_only) then
1434 [ + - ]: 32576 : if (.not.gpu) then
1435 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(idir,j,jb,icol)
1436 [ + + ]: 130304 : do idir=1,3
1437 [ + + ]: 833792 : do j=1,nbj
1438 : 703488 : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1439 [ + + ]: 2208192 : psinablapsi(1:2,idir,jb)=psinablapsi(1:2,idir,jb)+zblk(1:2,jb,icol)
1440 : : end do
1441 : : end do
1442 : : !$OMP END PARALLEL DO
1443 : : else
1444 : : #ifdef HAVE_OPENMP_OFFLOAD
1445 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) MAP(to:zblk,psinablapsi) PRIVATE(idir,j,jb,icol)
1446 : : do idir=1,3
1447 : : do j=1,nbj
1448 : : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1449 : : psinablapsi(1,idir,jb)=psinablapsi(1,idir,jb)+zblk(1,jb,icol)
1450 : : psinablapsi(2,idir,jb)=psinablapsi(2,idir,jb)+zblk(2,jb,icol)
1451 : : end do
1452 : : end do
1453 : : #endif
1454 : : end if
1455 : 32576 : return
1456 : : end if
1457 : :
1458 [ + - ]: 316 : if (.not.gpu) then
1459 : : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(idir,j,ib,jb,icol,ii)
1460 [ + + ]: 1264 : do idir=1,3
1461 [ + + ]: 11338 : do j=1,nbj
1462 [ + + ]: 133191 : do ib=1,nrows
1463 : 122169 : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1464 [ - + ]: 132243 : if (all_rows) then
1465 : 0 : ii=(jb-1)*mband+ib-ishift
1466 [ # # ]: 0 : psinablapsi(1:2,idir,ii)=psinablapsi(1:2,idir,ii)+zblk(1:2,ib,icol)
1467 [ + + ]: 122169 : else if (ib<=jb) then
1468 [ - + ]: 67659 : ii=merge((jb*(jb-1))/2,(jb-1)*mband,store_half)+ib-ishift
1469 [ + + ]: 202977 : psinablapsi(1:2,idir,ii)=psinablapsi(1:2,idir,ii)+zblk(1:2,ib,icol)
1470 : : ! Second half of the (n,m) matrix
1471 [ - + - - ]: 67659 : if ((.not.store_half).and.ib<jb) then
1472 : 0 : ii=(ib-1)*mband+jb
1473 : 0 : psinablapsi(1,idir,ii)=psinablapsi(1,idir,ii)+zblk(1,ib,icol)
1474 : 0 : psinablapsi(2,idir,ii)=psinablapsi(2,idir,ii)-zblk(2,ib,icol)
1475 : : end if
1476 : : end if
1477 : : end do
1478 : : end do
1479 : : end do
1480 : : !$OMP END PARALLEL DO
1481 : : else
1482 : : #ifdef HAVE_OPENMP_OFFLOAD
1483 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) MAP(to:zblk,psinablapsi) PRIVATE(idir,j,ib,jb,icol,ii)
1484 : : do idir=1,3
1485 : : do j=1,nbj
1486 : : do ib=1,nrows
1487 : : jb=j0+j-1 ; icol=(idir-1)*nbj+j
1488 : : if (all_rows) then
1489 : : ii=(jb-1)*mband+ib-ishift
1490 : : psinablapsi(1,idir,ii)=psinablapsi(1,idir,ii)+zblk(1,ib,icol)
1491 : : psinablapsi(2,idir,ii)=psinablapsi(2,idir,ii)+zblk(2,ib,icol)
1492 : : else if (ib<=jb) then
1493 : : ii=merge((jb*(jb-1))/2,(jb-1)*mband,store_half)+ib-ishift
1494 : : psinablapsi(1,idir,ii)=psinablapsi(1,idir,ii)+zblk(1,ib,icol)
1495 : : psinablapsi(2,idir,ii)=psinablapsi(2,idir,ii)+zblk(2,ib,icol)
1496 : : ! Second half of the (n,m) matrix
1497 : : if ((.not.store_half).and.ib<jb) then
1498 : : ii=(ib-1)*mband+jb
1499 : : psinablapsi(1,idir,ii)=psinablapsi(1,idir,ii)+zblk(1,ib,icol)
1500 : : psinablapsi(2,idir,ii)=psinablapsi(2,idir,ii)-zblk(2,ib,icol)
1501 : : end if
1502 : : end if
1503 : : end do
1504 : : end do
1505 : : end do
1506 : : #endif
1507 : : end if
1508 : :
1509 : : end subroutine optics_paw_scatter
1510 : : !!***
1511 : :
1512 : : !----------------------------------------------------------------------
1513 : :
1514 : : !!****f* m_paw_optics/optics_paw_core
1515 : : !! NAME
1516 : : !! optics_paw_core
1517 : : !!
1518 : : !! FUNCTION
1519 : : !! Compute matrix elements need for X spectr. (in the PAW context) and store them in a file
1520 : : !! Matrix elements = <Phi_core|Nabla|Phi_j>
1521 : : !!
1522 : : !! COPYRIGHT
1523 : : !! Copyright (C) 2005-2026 ABINIT group (SM,MT,NB)
1524 : : !! This file is distributed under the terms of the
1525 : : !! GNU General Public License, see ~ABINIT/COPYING
1526 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1527 : : !!
1528 : : !! INPUTS
1529 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
1530 : : !! cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
1531 : : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
1532 : : !! dtfil <type(datafiles_type)>=variables related to files
1533 : : !! dtset <type(dataset_type)>=all input variables for this dataset
1534 : : !! filpsp(ntypat)=name(s) of the pseudopotential file(s)
1535 : : !! mband=maximum number of bands
1536 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
1537 : : !! mkmem =number of k points treated by this node.
1538 : : !! mpi_enreg=information about MPI parallelization
1539 : : !! mpsang =1+maximum angular momentum for nonlocal pseudopotentials
1540 : : !! natom=number of atoms in cell.
1541 : : !! nkpt=number of k points.
1542 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
1543 : : !! pawang <type(pawang_type)>= PAW ANGular mesh discretization and related data
1544 : : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
1545 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
1546 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
1547 : : !! znucl(ntypat)=atomic number of atom type
1548 : : !!
1549 : : !! OUTPUT
1550 : : !! (only writing in a file)
1551 : : !!
1552 : : !! SOURCE
1553 : :
1554 : 14 : subroutine optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen0,filpsp,hdr,&
1555 : : & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,&
1556 : 14 : & pawang,pawrad,pawrhoij,pawtab,znucl,rcpaw)
1557 : :
1558 : : !Arguments ------------------------------------
1559 : : !scalars
1560 : : integer,intent(in) :: mband,mcprj,mkmem,mpsang,natom,nkpt,nsppol
1561 : : type(MPI_type),intent(in) :: mpi_enreg
1562 : : type(datafiles_type),intent(in) :: dtfil
1563 : : type(dataset_type),intent(in) :: dtset
1564 : : type(hdr_type),intent(inout) :: hdr
1565 : : type(pawang_type),intent(in) :: pawang
1566 : : type(rcpaw_type), intent(in), pointer,optional :: rcpaw
1567 : : !arrays
1568 : : integer,intent(in) :: atindx1(natom),dimcprj(natom)
1569 : : character(len=fnlen),intent(in) :: filpsp(dtset%ntypat)
1570 : : real(dp),intent(in) :: eigen0(mband*nkpt*nsppol),znucl(dtset%ntypat)
1571 : : type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
1572 : : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
1573 : : type(pawrhoij_type),intent(inout) :: pawrhoij(mpi_enreg%my_natom)
1574 : : type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
1575 : :
1576 : : !Local variables-------------------------------
1577 : : !scalars
1578 : : integer,parameter :: master=0
1579 : : integer :: bdtot_index,cplex,etiq,iatom,ic,ibg,idir,ilmn2
1580 : : integer :: ierr,ikpt,ilmn,iln,ount,is,my_jb,nphicor_lm
1581 : : integer :: iorder_cprj,ispinor,isppol,istwf_k,itypat,itypat2
1582 : : integer :: jb,jbsp,jlmn,lmn_size,lmncmax,mband_cprj,ncid,varid
1583 : : integer :: me,my_nspinor,nband_cprj_k,option_core,pnp_size
1584 : : integer :: nband_k,nphicor,ncorespinor,sender,iomode,fformopt,master_spfftband
1585 : : integer :: spaceComm_band,spaceComm_bandspinorfft,spaceComm_fft,spaceComm_kpt
1586 : : integer :: spaceComm_spinor,spaceComm_bandspinor,spaceComm_spinorfft,spaceComm_w
1587 : : logical :: already_has_nabla,cprj_paral_band,mykpt,myband,use_rcpaw_data
1588 : : logical :: iomode_etsf_mpiio,use_spinorbit
1589 : : logical :: i_am_master,i_am_master_band,i_am_master_spfft
1590 : : real(dp) :: cpnm1,cpnm2,el_temp
1591 : : character(len=500) :: msg
1592 : : !arrays
1593 : : integer :: nc_count(6),nc_start(6),nc_stride(6),tmp_shape(3)
1594 : 14 : integer,allocatable :: lcor(:,:),ncor(:,:),kappacor(:,:),nphicor_arr(:),lmn2ln(:,:)
1595 : : real(dp) :: tsec(2)
1596 : 14 : real(dp),allocatable :: energy_cor(:,:),occ_cor(:,:),maxocc_cor(:,:)
1597 : 14 : real(dp),allocatable :: psinablapsi(:,:,:,:,:),psinablapsi_soc(:,:,:,:,:),psinablapsi2(:,:,:,:)
1598 : 14 : real(dp),pointer :: soc_ij(:,:,:)
1599 : 14 : type(coeff5_type),allocatable,target :: phisocphj(:)
1600 : 14 : type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
1601 : : type(nctkdim_t) :: ncdims(3)
1602 : : type(nctkarr_t) :: nctk_arrays(9)
1603 : 14 : type(atomorb_type),allocatable :: atm(:)
1604 : :
1605 : : ! ************************************************************************
1606 : :
1607 : : DBG_ENTER("COLL")
1608 : 14 : call timab(1185,1,tsec)
1609 : : ABI_NVTX_START_RANGE(NVTX_OPTICS_PAW_CORE)
1610 : :
1611 : : !Compatibility tests
1612 : 14 : msg="mkmem==0 not supported anymore!"
1613 [ - + ]: 14 : ABI_CHECK(mkmem/=0,msg)
1614 : : !Probably should check for spinor parallelism, because thats likely to not work correctly
1615 : 14 : msg="Spinor parallelism not implemented for optics_paw_core!"
1616 [ - + ]: 14 : ABI_CHECK(dtset%npspinor==1,msg)
1617 : : !Is mpi_enreg initialized?
1618 : : if (xmpi_paral==1) then
1619 [ + + ]: 56 : tmp_shape = shape(mpi_enreg%proc_distrb)
1620 [ - + ]: 14 : if (nkpt > tmp_shape(1)) then
1621 : 0 : ABI_BUG('problem with proc_distrb!')
1622 : : end if
1623 : : end if
1624 : :
1625 : : !Init parallelism
1626 : 14 : spaceComm_w=mpi_enreg%comm_cell
1627 [ + + ]: 14 : if (mpi_enreg%paral_kgb==1) then
1628 : 8 : spaceComm_kpt=mpi_enreg%comm_kpt
1629 : 8 : spaceComm_fft=mpi_enreg%comm_fft
1630 : 8 : spaceComm_band=mpi_enreg%comm_band
1631 : 8 : spaceComm_spinor=mpi_enreg%comm_spinor
1632 : 8 : spaceComm_bandspinor=mpi_enreg%comm_bandspinor
1633 : 8 : spaceComm_spinorfft=mpi_enreg%comm_spinorfft
1634 : 8 : spaceComm_bandspinorfft=mpi_enreg%comm_bandspinorfft
1635 : : else
1636 : 6 : spaceComm_kpt=mpi_enreg%comm_kpt
1637 : 6 : spaceComm_fft=xmpi_comm_self
1638 : 6 : spaceComm_band=mpi_enreg%comm_band
1639 : 6 : spaceComm_spinor=xmpi_comm_self
1640 : 6 : spaceComm_bandspinor=spaceComm_band
1641 : 6 : spaceComm_spinorfft=xmpi_comm_self
1642 : 6 : spaceComm_bandspinorfft=xmpi_comm_self
1643 : : end if
1644 : 14 : me=xmpi_comm_rank(spaceComm_w)
1645 : 14 : i_am_master=(me==master)
1646 : 14 : i_am_master_band=(xmpi_comm_rank(spaceComm_band)==master)
1647 : 14 : i_am_master_spfft=(xmpi_comm_rank(spaceComm_spinorfft)==master)
1648 : 14 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1649 : :
1650 : : !------------------------------------------------------------------------------------------------
1651 : : !1- Reading of core wavefunctions
1652 : : !------------------------------------------------------------------------------------------------
1653 : : #if defined(HAVE_GPU_MARKERS)
1654 : : call nvtxStartRange("read")
1655 : : #endif
1656 : : !TODO At present, impose 2-spinor simulataneously for core and valence WF
1657 [ - + + - : 58 : ABI_MALLOC(atm,(dtset%ntypat))
+ - + - -
+ - + +
+ ]
1658 [ + + ]: 30 : do itypat=1,dtset%ntypat
1659 : 16 : use_rcpaw_data=.false.
1660 [ + + ]: 16 : if(present(rcpaw)) then
1661 [ + - ]: 1 : if(associated(rcpaw)) then
1662 [ + - ]: 1 : if(allocated(rcpaw%atm)) then
1663 [ + - ]: 1 : if(dtset%ntypat==size(rcpaw%atm)) then
1664 : 1 : call copy_atomorb(rcpaw%atm(itypat),atm(itypat))
1665 : : use_rcpaw_data=.true.
1666 : : endif
1667 : : endif
1668 : : endif
1669 : : endif
1670 : 14 : if(.not.use_rcpaw_data) then
1671 : 15 : call pawpsp_init_core(atm(itypat),filpsp(itypat),radmesh=pawrad(itypat))
1672 : : endif
1673 : : enddo
1674 : :
1675 [ - + ]: 14 : if(dtset%cwfs_wouth==1) then
1676 [ # # ]: 0 : do itypat=1,dtset%ntypat
1677 [ # # ]: 0 : if(atm(itypat)%mult==1) then
1678 [ # # ]: 0 : do itypat2=1,dtset%ntypat
1679 : : if(atm(itypat2)%mult>1.and.atm(itypat2)%znucl==atm(itypat)%znucl.and.&
1680 : : & atm(itypat2)%mesh_size==atm(itypat)%mesh_size.and.&
1681 : : & atm(itypat2)%ln_size==atm(itypat)%ln_size.and.&
1682 [ # # # # : 0 : & atm(itypat2)%nsppol==atm(itypat)%nsppol.and.&
# # # # #
# # # ]
1683 : 0 : & atm(itypat2)%zcore_orig>atm(itypat)%zcore_orig) then
1684 : 0 : write(std_out,*) 'Core wfs of typat ',itypat,' replaced by those of typat ',itypat2
1685 [ # # # # : 0 : atm(itypat)%phi=atm(itypat2)%phi
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
1686 : : exit
1687 : : endif
1688 : : enddo
1689 : : endif
1690 : : enddo
1691 : : endif
1692 : :
1693 : :
1694 : : !nphicor=max. number of core (l,n) states over atom types
1695 : : !nphicor_lm=max. number of core (l,m,n) channels over atom types
1696 : : ! (spinor components are removed below if nspinor=2)
1697 : 14 : nphicor=0
1698 : 14 : ncorespinor=0
1699 : 14 : nphicor_lm=0
1700 [ + + ]: 30 : do itypat=1,dtset%ntypat
1701 : 16 : nphicor=max(nphicor,atm(itypat)%ln_size)
1702 : 16 : nphicor_lm=max(nphicor_lm,atm(itypat)%lmn_size)
1703 [ - + ]: 16 : if(atm(itypat)%nsppol>1) ABI_ERROR("nsppol>1 is work in progress for optics_paw_core")
1704 [ - + ]: 16 : if(atm(itypat)%nspinor/=dtset%nspinor) ABI_ERROR("Core and valence not same number of spinors")
1705 : 30 : ncorespinor=max(ncorespinor,atm(itypat)%nspinor)
1706 : : enddo
1707 [ + + ]: 14 : if(dtset%nspinor==2) nphicor_lm=nphicor_lm/2
1708 : :
1709 [ - + + - : 56 : ABI_MALLOC(lmn2ln,(nphicor_lm,dtset%ntypat))
- + + - +
- + - - +
- + ]
1710 [ + + + + ]: 140 : lmn2ln=0
1711 [ + + ]: 30 : do itypat=1,dtset%ntypat
1712 : 16 : lmncmax=atm(itypat)%lmn_size
1713 [ + + ]: 16 : if(dtset%nspinor==2) lmncmax=lmncmax/2
1714 [ + + ]: 132 : do ilmn=1,lmncmax
1715 [ + + ]: 118 : if(dtset%nspinor==1) then
1716 : 42 : lmn2ln(ilmn,itypat)=atm(itypat)%indlmn(5,ilmn)
1717 : : else
1718 [ + - ]: 594 : do ilmn2=1,atm(itypat)%lmn_size
1719 [ + + ]: 594 : if(atm(itypat)%indlmn(9,ilmn2)==ilmn) then
1720 : 60 : lmn2ln(ilmn,itypat)=atm(itypat)%indlmn(5,ilmn2)
1721 : 60 : exit
1722 : : endif
1723 : : enddo
1724 : : endif
1725 : : enddo
1726 : : enddo
1727 : : #if defined(HAVE_GPU_MARKERS)
1728 : : call nvtxEndRange()
1729 : : #endif
1730 : :
1731 : : !----------------------------------------------------------------------------------
1732 : : !2- Computation of phipphj=<phi_i|nabla|phi_core>
1733 : : !----------------------------------------------------------------------------------
1734 : :
1735 : : #if defined(HAVE_GPU_MARKERS)
1736 : : call nvtxStartRange("pawnabla_core_init")
1737 : : #endif
1738 [ + - - + ]: 14 : already_has_nabla=all(pawtab(:)%has_nabla==3)
1739 [ + + + - : 14 : if (ncorespinor==2) already_has_nabla=all(pawtab(:)%has_nabla==4)
- + ]
1740 : : ! Should check whether this would work with spinor parallelism
1741 : 14 : call pawnabla_core_init(mpsang,dtset%ntypat,pawrad,pawtab,atm)
1742 : :
1743 : : !Get electronic temperature from dtset
1744 [ + - - - : 14 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
- - ]
1745 : : #if defined(HAVE_GPU_MARKERS)
1746 : : call nvtxEndRange()
1747 : : #endif
1748 : :
1749 : : !Compute spin-orbit contributions if necessary
1750 [ + + + - ]: 14 : use_spinorbit=(dtset%pawspnorb==1.and.dtset%userie/=111) ! For testing purpose
1751 : : if (use_spinorbit) then
1752 : : #if defined(HAVE_GPU_MARKERS)
1753 : : call nvtxStartRange("pawnabla_soc_init")
1754 : : #endif
1755 : 6 : option_core=1
1756 : : call pawnabla_soc_init(el_temp,phisocphj,option_core,dtset%ixc,mpi_enreg%my_natom,natom,&
1757 : : & dtset%nspden,dtset%ntypat,pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,&
1758 : : & dtset%spnorbscl,dtset%typat,dtset%xc_denpos,dtset%xc_taupos,znucl,&
1759 : 6 : & atm=atm,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1760 : : #if defined(HAVE_GPU_MARKERS)
1761 : : call nvtxEndRange()
1762 : : #endif
1763 : : end if
1764 : :
1765 : : !----------------------------------------------------------------------------------
1766 : : !3- Opening of OPT2 file and header writing
1767 : : !----------------------------------------------------------------------------------
1768 : :
1769 : : #if defined(HAVE_GPU_MARKERS)
1770 : : call nvtxStartRange("header_write")
1771 : : #endif
1772 : : !I/O mode is netCDF or Fortran
1773 : : iomode=merge(IO_MODE_ETSF,IO_MODE_FORTRAN_MASTER,dtset%iomode==IO_MODE_ETSF)
1774 : 14 : if (use_netcdf_forced) iomode=IO_MODE_ETSF
1775 : :
1776 [ - + + - : 56 : ABI_MALLOC(energy_cor,(nphicor,dtset%ntypat))
- + + - +
- + - - +
- + ]
1777 [ - + + - : 56 : ABI_MALLOC(ncor,(nphicor,dtset%ntypat))
- + + - +
- + - - +
- + ]
1778 [ - + + - : 42 : ABI_MALLOC(lcor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1779 [ - + + - : 42 : ABI_MALLOC(kappacor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1780 [ - + + - : 42 : ABI_MALLOC(occ_cor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1781 [ - + + - : 42 : ABI_MALLOC(maxocc_cor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1782 [ + + + + ]: 84 : energy_cor=zero
1783 [ + + + + ]: 84 : ncor=0
1784 [ + + + + ]: 84 : lcor=0
1785 [ + + + + ]: 84 : kappacor=0
1786 [ + + + + ]: 84 : occ_cor=zero
1787 [ + + + + ]: 84 : maxocc_cor=zero
1788 [ + + ]: 30 : do itypat=1,dtset%ntypat
1789 [ + + ]: 80 : do iln=1,atm(itypat)%ln_size
1790 : 50 : energy_cor(iln,itypat)=atm(itypat)%eig(iln,1)
1791 : 50 : ncor(iln,itypat)=atm(itypat)%indln(2,iln)
1792 : 50 : lcor(iln,itypat)=atm(itypat)%indln(1,iln)
1793 : 50 : occ_cor(iln,itypat)=atm(itypat)%occ(iln,1)
1794 : 50 : maxocc_cor(iln,itypat)=atm(itypat)%max_occ(iln,1)
1795 [ + + ]: 66 : if(atm(itypat)%dirac) then
1796 : 24 : kappacor(iln,itypat)=atm(itypat)%kappa(iln)
1797 : : else
1798 : 26 : kappacor(iln,itypat)=0
1799 : : endif
1800 : : enddo
1801 : : enddo
1802 : :
1803 : : !(master proc only)
1804 [ + + ]: 14 : if (i_am_master) then
1805 : : ! ====> NETCDF format
1806 : : if (iomode==IO_MODE_ETSF) then
1807 : 8 : fformopt=611
1808 [ - + + - : 24 : ABI_MALLOC(nphicor_arr,(dtset%ntypat))
+ - + - -
+ - + ]
1809 [ + + ]: 18 : do itypat=1,dtset%ntypat
1810 : 18 : nphicor_arr(itypat)=atm(itypat)%ln_size
1811 : : enddo
1812 : : ! Open/create nc file
1813 [ - + ]: 8 : NCF_CHECK(nctk_open_create(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),xmpi_comm_self))
1814 : : ! Write header data
1815 [ - + ]: 8 : NCF_CHECK(hdr%ncwrite(ncid,fformopt,nc_define=.true.))
1816 : : ! Define additional dimensions
1817 : 8 : ncdims(1)%name="number_of_atom_types"
1818 : 8 : ncdims(1)%value=dtset%ntypat
1819 : 8 : ncdims(2)%name="max_number_of_core_states"
1820 : 8 : ncdims(2)%value=nphicor
1821 : 8 : ncdims(3)%name="number_of_core_spinor_components"
1822 : 8 : ncdims(3)%value=ncorespinor
1823 [ - + ]: 8 : NCF_CHECK(nctk_def_dims(ncid,ncdims))
1824 : 8 : nctk_arrays(1)%name="number_of_core_states"
1825 : 8 : nctk_arrays(1)%dtype="int"
1826 : 8 : nctk_arrays(1)%shape_str="number_of_atom_types"
1827 : 8 : nctk_arrays(2)%name="eigenvalues_core"
1828 : 8 : nctk_arrays(2)%dtype="dp"
1829 : 8 : nctk_arrays(2)%shape_str="max_number_of_core_states,number_of_atom_types"
1830 : 8 : nctk_arrays(3)%name="dipole_core_valence"
1831 : 8 : nctk_arrays(3)%dtype="dp"
1832 : : nctk_arrays(3)%shape_str=&
1833 : : & "number_of_cartesian_directions,max_number_of_core_states,"// &
1834 : 8 : & "number_of_atoms,max_number_of_states,number_of_kpoints,number_of_spins"
1835 : 8 : nctk_arrays(4)%name="n_quantum_number_core"
1836 : 8 : nctk_arrays(4)%dtype="int"
1837 : 8 : nctk_arrays(4)%shape_str="max_number_of_core_states,number_of_atom_types"
1838 : 8 : nctk_arrays(5)%name="l_quantum_number_core"
1839 : 8 : nctk_arrays(5)%dtype="int"
1840 : 8 : nctk_arrays(5)%shape_str="max_number_of_core_states,number_of_atom_types"
1841 : 8 : nctk_arrays(6)%name="kappa_core"
1842 : 8 : nctk_arrays(6)%dtype="int"
1843 : 8 : nctk_arrays(6)%shape_str="max_number_of_core_states,number_of_atom_types"
1844 : 8 : nctk_arrays(7)%name="number_of_core_states"
1845 : 8 : nctk_arrays(7)%dtype="int"
1846 : 8 : nctk_arrays(7)%shape_str="number_of_atom_types"
1847 : 8 : nctk_arrays(8)%name="occupation_core"
1848 : 8 : nctk_arrays(8)%dtype="dp"
1849 : 8 : nctk_arrays(8)%shape_str="max_number_of_core_states,number_of_atom_types"
1850 : 8 : nctk_arrays(9)%name="max_occupation_core"
1851 : 8 : nctk_arrays(9)%dtype="dp"
1852 : 8 : nctk_arrays(9)%shape_str="max_number_of_core_states,number_of_atom_types"
1853 [ - + ]: 8 : NCF_CHECK(nctk_def_arrays(ncid, nctk_arrays))
1854 [ - + ]: 8 : NCF_CHECK(nctk_set_atomic_units(ncid, "eigenvalues_core"))
1855 [ - + ]: 8 : NCF_CHECK(nctk_set_atomic_units(ncid, "dipole_core_valence"))
1856 : : ! Write core states
1857 [ - + ]: 8 : NCF_CHECK(nctk_set_datamode(ncid))
1858 : 8 : varid=nctk_idname(ncid,"eigenvalues_core")
1859 [ - + ]: 8 : NCF_CHECK(nf90_put_var(ncid,varid,energy_cor))
1860 : 8 : varid=nctk_idname(ncid,"n_quantum_number_core")
1861 [ - + ]: 8 : NCF_CHECK(nf90_put_var(ncid,varid,ncor))
1862 : 8 : varid=nctk_idname(ncid,"l_quantum_number_core")
1863 [ - + ]: 8 : NCF_CHECK(nf90_put_var(ncid,varid,lcor))
1864 : 8 : varid=nctk_idname(ncid,"kappa_core")
1865 [ - + ]: 8 : NCF_CHECK(nf90_put_var(ncid,varid,kappacor))
1866 : 8 : varid=nctk_idname(ncid,"occupation_core")
1867 [ - + ]: 8 : NCF_CHECK(nf90_put_var(ncid,varid,occ_cor))
1868 : 8 : varid=nctk_idname(ncid,"max_occupation_core")
1869 [ - + ]: 8 : NCF_CHECK(nf90_put_var(ncid,varid,maxocc_cor))
1870 : 8 : varid=nctk_idname(ncid,"number_of_core_states")
1871 [ - + ]: 8 : NCF_CHECK(nf90_put_var(ncid,varid,nphicor_arr))
1872 : : ! Write eigenvalues
1873 : 8 : varid=nctk_idname(ncid,"eigenvalues")
1874 [ + + - + : 32 : NCF_CHECK(nf90_put_var(ncid,varid,reshape(eigen0,[mband,nkpt,nsppol])))
- - ]
1875 : : !Close file here because the rest has possibly to be written with collective I/O
1876 [ - + ]: 8 : NCF_CHECK(nf90_close(ncid))
1877 [ - + ]: 8 : ABI_FREE(nphicor_arr)
1878 : : ! ====> Standard FORTRAN binary file format
1879 : : else if (iomode==IO_MODE_FORTRAN_MASTER) then
1880 : : fformopt=612 ! MT 12sept21: change the OPT2 Fortran file format
1881 : : if (2*nphicor*natom*mband>2**30) fformopt=613 ! Format for large file records
1882 : : if (open_file(dtfil%fnameabo_app_opt2,msg,newunit=ount,form="unformatted",status="unknown")/= 0) then
1883 : : ABI_ERROR(msg)
1884 : : end if
1885 : : call hdr%fort_write(ount,fformopt,ierr,rewind=.true.)
1886 : : write(ount)(eigen0(jb),jb=1,mband*nkpt*nsppol)
1887 : : write(ount) nphicor
1888 : : do itypat=1,dtset%ntypat
1889 : : write(ount) atm(itypat)%ln_size
1890 : : do iln=1,nphicor
1891 : : write(ount) ncor(iln,itypat),lcor(iln,itypat),kappacor(iln,itypat),occ_cor(iln,itypat),maxocc_cor(iln,itypat),&
1892 : : energy_cor(iln,itypat)
1893 : : end do
1894 : : enddo
1895 : : else
1896 : : msg = "Wrong OPT2 file format!"
1897 : : ABI_BUG(msg)
1898 : : end if ! File format
1899 : : end if ! master node
1900 : 14 : call xmpi_bcast(iomode,master,spaceComm_w,ierr) ! Seems mandatory; why ?
1901 : 14 : call xmpi_bcast(fformopt,master,spaceComm_w,ierr)
1902 [ + - - + ]: 14 : iomode_etsf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
1903 : :
1904 [ - + ]: 14 : ABI_FREE(ncor)
1905 [ - + ]: 14 : ABI_FREE(lcor)
1906 [ - + ]: 14 : ABI_FREE(kappacor)
1907 [ - + ]: 14 : ABI_FREE(occ_cor)
1908 [ - + ]: 14 : ABI_FREE(maxocc_cor)
1909 [ - + ]: 14 : ABI_FREE(energy_cor)
1910 : : #if defined(HAVE_GPU_MARKERS)
1911 : : call nvtxEndRange()
1912 : : #endif
1913 : :
1914 : : !----------------------------------------------------------------------------------
1915 : : !4- Computation of <psi_n|p_i>(<phi_i|-i.nabla|phi_core>)
1916 : : !----------------------------------------------------------------------------------
1917 : :
1918 : : #if defined(HAVE_GPU_MARKERS)
1919 : : call nvtxStartRange("pre_spinkpt_loop")
1920 : : #endif
1921 : : !Prepare core-valence dipoles writing
1922 : : !In case of netCDF access to OPT2 file, prepare collective I/O
1923 [ + - ]: 14 : if (iomode == IO_MODE_ETSF) then
1924 [ + - ]: 14 : if (iomode_etsf_mpiio) then
1925 [ + + ]: 14 : if (i_am_master_spfft) then
1926 [ - + ]: 10 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),spaceComm_band))
1927 : 10 : varid=nctk_idname(ncid,"dipole_core_valence")
1928 [ + + ]: 10 : if (xmpi_comm_size(spaceComm_w)>1) then
1929 [ - + ]: 4 : NCF_CHECK(nctk_set_collective(ncid,varid))
1930 : : end if
1931 [ - + ]: 10 : NCF_CHECK(nctk_set_datamode(ncid))
1932 : : end if
1933 [ # # ]: 0 : else if (i_am_master) then
1934 [ # # ]: 0 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),xmpi_comm_self))
1935 : 0 : varid=nctk_idname(ncid,"dipole_core_valence")
1936 : : if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
1937 : : NCF_CHECK(nctk_set_collective(ncid,varid))
1938 : : end if
1939 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
1940 : : end if
1941 : : end if
1942 [ + - ]: 14 : if (iomode_etsf_mpiio) then
1943 : : !If MPI-IO, store only elements for one band
1944 [ - + + - : 56 : ABI_MALLOC(psinablapsi,(2,3,nphicor_lm,natom,1))
- + + - +
- + - - +
- + ]
1945 [ - + + - : 56 : ABI_MALLOC(psinablapsi2,(3,nphicor,natom,1))
- + + - +
- + - - +
- + ]
1946 [ + + ]: 14 : if (use_spinorbit) then
1947 [ - + + - : 18 : ABI_MALLOC(psinablapsi_soc,(2,3,nphicor_lm,natom,1))
- + + - -
+ - + ]
1948 : : end if
1949 : : else
1950 : : !If not, store the elements for all bands
1951 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,nphicor_lm,natom,mband))
# # # # #
# # # # #
# # # # #
# ]
1952 [ # # # # : 0 : ABI_MALLOC(psinablapsi2,(3,nphicor,natom,mband))
# # # # #
# # # # #
# # # # #
# ]
1953 [ # # ]: 0 : if (use_spinorbit) then
1954 [ # # # # : 0 : ABI_MALLOC(psinablapsi_soc,(2,3,nphicor_lm,natom,mband))
# # # # #
# # # # #
# # ]
1955 : : end if
1956 : : end if
1957 [ + + ]: 70 : pnp_size=size(psinablapsi2)
1958 : :
1959 : : !Determine if cprj datastructure is distributed over bands
1960 : 14 : mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
1961 : 14 : cprj_paral_band=(mband_cprj<mband)
1962 : : #if defined(HAVE_GPU_MARKERS)
1963 : : call nvtxEndRange()
1964 : : #endif
1965 : :
1966 : : !LOOP OVER SPINS
1967 : 14 : ibg=0
1968 : 14 : bdtot_index=0
1969 [ + + ]: 28 : do isppol=1,nsppol
1970 : :
1971 : : ! LOOP OVER k POINTS
1972 [ + + ]: 44 : do ikpt=1,nkpt
1973 : : #if defined(HAVE_GPU_MARKERS)
1974 : : call nvtxStartRange("kpoint step")
1975 : : #endif
1976 : :
1977 : 16 : etiq=ikpt+(isppol-1)*nkpt
1978 : 16 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1979 [ + + ]: 260 : master_spfftband=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
1980 : :
1981 : : ! Select k-points for current proc
1982 : 16 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt))
1983 [ + - ]: 16 : if (mykpt) then
1984 : : #if defined(HAVE_GPU_MARKERS)
1985 : : call nvtxStartRange("pawcprj_alloc")
1986 : : #endif
1987 : :
1988 : : ! Data depending on k-point
1989 : 16 : istwf_k=dtset%istwfk(ikpt)
1990 [ + + ]: 16 : cplex=2;if (istwf_k>1) cplex=1
1991 : :
1992 : : ! Extract cprj for this k-point
1993 [ + + ]: 16 : nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
1994 [ + + ]: 16 : if (mkmem*nsppol/=1) then
1995 : 4 : iorder_cprj=0
1996 [ - + + - : 208 : ABI_MALLOC(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
- + + - +
- + - - +
- + + + +
+ ]
1997 : 4 : call pawcprj_alloc(cprj_k_loc,0,dimcprj)
1998 : : call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
1999 : : & mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
2000 [ - + ]: 4 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2001 : : else
2002 : 12 : cprj_k_loc => cprj
2003 : : end if
2004 : :
2005 : : ! if cprj are distributed over bands, gather them (because we need to mix bands)
2006 [ + + ]: 16 : if (cprj_paral_band) then
2007 [ - + + - : 1032 : ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k))
- + + - +
- + - - +
- + + + +
+ ]
2008 : 8 : call pawcprj_alloc(cprj_k,0,dimcprj)
2009 : : call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom,my_nspinor*nband_cprj_k, &
2010 : : & my_nspinor*mpi_enreg%bandpp,&
2011 : 8 : & dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
2012 : : else
2013 : 8 : cprj_k => cprj_k_loc
2014 : : end if
2015 : : #if defined(HAVE_GPU_MARKERS)
2016 : : call nvtxEndRange()
2017 : : #endif
2018 : :
2019 : : #if defined(HAVE_GPU_MARKERS)
2020 : : call nvtxStartRange("jband_loop",5)
2021 : : #endif
2022 : : ! Loops over bands
2023 [ + + ]: 260 : do jb=1,nband_k
2024 : : !If MPI-IO, store only ib elements for each jb ; if not, store all (ib,jb) pairs
2025 [ - + ]: 244 : my_jb=merge(1,jb,iomode_etsf_mpiio)
2026 : :
2027 [ + + + + : 61432 : psinablapsi(:,:,:,:,my_jb)=zero
+ + + + ]
2028 [ + + + + : 12168 : psinablapsi2(:,:,:,my_jb)=zero
+ + ]
2029 [ + + + + : 42664 : if (use_spinorbit) psinablapsi_soc(:,:,:,:,my_jb)=zero
+ + + + +
+ ]
2030 : :
2031 : : ! Computation of <psi_n|p_i><phi_i|-i.nabla|phi_core>
2032 : : ! ----------------------------------------------------------------------------------
2033 : :
2034 : : ! Select bands for current proc
2035 : 244 : myband=.true.
2036 [ + + ]: 244 : if (mpi_enreg%paral_kgb==1) then
2037 : 120 : myband=(mod(jb-1,mpi_enreg%nproc_band)==mpi_enreg%me_band)
2038 : 124 : else if (xmpi_paral==1) then
2039 : 124 : myband=(abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-mpi_enreg%me_kpt)==0)
2040 : : end if
2041 [ + + ]: 244 : if (myband) then
2042 : :
2043 : 184 : jbsp=(jb-1)*my_nspinor
2044 : :
2045 : : ! 1-spinor case
2046 [ + + + - ]: 184 : if (dtset%nspinor==1.and.ncorespinor==1) then
2047 : 104 : jbsp=jbsp+1
2048 [ + + ]: 104 : if (cplex==1) then ! Real WF case
2049 [ + + ]: 146 : do iatom=1,natom
2050 : 104 : itypat=dtset%typat(iatom)
2051 : 104 : lmn_size=pawtab(itypat)%lmn_size
2052 : 104 : lmncmax=atm(itypat)%lmn_size
2053 [ + + ]: 1106 : do jlmn=1,lmn_size
2054 [ + + ]: 4328 : do ilmn=1,lmncmax
2055 : 3264 : ic=ilmn !atm(itypat)%indlmn(5,ilmn)
2056 : 3264 : cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
2057 : : psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
2058 [ + + ]: 14016 : & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
2059 : : end do !ilmn
2060 : : end do !jlmn
2061 : : end do !iatom
2062 : : else ! Complex WF case
2063 [ + + ]: 246 : do iatom=1,natom
2064 : 184 : itypat=dtset%typat(iatom)
2065 : 184 : lmn_size=pawtab(itypat)%lmn_size
2066 : 184 : lmncmax=atm(itypat)%lmn_size
2067 [ + + ]: 1846 : do jlmn=1,lmn_size
2068 [ + + ]: 8248 : do ilmn=1,lmncmax
2069 : 6464 : ic=ilmn !atm(itypat)%indlmn(5,ilmn)
2070 : 6464 : cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
2071 : 6464 : cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn)
2072 : : psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
2073 [ + + ]: 25856 : & -cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
2074 : : psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
2075 [ + + ]: 27456 : & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
2076 : : end do !ilmn
2077 : : end do !jlmn
2078 : : end do !iatom
2079 : : end if
2080 : :
2081 : : ! 2-spinor case
2082 [ + - + - ]: 80 : else if (dtset%nspinor==2.and.ncorespinor==2) then
2083 [ + + ]: 240 : do ispinor=1,my_nspinor
2084 : 160 : jbsp=jbsp+1
2085 [ + + ]: 760 : do iatom=1,natom
2086 : 520 : itypat=dtset%typat(iatom)
2087 : 520 : lmn_size=pawtab(itypat)%lmn_size
2088 : 520 : lmncmax=atm(itypat)%lmn_size
2089 [ + + ]: 5240 : do jlmn=1,lmn_size
2090 [ + + ]: 96280 : do ilmn=1,lmncmax
2091 : 91200 : is=atm(itypat)%indlmn(6,ilmn)
2092 [ + + ]: 95760 : if (modulo(jbsp,2)==modulo(is,2)) then ! Nabla is a spin-diagonal operator
2093 : 45600 : ic=atm(itypat)%indlmn(9,ilmn)
2094 [ + + ]: 45600 : if (ic>0) then
2095 : 31920 : cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
2096 : 31920 : cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn)
2097 : : psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
2098 : : & +cpnm1*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
2099 [ + + ]: 127680 : & -cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
2100 : : psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
2101 : : & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn) &
2102 [ + + ]: 127680 : & -cpnm2*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn)
2103 : : end if
2104 : : end if
2105 : : end do ! ilmn
2106 : : end do !jlmn
2107 : : end do !iatom
2108 : : end do !ispinor
2109 : : else
2110 : 0 : msg="Core and valence WF should have the same spinor representation!"
2111 : 0 : ABI_BUG(msg)
2112 : : !N. Brouwer initial coding: should be justified!
2113 : : !if (dtset%nspinor==1.and.ncorespinor==2) then ! Average core spinors
2114 : : ! psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
2115 : : ! & +half_sqrt2*(cpnm1*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
2116 : : ! & +cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn))
2117 : : ! psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
2118 : : ! & +half_sqrt2*(cpnm2*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
2119 : : ! & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn))
2120 : : !else if (dtset%nspinor==2.and.ncorespinor==1) then ! Average valence spinors
2121 : : ! psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
2122 : : ! & +half_sqrt2*cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
2123 : : ! psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
2124 : : ! & -half_sqrt2*cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
2125 : : !endif
2126 : : endif
2127 : :
2128 : : ! Spin-orbit coupling contribution:
2129 : : ! Sum_i,ss'[<psi^s_n|p_i><phi_i|1/4 Alpha^2 (Sigma^ss' X dV/dr)|phj_core^s'>]
2130 [ + + ]: 184 : if (use_spinorbit) then
2131 : : ! Add: Sum_i,ss'[<Psi^s_n|p_i> (Sigma^ss' X g_ij_core^s']
2132 : : ! where: g_ij_core^s = <Phi_i| 1/4 Alpha^2 dV/dr vec(r)/r) |Phj_core^s>
2133 : : ! Note that:
2134 : : ! if phi_cor_jlmn_cor(:) is up:
2135 : : ! phisocphj(:)%value(re:im,1,idir,ilmn,jlmn_cor) is (Sigma^up-up X g_ij_core^up)
2136 : : ! phisocphj(:)%value(re:im,2,idir,ilmn,jlmn_cor) is (Sigma^dn-up X g_ij_core^up)
2137 : : ! if phi_cor_jlmn_cor(:) is down:
2138 : : ! phisocphj(:)%value(re:im,1,idir,ilmn,jlmn_cor) is (Sigma^up-dn X g_ij_core^dn)
2139 : : ! phisocphj(:)%value(re:im,2,idir,ilmn,jlmn_cor) is (Sigma^dn-dn X g_ij_core^dn)
2140 : : ! Not compatible with parallelization over spinors
2141 : 80 : jbsp=1+(jb-1)*dtset%nspinor
2142 [ + + ]: 340 : do iatom=1,natom
2143 : 260 : itypat=dtset%typat(iatom)
2144 : 260 : lmn_size=pawtab(itypat)%lmn_size
2145 : 260 : lmncmax=atm(itypat)%lmn_size
2146 [ + + ]: 2620 : do jlmn=1,lmn_size
2147 [ + + ]: 48140 : do ilmn=1,lmncmax
2148 [ + - ]: 45600 : if(dtset%nspinor==1) then
2149 : : ic=ilmn !atm(itypat)%indlmn(5,ilmn)
2150 : : else
2151 : 45600 : ic=atm(itypat)%indlmn(9,ilmn)
2152 : : endif
2153 [ + + ]: 47880 : if (ic>0) then
2154 : 31920 : soc_ij => phisocphj(iatom)%value(:,:,:,jlmn,ilmn)
2155 : : !Contribution from real part of <Psi^s_n|p_i>
2156 : 31920 : cpnm1=cprj_k(iatom,jbsp )%cp(1,jlmn)
2157 : 31920 : cpnm2=cprj_k(iatom,jbsp+1)%cp(1,jlmn)
2158 [ + + ]: 127680 : do idir=1,3
2159 : : psinablapsi_soc(1,idir,ic,iatom,my_jb)=psinablapsi_soc(1,idir,ic,iatom,my_jb) &
2160 : 95760 : & +soc_ij(1,1,idir)*cpnm1+soc_ij(1,2,idir)*cpnm2
2161 : : psinablapsi_soc(2,idir,ic,iatom,my_jb)=psinablapsi_soc(2,idir,ic,iatom,my_jb) &
2162 : 127680 : & +soc_ij(2,1,idir)*cpnm1+soc_ij(2,2,idir)*cpnm2
2163 : : end do
2164 : : !Contribution from imaginary part of <Psi^s_n|p_i>
2165 : 31920 : cpnm1=cprj_k(iatom,jbsp )%cp(2,jlmn)
2166 : 31920 : cpnm2=cprj_k(iatom,jbsp+1)%cp(2,jlmn)
2167 [ + + ]: 127680 : do idir=1,3
2168 : : psinablapsi_soc(1,idir,ic,iatom,my_jb)=psinablapsi_soc(1,idir,ic,iatom,my_jb) &
2169 : 95760 : & +soc_ij(2,1,idir)*cpnm1+soc_ij(2,2,idir)*cpnm2
2170 : : psinablapsi_soc(2,idir,ic,iatom,my_jb)=psinablapsi_soc(2,idir,ic,iatom,my_jb) &
2171 : 127680 : & -soc_ij(1,1,idir)*cpnm1-soc_ij(1,2,idir)*cpnm2
2172 : : end do
2173 : : end if
2174 : : end do ! ilmn
2175 : : end do ! jlmn
2176 : : end do ! iatom
2177 : : end if ! use_spinorbit
2178 : :
2179 [ + - - + ]: 184 : if (iomode_etsf_mpiio.and.mpi_enreg%paral_spinor==1) then
2180 : : #if defined(HAVE_GPU_MARKERS)
2181 : : call nvtxStartRange("xmpi_sum",17)
2182 : : #endif
2183 : 0 : call timab(48,1,tsec)
2184 : 0 : call xmpi_sum_master(psinablapsi,master,spaceComm_spinor,ierr)
2185 : 0 : call timab(48,2,tsec)
2186 : : #if defined(HAVE_GPU_MARKERS)
2187 : : call nvtxEndRange()
2188 : : #endif
2189 : : end if
2190 : :
2191 : : end if ! myband
2192 : :
2193 : :
2194 : : ! Write to OPT2 file in case of MPI-IO
2195 [ + + ]: 260 : if (iomode_etsf_mpiio.and.i_am_master_spfft) then
2196 : : #if defined(HAVE_GPU_MARKERS)
2197 : : call nvtxStartRange("write_mpiio")
2198 : : #endif
2199 [ + + ]: 1288 : nc_start=[1,1,1,jb,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
2200 [ + + ]: 184 : if (myband) then
2201 [ + + - + : 18454 : if (use_spinorbit) psinablapsi=psinablapsi+psinablapsi_soc
+ + + + +
+ + + +
+ ]
2202 : : ! Take the square and sum over m
2203 [ + + ]: 582 : do iatom=1,natom
2204 : 428 : itypat=dtset%typat(iatom)
2205 : 428 : lmncmax=atm(itypat)%lmn_size
2206 [ + + ]: 428 : if(dtset%nspinor==2)lmncmax=lmncmax/2
2207 [ + + ]: 3366 : do ilmn=1,lmncmax
2208 : 2784 : ic=lmn2ln(ilmn,itypat)
2209 : : psinablapsi2(:,ic,iatom,my_jb)=psinablapsi2(:,ic,iatom,my_jb)+psinablapsi(1,:,ilmn,iatom,my_jb)**2+&
2210 [ + + ]: 11564 : & psinablapsi(2,:,ilmn,iatom,my_jb)**2
2211 : : enddo
2212 : : enddo
2213 [ + + ]: 1078 : nc_count=[3,nphicor,natom,1,1,1]
2214 : : else
2215 : 30 : nc_count=[0,0,0,0,0,0]
2216 : : end if
2217 [ - + ]: 184 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
2218 : : #if defined(HAVE_GPU_MARKERS)
2219 : : call nvtxEndRange()
2220 : : #endif
2221 : : end if
2222 : :
2223 : : end do ! jb
2224 : : #if defined(HAVE_GPU_MARKERS)
2225 : : call nvtxEndRange()
2226 : : #endif
2227 : :
2228 : : #if defined(HAVE_GPU_MARKERS)
2229 : : call nvtxStartRange("pawcprj_free")
2230 : : #endif
2231 [ + - ]: 16 : if (mkmem/=0) then
2232 : 16 : ibg = ibg + my_nspinor*nband_cprj_k
2233 : : end if
2234 : :
2235 [ + + ]: 16 : if (cprj_paral_band) then
2236 : 8 : call pawcprj_free(cprj_k)
2237 [ - + ]: 8 : ABI_FREE(cprj_k)
2238 : : end if
2239 [ + + ]: 16 : if (mkmem*nsppol/=1) then
2240 : 4 : call pawcprj_free(cprj_k_loc)
2241 [ - + ]: 4 : ABI_FREE(cprj_k_loc)
2242 : : end if
2243 : : #if defined(HAVE_GPU_MARKERS)
2244 : : call nvtxEndRange()
2245 : : #endif
2246 : :
2247 : : ! Write to OPT2 file if not MPI-IO
2248 : :
2249 : : #if defined(HAVE_GPU_MARKERS)
2250 : : call nvtxStartRange("write")
2251 : : #endif
2252 : : ! >>> Reduction in case of parallelism
2253 [ - + ]: 16 : if (.not.iomode_etsf_mpiio) then
2254 : 0 : call timab(48,1,tsec)
2255 : 0 : call xmpi_sum_master(psinablapsi,master,spaceComm_bandspinor,ierr)
2256 : 0 : call timab(48,2,tsec)
2257 [ # # ]: 0 : if (use_spinorbit) then
2258 : 0 : call xmpi_sum_master(psinablapsi_soc,master,spaceComm_band,ierr)
2259 [ # # # # : 0 : psinablapsi=psinablapsi+psinablapsi_soc
# # # # #
# # # ]
2260 : : end if
2261 : : ! Take the square and sum over m
2262 [ # # ]: 0 : do jb=1,nband_k
2263 : 0 : my_jb=merge(1,jb,iomode_etsf_mpiio)
2264 [ # # ]: 0 : do iatom=1,natom
2265 : 0 : itypat=dtset%typat(iatom)
2266 : 0 : lmncmax=atm(itypat)%lmn_size
2267 [ # # ]: 0 : if(dtset%nspinor==2)lmncmax=lmncmax/2
2268 [ # # ]: 0 : do ilmn=1,lmncmax
2269 : 0 : ic=lmn2ln(ilmn,itypat)
2270 : : psinablapsi2(:,ic,iatom,my_jb)=psinablapsi2(:,ic,iatom,my_jb)+psinablapsi(1,:,ilmn,iatom,my_jb)**2+&
2271 [ # # ]: 0 : & psinablapsi(2,:,ilmn,iatom,my_jb)**2
2272 : : enddo
2273 : : enddo
2274 : : enddo
2275 : : end if
2276 : :
2277 : : ! >>> This my kpt and I am the master node: I write the data
2278 : : if (.not.iomode_etsf_mpiio) then
2279 [ # # ]: 0 : if (i_am_master) then
2280 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
2281 [ # # ]: 0 : nc_start=[1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
2282 [ # # ]: 0 : nc_count=[3,nphicor,natom,mband,1,1]
2283 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
2284 : : else
2285 [ # # ]: 0 : if (fformopt==612) then ! New OPT2 file format
2286 : 0 : write(ount) (((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
2287 : 0 : write(ount) (((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
2288 : 0 : write(ount) (((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
2289 [ # # ]: 0 : else if (fformopt==613) then ! Large OPT2 file format
2290 [ # # ]: 0 : do jb=1,nband_k
2291 : 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
2292 : 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
2293 : 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
2294 : : end do
2295 : : else ! Old OPT2 file format
2296 : : !The old writing was not efficient (indexes order is bad)
2297 [ # # ]: 0 : do iatom=1,natom
2298 [ # # # # ]: 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
2299 [ # # # # ]: 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
2300 [ # # # # ]: 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
2301 : : end do
2302 : : end if
2303 : : end if
2304 : :
2305 : : ! >>> This my kpt and I am not the master node: I send the data
2306 [ # # ]: 0 : else if (i_am_master_band.and.i_am_master_spfft) then
2307 [ # # ]: 0 : if (mpi_enreg%me_kpt/=master_spfftband) then
2308 : 0 : ABI_BUG('Problem with band communicator!')
2309 : : end if
2310 : 0 : call xmpi_exch(psinablapsi2,pnp_size,mpi_enreg%me_kpt,psinablapsi2,master,spaceComm_kpt,etiq,ierr)
2311 : : end if
2312 : : end if
2313 : : #if defined(HAVE_GPU_MARKERS)
2314 : : call nvtxEndRange()
2315 : : #endif
2316 : :
2317 : : ! >>> This is not my kpt and I am the master node: I receive the data and I write
2318 [ # # ]: 0 : elseif ((.not.iomode_etsf_mpiio).and.i_am_master) then ! mykpt
2319 : 0 : sender=master_spfftband
2320 : 0 : call xmpi_exch(psinablapsi2,pnp_size,sender,psinablapsi2,master,spaceComm_kpt,etiq,ierr)
2321 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
2322 [ # # ]: 0 : nc_start=[1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
2323 [ # # ]: 0 : nc_count=[3,nphicor,natom,mband,1,1]
2324 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
2325 : : else
2326 [ # # ]: 0 : if (fformopt==612) then ! New OPT2 file format
2327 : 0 : write(ount) (((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
2328 : 0 : write(ount) (((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
2329 : 0 : write(ount) (((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
2330 [ # # ]: 0 : else if (fformopt==613) then ! Large OPT2 file format
2331 [ # # ]: 0 : do jb=1,nband_k
2332 : 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
2333 : 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
2334 : 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
2335 : : end do
2336 : : else ! Old OPT2 file format
2337 : : !The old writing was not efficient (indexes order is bad)
2338 [ # # ]: 0 : do iatom=1,natom
2339 [ # # # # ]: 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
2340 [ # # # # ]: 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
2341 [ # # # # ]: 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
2342 : : end do
2343 : : end if
2344 : : end if
2345 : : end if ! mykpt
2346 : :
2347 : 30 : bdtot_index=bdtot_index+nband_k
2348 : :
2349 : : ! End loop on spin,kpt
2350 : : #if defined(HAVE_GPU_MARKERS)
2351 : : call nvtxEndRange()
2352 : : #endif
2353 : : end do ! ikpt
2354 : : end do !isppol
2355 : :
2356 : : !Close file
2357 [ + + + + ]: 14 : if (i_am_master.or.(iomode_etsf_mpiio.and.i_am_master_spfft)) then
2358 [ + - ]: 10 : if (iomode==IO_MODE_ETSF) then
2359 [ - + ]: 10 : NCF_CHECK(nf90_close(ncid))
2360 : : else
2361 : 0 : ierr=close_unit(ount,msg)
2362 [ # # ]: 0 : ABI_CHECK(ierr==0,"Error while closing OPT2 file")
2363 : : end if
2364 : : end if
2365 : :
2366 : : !Datastructures deallocations
2367 [ + + ]: 30 : do itypat=1,dtset%ntypat
2368 : 30 : call destroy_atomorb(atm(itypat))
2369 : : enddo
2370 [ + - + + : 30 : ABI_FREE(atm)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
2371 [ - + ]: 14 : ABI_FREE(psinablapsi)
2372 [ - + ]: 14 : ABI_FREE(psinablapsi2)
2373 [ - + ]: 14 : ABI_FREE(lmn2ln)
2374 [ + + ]: 14 : if (use_spinorbit) then
2375 [ - + ]: 6 : ABI_FREE(psinablapsi_soc)
2376 [ + + ]: 27 : do iatom=1,natom
2377 [ - + ]: 27 : ABI_FREE(phisocphj(iatom)%value)
2378 : : end do
2379 [ + - + + : 27 : ABI_FREE(phisocphj)
- + - + ]
2380 : : end if
2381 [ + - ]: 14 : if (.not.already_has_nabla) then
2382 [ + + ]: 30 : do itypat=1,dtset%ntypat
2383 [ + - ]: 16 : if (allocated(pawtab(itypat)%nabla_ij)) then
2384 : 16 : ABI_FREE(pawtab(itypat)%nabla_ij)
2385 : 16 : pawtab(itypat)%has_nabla=0
2386 : : end if
2387 [ + + ]: 30 : if (allocated(pawtab(itypat)%nabla_im_ij)) then
2388 : 6 : ABI_FREE(pawtab(itypat)%nabla_im_ij)
2389 : : end if
2390 : : end do
2391 : : end if
2392 : :
2393 : : ABI_NVTX_END_RANGE()
2394 : 14 : call timab(1185,2,tsec)
2395 : : DBG_EXIT("COLL")
2396 : :
2397 [ - + - + : 28 : end subroutine optics_paw_core
- - - - -
+ - + - +
- + - + -
+ - + ]
2398 : : !!***
2399 : :
2400 : : !----------------------------------------------------------------------
2401 : :
2402 : : !!****f* m_paw_optics/linear_optics_paw
2403 : : !! NAME
2404 : : !! linear_optics_paw
2405 : : !!
2406 : : !! FUNCTION
2407 : : !! This program computes the elements of the optical frequency dependent
2408 : : !! linear susceptiblity using matrix elements <-i Nabla> obtained from a
2409 : : !! PAW ground state calculation. It uses formula 17 from Gadoc et al,
2410 : : !! Phys. Rev. B 73, 045112 (2006) [[cite:Gajdo2006]] together with a scissors correction. It uses
2411 : : !! a Kramers-Kronig transform to compute the real part from the imaginary part, and
2412 : : !! it will work on all types of unit cells. It outputs all tensor elements of
2413 : : !! both the real and imaginary parts.
2414 : : !!
2415 : : !! INPUTS
2416 : : !! filnam: base of file names to read data from
2417 : : !! mpi_enreg: mpi set up variable, not used in this code
2418 : : !!
2419 : : !! OUTPUT
2420 : : !! _real and _imag output files
2421 : : !!
2422 : : !! NOTES
2423 : : !! This routine is not tested
2424 : : !!
2425 : : !! SOURCE
2426 : :
2427 : 0 : subroutine linear_optics_paw(filnam,filnam_out)
2428 : :
2429 : : !Arguments -----------------------------------
2430 : : !scalars
2431 : : character(len=fnlen),intent(in) :: filnam,filnam_out
2432 : :
2433 : : !Local variables-------------------------------
2434 : : integer,parameter :: master=0
2435 : : integer :: iomode,bantot,bdtot_index,fform1,headform
2436 : : integer :: iband,ierr,ii,ikpt,iom,iout,isppol,isym,jband,jj,me,mband
2437 : : integer :: method,mom,nband_k,nkpt,nspinor,nsppol,nsym,occopt,only_check
2438 : : integer :: rdwr,spaceComm,inpunt,reunt,imunt,wfunt
2439 : 0 : integer,allocatable :: nband(:),symrel(:,:,:)
2440 : : real(dp) :: del,dom,fij,gdelta,omin,omax,paijpbij(2),mbpt_sciss,wij,ucvol
2441 : : real(dp) :: diffwp, diffwm
2442 : : real(dp) :: e2rot(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),rprimdinv(3,3),symd(3,3),symdinv(3,3)
2443 : 0 : real(dp),allocatable :: e1(:,:,:),e2(:,:,:,:),epsilon_tot(:,:,:,:),eigen0(:),eig0_k(:)
2444 : 0 : real(dp),allocatable :: kpts(:,:),occ(:),occ_k(:),oml1(:),wtk(:)
2445 : 0 : complex(dp),allocatable :: eps_work(:)
2446 : : character(len=fnlen) :: filnam1,filnam_gen
2447 : : character(len=500) :: msg
2448 : 0 : type(hdr_type) :: hdr
2449 : : type(wffile_type) :: wff1
2450 : : !arrays
2451 : 0 : real(dp),allocatable :: psinablapsi(:,:,:,:)
2452 : : ! *********************************************************************************
2453 : :
2454 : : DBG_ENTER("COLL")
2455 : : ABI_NVTX_START_RANGE(NVTX_LINEAR_OPTICS_PAW)
2456 : :
2457 : : !write(std_out,'(a)')' Give the name of the output file ...'
2458 : : !read(std_in, '(a)') filnam_out
2459 : : !write(std_out,'(a)')' The name of the output file is :',filnam_out
2460 : :
2461 : : !Read data file
2462 [ # # ]: 0 : if (open_file(filnam,msg,newunit=inpunt,form='formatted') /= 0 ) then
2463 : 0 : ABI_ERROR(msg)
2464 : : end if
2465 : :
2466 : 0 : rewind(inpunt)
2467 : 0 : read(inpunt,*)
2468 : 0 : read(inpunt,'(a)')filnam_gen ! generic name for the files
2469 [ # # ]: 0 : filnam1=trim(filnam_gen)//'_OPT' ! nabla matrix elements file
2470 : :
2471 : : !Open the Wavefunction and optic files
2472 : : !These default values are typical of sequential use
2473 : 0 : iomode=IO_MODE_FORTRAN ; spaceComm=xmpi_comm_self; me=0
2474 : :
2475 : : ! Read the header of the optic files
2476 : 0 : call hdr%from_fname(filnam1, fform1, spaceComm)
2477 : 0 : call hdr%free()
2478 [ # # ]: 0 : if (fform1 /= 610) then
2479 : 0 : ABI_ERROR("Abinit8 requires an OPT file with fform = 610")
2480 : : end if
2481 : :
2482 : : !Open the conducti optic files
2483 : 0 : wfunt = get_unit()
2484 : 0 : call WffOpen(iomode,spaceComm,filnam1,ierr,wff1,master,me,wfunt)
2485 : :
2486 : : !Read the header from Ground state file
2487 : 0 : rdwr=1
2488 : 0 : call hdr_io(fform1,hdr,rdwr,wff1)
2489 : :
2490 : : !Extract info from the header
2491 : 0 : headform=hdr%headform
2492 : 0 : bantot=hdr%bantot
2493 : 0 : nkpt=hdr%nkpt
2494 [ # # # # : 0 : ABI_MALLOC(kpts,(3,nkpt))
# # # # #
# # # ]
2495 [ # # # # : 0 : ABI_MALLOC(wtk,(nkpt))
# # # # #
# # # ]
2496 [ # # # # ]: 0 : kpts(:,:)=hdr%kptns(:,:)
2497 [ # # ]: 0 : wtk(:)=hdr%wtk(:)
2498 : 0 : nspinor=hdr%nspinor
2499 : 0 : nsppol=hdr%nsppol
2500 : 0 : occopt=hdr%occopt
2501 [ # # # # ]: 0 : rprimd(:,:)=hdr%rprimd(:,:)
2502 : 0 : rprimdinv(:,:) = rprimd(:,:)
2503 : 0 : call matrginv(rprimdinv,3,3) ! need the inverse of rprimd to symmetrize the tensors
2504 [ # # # # : 0 : ABI_MALLOC(nband,(nkpt*nsppol))
# # # # #
# # # ]
2505 [ # # # # : 0 : ABI_MALLOC(occ,(bantot))
# # # # #
# # # ]
2506 [ # # ]: 0 : occ(1:bantot)=hdr%occ(1:bantot)
2507 [ # # ]: 0 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
2508 : 0 : nsym=hdr%nsym
2509 [ # # # # : 0 : ABI_MALLOC(symrel,(3,3,nsym))
# # # # #
# # # ]
2510 [ # # # # : 0 : symrel(:,:,:)=hdr%symrel(:,:,:)
# # ]
2511 : :
2512 : : !Get mband, as the maximum value of nband(nkpt)
2513 [ # # ]: 0 : mband=maxval(nband(:))
2514 : :
2515 : : !get ucvol etc.
2516 : 0 : iout = -1
2517 : 0 : call metric(gmet,gprimd,iout,rmet,rprimd,ucvol)
2518 : :
2519 : 0 : write(std_out,*)
2520 : 0 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1)
2521 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2)
2522 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3)
2523 : 0 : write(std_out,*)
2524 : 0 : write(std_out,'(a,3f10.5,a)' )' rprimdinv =',rprimdinv(1:3,1)
2525 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimdinv(1:3,2)
2526 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimdinv(1:3,3)
2527 : 0 : write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband
2528 : :
2529 : : !get eigen0
2530 [ # # # # : 0 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
# # # # #
# # # ]
2531 : 0 : read(wfunt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
2532 : :
2533 : 0 : read(inpunt,*)mbpt_sciss
2534 : 0 : read(inpunt,*)dom,omin,omax,mom
2535 : 0 : close(inpunt)
2536 : :
2537 [ # # # # : 0 : ABI_MALLOC(oml1,(mom))
# # # # #
# # # ]
2538 [ # # # # : 0 : ABI_MALLOC(e1,(3,3,mom))
# # # # #
# # # ]
2539 [ # # # # : 0 : ABI_MALLOC(e2,(2,3,3,mom))
# # # # #
# # # ]
2540 [ # # # # : 0 : ABI_MALLOC(epsilon_tot,(2,3,3,mom))
# # # # ]
2541 [ # # # # : 0 : ABI_MALLOC(eps_work,(mom))
# # # # #
# # # ]
2542 : 0 : del=(omax-omin)/(mom-1)
2543 [ # # ]: 0 : do iom=1,mom
2544 : 0 : oml1(iom)=omin+dble(iom-1)*del
2545 : : end do
2546 : 0 : write(std_out,'(a,i8,4f10.5,a)')' npts,omin,omax,width,mbpt_sciss =',mom,omin,omax,dom,mbpt_sciss,' Ha'
2547 : :
2548 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,mband,mband))
# # # # #
# # # # #
# # ]
2549 : :
2550 : : !loop over spin components
2551 [ # # ]: 0 : do isppol=1,nsppol
2552 : : bdtot_index = 0
2553 : : ! loop over k points
2554 [ # # ]: 0 : do ikpt=1,nkpt
2555 : : !
2556 : : ! number of bands for this k point
2557 : 0 : nband_k=nband(ikpt+(isppol-1)*nkpt)
2558 [ # # # # : 0 : ABI_MALLOC(eig0_k,(nband_k))
# # # # #
# # # ]
2559 [ # # # # : 0 : ABI_MALLOC(occ_k,(nband_k))
# # # # ]
2560 : : ! eigenvalues for this k-point
2561 [ # # ]: 0 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
2562 : : ! occupation numbers for this k-point
2563 [ # # ]: 0 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
2564 : : ! values of -i*nabla matrix elements for this k point
2565 [ # # # # : 0 : psinablapsi=zero
# # # # ]
2566 [ # # # # : 0 : read(wfunt)((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k)
# # # # ]
2567 [ # # # # : 0 : read(wfunt)((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k)
# # # # ]
2568 [ # # # # : 0 : read(wfunt)((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k)
# # # # ]
2569 : :
2570 : : ! occupation numbers for k-point
2571 [ # # ]: 0 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
2572 : : ! accumulate e2 for this k point, Eq. 17 from PRB 73, 045112 (2006) [[cite:Gajdo2006]]
2573 [ # # ]: 0 : do iband = 1, nband_k
2574 [ # # ]: 0 : do jband = 1, nband_k
2575 : 0 : fij = occ_k(iband) - occ_k(jband) !occ number difference
2576 : 0 : wij = eig0_k(iband) - eig0_k(jband) !energy difference
2577 [ # # ]: 0 : if (abs(fij) > zero) then ! only consider states of differing occupation numbers
2578 [ # # ]: 0 : do ii = 1, 3
2579 [ # # ]: 0 : do jj = 1, 3
2580 : : paijpbij(1) = psinablapsi(1,ii,iband,jband)*psinablapsi(1,jj,iband,jband) + &
2581 : 0 : & psinablapsi(2,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
2582 : : paijpbij(2) = psinablapsi(2,ii,iband,jband)*psinablapsi(1,jj,iband,jband) - &
2583 : 0 : & psinablapsi(1,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
2584 [ # # ]: 0 : do iom = 1, mom
2585 : : ! original version
2586 : : ! diffw = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
2587 : : ! gdelta = exp(-diffw*diffw/(4.0*dom*dom))/(2.0*dom*sqrt(pi)) ! delta fnc resolved as Gaussian
2588 : : ! e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(oml1(iom)*oml1(iom))
2589 : : ! e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(oml1(iom)*oml1(iom))
2590 : 0 : diffwm = wij - mbpt_sciss + oml1(iom) ! apply scissors term here
2591 : 0 : diffwp = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
2592 : 0 : gdelta = exp(-diffwp*diffwp/(4.0*dom*dom))/(2.0*dom*sqrt(pi))
2593 : 0 : e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(wij*wij)
2594 : 0 : e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(wij*wij)
2595 : : end do ! end loop over spectral points
2596 : : end do ! end loop over jj = 1, 3
2597 : : end do ! end loop over ii = 1, 3
2598 : : end if ! end selection on fij /= 0
2599 : : end do ! end loop over jband
2600 : : end do ! end loop over iband
2601 : :
2602 [ # # ]: 0 : ABI_FREE(eig0_k)
2603 [ # # ]: 0 : ABI_FREE(occ_k)
2604 : 0 : bdtot_index=bdtot_index+nband_k
2605 : : end do ! end loop over k points
2606 : : end do ! end loop over spin polarizations
2607 : :
2608 : : !here apply nsym symrel transformations to reconstruct full tensor from IBZ part
2609 [ # # # # : 0 : epsilon_tot(:,:,:,:) = zero
# # # # ]
2610 [ # # ]: 0 : do isym = 1, nsym
2611 [ # # # # : 0 : symd(:,:)=matmul(rprimd(:,:),matmul(symrel(:,:,isym),rprimdinv(:,:)))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
2612 : 0 : symdinv(:,:)=symd(:,:)
2613 : 0 : call matrginv(symdinv,3,3)
2614 [ # # ]: 0 : do iom = 1, mom
2615 [ # # # # : 0 : e2rot(:,:)=matmul(symdinv(:,:),matmul(e2(1,:,:,iom),symd(:,:)))
# # # # #
# # # # #
# # # # #
# ]
2616 [ # # # # ]: 0 : epsilon_tot(2,:,:,iom) = epsilon_tot(2,:,:,iom)+e2rot(:,:)/nsym
2617 : : end do
2618 : : end do
2619 : :
2620 : : !generate e1 from e2 via KK transforma
2621 : 0 : method=0 ! use naive integration ( = 1 for simpson)
2622 : 0 : only_check=0 ! compute real part of eps in kk routine
2623 [ # # ]: 0 : do ii = 1, 3
2624 [ # # ]: 0 : do jj = 1, 3
2625 [ # # ]: 0 : eps_work(:) = cmplx(0.0,epsilon_tot(2,ii,jj,:), kind=dp)
2626 : 0 : call kramerskronig(mom,oml1,eps_work,method,only_check)
2627 [ # # ]: 0 : epsilon_tot(1,ii,jj,:) = real(eps_work(:))
2628 [ # # # # ]: 0 : if (ii /= jj) epsilon_tot(1,ii,jj,:) = epsilon_tot(1,ii,jj,:)- 1.0
2629 : : end do ! end loop over jj
2630 : : end do ! end loop over ii
2631 : :
2632 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_imag',msg,newunit=reunt,form='formatted') /= 0) then
2633 : 0 : ABI_ERROR(msg)
2634 : : end if
2635 : :
2636 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_real',msg,unit=imunt,form='formatted') /= 0) then
2637 : 0 : ABI_ERROR(msg)
2638 : : end if
2639 : :
2640 : 0 : write(reunt,'(a12,6a13)')' # Energy/Ha ','eps_2_xx','eps_2_yy','eps_2_zz',&
2641 : 0 : & 'eps_2_yz','eps_2_xz','eps_2_xy'
2642 : 0 : write(imunt,'(a12,6a13)')' # Energy/Ha ','eps_1_xx','eps_1_yy','eps_1_zz',&
2643 : 0 : & 'eps_1_yz','eps_1_xz','eps_1_xy'
2644 : :
2645 [ # # ]: 0 : do iom = 1, mom
2646 : 0 : write(reunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',&
2647 : 0 : & epsilon_tot(2,1,1,iom),' ',epsilon_tot(2,2,2,iom),' ',epsilon_tot(2,3,3,iom),' ',&
2648 : 0 : & epsilon_tot(2,2,3,iom),' ',epsilon_tot(2,1,3,iom),' ',epsilon_tot(2,1,2,iom)
2649 : 0 : write(imunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',&
2650 : 0 : & epsilon_tot(1,1,1,iom),' ',epsilon_tot(1,2,2,iom),' ',epsilon_tot(1,3,3,iom),' ',&
2651 : 0 : & epsilon_tot(1,2,3,iom),' ',epsilon_tot(1,1,3,iom),' ',epsilon_tot(1,1,2,iom)
2652 : : end do
2653 : :
2654 : 0 : close(reunt)
2655 : 0 : close(imunt)
2656 : :
2657 [ # # ]: 0 : ABI_FREE(nband)
2658 [ # # ]: 0 : ABI_FREE(oml1)
2659 [ # # ]: 0 : ABI_FREE(e2)
2660 [ # # ]: 0 : ABI_FREE(e1)
2661 [ # # ]: 0 : ABI_FREE(occ)
2662 [ # # ]: 0 : ABI_FREE(psinablapsi)
2663 [ # # ]: 0 : ABI_FREE(eigen0)
2664 [ # # ]: 0 : ABI_FREE(wtk)
2665 [ # # ]: 0 : ABI_FREE(kpts)
2666 : :
2667 : 0 : call hdr%free()
2668 : :
2669 : : ABI_NVTX_END_RANGE()
2670 : : DBG_EXIT("COLL")
2671 : :
2672 [ # # # # : 0 : end subroutine linear_optics_paw
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
2673 : : !!***
2674 : :
2675 : : !----------------------------------------------------------------------
2676 : :
2677 : : !!****f* m_paw_optics/pawnabla_soc_init
2678 : : !! NAME
2679 : : !! pawnabla_soc_init
2680 : : !!
2681 : : !! FUNCTION
2682 : : !! Compute the PAW SOC contribution(s) to the momentum PAW matrix elements,
2683 : : !! i.e. <Phi_i|1/4 Alpha^2 dV/dr (Sigma X vec(r)/r) |Phi_j>
2684 : : !! where:
2685 : : !! {Phi_i}= AE partial waves
2686 : : !! Alpha = inverse of fine structure constant
2687 : : !! Sigma^alpha,beta= Pauli matrices
2688 : : !! X = cross product
2689 : : !!
2690 : : !! There are 2 typical uses:
2691 : : !! - Valence-valence terms: Phi_i and Phi_j are PAW AE partial waves (unpolarized)
2692 : : !! - Core-valence terms: Phi_i is are AE partial waves and Phi_j are core AE wave-functions (spinors)
2693 : : !!
2694 : : !! In practice we compute:
2695 : : !! (Sigma^up-up X g_ij) and (Sigma^up-dn X g_ij) (X = vector cross product)
2696 : : !! (Sigma^dn-up X g_ij) and (Sigma^dn-dn X g_ij)
2697 : : !! where:
2698 : : !! g_ij= 1/4 Alpha^2 Int_[Phi_i(r)/r Phi_j(r)/r dV(r)/dr r^2 dr] . Gvec_ij
2699 : : !! and Gvec_ij= Int[S_limi S_ljmj vec(r)/r dOmega] (Gaunt coefficients)
2700 : : !!
2701 : : !! COPYRIGHT
2702 : : !! Copyright (C) 2021-2026 ABINIT group (NBrouwer,MT)
2703 : : !! This file is distributed under the terms of the
2704 : : !! GNU General Public License, see ~ABINIT/COPYING
2705 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2706 : : !!
2707 : : !! INPUTS
2708 : : !! el_temp=electronic temperature (hartree)
2709 : : !! ixc= choice of exchange-correlation scheme (see above, and below)
2710 : : !! my_natom=number of atoms treated by current processor
2711 : : !! natom=total number of atoms in cell
2712 : : !! nspden=number of spin-density components
2713 : : !! ntypat=number of types of atoms in unit cell.
2714 : : !! option_core=Type of calculation: 0=valence-valence, 1=core-valence
2715 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
2716 : : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
2717 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
2718 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
2719 : : !! pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
2720 : : !! spnorbscl=scaling factor for spin-orbit coupling
2721 : : !! typat(natom) =Type of each atoms
2722 : : !! xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
2723 : : !! xc_taupos= lowest allowed kinetic energy density (for mGGA XC functionals)
2724 : : !! znucl(ntypat)=gives the nuclear charge for all types of atoms
2725 : : !! [atm <type(paw_atomorb_type)>]=--optional-- structure containing core info
2726 : : !! [mpi_atmtab(:)]=--optional-- indexes of the atoms treated by current proc
2727 : : !! [comm_atom]=--optional-- MPI communicator over atoms
2728 : : !!
2729 : : !! OUTPUT
2730 : : ! phisocphj(dtset%natom) <type(coeff4_type)>= stores soc coefficients:
2731 : : !!! If option_core==0 or nspinor_cor==1:
2732 : : !! phisocphj(iat)%value(1,1,idir,ilmn,jlmn) is real part of (Sigma^up-up X g_ij)
2733 : : !! phisocphj(iat)%value(2,1,idir,ilmn,jlmn) is imaginary part of (Sigma^up-up X g_ij)
2734 : : !! phisocphj(iat)%value(1,2,idir,ilmn,jlmn) is real part of (Sigma^up-dn X g_ij)
2735 : : !! phisocphj(iat)%value(2,2,idir,ilmn,jlmn) is imaginary part of (Sigma^up-dn X g_ij)
2736 : : !! If option_core==1 and nspinor_cor==2 (core-valence with spinorial core WF):
2737 : : !! phisocphj(iat)%value(1,1,idir,ilmn,2*jlmn-1) is real part of (Sigma^up-up X g_ij^up)
2738 : : !! phisocphj(iat)%value(2,1,idir,ilmn,2*jlmn-1) is imaginary part of (Sigma^up-up X g_ij^up)
2739 : : !! phisocphj(iat)%value(1,2,idir,ilmn,2*jlmn-1) is real part of (Sigma^dn-up X g_ij^dn)
2740 : : !! phisocphj(iat)%value(2,2,idir,ilmn,2*jlmn-1) is imaginary part of (Sigma^dn-up X g_ij^dn)
2741 : : !! phisocphj(iat)%value(1,1,idir,ilmn,2*jlmn ) is real part of (Sigma^up-dn X g_ij^up)
2742 : : !! phisocphj(iat)%value(2,1,idir,ilmn,2*jlmn ) is imaginary part of (Sigma^up-dn X g_ij^up)
2743 : : !! phisocphj(iat)%value(1,2,idir,ilmn,2*jlmn ) is real part of (Sigma^dn-dn X g_ij^dn)
2744 : : !! phisocphj(iat)%value(2,2,idir,ilmn,2*jlmn ) is imaginary part of (Sigma^dn-dn X g_ij^dn)
2745 : : !! (idir=cartesian direction)
2746 : : !!
2747 : : !! SIDE EFFECTS
2748 : : !!
2749 : : !! NOTES
2750 : : !! If Phi_j is not polarized,
2751 : : !! (Sigma^dn-dn X g_ij)=-(Sigma^up-up X g_ij)
2752 : : !! (Sigma^dn-up X g_ij)= (Sigma^up-dn X g_ij)^*
2753 : : !! So, we store only 2 components.
2754 : : !! If Phi_j is polarized, the spin component is included in the last dimension of
2755 : : !! phisocphj(iat)%value, i.e. lmn_size_cor=2*lmn_size
2756 : : !!
2757 : : !! SOURCE
2758 : :
2759 : 11 : subroutine pawnabla_soc_init(el_temp,phisocphj,option_core,ixc,my_natom,natom,nspden,ntypat,pawang, &
2760 : 22 : & pawrad,pawrhoij,pawtab,pawxcdev,spnorbscl,typat,xc_denpos,xc_taupos,znucl, &
2761 [ + - + + : 11 : & atm,mpi_atmtab,comm_atom) ! Optional arguments
- + ]
2762 : :
2763 : : !Arguments ------------------------------------
2764 : : !scalars
2765 : : integer,intent(in) :: ixc,my_natom,natom,nspden,ntypat,option_core,pawxcdev
2766 : : integer,optional,intent(in) :: comm_atom
2767 : : real(dp),intent(in) :: el_temp,spnorbscl,xc_denpos,xc_taupos
2768 : : type(pawang_type),intent(in) :: pawang
2769 : : !arrays
2770 : : integer,intent(in) :: typat(natom)
2771 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
2772 : : real(dp),intent(in) :: znucl(ntypat)
2773 : : type(coeff5_type),allocatable,target,intent(inout) :: phisocphj(:)
2774 : : type(pawrad_type),target,intent(in) :: pawrad(ntypat)
2775 : : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
2776 : : type(pawtab_type),target,intent(in) :: pawtab(ntypat)
2777 : : type(atomorb_type), intent(in), target, optional :: atm(ntypat)
2778 : :
2779 : : !Local variables-------------------------------
2780 : : !scalars
2781 : : real(dp),parameter :: one_over_fourpi = one/sqrt(four_pi)
2782 : : real(dp),parameter :: sqr_fourpi_over_3 = sqrt(four_pi/3)
2783 : : real(dp),parameter :: QuarterFineStruct2=(half/InvFineStruct)**2
2784 : : real(dp),parameter :: hyb_mixing_ = 0.0_dp ! Fake value to be updated in the future
2785 : : integer :: iatom,iatom_tot,itypat,ii,jj,ierr,ipts,ignt,sgnkappa
2786 : : integer :: idum,option,usenhat,usekden,usecore,xclevel,nkxc,my_comm_atom
2787 : : integer :: mesh_size,lmn_size,lmn2_size,lmn_size_j,lmn_size_cor
2788 : : integer :: lm_size,ln_size,ln_size_j
2789 : : integer :: ilmn,ilm,iln,jl,jm,jm_re,jm_im,jlmn,jlm,jlm_re,jlm_im,jln,js,klm_re,klm_im
2790 : : logical :: my_atmtab_allocated,paral_atom
2791 : : real(dp) :: avg,cgc,compch_sph_dum,eexc_dum,ssxc_dum,eexcdc_dum,jmj
2792 : : real(dp) :: fact_re,fact_im,gx_re,gx_im,gy_re,gy_im,gz_re,gz_im,if3
2793 : : character(len=500) :: msg
2794 : : !arrays
2795 : 11 : integer,pointer :: my_atmtab(:),indlmn(:,:),indlmn_j(:,:)
2796 : 11 : logical,allocatable :: lmselect(:)
2797 : : real(dp) :: nhat_dum(1,1,1),trho_dum(1,1,1),kxc_dum(1,1,1),k3xc_dum(1,1,1)
2798 : 11 : real(dp),allocatable :: rho1(:,:,:),tau1(:,:,:),rhosph(:),vhartree(:),vxc(:,:,:)
2799 : 11 : real(dp),allocatable :: intf3(:,:),potsph(:),dVdr(:),func(:)
2800 : 11 : real(dp),pointer :: phi_j(:,:),soc_ij(:,:,:,:,:)
2801 : : type(pawrad_type),pointer :: pawrd
2802 : : type(pawtab_type),pointer :: pawtb
2803 : :
2804 : : ! ************************************************************************
2805 : : ABI_NVTX_START_RANGE(NVTX_PAWNABLA_SOC_INIT)
2806 : :
2807 : : !Some checks in case of core-valence (option_core=1)
2808 : 11 : if (option_core/=0.and.option_core/=1) then
2809 : 0 : msg='Wrong option_core value!'
2810 : 0 : ABI_BUG(msg)
2811 : : end if
2812 [ + + ]: 11 : if (option_core==1) then
2813 : : ! Check if we have the optional arguments
2814 [ - + ]: 6 : if (.not.present(atm)) then
2815 : 0 : msg='For core-valence calculation, need atm data!'
2816 : 0 : ABI_BUG(msg)
2817 : : end if
2818 [ + + ]: 12 : do itypat=1,ntypat
2819 : : ! Check if we have relativistic core wave functions
2820 [ - + ]: 12 : if (size(atm(itypat)%indlmn,1)<8) then
2821 : 0 : write(msg,'(a)') 'Wrong 1st dim. of indlmn_cor in pawnabla_soc_init (need spinors)!'
2822 : 0 : ABI_BUG(msg)
2823 : : end if
2824 : : enddo
2825 : : endif
2826 : :
2827 : : !Some useful variables
2828 : 11 : usekden=pawxc_get_usekden(ixc)
2829 : 11 : usecore=1 ; nkxc=0 ; usenhat=0
2830 : 11 : xclevel=pawxc_get_xclevel(ixc)
2831 [ + + ]: 11 : if (option_core==1) then
2832 : : lmn_size_cor=0
2833 [ + + ]: 12 : do itypat=1,ntypat
2834 : 12 : lmn_size_cor=max(lmn_size_cor,atm(itypat)%lmn_size) !Includes spinors
2835 : : enddo
2836 : : end if
2837 : :
2838 : : !Prepare output arrays
2839 [ - + ]: 11 : if (allocated(phisocphj)) then
2840 [ # # ]: 0 : do iatom=1,natom
2841 [ # # ]: 0 : if (allocated(phisocphj(iatom)%value)) then
2842 : 0 : ABI_FREE(phisocphj(iatom)%value)
2843 : : end if
2844 : : end do
2845 [ # # # # : 0 : ABI_FREE(phisocphj)
# # # # ]
2846 : : end if
2847 [ - + + - : 74 : ABI_MALLOC(phisocphj,(natom))
+ - + - -
+ - + - +
+ + ]
2848 [ + + ]: 52 : do iatom=1,natom
2849 : 41 : lmn_size=pawtab(typat(iatom))%lmn_size
2850 [ + + ]: 41 : if (option_core==0) then
2851 [ - + + - : 80 : ABI_MALLOC(phisocphj(iatom)%value,(2,2,3,lmn_size,lmn_size))
- + + - +
- + - - +
- + - + ]
2852 : : else
2853 : : !lmn_size_cor is double because it contains the spin component
2854 [ - + + - : 84 : ABI_MALLOC(phisocphj(iatom)%value,(2,2,3,lmn_size,lmn_size_cor))
- + + - +
- + - - +
- + - + ]
2855 : : end if
2856 [ + + + + : 107112 : phisocphj(iatom)%value=zero
+ + + + +
+ ]
2857 : : end do
2858 : :
2859 : : !Set up parallelism over atoms
2860 [ + - + + ]: 11 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2861 [ + + ]: 11 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2862 [ + - ]: 11 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2863 : : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
2864 : 11 : & my_natom_ref=my_natom)
2865 : :
2866 : : !-------------------------------------------------------------------
2867 : : !Loop over atoms
2868 : :
2869 [ + + ]: 36 : do iatom=1,my_natom
2870 : :
2871 : : ! Atom-dependent data
2872 : 25 : itypat=typat(iatom)
2873 : 25 : pawrd => pawrad(itypat)
2874 : 25 : pawtb => pawtab(itypat)
2875 : 25 : lmn_size=pawtb%lmn_size
2876 : 25 : lmn2_size=pawtb%lmn2_size
2877 : 25 : lm_size=pawtb%lcut_size**2
2878 : 25 : ln_size=pawtb%basis_size
2879 : 25 : mesh_size=pawtb%mesh_size
2880 : 25 : indlmn => pawtb%indlmn
2881 [ - + + - : 75 : ABI_MALLOC(lmselect,(lm_size))
+ - + - -
+ - + ]
2882 [ + + ]: 266 : lmselect(:)=.true.
2883 : :
2884 : : ! Distinguish valence-valence and core-valence cases
2885 [ + + ]: 25 : if (option_core==0) then
2886 : 12 : ln_size_j=ln_size
2887 : 12 : lmn_size_j=lmn_size
2888 : 12 : indlmn_j => pawtb%indlmn
2889 : 12 : phi_j => pawtb%phi
2890 : : else
2891 : 13 : ln_size_j=atm(itypat)%ln_size
2892 : 13 : lmn_size_j=atm(itypat)%lmn_size
2893 : 13 : indlmn_j => atm(itypat)%indlmn
2894 : 13 : phi_j => atm(itypat)%phi(:,:,1)
2895 [ + - ]: 13 : if (atm(itypat)%mesh_size<mesh_size) then
2896 : 0 : msg='mesh_size and mesh_size_cor not compatible!'
2897 : 0 : ABI_BUG(msg)
2898 : : end if
2899 : : endif
2900 : :
2901 : : ! Manage parallelism
2902 [ + + ]: 25 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
2903 : 25 : soc_ij => phisocphj(iatom_tot)%value(:,:,:,:,:)
2904 : :
2905 : : !-------------------------------------------------------------------
2906 : : !Compute all-electron density
2907 : :
2908 [ - + + - : 125 : ABI_MALLOC(rho1,(mesh_size,lm_size,nspden))
- + + - -
+ + - + -
+ - - + -
+ ]
2909 [ + + + + : 307265 : option=2 ; rho1=zero
+ + ]
2910 : : call pawdensities(compch_sph_dum,1,iatom_tot,lmselect,lmselect,lm_size,&
2911 : : & nhat_dum,nspden,-1,0,option,-1,0,pawang,0,pawrd,&
2912 : 25 : & pawrhoij(iatom),pawtb,rho1,trho_dum)
2913 [ - + ]: 25 : if (usekden==1) then
2914 [ # # # # : 0 : ABI_MALLOC(tau1,(mesh_size,lm_size,nspden))
# # # # #
# # # # #
# # # # ]
2915 [ # # # # : 0 : tau1=zero
# # ]
2916 : : call pawkindensities(1,lmselect,lm_size,nspden,-1,option,-1,&
2917 : 0 : & pawang,pawrd,pawrhoij(iatom),pawtb,tau1,trho_dum)
2918 : : end if
2919 : :
2920 : : !-------------------------------------------------------------------
2921 : : !Compute spherical potential and compute its first derivative dV/dr
2922 : :
2923 [ - + + - : 75 : ABI_MALLOC(potsph,(mesh_size))
+ - + - -
+ - + -
+ ]
2924 [ - + + - : 50 : ABI_MALLOC(dVdr,(mesh_size))
- + - + -
+ ]
2925 [ + + + + ]: 63157 : potsph=zero ; dVdr=zero
2926 : :
2927 : : ! Compute XC potential
2928 : 25 : option=1
2929 [ + - ]: 25 : if (pawxcdev/=0) then
2930 [ - + + - : 100 : ABI_MALLOC(vxc,(mesh_size,lm_size,nspden))
- + + - -
+ + - - +
- + ]
2931 [ + + + + : 307265 : vxc=zero
+ + ]
2932 : : call pawxcm(pawtb%coredens,eexc_dum,eexcdc_dum,ssxc_dum,idum,hyb_mixing_,ixc,kxc_dum,lm_size,&
2933 : : & lmselect,nhat_dum,nkxc,.false.,mesh_size,nspden,option,pawang,pawrd,&
2934 : 25 : & pawxcdev,rho1,usecore,usenhat,vxc,xclevel,xc_denpos,el_temp)
2935 [ + + ]: 31591 : potsph(1:mesh_size)=half*(vxc(1:mesh_size,1,1)+vxc(1:mesh_size,1,nspden))
2936 : : else
2937 [ # # # # : 0 : ABI_MALLOC(vxc,(mesh_size,pawang%angl_size,nspden))
# # # # #
# # # # #
# # # # #
# ]
2938 [ # # # # : 0 : vxc=zero
# # ]
2939 : : call pawxc(pawtb%coredens,eexc_dum,eexcdc_dum,ssxc_dum,hyb_mixing_,ixc,kxc_dum,k3xc_dum,&
2940 : : & lm_size,lmselect,nhat_dum,nkxc,nkxc,.false.,mesh_size,nspden,option,pawang,&
2941 : : & pawrd,rho1,usecore,usenhat,vxc,xclevel,xc_denpos,el_temp,&
2942 : 0 : & coretau=pawtb%coretau,taur=tau1,xc_taupos=xc_taupos)
2943 [ # # ]: 0 : potsph(1:mesh_size)=zero
2944 [ # # ]: 0 : do ipts=1,pawang%angl_size
2945 : : potsph(1:mesh_size)=potsph(1:mesh_size) &
2946 [ # # ]: 0 : & +half*(vxc(1:mesh_size,ipts,1)+vxc(1:mesh_size,ipts,nspden))*pawang%angwgth(ipts)
2947 : : end do
2948 [ # # ]: 0 : potsph(1:mesh_size)=sqrt(four_pi)*potsph(1:mesh_size)
2949 : : end if
2950 : :
2951 : : ! Compute Hartree potentialHalfFineStruct2
2952 [ - + + - : 50 : ABI_MALLOC(vhartree,(mesh_size))
- + - + -
+ ]
2953 [ - + + - : 50 : ABI_MALLOC(rhosph,(mesh_size))
- + - + -
+ ]
2954 [ + + + + ]: 63157 : vhartree=zero ; rhosph=zero
2955 [ + + ]: 31591 : rhosph(1:mesh_size)=rho1(1:mesh_size,1,1)
2956 [ + + ]: 31591 : if (usecore==1) rhosph(1:mesh_size)=rhosph(1:mesh_size)+sqrt(four_pi)*pawtb%coredens(1:mesh_size)
2957 [ + + ]: 31591 : rhosph(1:mesh_size)=rhosph(1:mesh_size)*four_pi*pawrd%rad(1:mesh_size)**2
2958 : 25 : call poisson(rhosph,0,pawrd,vhartree)
2959 [ + + ]: 31566 : vhartree(2:mesh_size)=(vhartree(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrd%rad(2:mesh_size)
2960 : 25 : call pawrad_deducer0(vhartree,mesh_size,pawrd)
2961 [ + + ]: 31591 : potsph(1:mesh_size)=potsph(1:mesh_size)+vhartree(1:mesh_size)
2962 : :
2963 : : ! Apply angular scaling factor
2964 [ + + ]: 31591 : potsph(1:mesh_size)=one_over_fourpi*potsph(1:mesh_size)
2965 : :
2966 : : ! Compute 1st derivative of potential
2967 : 25 : call nderiv_gen(dVdr,potsph,pawrd)
2968 : :
2969 : : ! Multiply by relativistic factor
2970 [ + + ]: 31591 : dVdr(1:mesh_size)=dVdr(1:mesh_size)*(one/(one-potsph(1:mesh_size)*half/InvFineStruct**2)**2)
2971 : :
2972 [ - + ]: 25 : ABI_FREE(vxc)
2973 [ - + ]: 25 : ABI_FREE(vhartree)
2974 [ - + ]: 25 : ABI_FREE(potsph)
2975 [ - + ]: 25 : ABI_FREE(rhosph)
2976 [ - + ]: 25 : ABI_FREE(lmselect)
2977 [ - + ]: 25 : ABI_FREE(rho1)
2978 [ - + ]: 25 : if (usekden==1) then
2979 [ # # ]: 0 : ABI_FREE(tau1)
2980 : : end if
2981 : :
2982 : : !-------------------------------------------------------------------
2983 : : !Compute radial and angular contributions
2984 : :
2985 [ - + + - : 100 : ABI_MALLOC(intf3,(ln_size,ln_size_j))
- + + - +
- + - - +
- + ]
2986 [ + + + + ]: 533 : intf3=zero
2987 : :
2988 : : ! >>>> Calculate f_3= alpha^2/4 int[dr ui*(r) uj(r) dV/dr]
2989 [ - + + - : 50 : ABI_MALLOC(func,(mesh_size))
- + - + -
+ ]
2990 [ + + ]: 125 : do jln=1,ln_size_j
2991 [ + + ]: 533 : do iln=1,ln_size
2992 [ + + ]: 516904 : func(1:mesh_size)=dVdr(1:mesh_size)*pawtb%phi(1:mesh_size,iln)*phi_j(1:mesh_size,jln)
2993 : 508 : call simp_gen(intf3(iln,jln),func,pawrd)
2994 : : end do
2995 : : end do
2996 [ + + + + ]: 533 : intf3(:,:)=QuarterFineStruct2*spnorbscl*intf3(:,:)
2997 [ - + ]: 25 : ABI_FREE(func)
2998 : :
2999 : : ! Loop over initial states (valence or core, according to option_core)
3000 [ + + ]: 381 : do jlmn=1,lmn_size_j
3001 : 356 : jl=indlmn_j(1,jlmn)
3002 : 356 : jm=indlmn_j(2,jlmn)
3003 : 356 : jlm=indlmn_j(4,jlmn)
3004 : 356 : jln=indlmn_j(5,jlmn)
3005 : :
3006 : : ! In case of spinorial core wave function, we have to handle imaginary
3007 : : ! spherical harmonics as linear combination of real spherical harmonics
3008 : : ! See Brouwer et al, CPC 266, 108029 (2021), equation (21)
3009 [ + + ]: 356 : if (option_core==1) then
3010 : 260 : jm_re= abs(jm)
3011 : 260 : jm_im=-abs(jm)
3012 : 260 : jlm_re=jl*(jl+1)+jm_re+1
3013 : 260 : jlm_im=jl*(jl+1)+jm_im+1
3014 : : ! Calculate spinor dependent coefficients
3015 : 260 : sgnkappa=indlmn_j(3,jlmn) !sign of kappa
3016 : 260 : jmj=half*indlmn_j(8,jlmn) !2mj is stored in indlmn_cor
3017 : 260 : js=indlmn_j(6,jlmn) !1 is up, 2 is down
3018 [ + + ]: 260 : if (sgnkappa==1) then
3019 [ + + ]: 52 : if(js==1) then
3020 : 26 : cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
3021 : : else
3022 : 26 : cgc=-sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
3023 : : endif
3024 : : else
3025 [ + + ]: 208 : if(js==1) then
3026 : 65 : cgc= sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
3027 : : else
3028 : 143 : cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
3029 : : endif
3030 : : endif
3031 : :
3032 : : ! Calculate factors to convert from complex to real sph. harm.
3033 [ + + ]: 260 : if (jm<0) then
3034 : 117 : fact_re=sqr_fourpi_over_3*half_sqrt2*cgc
3035 : 117 : fact_im=-fact_re
3036 [ + + ]: 143 : else if (jm>0) then
3037 : 39 : fact_re=sqr_fourpi_over_3*half_sqrt2*cgc*(-1)**jm
3038 : 39 : fact_im=fact_re
3039 : : else
3040 : : fact_re=sqr_fourpi_over_3
3041 : : fact_im=0
3042 : : end if
3043 : : else ! valence-valence case (real)
3044 : : js=1
3045 : : jlm_re=jlm ; jlm_im=jlm
3046 : : fact_re=sqr_fourpi_over_3
3047 : : fact_im=0
3048 : : end if
3049 : :
3050 : : ! Loop over final states
3051 [ + + ]: 3429 : do ilmn=1,lmn_size
3052 : 3048 : ilm=indlmn(4,ilmn)
3053 : 3048 : iln=indlmn(5,ilmn)
3054 : :
3055 : : ! >>>> Calculate g_ij=(g_x,g_y,g_z) = sqrt(4pi/3) int dOmega Ylm Ylm' S1-1,0,1
3056 : : ! using real Gaunt coefficients
3057 : 3048 : gx_re=zero;gy_re=zero;gz_re=zero
3058 : 3048 : gx_im=zero;gy_im=zero;gz_im=zero
3059 : 3048 : if3=zero
3060 : :
3061 : : ! jl was set as a flag for invalid combinations
3062 : : ! i.e. m=-(l+1) or m=(l+1)
3063 : : ! In these cases, cgc=0 ; so gx=gy=gz=0
3064 [ + + ]: 3048 : if (jl/=-1) then
3065 : 2364 : if3=intf3(iln,jln)
3066 [ + + ]: 2364 : klm_re=merge((jlm_re*(jlm_re-1))/2+ilm,(ilm*(ilm-1))/2+jlm_re,ilm<=jlm_re)
3067 : :
3068 : : ! Real parts
3069 : : ! M=-1
3070 : 2364 : ignt=pawang%gntselect(2,klm_re) !get index for L=1 M =-1 ilm jlm_re
3071 [ + + ]: 2364 : if (ignt/=0) gy_re=fact_re*pawang%realgnt(ignt)
3072 : : ! M=0
3073 : 2364 : ignt=pawang%gntselect(3,klm_re) !get index for L=1 M = 0 ilm jlm_re
3074 [ + + ]: 2364 : if (ignt/=0) gz_re=fact_re*pawang%realgnt(ignt)
3075 : : ! M=1
3076 : 2364 : ignt=pawang%gntselect(4,klm_re) !get index for L=1 M = 1 ilm jlm_re
3077 [ + + ]: 2364 : if (ignt/=0) gx_re=fact_re*pawang%realgnt(ignt)
3078 : :
3079 : : ! Imaginary parts
3080 [ + + ]: 2364 : if (option_core==1) then
3081 [ + + ]: 1596 : klm_im=merge((jlm_im*(jlm_im-1))/2+ilm,(ilm*(ilm-1))/2+jlm_im,ilm<=jlm_im)
3082 : : ! M=-1
3083 : 1596 : ignt=pawang%gntselect(2,klm_im) !get index for L=1 M =-1 ilm jlm_im
3084 [ + + ]: 1596 : if (ignt/=0) gy_im=fact_im*pawang%realgnt(ignt)
3085 : : ! M=0
3086 : 1596 : ignt=pawang%gntselect(3,klm_im) !get index for L=1 M = 0 ilm jlm_im
3087 [ + + ]: 1596 : if (ignt/=0) gz_im=fact_im*pawang%realgnt(ignt)
3088 : : ! M=1
3089 : 1596 : ignt=pawang%gntselect(4,klm_im) !get index for L=1 M = 1 ilm jlm_im
3090 [ + + ]: 1596 : if (ignt/=0) gx_im=fact_im*pawang%realgnt(ignt)
3091 : : end if
3092 : : end if
3093 : :
3094 : : ! >>>> Calculate Sigma X g_ij
3095 : :
3096 [ + + + + ]: 3404 : if (option_core==0.or.js==1) then
3097 : : !(Sigma^up-up X gij)_x = -gy*f_3
3098 : 1566 : soc_ij(1,1,1,ilmn,jlmn)=-if3*gy_re ! real part
3099 : 1566 : soc_ij(2,1,1,ilmn,jlmn)=-if3*gy_im ! imag part
3100 : : !(Sigma^up-up X gij)_y = gx*f_3
3101 : 1566 : soc_ij(1,1,2,ilmn,jlmn)= if3*gx_re ! real part
3102 : 1566 : soc_ij(2,1,2,ilmn,jlmn)= if3*gx_im ! imag part
3103 : : !(Sigma^up-up X gij)_z = 0
3104 : 1566 : soc_ij(1,1,3,ilmn,jlmn)= zero ! real part
3105 : 1566 : soc_ij(2,1,3,ilmn,jlmn)= zero ! imag part
3106 : :
3107 : : !(Sigma^dn-up X gij)_x = i.gz*f_3
3108 : 1566 : soc_ij(1,2,1,ilmn,jlmn)=-if3*gz_im ! real part
3109 : 1566 : soc_ij(2,2,1,ilmn,jlmn)= if3*gz_re ! imag part
3110 : : !(Sigma^dn-up X gij)_y = -gz*f_3
3111 : 1566 : soc_ij(1,2,2,ilmn,jlmn)=-if3*gz_re ! real part
3112 : 1566 : soc_ij(2,2,2,ilmn,jlmn)=-if3*gz_im ! imag part
3113 : : !(Sigma^dn-up X gij)_z = (gy-i.gx)*f_3
3114 : 1566 : soc_ij(1,2,3,ilmn,jlmn)= if3*(gy_re+gx_im) ! real part
3115 : 1566 : soc_ij(2,2,3,ilmn,jlmn)= if3*(gy_im-gx_re) ! imag part
3116 : :
3117 [ + - + + ]: 1482 : else if (option_core==1.and.js==2) then
3118 : : !(Sigma^up-dn X gij^dn)_x = -i.gz^dn*f_3
3119 : 798 : soc_ij(1,1,1,ilmn,jlmn)= if3*gz_im ! real part
3120 : 798 : soc_ij(2,1,1,ilmn,jlmn)=-if3*gz_re ! imag part
3121 : : !(Sigma^up-dn X gij^dn)_y = -gz^dn*f_3
3122 : 798 : soc_ij(1,1,2,ilmn,jlmn)=-if3*gz_re ! real part
3123 : 798 : soc_ij(2,1,2,ilmn,jlmn)=-if3*gz_im ! imag part
3124 : : !(Sigma^up-dn X gij^dn)_z = (gy^dn+i.gx^dn)*f_3
3125 : 798 : soc_ij(1,1,3,ilmn,jlmn)= if3*(gy_re-gx_im) ! real part
3126 : 798 : soc_ij(2,1,3,ilmn,jlmn)= if3*(gy_im+gx_re) ! imag part
3127 : :
3128 : : !(Sigma^dn-dn X gij^dn)_x = gy^dn*f_3
3129 : 798 : soc_ij(1,2,1,ilmn,jlmn)= if3*gy_re ! real part
3130 : 798 : soc_ij(2,2,1,ilmn,jlmn)= if3*gy_im ! imag part
3131 : : !(Sigma^dn-dn X gij^dn)_y = -gx^dn*f_3
3132 : 798 : soc_ij(1,2,2,ilmn,jlmn)=-if3*gx_re ! real part
3133 : 798 : soc_ij(2,2,2,ilmn,jlmn)=-if3*gx_im ! imag part
3134 : : !(Sigma^dn-dn X gij^dn)_z = 0
3135 : 798 : soc_ij(1,2,3,ilmn,jlmn)= zero ! real part
3136 : 798 : soc_ij(2,2,3,ilmn,jlmn)= zero ! imag part
3137 : : end if
3138 : :
3139 : : end do ! ilmn
3140 : : end do ! jlmn
3141 : :
3142 : : ! Symetrization
3143 [ + + + - ]: 25 : if (option_core==0.and.lmn_size>1) then
3144 [ + + ]: 96 : do jlmn=2,lmn_size
3145 [ + + ]: 432 : do ilmn=1,jlmn-1
3146 [ + + ]: 1428 : do ii=1,3
3147 [ + + ]: 3360 : do jj=1,2
3148 : 2016 : avg=half*(soc_ij(1,jj,ii,ilmn,jlmn)+soc_ij(1,jj,ii,jlmn,ilmn))
3149 : 2016 : soc_ij(1,jj,ii,ilmn,jlmn)=avg ; soc_ij(1,jj,ii,jlmn,ilmn)=avg
3150 : 2016 : avg=half*(soc_ij(2,jj,ii,ilmn,jlmn)+soc_ij(2,jj,ii,ilmn,jlmn))
3151 : 3024 : soc_ij(2,jj,ii,ilmn,jlmn)=avg ; soc_ij(2,jj,ii,jlmn,ilmn)=avg
3152 : : end do
3153 : : end do
3154 : : end do
3155 : : end do
3156 : : end if
3157 : :
3158 [ - + ]: 25 : ABI_FREE(dVdr)
3159 [ - + ]: 36 : ABI_FREE(intf3)
3160 : :
3161 : : end do ! iatom
3162 : :
3163 : : !Reduction in case of parallelism
3164 [ + + ]: 11 : if (paral_atom) then
3165 : 8 : call xmpi_sum(phisocphj,my_comm_atom,ierr)
3166 : : end if
3167 : :
3168 : : !Destroy atom table(s) used for parallelism
3169 : 11 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
3170 : : ABI_NVTX_END_RANGE()
3171 : :
3172 [ - + - + : 22 : end subroutine pawnabla_soc_init
- + - + -
+ - + -
+ ]
3173 : : !!***
3174 : :
3175 : : !----------------------------------------------------------------------
3176 : :
3177 : 0 : END MODULE m_paw_optics
3178 : : !!***
|
