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
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_paral_atom, only : get_my_atmtab,free_my_atmtab
63 : : #if defined(HAVE_GPU_MARKERS)
64 : : use m_nvtx_data
65 : : #endif
66 : :
67 : : implicit none
68 : :
69 : : private
70 : :
71 : : !public procedures.
72 : : public :: optics_paw
73 : : public :: optics_paw_core
74 : : public :: linear_optics_paw
75 : :
76 : : !I/O parameters
77 : : !Set to true to force the use of netCDF when available
78 : : ! overriding the value of dtset%iomode
79 : : logical,parameter :: use_netcdf_forced=.true.
80 : : !Set to true to use netcdf-MPIIO when available
81 : : logical,parameter :: use_netcdf_mpiio=.true.
82 : : !Set to true to compute/write only half of the (n,m) dipoles
83 : : logical,parameter :: compute_half_dipoles=.true.
84 : : !Set to true to use unlimited dimensions in netCDF file (experimental)
85 : : ! This is not mandatory because we know exactly the amount of data to write
86 : : ! and seems to impact performances negatively...
87 : : ! Not compatible with compute_half_dipoles=.true.
88 : : logical,parameter :: use_netcdf_unlimited=.false.
89 : :
90 : : CONTAINS !========================================================================================
91 : : !!***
92 : :
93 : : !----------------------------------------------------------------------
94 : :
95 : : !!****f* m_paw_optics/optics_paw
96 : : !! NAME
97 : : !! optics_paw
98 : : !!
99 : : !! FUNCTION
100 : : !! Compute matrix elements need for optical conductivity (in the PAW context) and store them in a file
101 : : !! Matrix elements = <Phi_i|Nabla|Phi_j>
102 : : !!
103 : : !! INPUTS
104 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
105 : : !! cg(2,mcg)=planewave coefficients of wavefunctions.
106 : : !! cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk>
107 : : !! and each |p_lmn> non-local projector
108 : : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
109 : : !! dtfil <type(datafiles_type)>=variables related to files
110 : : !! dtset <type(dataset_type)>=all input variables for this dataset
111 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
112 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
113 : : !! mband=maximum number of bands
114 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
115 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
116 : : !! mkmem =number of k points treated by this node.
117 : : !! mpi_enreg=information about MPI parallelization
118 : : !! mpsang =1+maximum angular momentum for nonlocal pseudopotentials
119 : : !! mpw=maximum dimensioned size of npw.
120 : : !! natom=number of atoms in cell.
121 : : !! nkpt=number of k points.
122 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
123 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
124 : : !! pawang <type(pawang_type)>= PAW ANGular mesh discretization and related data
125 : : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
126 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
127 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
128 : : !! znucl(ntypat)=atomic number of atom type
129 : : !!
130 : : !! OUTPUT
131 : : !! psinablapsi_out=contains the matrix elements
132 : : !! (The size of the psinablapsi_out decide whether we compute the full matrix or
133 : : !! only the diagonal part)
134 : : !! (if not present only writing in a file)
135 : : !!
136 : : !! SIDE EFFECTS
137 : : !!
138 : : !! NOTES
139 : : !!
140 : : !! SOURCE
141 : :
142 : 2294 : subroutine optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen0,gprimd,hdr,kg,&
143 : 1147 : & mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,&
144 [ + + + - : 1147 : & pawang,pawrad,pawrhoij,pawtab,znucl,psinablapsi_out)
- + ]
145 : :
146 : : !Arguments ------------------------------------
147 : : !scalars
148 : : integer,intent(in) :: mband,mcg,mcprj,mkmem,mpsang,mpw,natom,nkpt,nsppol
149 : : type(MPI_type),intent(in) :: mpi_enreg
150 : : type(datafiles_type),intent(in) :: dtfil
151 : : type(dataset_type),intent(in) :: dtset
152 : : type(hdr_type),intent(inout) :: hdr
153 : : type(pawang_type),intent(in) :: pawang
154 : : real(dp),optional,target,intent(out) :: psinablapsi_out(:,:,:,:)
155 : : !arrays
156 : : integer,intent(in) :: atindx1(natom),dimcprj(natom),npwarr(nkpt)
157 : : integer,intent(in),target :: kg(3,mpw*mkmem)
158 : : real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
159 : : real(dp),intent(in) :: gprimd(3,3),znucl(dtset%ntypat)
160 : : real(dp),intent(inout) :: cg(2,mcg)
161 : : type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
162 : : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
163 : : type(pawrhoij_type),intent(inout) :: pawrhoij(mpi_enreg%my_natom)
164 : : type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
165 : :
166 : : !Local variables-------------------------------
167 : : !scalars
168 : : integer,parameter :: master=0
169 : : integer :: bsize,iomode,bdtot_index,cplex,etiq,fformopt,iatom,ib,ibmax,ibmin,ibg,ibsp
170 : : integer :: ibshift,icg,ierr,ikg,ikpt,ilmn,ount,ncid,varid,idir
171 : : integer :: iorder_cprj,ipw,ispinor,isppol,istwf_k,itypat,iwavef
172 : : integer :: jb,jbshift,jbsp,my_jb,jlmn,jwavef,lmn_size,mband_cprj,option_core
173 : : integer :: my_nspinor,nband_k,nband_cprj_k,npw_k,sender,me,master_spfftband,pnp_size
174 : : integer :: spaceComm_band,spaceComm_bandspinorfft,spaceComm_fft,spaceComm_kpt
175 : : integer :: spaceComm_spinor,spaceComm_bandspinor,spaceComm_spinorfft,spaceComm_w
176 : : integer, parameter :: NO_FILE_OUT=-1
177 : : logical :: already_has_nabla,cprj_paral_band,myband,mykpt,iomode_etsf_mpiio
178 : : logical :: i_am_master,i_am_master_kpt,i_am_master_band,i_am_master_spfft,nc_unlimited,store_half_dipoles
179 : : logical :: diag_only
180 : : real(dp) :: cgnm1,cgnm2,cpnm1,cpnm2,cpnm11,cpnm22,cpnm12,cpnm21,cpnm_11m22,cpnm_21p12,cpnm_21m12,el_temp
181 : : character(len=500) :: msg
182 : : type(nctkdim_t) :: nctkdim
183 : : !arrays
184 : : 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)
185 : 1147 : integer, contiguous, pointer :: kg_k(:,:)
186 : : real(dp) :: kpoint(3),tsec(2),nabla_ij(3)
187 : 1147 : real(dp),allocatable :: kpg_k(:,:)
188 : 1147 : real(dp),pointer :: psinablapsi(:,:,:),psinablapsi_paw(:,:,:),psinablapsi_soc(:,:,:)
189 : 1147 : real(dp),pointer :: soc_ij(:,:,:)
190 : 1147 : type(coeff5_type),allocatable,target :: phisocphj(:)
191 : 1147 : type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
192 : : type(nctkarr_t) :: nctk_arrays(1)
193 : :
194 : : ! ************************************************************************
195 : :
196 : : DBG_ENTER("COLL")
197 : : ABI_NVTX_START_RANGE(NVTX_OPTICS_PAW)
198 : :
199 : : !Compatibility tests
200 : 0 : ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
201 [ - + - - ]: 1147 : ABI_CHECK(mpi_enreg%paral_spinor==0.or.dtset%pawspnorb==0,"spinor parallelization not supported with SOC!")
202 : : ! MJV 6/12/2008: looks like mpi_enreg may not be completely initialized here
203 : : if (xmpi_paral==1) then
204 [ + + ]: 4588 : tmp_shape = shape(mpi_enreg%proc_distrb)
205 [ - + ]: 1147 : if (nkpt > tmp_shape(1)) then
206 : 0 : ABI_BUG('problem with proc_distrb!')
207 : : end if
208 : : end if
209 : :
210 : : !Init parallelism
211 : 1147 : spaceComm_w=mpi_enreg%comm_cell
212 [ + + ]: 1147 : if (mpi_enreg%paral_kgb==1) then
213 : 836 : spaceComm_kpt=mpi_enreg%comm_kpt
214 : 836 : spaceComm_fft=mpi_enreg%comm_fft
215 : 836 : spaceComm_band=mpi_enreg%comm_band
216 : 836 : spaceComm_spinor=mpi_enreg%comm_spinor
217 : 836 : spaceComm_bandspinor=mpi_enreg%comm_bandspinor
218 : 836 : spaceComm_spinorfft=mpi_enreg%comm_spinorfft
219 : 836 : spaceComm_bandspinorfft=mpi_enreg%comm_bandspinorfft
220 : : else
221 : 311 : spaceComm_kpt=mpi_enreg%comm_kpt
222 : 311 : spaceComm_fft=xmpi_comm_self
223 : 311 : spaceComm_band=mpi_enreg%comm_band
224 : 311 : spaceComm_spinor=xmpi_comm_self
225 : 311 : spaceComm_bandspinor=spaceComm_band
226 : 311 : spaceComm_spinorfft=xmpi_comm_self
227 : 311 : spaceComm_bandspinorfft=xmpi_comm_self
228 : : end if
229 : 1147 : me=xmpi_comm_rank(spaceComm_w)
230 : 1147 : i_am_master=(me==master)
231 : 1147 : i_am_master_kpt=(xmpi_comm_rank(spaceComm_kpt)==master)
232 : 1147 : i_am_master_band=(xmpi_comm_rank(spaceComm_band)==master)
233 : 1147 : i_am_master_spfft=(xmpi_comm_rank(spaceComm_spinorfft)==master)
234 : 1147 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
235 : :
236 : : !Check wether we write in file or save the matrix elements in psinablapsi_out
237 : : !and if we need to compute the full matrix or only the diagonal part
238 : 1147 : iomode=dtset%iomode
239 : 1147 : diag_only = .false.
240 : 1147 : store_half_dipoles = .false.
241 [ + + ]: 1147 : if (present(psinablapsi_out)) then
242 : 1121 : iomode = NO_FILE_OUT
243 : 1121 : iomode_etsf_mpiio = .false.
244 : 1121 : nc_unlimited = .false.
245 : 1121 : bsize = size(psinablapsi_out(1,1,:,1))
246 [ - + ]: 1121 : if (bsize == mband) then
247 : : diag_only = .true.
248 : : store_half_dipoles = .true. !FB ??
249 [ # # ]: 0 : else if (bsize == mband*(mband+1)/2) then
250 : : diag_only = .false.
251 : : store_half_dipoles = .true.
252 [ # # ]: 0 : else if (bsize == mband**2) then
253 : : diag_only = .false.
254 : : store_half_dipoles = .false.
255 : : else
256 : 0 : msg = "Wrong dimensions of psinablapsi_out!"
257 : 0 : ABI_ERROR(msg)
258 : : end if
259 [ + + + + : 2906753 : psinablapsi_out = zero
+ + + + ]
260 : : end if
261 : :
262 : : !----------------------------------------------------------------------------------
263 : : !1- Opening of OPT file and header writing
264 : : !----------------------------------------------------------------------------------
265 : :
266 [ + + ]: 1147 : if (iomode /= NO_FILE_OUT) then
267 : : ! I/O mode is netCDF or Fortran
268 : : iomode=merge(IO_MODE_ETSF,IO_MODE_FORTRAN_MASTER,dtset%iomode==IO_MODE_ETSF)
269 : 26 : if (use_netcdf_forced) iomode=IO_MODE_ETSF
270 : :
271 : : ! (master proc only)
272 [ + + ]: 26 : if (i_am_master) then
273 : 8 : fformopt=610 ; if (compute_half_dipoles) fformopt=620
274 : : ! ====> NETCDF format
275 : : if (iomode==IO_MODE_ETSF) then
276 : : ! Open/create nc file
277 [ - + ]: 8 : NCF_CHECK(nctk_open_create(ncid,nctk_ncify(dtfil%fnameabo_app_opt),xmpi_comm_self))
278 : : ! Write header data
279 [ - + ]: 8 : NCF_CHECK(hdr%ncwrite(ncid,fformopt,nc_define=.true.))
280 : : ! Define dims and array for dipole variables
281 : 8 : nctk_arrays(1)%name="dipole_valence_valence"
282 : 8 : nctk_arrays(1)%dtype="dp"
283 : 8 : nc_unlimited=(use_netcdf_unlimited.and.(.not.(nctk_has_mpiio.and.use_netcdf_mpiio)))
284 : : if (nc_unlimited) then
285 : : nctkdim%name="unlimited_bands"
286 : : nctkdim%value=NF90_UNLIMITED
287 : : NCF_CHECK(nctk_def_dims(ncid,nctkdim))
288 : : nctk_arrays(1)%shape_str=&
289 : : & "complex,number_of_cartesian_directions,max_number_of_states,number_of_kpoints,number_of_spins,unlimited_bands"
290 : : else if (compute_half_dipoles) then
291 : 8 : nctkdim%name="max_number_of_state_pairs"
292 : 8 : nctkdim%value=(mband*(mband+1))/2
293 [ - + ]: 8 : NCF_CHECK(nctk_def_dims(ncid,nctkdim))
294 : : nctk_arrays(1)%shape_str=&
295 : 8 : & "complex,number_of_cartesian_directions,max_number_of_state_pairs,number_of_kpoints,number_of_spins"
296 : : else
297 : : nctk_arrays(1)%shape_str=&
298 : : & "complex,number_of_cartesian_directions,max_number_of_states,max_number_of_states,number_of_kpoints,number_of_spins"
299 : : end if
300 [ - + ]: 8 : NCF_CHECK(nctk_def_arrays(ncid, nctk_arrays))
301 [ - + ]: 8 : NCF_CHECK(nctk_set_atomic_units(ncid, "dipole_valence_valence"))
302 : : ! Write eigenvalues
303 [ - + ]: 8 : NCF_CHECK(nctk_set_datamode(ncid))
304 : 8 : varid=nctk_idname(ncid,"eigenvalues")
305 [ + + - + : 32 : NCF_CHECK(nf90_put_var(ncid,varid,reshape(eigen0,[mband,nkpt,nsppol])))
- - ]
306 : : !Close file here because the rest has possibly to be written with collective I/O
307 [ - + ]: 8 : NCF_CHECK(nf90_close(ncid))
308 : : ! ====> Standard FORTRAN binary format
309 : : else if (iomode==IO_MODE_FORTRAN_MASTER) then
310 : : if (open_file(dtfil%fnameabo_app_opt,msg,newunit=ount,form="unformatted",status="unknown")/= 0) then
311 : : ABI_ERROR(msg)
312 : : end if
313 : : call hdr%fort_write(ount,fformopt,ierr,rewind=.true.)
314 : : write(ount)(eigen0(ib),ib=1,mband*nkpt*nsppol)
315 : : else
316 : : msg = "Wrong OPT file format!"
317 : : ABI_BUG(msg)
318 : : end if ! File format
319 : : end if ! master node
320 : 26 : call xmpi_bcast(iomode,master,spaceComm_w,ierr) ! Seems mandatory; why ?
321 [ + - - + ]: 26 : iomode_etsf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
322 : : nc_unlimited=(iomode==IO_MODE_ETSF.and.use_netcdf_unlimited.and.(.not.iomode_etsf_mpiio)) ! UNLIMITED not compatible with mpi-io
323 : : store_half_dipoles=(compute_half_dipoles.and.(.not.nc_unlimited))
324 : : end if
325 : :
326 : : !----------------------------------------------------------------------------------
327 : : !2- Computation of on-site contribution: <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>
328 : : !----------------------------------------------------------------------------------
329 : :
330 [ + - - + ]: 1147 : already_has_nabla=all(pawtab(:)%has_nabla==2)
331 : 1147 : call pawnabla_init(mpsang,dtset%ntypat,pawrad,pawtab)
332 : :
333 : : !Get electronic temperature from dtset
334 [ - + - - : 1147 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
- - ]
335 : :
336 : : !Compute spin-orbit contributions if necessary
337 [ + + ]: 1147 : if (dtset%pawspnorb==1) then
338 : 1 : option_core=0
339 : : call pawnabla_soc_init(el_temp,phisocphj,option_core,dtset%ixc,mpi_enreg%my_natom,natom,&
340 : : & dtset%nspden,dtset%ntypat,pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,&
341 : : & dtset%spnorbscl,dtset%typat,dtset%xc_denpos,dtset%xc_taupos,znucl,&
342 : 1 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
343 : : end if
344 : :
345 : : !----------------------------------------------------------------------------------
346 : : !3- Computation of <psi_n|-i.nabla|psi_m> for each k
347 : : !----------------------------------------------------------------------------------
348 : :
349 : : !Prepare valence-valence dipoles writing
350 : : !In case of netCDF access to OPT file, prepare collective I/O
351 [ + + ]: 1147 : if (iomode == IO_MODE_ETSF) then
352 [ + - ]: 26 : if (iomode_etsf_mpiio) then
353 [ + + ]: 26 : if (i_am_master_spfft) then
354 [ - + ]: 16 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt),spaceComm_band))
355 : 16 : varid=nctk_idname(ncid,"dipole_valence_valence")
356 [ + + ]: 16 : if (xmpi_comm_size(spaceComm_w)>1) then
357 [ - + ]: 14 : NCF_CHECK(nctk_set_collective(ncid,varid))
358 : : end if
359 [ - + ]: 16 : NCF_CHECK(nctk_set_datamode(ncid))
360 : : end if
361 [ # # ]: 0 : else if (i_am_master) then
362 [ # # ]: 0 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt),xmpi_comm_self))
363 : 0 : varid=nctk_idname(ncid,"dipole_valence_valence")
364 : : if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
365 : : NCF_CHECK(nctk_set_collective(ncid,varid))
366 : : end if
367 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
368 : : end if
369 : : end if
370 [ + + ]: 1147 : if (iomode /= NO_FILE_OUT) then
371 [ + - ]: 26 : if (iomode_etsf_mpiio) then
372 : : !If MPI-IO, store only ib elements for each jb
373 [ - + + - : 78 : ABI_MALLOC(psinablapsi,(2,3,mband))
+ - + - -
+ - + ]
374 [ - + + - : 52 : ABI_MALLOC(psinablapsi_paw,(2,3,mband))
- + - + ]
375 [ + + ]: 26 : if (dtset%pawspnorb==1) then
376 [ - + + - : 2 : ABI_MALLOC(psinablapsi_soc,(2,3,mband))
- + - + ]
377 : : end if
378 : : else
379 : : !If not, store all (ib,jb) pairs (or half)
380 [ # # ]: 0 : bsize=mband**2 ; if (store_half_dipoles) bsize=(mband*(mband+1))/2
381 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,bsize))
# # # # #
# # # ]
382 [ # # # # : 0 : ABI_MALLOC(psinablapsi_paw,(2,3,bsize))
# # # # ]
383 [ # # ]: 0 : if (dtset%pawspnorb==1) then
384 [ # # # # : 0 : ABI_MALLOC(psinablapsi_soc,(2,3,bsize))
# # # # ]
385 : : end if
386 [ # # # # : 0 : psinablapsi=zero
# # ]
387 : : end if
388 [ + + ]: 104 : pnp_size=size(psinablapsi)
389 : : else
390 [ - + ]: 1121 : if (diag_only) then
391 : : bsize = mband
392 [ # # ]: 0 : else if (store_half_dipoles) then
393 : 0 : bsize=(mband*(mband+1))/2
394 : : else
395 : 0 : bsize=mband**2
396 : : end if
397 [ - + + - : 3363 : ABI_MALLOC(psinablapsi_paw,(2,3,bsize))
+ - + - -
+ - + ]
398 [ - + ]: 1121 : if (dtset%pawspnorb==1) then
399 [ # # # # : 0 : ABI_MALLOC(psinablapsi_soc,(2,3,bsize))
# # # # ]
400 : : end if
401 [ + + ]: 4484 : pnp_size=size(psinablapsi_paw)
402 : : end if
403 : :
404 : : !Determine if cprj datastructure is distributed over bands
405 : 1147 : mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
406 : 1147 : cprj_paral_band=(mband_cprj<mband)
407 : :
408 : : !LOOP OVER SPINS
409 : 1147 : ibg=0;icg=0
410 : 1147 : bdtot_index=0
411 [ + + ]: 2294 : do isppol=1,nsppol
412 : :
413 : : ! LOOP OVER k POINTS
414 : 1147 : ikg=0
415 [ + + ]: 38408 : do ikpt=1,nkpt
416 : :
417 [ + + + - ]: 36114 : if (iomode == NO_FILE_OUT) psinablapsi => psinablapsi_out(:,:,:,ikpt)
418 : :
419 : 36114 : etiq=ikpt+(isppol-1)*nkpt
420 : 36114 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
421 [ + + ]: 326000 : master_spfftband=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
422 : :
423 : : ! Select k-points for current proc
424 : 36114 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt))
425 [ + + ]: 36114 : if (mykpt) then
426 : :
427 : : ! Data depending on k-point
428 : 16480 : npw_k=npwarr(ikpt)
429 : 16480 : istwf_k=dtset%istwfk(ikpt)
430 [ - + ]: 16480 : cplex=2;if (istwf_k>1) cplex=1
431 [ + + ]: 65920 : kpoint(:)=dtset%kptns(:,ikpt)
432 : :
433 : : ! Extract cprj for this k-point
434 [ + + ]: 16480 : nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
435 [ + + ]: 16480 : if (mkmem*nsppol/=1) then
436 : 16477 : iorder_cprj=0
437 [ - + + - : 383488 : ABI_MALLOC(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
- + + - +
- + - - +
- + + + +
+ ]
438 : 16477 : call pawcprj_alloc(cprj_k_loc,0,dimcprj)
439 : : call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
440 : : & mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
441 [ - + ]: 16477 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
442 : : else
443 : 3 : cprj_k_loc => cprj
444 : : end if
445 : :
446 : : ! if cprj are distributed over bands, gather them (because we need to mix bands)
447 [ + + ]: 16480 : if (cprj_paral_band) then
448 [ - + + - : 185984 : ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k))
- + + - +
- + - - +
- + + + +
+ ]
449 : 6608 : call pawcprj_alloc(cprj_k,0,dimcprj)
450 : : call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom, &
451 : : & my_nspinor*nband_cprj_k,my_nspinor*mpi_enreg%bandpp,&
452 : 6608 : & dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
453 : : else
454 : 9872 : cprj_k => cprj_k_loc
455 : : end if
456 : :
457 : : ! Compute k+G in cartesian coordinates
458 [ - + + - : 49440 : ABI_MALLOC(kpg_k,(3,npw_k*dtset%nspinor))
+ - + - -
+ - + ]
459 : 16480 : kg_k => kg(:,1+ikg:npw_k+ikg)
460 [ + + ]: 2327848 : do ipw=1,npw_k
461 : : kpg_k(1,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1) &
462 : : & +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2) &
463 : 2311368 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3)
464 : : kpg_k(2,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1) &
465 : : & +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2) &
466 : 2311368 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3)
467 : : kpg_k(3,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1) &
468 : : & +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2) &
469 : 2327848 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3)
470 : : end do
471 [ + + + + ]: 9261952 : kpg_k(:,1:npw_k)=two_pi*kpg_k(:,1:npw_k)
472 [ + + + + : 18956 : if (dtset%nspinor==2) kpg_k(1:3,npw_k+1:2*npw_k)=kpg_k(1:3,1:npw_k)
+ + ]
473 : :
474 : : ! Loops over bands
475 [ + + ]: 149094 : do jb=1,nband_k
476 : 132614 : jwavef=(jb-1)*npw_k*my_nspinor+icg
477 : : !If MPI-IO, compute all (ib,jb) pairs ; if not, compute only ib<=jb
478 [ + + ]: 132614 : ibmax=merge(nband_k,jb,iomode_etsf_mpiio)
479 : : !If MPI-IO, store only ib elements for each jb ; if not, store all (ib,jb) pairs
480 [ + - ]: 130304 : my_jb=merge(1,jb,(iomode_etsf_mpiio .or. diag_only))
481 : : !If diag_only then ibmin = ibmax = jb
482 [ + + ]: 132614 : ibmax=merge(jb,ibmax,diag_only)
483 : 132614 : ibmin=merge(jb,1,diag_only)
484 : :
485 : : ! Fill output arrays with zeros
486 [ + - ]: 132614 : if (store_half_dipoles) then
487 : 132614 : jbshift=(my_jb*(my_jb-1))/2 ; bsize=nband_k
488 : : else
489 : 0 : jbshift=(my_jb-1)*mband ; bsize=mband
490 : : end if
491 [ + + ]: 132614 : if ((.not. diag_only) .or. jb == 1) then
492 [ + + + + : 1600238 : psinablapsi(:,:,jbshift+1:jbshift+bsize)=zero
+ + ]
493 [ + + + + : 1600238 : psinablapsi_paw(:,:,jbshift+1:jbshift+bsize)=zero
+ + ]
494 [ + + + + : 22598 : if (dtset%pawspnorb==1) psinablapsi_soc(:,:,jbshift+1:jbshift+bsize)=zero
+ + + + ]
495 : : end if
496 : :
497 : : ! 2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m>
498 : : ! ----------------------------------------------------------------------------------
499 : : ! Note: <psi_nk|-i.nabla|psi_mk> => Sum_g[ <G|Psi_nk>^* <G|Psi_mk> G ]
500 : :
501 : : ! Select bands for current proc
502 : 132614 : myband=.true.
503 [ + + + + ]: 132614 : if (xmpi_paral==1.and.mpi_enreg%paral_kgb/=1.and.(.not.iomode_etsf_mpiio)) then
504 : 78080 : myband=(abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-mpi_enreg%me_kpt)==0)
505 : : end if
506 [ + - ]: 132614 : if (myband) then
507 : :
508 : : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(ib,iwavef,ipw,cgnm1,cgnm2,ibshift) &
509 : : !$OMP SHARED(ibmin,ibmax,cg,kpg_k,jwavef,jbshift,npw_k,my_nspinor,icg,istwf_k,psinablapsi, &
510 : : !$OMP jb,store_half_dipoles,iomode_etsf_mpiio,mband) &
511 : : !$OMP SCHEDULE(STATIC)
512 [ + + ]: 290778 : do ib=ibmin,ibmax
513 : 158164 : iwavef=(ib-1)*npw_k*my_nspinor+icg
514 : :
515 : : ! (C_nk^*)*C_mk*(k+g) is expressed in cartesian coordinates
516 [ - + ]: 158164 : if (istwf_k>1) then
517 : : !In this case (istwfk>1): Sum_g>=g0[ 2i.Im(<G|Psi_nk>^* <G|Psi_mk> G) ]
518 : : !G=k+g=0 term is included but do not contribute
519 [ # # ]: 0 : do ipw=1,npw_k*my_nspinor
520 : 0 : cgnm2=two*(cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef))
521 [ # # ]: 0 : psinablapsi(2,1:3,jbshift+ib)=psinablapsi(2,1:3,jbshift+ib)+cgnm2*kpg_k(1:3,ipw)
522 : : end do
523 : : else
524 [ + + ]: 21140832 : do ipw=1,npw_k*my_nspinor
525 : 20982668 : cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef)
526 : 20982668 : cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef)
527 [ + + ]: 83930672 : psinablapsi(1,1:3,jbshift+ib)=psinablapsi(1,1:3,jbshift+ib)+cgnm1*kpg_k(1:3,ipw)
528 [ + + ]: 84088836 : psinablapsi(2,1:3,jbshift+ib)=psinablapsi(2,1:3,jbshift+ib)+cgnm2*kpg_k(1:3,ipw)
529 : : end do
530 : : end if
531 : :
532 : : ! Second half of the (n,m) matrix
533 [ - + - - ]: 290778 : if ((ib/=jb).and.(.not.store_half_dipoles).and.(.not.iomode_etsf_mpiio)) then
534 : 0 : ibshift=(ib-1)*mband
535 [ # # ]: 0 : psinablapsi(1,1:3,ibshift+jb)= psinablapsi(1,1:3,jbshift+ib)
536 [ # # ]: 0 : psinablapsi(2,1:3,ibshift+jb)=-psinablapsi(2,1:3,jbshift+ib)
537 : : end if
538 : :
539 : : end do ! ib
540 : : !$OMP END PARALLEL DO
541 : :
542 : : ! Reduction in case of parallelism
543 [ + + ]: 132614 : if (iomode_etsf_mpiio) then
544 : 2310 : call timab(48,1,tsec)
545 [ + + ]: 2310 : if (mpi_enreg%paral_kgb==1) then
546 : 2160 : call xmpi_sum_master(psinablapsi,master,spaceComm_spinorfft,ierr)
547 : : end if
548 [ + + ]: 2310 : if (i_am_master_spfft) then
549 : 1230 : call xmpi_sum(psinablapsi,spaceComm_band,ierr)
550 : : end if
551 : 2310 : call timab(48,2,tsec)
552 : : end if
553 : :
554 : : end if ! myband
555 : :
556 : :
557 : : ! 2-B Computation of <psi_n|p_i><p_j|psi_m>(<phi_i|-i.nabla|phi_j>-<tphi_i|-i.nabla|tphi_j>)
558 : : ! ----------------------------------------------------------------------------------
559 : : ! Non relativistic contribution
560 : : ! Note: <psi|-i.nabla|psi_mk>
561 : : ! => -i Sum_ij[ <p_i|Psi_nk>^* <p_j|Psi_mk> (<Phi_i|Nabla|Phi_j>-<tPhi_i|Nabla|tPhi_j>) ]
562 : :
563 : : ! Select bands for current proc
564 : 132614 : myband=.true.
565 [ + + ]: 132614 : if (mpi_enreg%paral_kgb==1) then
566 : 54384 : myband=(mod(jb-1,mpi_enreg%nproc_band)==mpi_enreg%me_band)
567 : 78230 : else if (xmpi_paral==1) then
568 : 78230 : myband=(abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-mpi_enreg%me_kpt)==0)
569 : : end if
570 [ + + ]: 132614 : if (myband) then
571 : : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
572 : : !$OMP PRIVATE(ib,ibsp,jbsp,ispinor,iatom,itypat,lmn_size,jlmn,ilmn,idir, &
573 : : !$OMP cpnm1,cpnm2,cpnm11,cpnm22,cpnm12,cpnm21,cpnm_11m22,cpnm_21p12,cpnm_21m12, &
574 : : !$OMP nabla_ij,soc_ij,ibshift) &
575 : : !$OMP SHARED(ibmin,ibmax,jb,jbshift,my_nspinor,cplex,natom,dtset,pawtab,cprj_k, &
576 : : !$OMP phisocphj,psinablapsi_paw,psinablapsi_soc,store_half_dipoles,iomode_etsf_mpiio,mband)
577 [ + + ]: 230754 : do ib=ibmin,ibmax
578 : :
579 : 124852 : ibsp=(ib-1)*my_nspinor ; jbsp=(jb-1)*my_nspinor
580 [ + + ]: 250104 : do ispinor=1,my_nspinor
581 : 125252 : ibsp=ibsp+1;jbsp=jbsp+1
582 [ - + ]: 250104 : if (cplex==1) then
583 [ # # ]: 0 : do iatom=1,natom
584 : 0 : itypat=dtset%typat(iatom)
585 : 0 : lmn_size=pawtab(itypat)%lmn_size
586 [ # # ]: 0 : do jlmn=1,lmn_size
587 [ # # ]: 0 : do ilmn=1,lmn_size
588 [ # # ]: 0 : nabla_ij(:)=pawtab(itypat)%nabla_ij(:,ilmn,jlmn)
589 [ # # # # ]: 0 : if (ib>jb) nabla_ij(:)=-pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
590 : 0 : cpnm1=cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn)
591 [ # # ]: 0 : if (dtset%nspinor==2) cpnm1=cpnm1+cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn)
592 [ # # ]: 0 : psinablapsi_paw(2,:,jbshift+ib)=psinablapsi_paw(2,:,jbshift+ib)-cpnm1*nabla_ij(:)
593 : : end do !ilmn
594 : : end do !jlmn
595 : : end do !iatom
596 : : else
597 [ + + ]: 377556 : do iatom=1,natom
598 : 252304 : itypat=dtset%typat(iatom)
599 : 252304 : lmn_size=pawtab(itypat)%lmn_size
600 [ + + ]: 2395988 : do jlmn=1,lmn_size
601 [ + + ]: 18418192 : do ilmn=1,lmn_size
602 [ + + ]: 64589824 : nabla_ij(:)=pawtab(itypat)%nabla_ij(:,ilmn,jlmn)
603 [ + + + + ]: 20022016 : if (ib>jb) nabla_ij(:)=-pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
604 : : cpnm1=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn) &
605 : 16147456 : & +cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn))
606 : : cpnm2=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn) &
607 : 16147456 : & -cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn))
608 [ + + ]: 64589824 : psinablapsi_paw(1,:,jbshift+ib)=psinablapsi_paw(1,:,jbshift+ib)+cpnm2*nabla_ij(:)
609 [ + + ]: 66608256 : psinablapsi_paw(2,:,jbshift+ib)=psinablapsi_paw(2,:,jbshift+ib)-cpnm1*nabla_ij(:)
610 : : end do !ilmn
611 : : end do !jlmn
612 : : end do !iatom
613 : : end if
614 : : end do !ispinor
615 : :
616 : : ! 2-C Computation of Spin-orbit coupling contribution:
617 : : ! 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>]
618 : : ! ----------------------------------------------------------------------------------
619 [ + + ]: 124852 : if (dtset%pawspnorb==1) then
620 : : ! Add: Sum_ij,ss'[<Psi^s_n|p_i><p_j|Psi^s'_m> (Sigma^ss' X g_ij)]
621 : : ! =Sum_ij[ (<Psi^up_n|p_i><p_j|Psi^up_m>-<Psi^dn_n|p_i><p_j|Psi^dn_m>) (Sigma^up-up X g_ij)
622 : : ! +(<Psi^dn_n|p_i><p_j|Psi^up_m>+<Psi^up_n|p_i><p_j|Psi^dn_m>) (Re{Sigma^dn-up} X g_ij)
623 : : ! +(<Psi^dn_n|p_i><p_j|Psi^up_m>-<Psi^up_n|p_i><p_j|Psi^dn_m>) (Im{Sigma^dn-up} X g_ij) ]
624 : : ! where: g_ij = <Phi_i| 1/4 Alpha^2 dV/dr vec(r)/r) |Phi_j>
625 : : ! Note that:
626 : : ! phisocphj(:)%value(re:im,1,idir,ilmn,jlmn) is (Sigma^up-up X g_ij)
627 : : ! phisocphj(:)%value(re:im,2,idir,ilmn,jlmn) is (Sigma^dn-up X g_ij)
628 : : ! Not compatible with parallelization over spinors
629 : 400 : ibsp=1+(ib-1)*dtset%nspinor ; jbsp=1+(jb-1)*dtset%nspinor
630 [ + + ]: 2000 : do iatom=1,natom
631 : 1600 : itypat=dtset%typat(iatom)
632 : 1600 : lmn_size=pawtab(itypat)%lmn_size
633 [ + + ]: 14800 : do jlmn=1,lmn_size
634 [ + + ]: 116800 : do ilmn=1,lmn_size
635 : 102400 : soc_ij => phisocphj(iatom)%value(:,:,:,ilmn,jlmn)
636 : : !Contribution from real part of <Psi^s_n|p_i><p_j|Psi^s'_m>
637 : : cpnm11=cprj_k(iatom,ibsp )%cp(1,ilmn)*cprj_k(iatom,jbsp )%cp(1,jlmn) &
638 : 102400 : & +cprj_k(iatom,ibsp )%cp(2,ilmn)*cprj_k(iatom,jbsp )%cp(2,jlmn)
639 : : cpnm22=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn) &
640 : 102400 : & +cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn)
641 : : cpnm12=cprj_k(iatom,ibsp )%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn) &
642 : 102400 : & +cprj_k(iatom,ibsp )%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn)
643 : : cpnm21=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp )%cp(1,jlmn) &
644 : 102400 : & +cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp )%cp(2,jlmn)
645 : 102400 : cpnm_11m22=cpnm11-cpnm22
646 : 102400 : cpnm_21p12=cpnm21+cpnm12
647 : 102400 : cpnm_21m12=cpnm21-cpnm12
648 [ + + ]: 409600 : do idir=1,3
649 : : psinablapsi_soc(1,idir,jbshift+ib)=psinablapsi_soc(1,idir,jbshift+ib) &
650 : 307200 : & +soc_ij(1,1,idir)*cpnm_11m22+soc_ij(1,2,idir)*cpnm_21p12
651 : : psinablapsi_soc(2,idir,jbshift+ib)=psinablapsi_soc(2,idir,jbshift+ib) &
652 : 409600 : & +soc_ij(2,2,idir)*cpnm_21m12
653 : : end do
654 : : !Contribution from imaginary part of <Psi^s_n|p_i><p_j|Psi^s'_m>
655 : : cpnm11=cprj_k(iatom,ibsp )%cp(1,ilmn)*cprj_k(iatom,jbsp )%cp(2,jlmn) &
656 : 102400 : & -cprj_k(iatom,ibsp )%cp(2,ilmn)*cprj_k(iatom,jbsp )%cp(1,jlmn)
657 : : cpnm22=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn) &
658 : 102400 : & -cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn)
659 : : cpnm12=cprj_k(iatom,ibsp )%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn) &
660 : 102400 : & -cprj_k(iatom,ibsp )%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn)
661 : : cpnm21=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp )%cp(2,jlmn) &
662 : 102400 : & -cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp )%cp(1,jlmn)
663 : 102400 : cpnm_11m22=cpnm11-cpnm22
664 : 102400 : cpnm_21p12=cpnm21+cpnm12
665 : 102400 : cpnm_21m12=cpnm21-cpnm12
666 [ + + ]: 422400 : do idir=1,3
667 : : psinablapsi_soc(1,idir,jbshift+ib)=psinablapsi_soc(1,idir,jbshift+ib) &
668 : 307200 : & -soc_ij(2,2,idir)*cpnm_21m12
669 : : psinablapsi_soc(2,idir,jbshift+ib)=psinablapsi_soc(2,idir,jbshift+ib) &
670 : 409600 : & +soc_ij(1,1,idir)*cpnm_11m22+soc_ij(1,2,idir)*cpnm_21p12
671 : : end do
672 : : end do ! ilmn
673 : : end do ! jlmn
674 : : end do ! iatom
675 : : end if ! pawspnorb
676 : :
677 : : ! Second half of the (n,m) matrix
678 [ - + - - ]: 230754 : if ((ib/=jb).and.(.not.store_half_dipoles).and.(.not.iomode_etsf_mpiio)) then
679 : 0 : ibshift=(ib-1)*mband
680 [ # # ]: 0 : psinablapsi_paw(1,1:3,ibshift+jb)= psinablapsi_paw(1,1:3,jbshift+ib)
681 [ # # ]: 0 : psinablapsi_paw(2,1:3,ibshift+jb)=-psinablapsi_paw(2,1:3,jbshift+ib)
682 [ # # ]: 0 : if (dtset%pawspnorb==1) then
683 [ # # ]: 0 : psinablapsi_soc(1,1:3,ibshift+jb)= psinablapsi_soc(1,1:3,jbshift+ib)
684 [ # # ]: 0 : psinablapsi_soc(2,1:3,ibshift+jb)=-psinablapsi_soc(2,1:3,jbshift+ib)
685 : : end if
686 : : end if
687 : :
688 : : end do ! ib loop
689 : : !$OMP END PARALLEL DO
690 : :
691 [ + + - + ]: 105902 : if (iomode_etsf_mpiio.and.mpi_enreg%paral_spinor==1) then
692 : 0 : call timab(48,1,tsec)
693 : 0 : call xmpi_sum_master(psinablapsi_paw,master,spaceComm_spinor,ierr)
694 : 0 : call timab(48,2,tsec)
695 : : end if
696 : :
697 : : end if ! myband
698 : :
699 : : ! Write to OPT file in case of MPI-IO
700 [ + + ]: 44902 : if (iomode_etsf_mpiio.and.i_am_master_spfft) then
701 [ + + ]: 1230 : if (myband) then
702 [ + - + + : 183550 : psinablapsi=psinablapsi+psinablapsi_paw
+ + + + +
+ + + +
+ ]
703 [ + + + + : 8750 : if (dtset%pawspnorb==1) psinablapsi=psinablapsi+psinablapsi_soc
+ + + + +
+ + + +
+ ]
704 : : end if
705 [ + - ]: 1230 : if (store_half_dipoles) then
706 [ + + ]: 7380 : nc_start_5=[1,1,(jb*(jb-1))/2+1,ikpt,isppol];nc_stride_5=[1,1,1,1,1]
707 [ + + + + ]: 4980 : nc_count_5=[0,0,0,0,0];if (myband) nc_count_5=[2,3,jb,1,1]
708 [ - + ]: 1230 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
709 : : else
710 [ # # ]: 0 : nc_start_6=[1,1,1,jb,ikpt,isppol];nc_stride_6=[1,1,1,1,1,1]
711 [ # # # # ]: 0 : nc_count_6=[0,0,0,0,0,0];if (myband) nc_count_6=[2,3,mband,1,1,1]
712 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
713 : : end if
714 : : end if
715 : :
716 : : end do ! jb
717 : :
718 [ + - ]: 16480 : if (mkmem/=0) then
719 : 16480 : ibg = ibg + my_nspinor*nband_cprj_k
720 : 16480 : icg = icg + npw_k*my_nspinor*nband_k
721 : 16480 : ikg = ikg + npw_k
722 : : end if
723 : :
724 [ + + ]: 16480 : if (cprj_paral_band) then
725 : 6608 : call pawcprj_free(cprj_k)
726 [ - + ]: 6608 : ABI_FREE(cprj_k)
727 : : end if
728 [ + + ]: 16480 : if (mkmem*nsppol/=1) then
729 : 16477 : call pawcprj_free(cprj_k_loc)
730 [ - + ]: 16477 : ABI_FREE(cprj_k_loc)
731 : : end if
732 [ - + ]: 16480 : ABI_FREE(kpg_k)
733 : :
734 : : ! Write to OPT file if not MPI-IO
735 [ + + ]: 16480 : if (iomode /= NO_FILE_OUT) then
736 : : ! >>> Reduction in case of parallelism
737 [ - + ]: 192 : if (.not.iomode_etsf_mpiio) then
738 : 0 : call timab(48,1,tsec)
739 : 0 : call xmpi_sum_master(psinablapsi,master,spaceComm_bandspinorfft,ierr)
740 : 0 : call xmpi_sum_master(psinablapsi_paw,master,spaceComm_bandspinor,ierr)
741 : 0 : call timab(48,2,tsec)
742 [ # # # # : 0 : psinablapsi=psinablapsi+psinablapsi_paw
# # # # #
# # # #
# ]
743 [ # # ]: 0 : if (dtset%pawspnorb==1) then
744 : 0 : call xmpi_sum_master(psinablapsi_soc,master,spaceComm_band,ierr)
745 [ # # # # : 0 : psinablapsi=psinablapsi+psinablapsi_soc
# # # # #
# # # ]
746 : : end if
747 : : end if
748 : :
749 : : ! >>> This my kpt and I am the master node: I write the data
750 : : if (.not.iomode_etsf_mpiio) then
751 [ # # ]: 0 : if (i_am_master) then
752 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
753 : : if (nc_unlimited) then
754 : : 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]
755 : : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
756 [ # # ]: 0 : else if (.not.store_half_dipoles) then
757 [ # # # # ]: 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]
758 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
759 : : else
760 [ # # # # ]: 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]
761 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
762 : : end if
763 : : else
764 [ # # ]: 0 : bsize=nband_k**2;if (store_half_dipoles) bsize=(nband_k*(nband_k+1))/2
765 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,1,ib),ib=1,bsize)
766 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,2,ib),ib=1,bsize)
767 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,3,ib),ib=1,bsize)
768 : : end if
769 : :
770 : : ! >>> This my kpt and I am not the master node: I send the data
771 [ # # ]: 0 : else if (i_am_master_band.and.i_am_master_spfft) then
772 [ # # ]: 0 : if (mpi_enreg%me_kpt/=master_spfftband) then
773 : 0 : ABI_BUG('Problem with band communicator!')
774 : : end if
775 : 0 : call xmpi_exch(psinablapsi,pnp_size,mpi_enreg%me_kpt,psinablapsi,master,spaceComm_kpt,etiq,ierr)
776 : : end if
777 : : end if
778 : : else
779 : : ! >>> Reduction in case of parallelism
780 : 16288 : call timab(48,1,tsec)
781 : 16288 : call xmpi_sum(psinablapsi,spaceComm_bandspinorfft,ierr)
782 : 16288 : call xmpi_sum(psinablapsi_paw,spaceComm_bandspinor,ierr)
783 : 16288 : call timab(48,2,tsec)
784 [ + - + + : 2622368 : psinablapsi=psinablapsi+psinablapsi_paw
+ + + + +
+ + + +
+ ]
785 [ - + ]: 16288 : if (dtset%pawspnorb==1) then
786 : 0 : call xmpi_sum(psinablapsi_soc,spaceComm_band,ierr)
787 [ # # # # : 0 : psinablapsi=psinablapsi+psinablapsi_soc
# # # # #
# # # ]
788 : : end if
789 : : end if ! no_file_out
790 : :
791 : : ! >>> This is not my kpt and I am the master node: I receive the data and I write
792 [ + + + - ]: 19634 : elseif ((.not.iomode_etsf_mpiio).and.i_am_master.and.(iomode/=NO_FILE_OUT)) then ! mykpt
793 : 0 : sender=master_spfftband
794 : 0 : call xmpi_exch(psinablapsi,pnp_size,sender,psinablapsi,master,spaceComm_kpt,etiq,ierr)
795 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
796 : : if (nc_unlimited) then
797 : : 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]
798 : : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
799 [ # # ]: 0 : else if (.not.store_half_dipoles) then
800 [ # # # # ]: 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]
801 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
802 : : else
803 [ # # # # ]: 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]
804 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
805 : : end if
806 : : else
807 [ # # ]: 0 : bsize=nband_k**2;if (store_half_dipoles) bsize=(nband_k*(nband_k+1))/2
808 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,1,ib),ib=1,bsize)
809 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,2,ib),ib=1,bsize)
810 [ # # # # ]: 0 : write(ount)(psinablapsi(1:2,3,ib),ib=1,bsize)
811 : : end if
812 : : end if ! mykpt
813 : :
814 : 37261 : bdtot_index=bdtot_index+nband_k
815 : :
816 : : ! End loop on spin,kpt
817 : : end do ! ikpt
818 : : end do !isppol
819 : :
820 : : ! >>> Last reduction over k-points in case of parallelism
821 : : ! >>> if no output to file
822 [ + + ]: 1147 : if (iomode == NO_FILE_OUT) then
823 : 1121 : call xmpi_sum(psinablapsi_out,spaceComm_kpt,ierr)
824 : : end if
825 : :
826 : : !Close file
827 [ + + + + : 1147 : if ((i_am_master.or.(iomode_etsf_mpiio.and.i_am_master_spfft)).and.iomode /= NO_FILE_OUT) then
+ + ]
828 [ + - ]: 16 : if (iomode==IO_MODE_ETSF) then
829 [ - + ]: 16 : NCF_CHECK(nf90_close(ncid))
830 : : else
831 : 0 : ierr=close_unit(ount,msg)
832 [ # # ]: 0 : ABI_CHECK(ierr==0,"Error while closing OPT file")
833 : : end if
834 : : end if
835 : :
836 : : !Datastructures deallocations
837 [ + + ]: 1147 : if (iomode /= NO_FILE_OUT) then
838 [ - + ]: 26 : ABI_FREE(psinablapsi)
839 : : end if
840 [ - + ]: 1147 : ABI_FREE(psinablapsi_paw)
841 [ + + ]: 1147 : if (dtset%pawspnorb==1) then
842 [ - + ]: 1 : ABI_FREE(psinablapsi_soc)
843 [ + + ]: 5 : do iatom=1,natom
844 [ - + ]: 5 : ABI_FREE(phisocphj(iatom)%value)
845 : : end do
846 [ + - + + : 5 : ABI_FREE(phisocphj)
- + - + ]
847 : : end if
848 [ + - ]: 1147 : if (.not.already_has_nabla) then
849 [ + + ]: 2294 : do itypat=1,dtset%ntypat
850 [ + - ]: 2294 : if (allocated(pawtab(itypat)%nabla_ij)) then
851 : 1147 : ABI_FREE(pawtab(itypat)%nabla_ij)
852 : 1147 : pawtab(itypat)%has_nabla=0
853 : : end if
854 : : end do
855 : : end if
856 : :
857 : : ABI_NVTX_END_RANGE()
858 : : DBG_EXIT("COLL")
859 : :
860 [ - + - + : 2294 : end subroutine optics_paw
- - - - ]
861 : : !!***
862 : :
863 : : !----------------------------------------------------------------------
864 : :
865 : : !!****f* m_paw_optics/optics_paw_core
866 : : !! NAME
867 : : !! optics_paw_core
868 : : !!
869 : : !! FUNCTION
870 : : !! Compute matrix elements need for X spectr. (in the PAW context) and store them in a file
871 : : !! Matrix elements = <Phi_core|Nabla|Phi_j>
872 : : !!
873 : : !! COPYRIGHT
874 : : !! Copyright (C) 2005-2026 ABINIT group (SM,MT,NB)
875 : : !! This file is distributed under the terms of the
876 : : !! GNU General Public License, see ~ABINIT/COPYING
877 : : !! or http://www.gnu.org/copyleft/gpl.txt .
878 : : !!
879 : : !! INPUTS
880 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
881 : : !! cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
882 : : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
883 : : !! dtfil <type(datafiles_type)>=variables related to files
884 : : !! dtset <type(dataset_type)>=all input variables for this dataset
885 : : !! filpsp(ntypat)=name(s) of the pseudopotential file(s)
886 : : !! mband=maximum number of bands
887 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
888 : : !! mkmem =number of k points treated by this node.
889 : : !! mpi_enreg=information about MPI parallelization
890 : : !! mpsang =1+maximum angular momentum for nonlocal pseudopotentials
891 : : !! natom=number of atoms in cell.
892 : : !! nkpt=number of k points.
893 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
894 : : !! pawang <type(pawang_type)>= PAW ANGular mesh discretization and related data
895 : : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
896 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
897 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
898 : : !! znucl(ntypat)=atomic number of atom type
899 : : !!
900 : : !! OUTPUT
901 : : !! (only writing in a file)
902 : : !!
903 : : !! SOURCE
904 : :
905 : 3 : subroutine optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen0,filpsp,hdr,&
906 : : & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,&
907 : 3 : & pawang,pawrad,pawrhoij,pawtab,znucl,rcpaw)
908 : :
909 : : !Arguments ------------------------------------
910 : : !scalars
911 : : integer,intent(in) :: mband,mcprj,mkmem,mpsang,natom,nkpt,nsppol
912 : : type(MPI_type),intent(in) :: mpi_enreg
913 : : type(datafiles_type),intent(in) :: dtfil
914 : : type(dataset_type),intent(in) :: dtset
915 : : type(hdr_type),intent(inout) :: hdr
916 : : type(pawang_type),intent(in) :: pawang
917 : : type(rcpaw_type), intent(in), pointer,optional :: rcpaw
918 : : !arrays
919 : : integer,intent(in) :: atindx1(natom),dimcprj(natom)
920 : : character(len=fnlen),intent(in) :: filpsp(dtset%ntypat)
921 : : real(dp),intent(in) :: eigen0(mband*nkpt*nsppol),znucl(dtset%ntypat)
922 : : type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
923 : : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
924 : : type(pawrhoij_type),intent(inout) :: pawrhoij(mpi_enreg%my_natom)
925 : : type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
926 : :
927 : : !Local variables-------------------------------
928 : : !scalars
929 : : integer,parameter :: master=0
930 : : integer :: bdtot_index,cplex,etiq,iatom,ic,ibg,idir,ilmn2
931 : : integer :: ierr,ikpt,ilmn,iln,ount,is,my_jb,nphicor_lm
932 : : integer :: iorder_cprj,ispinor,isppol,istwf_k,itypat,itypat2
933 : : integer :: jb,jbsp,jlmn,lmn_size,lmncmax,mband_cprj,ncid,varid
934 : : integer :: me,my_nspinor,nband_cprj_k,option_core,pnp_size
935 : : integer :: nband_k,nphicor,ncorespinor,sender,iomode,fformopt,master_spfftband
936 : : integer :: spaceComm_band,spaceComm_bandspinorfft,spaceComm_fft,spaceComm_kpt
937 : : integer :: spaceComm_spinor,spaceComm_bandspinor,spaceComm_spinorfft,spaceComm_w
938 : : logical :: already_has_nabla,cprj_paral_band,mykpt,myband,use_rcpaw_data
939 : : logical :: iomode_etsf_mpiio,use_spinorbit
940 : : logical :: i_am_master,i_am_master_band,i_am_master_spfft
941 : : real(dp) :: cpnm1,cpnm2,el_temp
942 : : character(len=500) :: msg
943 : : !arrays
944 : : integer :: nc_count(6),nc_start(6),nc_stride(6),tmp_shape(3)
945 : 3 : integer,allocatable :: lcor(:,:),ncor(:,:),kappacor(:,:),nphicor_arr(:),lmn2ln(:,:)
946 : : real(dp) :: tsec(2)
947 : 3 : real(dp),allocatable :: energy_cor(:,:),occ_cor(:,:),maxocc_cor(:,:)
948 : 3 : real(dp),allocatable :: psinablapsi(:,:,:,:,:),psinablapsi_soc(:,:,:,:,:),psinablapsi2(:,:,:,:)
949 : 3 : real(dp),pointer :: soc_ij(:,:,:)
950 : 3 : type(coeff5_type),allocatable,target :: phisocphj(:)
951 : 3 : type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
952 : : type(nctkdim_t) :: ncdims(3)
953 : : type(nctkarr_t) :: nctk_arrays(9)
954 : 3 : type(atomorb_type),allocatable :: atm(:)
955 : :
956 : : ! ************************************************************************
957 : :
958 : : DBG_ENTER("COLL")
959 : : ABI_NVTX_START_RANGE(NVTX_OPTICS_PAW_CORE)
960 : :
961 : : !Compatibility tests
962 : 3 : msg="mkmem==0 not supported anymore!"
963 : 0 : ABI_CHECK(mkmem/=0,msg)
964 : : !Probably should check for spinor parallelism, because thats likely to not work correctly
965 : 3 : msg="Spinor parallelism not implemented for optics_paw_core!"
966 [ - + ]: 3 : ABI_CHECK(dtset%npspinor==1,msg)
967 : : !Is mpi_enreg initialized?
968 : : if (xmpi_paral==1) then
969 [ + + ]: 12 : tmp_shape = shape(mpi_enreg%proc_distrb)
970 [ - + ]: 3 : if (nkpt > tmp_shape(1)) then
971 : 0 : ABI_BUG('problem with proc_distrb!')
972 : : end if
973 : : end if
974 : :
975 : : !Init parallelism
976 : 3 : spaceComm_w=mpi_enreg%comm_cell
977 [ - + ]: 3 : if (mpi_enreg%paral_kgb==1) then
978 : 0 : spaceComm_kpt=mpi_enreg%comm_kpt
979 : 0 : spaceComm_fft=mpi_enreg%comm_fft
980 : 0 : spaceComm_band=mpi_enreg%comm_band
981 : 0 : spaceComm_spinor=mpi_enreg%comm_spinor
982 : 0 : spaceComm_bandspinor=mpi_enreg%comm_bandspinor
983 : 0 : spaceComm_spinorfft=mpi_enreg%comm_spinorfft
984 : 0 : spaceComm_bandspinorfft=mpi_enreg%comm_bandspinorfft
985 : : else
986 : 3 : spaceComm_kpt=mpi_enreg%comm_kpt
987 : 3 : spaceComm_fft=xmpi_comm_self
988 : 3 : spaceComm_band=mpi_enreg%comm_band
989 : 3 : spaceComm_spinor=xmpi_comm_self
990 : 3 : spaceComm_bandspinor=spaceComm_band
991 : 3 : spaceComm_spinorfft=xmpi_comm_self
992 : 3 : spaceComm_bandspinorfft=xmpi_comm_self
993 : : end if
994 : 3 : me=xmpi_comm_rank(spaceComm_w)
995 : 3 : i_am_master=(me==master)
996 : 3 : i_am_master_band=(xmpi_comm_rank(spaceComm_band)==master)
997 : 3 : i_am_master_spfft=(xmpi_comm_rank(spaceComm_spinorfft)==master)
998 : 3 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
999 : :
1000 : : !------------------------------------------------------------------------------------------------
1001 : : !1- Reading of core wavefunctions
1002 : : !------------------------------------------------------------------------------------------------
1003 : : !TODO At present, impose 2-spinor simulataneously for core and valence WF
1004 [ - + + - : 12 : ABI_MALLOC(atm,(dtset%ntypat))
+ - + - -
+ - + +
+ ]
1005 [ + + ]: 6 : do itypat=1,dtset%ntypat
1006 : 3 : use_rcpaw_data=.false.
1007 [ + + ]: 3 : if(present(rcpaw)) then
1008 [ + - ]: 1 : if(associated(rcpaw)) then
1009 [ + - ]: 1 : if(allocated(rcpaw%atm)) then
1010 [ + - ]: 1 : if(dtset%ntypat==size(rcpaw%atm)) then
1011 : 1 : call copy_atomorb(rcpaw%atm(itypat),atm(itypat))
1012 : : use_rcpaw_data=.true.
1013 : : endif
1014 : : endif
1015 : : endif
1016 : : endif
1017 : 3 : if(.not.use_rcpaw_data) then
1018 : 2 : call pawpsp_init_core(atm(itypat),filpsp(itypat),radmesh=pawrad(itypat))
1019 : : endif
1020 : : enddo
1021 : :
1022 [ - + ]: 3 : if(dtset%cwfs_wouth==1) then
1023 [ # # ]: 0 : do itypat=1,dtset%ntypat
1024 [ # # ]: 0 : if(atm(itypat)%mult==1) then
1025 [ # # ]: 0 : do itypat2=1,dtset%ntypat
1026 : : if(atm(itypat2)%mult>1.and.atm(itypat2)%znucl==atm(itypat)%znucl.and.&
1027 : : & atm(itypat2)%mesh_size==atm(itypat)%mesh_size.and.&
1028 : : & atm(itypat2)%ln_size==atm(itypat)%ln_size.and.&
1029 [ # # # # : 0 : & atm(itypat2)%nsppol==atm(itypat)%nsppol.and.&
# # # # #
# # # ]
1030 : 0 : & atm(itypat2)%zcore_orig>atm(itypat)%zcore_orig) then
1031 : 0 : write(std_out,*) 'Core wfs of typat ',itypat,' replaced by those of typat ',itypat2
1032 [ # # # # : 0 : atm(itypat)%phi=atm(itypat2)%phi
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
1033 : : exit
1034 : : endif
1035 : : enddo
1036 : : endif
1037 : : enddo
1038 : : endif
1039 : :
1040 : :
1041 : 3 : nphicor=0
1042 : 3 : ncorespinor=0
1043 : 3 : nphicor_lm=0
1044 [ + + ]: 6 : do itypat=1,dtset%ntypat
1045 : 3 : nphicor=max(nphicor,atm(itypat)%ln_size)
1046 : 3 : nphicor_lm=max(nphicor,atm(itypat)%lmn_size)
1047 [ - + ]: 3 : if(atm(itypat)%nsppol>1) ABI_ERROR("nsppol>1 is work in progress for optics_paw_core")
1048 [ - + ]: 3 : if(atm(itypat)%nspinor/=dtset%nspinor) ABI_ERROR("Core and valence not same number of spinors")
1049 : 6 : ncorespinor=max(ncorespinor,atm(itypat)%nspinor)
1050 : : enddo
1051 [ + + ]: 3 : if(dtset%nspinor==2) nphicor_lm=nphicor_lm/2
1052 : :
1053 [ - + + - : 12 : ABI_MALLOC(lmn2ln,(nphicor_lm,dtset%ntypat))
- + + - +
- + - - +
- + ]
1054 [ + + + + ]: 31 : lmn2ln=0
1055 [ + + ]: 6 : do itypat=1,dtset%ntypat
1056 : 3 : lmncmax=atm(itypat)%lmn_size
1057 [ + + ]: 3 : if(dtset%nspinor==2) lmncmax=lmncmax/2
1058 [ + + ]: 31 : do ilmn=1,lmncmax
1059 [ + + ]: 28 : if(dtset%nspinor==1) then
1060 : 5 : lmn2ln(ilmn,itypat)=atm(itypat)%indlmn(5,ilmn)
1061 : : else
1062 [ + - ]: 198 : do ilmn2=1,atm(itypat)%lmn_size
1063 [ + + ]: 198 : if(atm(itypat)%indlmn(9,ilmn2)==ilmn) then
1064 : 20 : lmn2ln(ilmn,itypat)=atm(itypat)%indlmn(5,ilmn2)
1065 : 20 : exit
1066 : : endif
1067 : : enddo
1068 : : endif
1069 : : enddo
1070 : : enddo
1071 : :
1072 : : !----------------------------------------------------------------------------------
1073 : : !2- Computation of phipphj=<phi_i|nabla|phi_core>
1074 : : !----------------------------------------------------------------------------------
1075 : :
1076 [ + - - + ]: 3 : already_has_nabla=all(pawtab(:)%has_nabla==3)
1077 [ + + + - : 3 : if (ncorespinor==2) already_has_nabla=all(pawtab(:)%has_nabla==4)
- + ]
1078 : : ! Should check whether this would work with spinor parallelism
1079 : 3 : call pawnabla_core_init(mpsang,dtset%ntypat,pawrad,pawtab,atm)
1080 : :
1081 : : !Get electronic temperature from dtset
1082 [ + - - - : 3 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
- - ]
1083 : :
1084 : : !Compute spin-orbit contributions if necessary
1085 [ + + + - ]: 3 : use_spinorbit=(dtset%pawspnorb==1.and.dtset%userie/=111) ! For testing purpose
1086 : : if (use_spinorbit) then
1087 : 2 : option_core=1
1088 : : call pawnabla_soc_init(el_temp,phisocphj,option_core,dtset%ixc,mpi_enreg%my_natom,natom,&
1089 : : & dtset%nspden,dtset%ntypat,pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,&
1090 : : & dtset%spnorbscl,dtset%typat,dtset%xc_denpos,dtset%xc_taupos,znucl,&
1091 : 2 : & atm=atm,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1092 : : end if
1093 : :
1094 : : !----------------------------------------------------------------------------------
1095 : : !3- Opening of OPT2 file and header writing
1096 : : !----------------------------------------------------------------------------------
1097 : :
1098 : : !I/O mode is netCDF or Fortran
1099 : : iomode=merge(IO_MODE_ETSF,IO_MODE_FORTRAN_MASTER,dtset%iomode==IO_MODE_ETSF)
1100 : 3 : if (use_netcdf_forced) iomode=IO_MODE_ETSF
1101 : :
1102 [ - + + - : 12 : ABI_MALLOC(energy_cor,(nphicor,dtset%ntypat))
- + + - +
- + - - +
- + ]
1103 [ - + + - : 12 : ABI_MALLOC(ncor,(nphicor,dtset%ntypat))
- + + - +
- + - - +
- + ]
1104 [ - + + - : 9 : ABI_MALLOC(lcor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1105 [ - + + - : 9 : ABI_MALLOC(kappacor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1106 [ - + + - : 9 : ABI_MALLOC(occ_cor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1107 [ - + + - : 9 : ABI_MALLOC(maxocc_cor,(nphicor,dtset%ntypat))
- + + - -
+ - + ]
1108 [ + + + + ]: 17 : energy_cor=zero
1109 [ + + + + ]: 17 : ncor=0
1110 [ + + + + ]: 17 : lcor=0
1111 [ + + + + ]: 17 : kappacor=0
1112 [ + + + + ]: 17 : occ_cor=zero
1113 [ + + + + ]: 17 : maxocc_cor=zero
1114 [ + + ]: 6 : do itypat=1,dtset%ntypat
1115 [ + + ]: 17 : do iln=1,atm(itypat)%ln_size
1116 : 11 : energy_cor(iln,itypat)=atm(itypat)%eig(iln,1)
1117 : 11 : ncor(iln,itypat)=atm(itypat)%indln(2,iln)
1118 : 11 : lcor(iln,itypat)=atm(itypat)%indln(1,iln)
1119 : 11 : occ_cor(iln,itypat)=atm(itypat)%occ(iln,1)
1120 : 11 : maxocc_cor(iln,itypat)=atm(itypat)%max_occ(iln,1)
1121 [ + + ]: 14 : if(atm(itypat)%dirac) then
1122 : 8 : kappacor(iln,itypat)=atm(itypat)%kappa(iln)
1123 : : else
1124 : 3 : kappacor(iln,itypat)=0
1125 : : endif
1126 : : enddo
1127 : : enddo
1128 : :
1129 : : !(master proc only)
1130 [ + - ]: 3 : if (i_am_master) then
1131 : : ! ====> NETCDF format
1132 : : if (iomode==IO_MODE_ETSF) then
1133 : 3 : fformopt=611
1134 [ - + + - : 9 : ABI_MALLOC(nphicor_arr,(dtset%ntypat))
+ - + - -
+ - + ]
1135 [ + + ]: 6 : do itypat=1,dtset%ntypat
1136 : 6 : nphicor_arr(itypat)=atm(itypat)%ln_size
1137 : : enddo
1138 : : ! Open/create nc file
1139 [ - + ]: 3 : NCF_CHECK(nctk_open_create(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),xmpi_comm_self))
1140 : : ! Write header data
1141 [ - + ]: 3 : NCF_CHECK(hdr%ncwrite(ncid,fformopt,nc_define=.true.))
1142 : : ! Define additional dimensions
1143 : 3 : ncdims(1)%name="number_of_atom_types"
1144 : 3 : ncdims(1)%value=dtset%ntypat
1145 : 3 : ncdims(2)%name="max_number_of_core_states"
1146 : 3 : ncdims(2)%value=nphicor
1147 : 3 : ncdims(3)%name="number_of_core_spinor_components"
1148 : 3 : ncdims(3)%value=ncorespinor
1149 [ - + ]: 3 : NCF_CHECK(nctk_def_dims(ncid,ncdims))
1150 : 3 : nctk_arrays(1)%name="number_of_core_states"
1151 : 3 : nctk_arrays(1)%dtype="int"
1152 : 3 : nctk_arrays(1)%shape_str="number_of_atom_types"
1153 : 3 : nctk_arrays(2)%name="eigenvalues_core"
1154 : 3 : nctk_arrays(2)%dtype="dp"
1155 : 3 : nctk_arrays(2)%shape_str="max_number_of_core_states,number_of_atom_types"
1156 : 3 : nctk_arrays(3)%name="dipole_core_valence"
1157 : 3 : nctk_arrays(3)%dtype="dp"
1158 : : nctk_arrays(3)%shape_str=&
1159 : : & "number_of_cartesian_directions,max_number_of_core_states,"// &
1160 : 3 : & "number_of_atoms,max_number_of_states,number_of_kpoints,number_of_spins"
1161 : 3 : nctk_arrays(4)%name="n_quantum_number_core"
1162 : 3 : nctk_arrays(4)%dtype="int"
1163 : 3 : nctk_arrays(4)%shape_str="max_number_of_core_states,number_of_atom_types"
1164 : 3 : nctk_arrays(5)%name="l_quantum_number_core"
1165 : 3 : nctk_arrays(5)%dtype="int"
1166 : 3 : nctk_arrays(5)%shape_str="max_number_of_core_states,number_of_atom_types"
1167 : 3 : nctk_arrays(6)%name="kappa_core"
1168 : 3 : nctk_arrays(6)%dtype="int"
1169 : 3 : nctk_arrays(6)%shape_str="max_number_of_core_states,number_of_atom_types"
1170 : 3 : nctk_arrays(7)%name="number_of_core_states"
1171 : 3 : nctk_arrays(7)%dtype="int"
1172 : 3 : nctk_arrays(7)%shape_str="number_of_atom_types"
1173 : 3 : nctk_arrays(8)%name="occupation_core"
1174 : 3 : nctk_arrays(8)%dtype="dp"
1175 : 3 : nctk_arrays(8)%shape_str="max_number_of_core_states,number_of_atom_types"
1176 : 3 : nctk_arrays(9)%name="max_occupation_core"
1177 : 3 : nctk_arrays(9)%dtype="dp"
1178 : 3 : nctk_arrays(9)%shape_str="max_number_of_core_states,number_of_atom_types"
1179 [ - + ]: 3 : NCF_CHECK(nctk_def_arrays(ncid, nctk_arrays))
1180 [ - + ]: 3 : NCF_CHECK(nctk_set_atomic_units(ncid, "eigenvalues_core"))
1181 [ - + ]: 3 : NCF_CHECK(nctk_set_atomic_units(ncid, "dipole_core_valence"))
1182 : : ! Write core states
1183 [ - + ]: 3 : NCF_CHECK(nctk_set_datamode(ncid))
1184 : 3 : varid=nctk_idname(ncid,"eigenvalues_core")
1185 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ncid,varid,energy_cor))
1186 : 3 : varid=nctk_idname(ncid,"n_quantum_number_core")
1187 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ncid,varid,ncor))
1188 : 3 : varid=nctk_idname(ncid,"l_quantum_number_core")
1189 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ncid,varid,lcor))
1190 : 3 : varid=nctk_idname(ncid,"kappa_core")
1191 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ncid,varid,kappacor))
1192 : 3 : varid=nctk_idname(ncid,"occupation_core")
1193 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ncid,varid,occ_cor))
1194 : 3 : varid=nctk_idname(ncid,"max_occupation_core")
1195 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ncid,varid,maxocc_cor))
1196 : 3 : varid=nctk_idname(ncid,"number_of_core_states")
1197 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ncid,varid,nphicor_arr))
1198 : : ! Write eigenvalues
1199 : 3 : varid=nctk_idname(ncid,"eigenvalues")
1200 [ + + - + : 12 : NCF_CHECK(nf90_put_var(ncid,varid,reshape(eigen0,[mband,nkpt,nsppol])))
- - ]
1201 : : !Close file here because the rest has possibly to be written with collective I/O
1202 [ - + ]: 3 : NCF_CHECK(nf90_close(ncid))
1203 [ - + ]: 3 : ABI_FREE(nphicor_arr)
1204 : : ! ====> Standard FORTRAN binary file format
1205 : : else if (iomode==IO_MODE_FORTRAN_MASTER) then
1206 : : fformopt=612 ! MT 12sept21: change the OPT2 Fortran file format
1207 : : if (2*nphicor*natom*mband>2**30) fformopt=613 ! Format for large file records
1208 : : if (open_file(dtfil%fnameabo_app_opt2,msg,newunit=ount,form="unformatted",status="unknown")/= 0) then
1209 : : ABI_ERROR(msg)
1210 : : end if
1211 : : call hdr%fort_write(ount,fformopt,ierr,rewind=.true.)
1212 : : write(ount)(eigen0(jb),jb=1,mband*nkpt*nsppol)
1213 : : write(ount) nphicor
1214 : : do itypat=1,dtset%ntypat
1215 : : write(ount) atm(itypat)%ln_size
1216 : : do iln=1,nphicor
1217 : : write(ount) ncor(iln,itypat),lcor(iln,itypat),kappacor(iln,itypat),occ_cor(iln,itypat),maxocc_cor(iln,itypat),&
1218 : : energy_cor(iln,itypat)
1219 : : end do
1220 : : enddo
1221 : : else
1222 : : msg = "Wrong OPT2 file format!"
1223 : : ABI_BUG(msg)
1224 : : end if ! File format
1225 : : end if ! master node
1226 : 3 : call xmpi_bcast(iomode,master,spaceComm_w,ierr) ! Seems mandatory; why ?
1227 : 3 : call xmpi_bcast(fformopt,master,spaceComm_w,ierr)
1228 [ + - - + ]: 3 : iomode_etsf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
1229 : :
1230 [ - + ]: 3 : ABI_FREE(ncor)
1231 [ - + ]: 3 : ABI_FREE(lcor)
1232 [ - + ]: 3 : ABI_FREE(kappacor)
1233 [ - + ]: 3 : ABI_FREE(occ_cor)
1234 [ - + ]: 3 : ABI_FREE(maxocc_cor)
1235 [ - + ]: 3 : ABI_FREE(energy_cor)
1236 : :
1237 : : !----------------------------------------------------------------------------------
1238 : : !4- Computation of <psi_n|p_i>(<phi_i|-i.nabla|phi_core>)
1239 : : !----------------------------------------------------------------------------------
1240 : :
1241 : : !Prepare core-valence dipoles writing
1242 : : !In case of netCDF access to OPT2 file, prepare collective I/O
1243 [ + - ]: 3 : if (iomode == IO_MODE_ETSF) then
1244 [ + - ]: 3 : if (iomode_etsf_mpiio) then
1245 [ + - ]: 3 : if (i_am_master_spfft) then
1246 [ - + ]: 3 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),spaceComm_band))
1247 : 3 : varid=nctk_idname(ncid,"dipole_core_valence")
1248 [ - + ]: 3 : if (xmpi_comm_size(spaceComm_w)>1) then
1249 [ # # ]: 0 : NCF_CHECK(nctk_set_collective(ncid,varid))
1250 : : end if
1251 [ - + ]: 3 : NCF_CHECK(nctk_set_datamode(ncid))
1252 : : end if
1253 [ # # ]: 0 : else if (i_am_master) then
1254 [ # # ]: 0 : NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),xmpi_comm_self))
1255 : 0 : varid=nctk_idname(ncid,"dipole_core_valence")
1256 : : if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
1257 : : NCF_CHECK(nctk_set_collective(ncid,varid))
1258 : : end if
1259 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
1260 : : end if
1261 : : end if
1262 [ + - ]: 3 : if (iomode_etsf_mpiio) then
1263 : : !If MPI-IO, store only elements for one band
1264 [ - + + - : 12 : ABI_MALLOC(psinablapsi,(2,3,nphicor_lm,natom,1))
- + + - +
- + - - +
- + ]
1265 [ - + + - : 12 : ABI_MALLOC(psinablapsi2,(3,nphicor,natom,1))
- + + - +
- + - - +
- + ]
1266 [ + + ]: 3 : if (use_spinorbit) then
1267 [ - + + - : 6 : ABI_MALLOC(psinablapsi_soc,(2,3,nphicor_lm,natom,1))
- + + - -
+ - + ]
1268 : : end if
1269 : : else
1270 : : !If not, store the elements for all bands
1271 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,nphicor_lm,natom,mband))
# # # # #
# # # # #
# # # # #
# ]
1272 [ # # # # : 0 : ABI_MALLOC(psinablapsi2,(3,nphicor,natom,mband))
# # # # #
# # # # #
# # # # #
# ]
1273 [ # # ]: 0 : if (use_spinorbit) then
1274 [ # # # # : 0 : ABI_MALLOC(psinablapsi_soc,(2,3,nphicor_lm,natom,mband))
# # # # #
# # # # #
# # ]
1275 : : end if
1276 : : end if
1277 [ + + ]: 15 : pnp_size=size(psinablapsi2)
1278 : :
1279 : : !Determine if cprj datastructure is distributed over bands
1280 : 3 : mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
1281 : 3 : cprj_paral_band=(mband_cprj<mband)
1282 : :
1283 : : !LOOP OVER SPINS
1284 : 3 : ibg=0
1285 : 3 : bdtot_index=0
1286 [ + + ]: 6 : do isppol=1,nsppol
1287 : :
1288 : : ! LOOP OVER k POINTS
1289 [ + + ]: 9 : do ikpt=1,nkpt
1290 : :
1291 : 3 : etiq=ikpt+(isppol-1)*nkpt
1292 : 3 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1293 [ + + ]: 53 : master_spfftband=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
1294 : :
1295 : : ! Select k-points for current proc
1296 : 3 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt))
1297 [ + - ]: 3 : if (mykpt) then
1298 : :
1299 : : ! Data depending on k-point
1300 : 3 : istwf_k=dtset%istwfk(ikpt)
1301 [ - + ]: 3 : cplex=2;if (istwf_k>1) cplex=1
1302 : :
1303 : : ! Extract cprj for this k-point
1304 [ - + ]: 3 : nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
1305 [ - + ]: 3 : if (mkmem*nsppol/=1) then
1306 : 0 : iorder_cprj=0
1307 [ # # # # : 0 : ABI_MALLOC(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
# # # # #
# # # # #
# # # # #
# ]
1308 : 0 : call pawcprj_alloc(cprj_k_loc,0,dimcprj)
1309 : : call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
1310 : : & mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
1311 [ # # ]: 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1312 : : else
1313 : 3 : cprj_k_loc => cprj
1314 : : end if
1315 : :
1316 : : ! if cprj are distributed over bands, gather them (because we need to mix bands)
1317 [ - + ]: 3 : if (cprj_paral_band) then
1318 [ # # # # : 0 : ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k))
# # # # #
# # # # #
# # # # #
# ]
1319 : 0 : call pawcprj_alloc(cprj_k,0,dimcprj)
1320 : : call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom,my_nspinor*nband_cprj_k, &
1321 : : & my_nspinor*mpi_enreg%bandpp,&
1322 : 0 : & dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
1323 : : else
1324 : 3 : cprj_k => cprj_k_loc
1325 : : end if
1326 : :
1327 : : ! Loops over bands
1328 [ + + ]: 53 : do jb=1,nband_k
1329 : : !If MPI-IO, store only ib elements for each jb ; if not, store all (ib,jb) pairs
1330 [ - + ]: 50 : my_jb=merge(1,jb,iomode_etsf_mpiio)
1331 : :
1332 [ + + + + : 12190 : psinablapsi(:,:,:,:,my_jb)=zero
+ + + + ]
1333 [ + + + + : 2270 : psinablapsi2(:,:,:,my_jb)=zero
+ + ]
1334 [ + + + + : 10150 : if (use_spinorbit) psinablapsi_soc(:,:,:,:,my_jb)=zero
+ + + + +
+ ]
1335 : :
1336 : : ! Computation of <psi_n|p_i><phi_i|-i.nabla|phi_core>
1337 : : ! ----------------------------------------------------------------------------------
1338 : :
1339 : : ! Select bands for current proc
1340 : 50 : myband=.true.
1341 [ - + ]: 50 : if (mpi_enreg%paral_kgb==1) then
1342 : 0 : myband=(mod(jb-1,mpi_enreg%nproc_band)==mpi_enreg%me_band)
1343 : 50 : else if (xmpi_paral==1) then
1344 : 50 : myband=(abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-mpi_enreg%me_kpt)==0)
1345 : : end if
1346 [ + - ]: 50 : if (myband) then
1347 : :
1348 : 50 : jbsp=(jb-1)*my_nspinor
1349 : :
1350 : : ! 1-spinor case
1351 [ + + + - ]: 50 : if (dtset%nspinor==1.and.ncorespinor==1) then
1352 : 10 : jbsp=jbsp+1
1353 [ - + ]: 10 : if (cplex==1) then ! Real WF case
1354 [ # # ]: 0 : do iatom=1,natom
1355 : 0 : itypat=dtset%typat(iatom)
1356 : 0 : lmn_size=pawtab(itypat)%lmn_size
1357 : 0 : lmncmax=atm(itypat)%lmn_size
1358 [ # # ]: 0 : do jlmn=1,lmn_size
1359 [ # # ]: 0 : do ilmn=1,lmncmax
1360 : 0 : ic=ilmn !atm(itypat)%indlmn(5,ilmn)
1361 : 0 : cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
1362 : : psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1363 [ # # ]: 0 : & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1364 : : end do !ilmn
1365 : : end do !jlmn
1366 : : end do !iatom
1367 : : else ! Complex WF case
1368 [ + + ]: 50 : do iatom=1,natom
1369 : 40 : itypat=dtset%typat(iatom)
1370 : 40 : lmn_size=pawtab(itypat)%lmn_size
1371 : 40 : lmncmax=atm(itypat)%lmn_size
1372 [ + + ]: 370 : do jlmn=1,lmn_size
1373 [ + + ]: 1960 : do ilmn=1,lmncmax
1374 : 1600 : ic=ilmn !atm(itypat)%indlmn(5,ilmn)
1375 : 1600 : cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
1376 : 1600 : cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn)
1377 : : psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1378 [ + + ]: 6400 : & -cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1379 : : psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1380 [ + + ]: 6720 : & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1381 : : end do !ilmn
1382 : : end do !jlmn
1383 : : end do !iatom
1384 : : end if
1385 : :
1386 : : ! 2-spinor case
1387 [ + - + - ]: 40 : else if (dtset%nspinor==2.and.ncorespinor==2) then
1388 [ + + ]: 120 : do ispinor=1,my_nspinor
1389 : 80 : jbsp=jbsp+1
1390 [ + + ]: 320 : do iatom=1,natom
1391 : 200 : itypat=dtset%typat(iatom)
1392 : 200 : lmn_size=pawtab(itypat)%lmn_size
1393 : 200 : lmncmax=atm(itypat)%lmn_size
1394 [ + + ]: 2280 : do jlmn=1,lmn_size
1395 [ + + ]: 42200 : do ilmn=1,lmncmax
1396 : 40000 : is=atm(itypat)%indlmn(6,ilmn)
1397 [ + + ]: 42000 : if (modulo(jbsp,2)==modulo(is,2)) then ! Nabla is a spin-diagonal operator
1398 : 20000 : ic=atm(itypat)%indlmn(9,ilmn)
1399 [ + + ]: 20000 : if (ic>0) then
1400 : 14000 : cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
1401 : 14000 : cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn)
1402 : : psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1403 : : & +cpnm1*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
1404 [ + + ]: 56000 : & -cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1405 : : psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1406 : : & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn) &
1407 [ + + ]: 56000 : & -cpnm2*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn)
1408 : : end if
1409 : : end if
1410 : : end do ! ilmn
1411 : : end do !jlmn
1412 : : end do !iatom
1413 : : end do !ispinor
1414 : : else
1415 : 0 : msg="Core and valence WF should have the same spinor representation!"
1416 : 0 : ABI_BUG(msg)
1417 : : !N. Brouwer initial coding: should be justified!
1418 : : !if (dtset%nspinor==1.and.ncorespinor==2) then ! Average core spinors
1419 : : ! psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1420 : : ! & +half_sqrt2*(cpnm1*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
1421 : : ! & +cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn))
1422 : : ! psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1423 : : ! & +half_sqrt2*(cpnm2*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
1424 : : ! & -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn))
1425 : : !else if (dtset%nspinor==2.and.ncorespinor==1) then ! Average valence spinors
1426 : : ! psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1427 : : ! & +half_sqrt2*cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1428 : : ! psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1429 : : ! & -half_sqrt2*cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1430 : : !endif
1431 : : endif
1432 : :
1433 : : ! Spin-orbit coupling contribution:
1434 : : ! Sum_i,ss'[<psi^s_n|p_i><phi_i|1/4 Alpha^2 (Sigma^ss' X dV/dr)|phj_core^s'>]
1435 [ + + ]: 50 : if (use_spinorbit) then
1436 : : ! Add: Sum_i,ss'[<Psi^s_n|p_i> (Sigma^ss' X g_ij_core^s']
1437 : : ! where: g_ij_core^s = <Phi_i| 1/4 Alpha^2 dV/dr vec(r)/r) |Phj_core^s>
1438 : : ! Note that:
1439 : : ! if phi_cor_jlmn_cor(:) is up:
1440 : : ! phisocphj(:)%value(re:im,1,idir,ilmn,jlmn_cor) is (Sigma^up-up X g_ij_core^up)
1441 : : ! phisocphj(:)%value(re:im,2,idir,ilmn,jlmn_cor) is (Sigma^dn-up X g_ij_core^up)
1442 : : ! if phi_cor_jlmn_cor(:) is down:
1443 : : ! phisocphj(:)%value(re:im,1,idir,ilmn,jlmn_cor) is (Sigma^up-dn X g_ij_core^dn)
1444 : : ! phisocphj(:)%value(re:im,2,idir,ilmn,jlmn_cor) is (Sigma^dn-dn X g_ij_core^dn)
1445 : : ! Not compatible with parallelization over spinors
1446 : 40 : jbsp=1+(jb-1)*dtset%nspinor
1447 [ + + ]: 140 : do iatom=1,natom
1448 : 100 : itypat=dtset%typat(iatom)
1449 : 100 : lmn_size=pawtab(itypat)%lmn_size
1450 : 100 : lmncmax=atm(itypat)%lmn_size
1451 [ + + ]: 1140 : do jlmn=1,lmn_size
1452 [ + + ]: 21100 : do ilmn=1,lmncmax
1453 [ + - ]: 20000 : if(dtset%nspinor==1) then
1454 : : ic=ilmn !atm(itypat)%indlmn(5,ilmn)
1455 : : else
1456 : 20000 : ic=atm(itypat)%indlmn(9,ilmn)
1457 : : endif
1458 [ + + ]: 21000 : if (ic>0) then
1459 : 14000 : soc_ij => phisocphj(iatom)%value(:,:,:,jlmn,ilmn)
1460 : : !Contribution from real part of <Psi^s_n|p_i>
1461 : 14000 : cpnm1=cprj_k(iatom,jbsp )%cp(1,jlmn)
1462 : 14000 : cpnm2=cprj_k(iatom,jbsp+1)%cp(1,jlmn)
1463 [ + + ]: 56000 : do idir=1,3
1464 : : psinablapsi_soc(1,idir,ic,iatom,my_jb)=psinablapsi_soc(1,idir,ic,iatom,my_jb) &
1465 : 42000 : & +soc_ij(1,1,idir)*cpnm1+soc_ij(1,2,idir)*cpnm2
1466 : : psinablapsi_soc(2,idir,ic,iatom,my_jb)=psinablapsi_soc(2,idir,ic,iatom,my_jb) &
1467 : 56000 : & +soc_ij(2,1,idir)*cpnm1+soc_ij(2,2,idir)*cpnm2
1468 : : end do
1469 : : !Contribution from imaginary part of <Psi^s_n|p_i>
1470 : 14000 : cpnm1=cprj_k(iatom,jbsp )%cp(2,jlmn)
1471 : 14000 : cpnm2=cprj_k(iatom,jbsp+1)%cp(2,jlmn)
1472 [ + + ]: 56000 : do idir=1,3
1473 : : psinablapsi_soc(1,idir,ic,iatom,my_jb)=psinablapsi_soc(1,idir,ic,iatom,my_jb) &
1474 : 42000 : & +soc_ij(2,1,idir)*cpnm1+soc_ij(2,2,idir)*cpnm2
1475 : : psinablapsi_soc(2,idir,ic,iatom,my_jb)=psinablapsi_soc(2,idir,ic,iatom,my_jb) &
1476 : 56000 : & -soc_ij(1,1,idir)*cpnm1-soc_ij(1,2,idir)*cpnm2
1477 : : end do
1478 : : end if
1479 : : end do ! ilmn
1480 : : end do ! jlmn
1481 : : end do ! iatom
1482 : : end if ! use_spinorbit
1483 : :
1484 [ + - - + ]: 50 : if (iomode_etsf_mpiio.and.mpi_enreg%paral_spinor==1) then
1485 : 0 : call timab(48,1,tsec)
1486 : 0 : call xmpi_sum_master(psinablapsi,master,spaceComm_spinor,ierr)
1487 : 0 : call timab(48,2,tsec)
1488 : : end if
1489 : :
1490 : : end if ! myband
1491 : :
1492 : :
1493 : : ! Write to OPT2 file in case of MPI-IO
1494 [ + - ]: 53 : if (iomode_etsf_mpiio.and.i_am_master_spfft) then
1495 [ + + ]: 350 : nc_start=[1,1,1,jb,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
1496 [ + - ]: 50 : if (myband) then
1497 [ + + - + : 10230 : if (use_spinorbit) psinablapsi=psinablapsi+psinablapsi_soc
+ + + + +
+ + + +
+ ]
1498 : : ! Take the square and sum over m
1499 [ + + ]: 190 : do iatom=1,natom
1500 : 140 : itypat=dtset%typat(iatom)
1501 : 140 : lmncmax=atm(itypat)%lmn_size
1502 [ + + ]: 140 : if(dtset%nspinor==2)lmncmax=lmncmax/2
1503 [ + + ]: 1390 : do ilmn=1,lmncmax
1504 : 1200 : ic=lmn2ln(ilmn,itypat)
1505 : : psinablapsi2(:,ic,iatom,my_jb)=psinablapsi2(:,ic,iatom,my_jb)+psinablapsi(1,:,ilmn,iatom,my_jb)**2+&
1506 [ + + ]: 4940 : & psinablapsi(2,:,ilmn,iatom,my_jb)**2
1507 : : enddo
1508 : : enddo
1509 [ + + ]: 350 : nc_count=[3,nphicor,natom,1,1,1]
1510 : : else
1511 : 0 : nc_count=[0,0,0,0,0,0]
1512 : : end if
1513 [ - + ]: 50 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
1514 : : end if
1515 : :
1516 : : end do ! jb
1517 : :
1518 [ + - ]: 3 : if (mkmem/=0) then
1519 : 3 : ibg = ibg + my_nspinor*nband_cprj_k
1520 : : end if
1521 : :
1522 [ - + ]: 3 : if (cprj_paral_band) then
1523 : 0 : call pawcprj_free(cprj_k)
1524 [ # # ]: 0 : ABI_FREE(cprj_k)
1525 : : end if
1526 [ - + ]: 3 : if (mkmem*nsppol/=1) then
1527 : 0 : call pawcprj_free(cprj_k_loc)
1528 [ # # ]: 0 : ABI_FREE(cprj_k_loc)
1529 : : end if
1530 : :
1531 : : ! Write to OPT2 file if not MPI-IO
1532 : :
1533 : : ! >>> Reduction in case of parallelism
1534 [ - + ]: 3 : if (.not.iomode_etsf_mpiio) then
1535 : 0 : call timab(48,1,tsec)
1536 : 0 : call xmpi_sum_master(psinablapsi,master,spaceComm_bandspinor,ierr)
1537 : 0 : call timab(48,2,tsec)
1538 [ # # ]: 0 : if (use_spinorbit) then
1539 : 0 : call xmpi_sum_master(psinablapsi_soc,master,spaceComm_band,ierr)
1540 [ # # # # : 0 : psinablapsi=psinablapsi+psinablapsi_soc
# # # # #
# # # ]
1541 : : end if
1542 : : ! Take the square and sum over m
1543 [ # # ]: 0 : do jb=1,nband_k
1544 : 0 : my_jb=merge(1,jb,iomode_etsf_mpiio)
1545 [ # # ]: 0 : do iatom=1,natom
1546 : 0 : itypat=dtset%typat(iatom)
1547 : 0 : lmncmax=atm(itypat)%lmn_size
1548 [ # # ]: 0 : if(dtset%nspinor==2)lmncmax=lmncmax/2
1549 [ # # ]: 0 : do ilmn=1,lmncmax
1550 : 0 : ic=lmn2ln(ilmn,itypat)
1551 : : psinablapsi2(:,ic,iatom,my_jb)=psinablapsi2(:,ic,iatom,my_jb)+psinablapsi(1,:,ilmn,iatom,my_jb)**2+&
1552 [ # # ]: 0 : & psinablapsi(2,:,ilmn,iatom,my_jb)**2
1553 : : enddo
1554 : : enddo
1555 : : enddo
1556 : : end if
1557 : :
1558 : : ! >>> This my kpt and I am the master node: I write the data
1559 : : if (.not.iomode_etsf_mpiio) then
1560 [ # # ]: 0 : if (i_am_master) then
1561 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
1562 [ # # ]: 0 : nc_start=[1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
1563 [ # # ]: 0 : nc_count=[3,nphicor,natom,mband,1,1]
1564 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
1565 : : else
1566 [ # # ]: 0 : if (fformopt==612) then ! New OPT2 file format
1567 : 0 : write(ount) (((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1568 : 0 : write(ount) (((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1569 : 0 : write(ount) (((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1570 [ # # ]: 0 : else if (fformopt==613) then ! Large OPT2 file format
1571 [ # # ]: 0 : do jb=1,nband_k
1572 : 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1573 : 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1574 : 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1575 : : end do
1576 : : else ! Old OPT2 file format
1577 : : !The old writing was not efficient (indexes order is bad)
1578 [ # # ]: 0 : do iatom=1,natom
1579 [ # # # # ]: 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1580 [ # # # # ]: 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1581 [ # # # # ]: 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1582 : : end do
1583 : : end if
1584 : : end if
1585 : :
1586 : : ! >>> This my kpt and I am not the master node: I send the data
1587 [ # # ]: 0 : else if (i_am_master_band.and.i_am_master_spfft) then
1588 [ # # ]: 0 : if (mpi_enreg%me_kpt/=master_spfftband) then
1589 : 0 : ABI_BUG('Problem with band communicator!')
1590 : : end if
1591 : 0 : call xmpi_exch(psinablapsi2,pnp_size,mpi_enreg%me_kpt,psinablapsi2,master,spaceComm_kpt,etiq,ierr)
1592 : : end if
1593 : : end if
1594 : :
1595 : : ! >>> This is not my kpt and I am the master node: I receive the data and I write
1596 [ # # ]: 0 : elseif ((.not.iomode_etsf_mpiio).and.i_am_master) then ! mykpt
1597 : 0 : sender=master_spfftband
1598 : 0 : call xmpi_exch(psinablapsi2,pnp_size,sender,psinablapsi2,master,spaceComm_kpt,etiq,ierr)
1599 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
1600 [ # # ]: 0 : nc_start=[1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
1601 [ # # ]: 0 : nc_count=[3,nphicor,natom,mband,1,1]
1602 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
1603 : : else
1604 [ # # ]: 0 : if (fformopt==612) then ! New OPT2 file format
1605 : 0 : write(ount) (((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1606 : 0 : write(ount) (((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1607 : 0 : write(ount) (((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1608 [ # # ]: 0 : else if (fformopt==613) then ! Large OPT2 file format
1609 [ # # ]: 0 : do jb=1,nband_k
1610 : 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1611 : 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1612 : 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1613 : : end do
1614 : : else ! Old OPT2 file format
1615 : : !The old writing was not efficient (indexes order is bad)
1616 [ # # ]: 0 : do iatom=1,natom
1617 [ # # # # ]: 0 : write(ount) ((psinablapsi2(1,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1618 [ # # # # ]: 0 : write(ount) ((psinablapsi2(2,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1619 [ # # # # ]: 0 : write(ount) ((psinablapsi2(3,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1620 : : end do
1621 : : end if
1622 : : end if
1623 : : end if ! mykpt
1624 : :
1625 : 6 : bdtot_index=bdtot_index+nband_k
1626 : :
1627 : : ! End loop on spin,kpt
1628 : : end do ! ikpt
1629 : : end do !isppol
1630 : :
1631 : : !Close file
1632 [ - + - - ]: 3 : if (i_am_master.or.(iomode_etsf_mpiio.and.i_am_master_spfft)) then
1633 [ + - ]: 3 : if (iomode==IO_MODE_ETSF) then
1634 [ - + ]: 3 : NCF_CHECK(nf90_close(ncid))
1635 : : else
1636 : 0 : ierr=close_unit(ount,msg)
1637 [ # # ]: 0 : ABI_CHECK(ierr==0,"Error while closing OPT2 file")
1638 : : end if
1639 : : end if
1640 : :
1641 : : !Datastructures deallocations
1642 [ + + ]: 6 : do itypat=1,dtset%ntypat
1643 : 6 : call destroy_atomorb(atm(itypat))
1644 : : enddo
1645 [ + - + + : 6 : ABI_FREE(atm)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
1646 [ - + ]: 3 : ABI_FREE(psinablapsi)
1647 [ - + ]: 3 : ABI_FREE(psinablapsi2)
1648 [ - + ]: 3 : ABI_FREE(lmn2ln)
1649 [ + + ]: 3 : if (use_spinorbit) then
1650 [ - + ]: 2 : ABI_FREE(psinablapsi_soc)
1651 [ + + ]: 7 : do iatom=1,natom
1652 [ - + ]: 7 : ABI_FREE(phisocphj(iatom)%value)
1653 : : end do
1654 [ + - + + : 7 : ABI_FREE(phisocphj)
- + - + ]
1655 : : end if
1656 [ + - ]: 3 : if (.not.already_has_nabla) then
1657 [ + + ]: 6 : do itypat=1,dtset%ntypat
1658 [ + - ]: 3 : if (allocated(pawtab(itypat)%nabla_ij)) then
1659 : 3 : ABI_FREE(pawtab(itypat)%nabla_ij)
1660 : 3 : pawtab(itypat)%has_nabla=0
1661 : : end if
1662 [ + + ]: 6 : if (allocated(pawtab(itypat)%nabla_im_ij)) then
1663 : 2 : ABI_FREE(pawtab(itypat)%nabla_im_ij)
1664 : : end if
1665 : : end do
1666 : : end if
1667 : :
1668 : : ABI_NVTX_END_RANGE()
1669 : : DBG_EXIT("COLL")
1670 : :
1671 [ - + - + : 6 : end subroutine optics_paw_core
- + - - -
- - + - +
- + - + -
+ - + -
+ ]
1672 : : !!***
1673 : :
1674 : : !----------------------------------------------------------------------
1675 : :
1676 : : !!****f* m_paw_optics/linear_optics_paw
1677 : : !! NAME
1678 : : !! linear_optics_paw
1679 : : !!
1680 : : !! FUNCTION
1681 : : !! This program computes the elements of the optical frequency dependent
1682 : : !! linear susceptiblity using matrix elements <-i Nabla> obtained from a
1683 : : !! PAW ground state calculation. It uses formula 17 from Gadoc et al,
1684 : : !! Phys. Rev. B 73, 045112 (2006) [[cite:Gajdo2006]] together with a scissors correction. It uses
1685 : : !! a Kramers-Kronig transform to compute the real part from the imaginary part, and
1686 : : !! it will work on all types of unit cells. It outputs all tensor elements of
1687 : : !! both the real and imaginary parts.
1688 : : !!
1689 : : !! INPUTS
1690 : : !! filnam: base of file names to read data from
1691 : : !! mpi_enreg: mpi set up variable, not used in this code
1692 : : !!
1693 : : !! OUTPUT
1694 : : !! _real and _imag output files
1695 : : !!
1696 : : !! NOTES
1697 : : !! This routine is not tested
1698 : : !!
1699 : : !! SOURCE
1700 : :
1701 : 0 : subroutine linear_optics_paw(filnam,filnam_out)
1702 : :
1703 : : !Arguments -----------------------------------
1704 : : !scalars
1705 : : character(len=fnlen),intent(in) :: filnam,filnam_out
1706 : :
1707 : : !Local variables-------------------------------
1708 : : integer,parameter :: master=0
1709 : : integer :: iomode,bantot,bdtot_index,fform1,headform
1710 : : integer :: iband,ierr,ii,ikpt,iom,iout,isppol,isym,jband,jj,me,mband
1711 : : integer :: method,mom,nband_k,nkpt,nspinor,nsppol,nsym,occopt,only_check
1712 : : integer :: rdwr,spaceComm,inpunt,reunt,imunt,wfunt
1713 : 0 : integer,allocatable :: nband(:),symrel(:,:,:)
1714 : : real(dp) :: del,dom,fij,gdelta,omin,omax,paijpbij(2),mbpt_sciss,wij,ucvol
1715 : : real(dp) :: diffwp, diffwm
1716 : : 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)
1717 : 0 : real(dp),allocatable :: e1(:,:,:),e2(:,:,:,:),epsilon_tot(:,:,:,:),eigen0(:),eig0_k(:)
1718 : 0 : real(dp),allocatable :: kpts(:,:),occ(:),occ_k(:),oml1(:),wtk(:)
1719 : 0 : complex(dp),allocatable :: eps_work(:)
1720 : : character(len=fnlen) :: filnam1,filnam_gen
1721 : : character(len=500) :: msg
1722 : 0 : type(hdr_type) :: hdr
1723 : : type(wffile_type) :: wff1
1724 : : !arrays
1725 : 0 : real(dp),allocatable :: psinablapsi(:,:,:,:)
1726 : : ! *********************************************************************************
1727 : :
1728 : : DBG_ENTER("COLL")
1729 : : ABI_NVTX_START_RANGE(NVTX_LINEAR_OPTICS_PAW)
1730 : :
1731 : : !write(std_out,'(a)')' Give the name of the output file ...'
1732 : : !read(std_in, '(a)') filnam_out
1733 : : !write(std_out,'(a)')' The name of the output file is :',filnam_out
1734 : :
1735 : : !Read data file
1736 [ # # ]: 0 : if (open_file(filnam,msg,newunit=inpunt,form='formatted') /= 0 ) then
1737 : 0 : ABI_ERROR(msg)
1738 : : end if
1739 : :
1740 : 0 : rewind(inpunt)
1741 : 0 : read(inpunt,*)
1742 : 0 : read(inpunt,'(a)')filnam_gen ! generic name for the files
1743 [ # # ]: 0 : filnam1=trim(filnam_gen)//'_OPT' ! nabla matrix elements file
1744 : :
1745 : : !Open the Wavefunction and optic files
1746 : : !These default values are typical of sequential use
1747 : 0 : iomode=IO_MODE_FORTRAN ; spaceComm=xmpi_comm_self; me=0
1748 : :
1749 : : ! Read the header of the optic files
1750 : 0 : call hdr%from_fname(filnam1, fform1, spaceComm)
1751 : 0 : call hdr%free()
1752 [ # # ]: 0 : if (fform1 /= 610) then
1753 : 0 : ABI_ERROR("Abinit8 requires an OPT file with fform = 610")
1754 : : end if
1755 : :
1756 : : !Open the conducti optic files
1757 : 0 : wfunt = get_unit()
1758 : 0 : call WffOpen(iomode,spaceComm,filnam1,ierr,wff1,master,me,wfunt)
1759 : :
1760 : : !Read the header from Ground state file
1761 : 0 : rdwr=1
1762 : 0 : call hdr_io(fform1,hdr,rdwr,wff1)
1763 : :
1764 : : !Extract info from the header
1765 : 0 : headform=hdr%headform
1766 : 0 : bantot=hdr%bantot
1767 : 0 : nkpt=hdr%nkpt
1768 [ # # # # : 0 : ABI_MALLOC(kpts,(3,nkpt))
# # # # #
# # # ]
1769 [ # # # # : 0 : ABI_MALLOC(wtk,(nkpt))
# # # # #
# # # ]
1770 [ # # # # ]: 0 : kpts(:,:)=hdr%kptns(:,:)
1771 [ # # ]: 0 : wtk(:)=hdr%wtk(:)
1772 : 0 : nspinor=hdr%nspinor
1773 : 0 : nsppol=hdr%nsppol
1774 : 0 : occopt=hdr%occopt
1775 [ # # # # ]: 0 : rprimd(:,:)=hdr%rprimd(:,:)
1776 : 0 : rprimdinv(:,:) = rprimd(:,:)
1777 : 0 : call matrginv(rprimdinv,3,3) ! need the inverse of rprimd to symmetrize the tensors
1778 [ # # # # : 0 : ABI_MALLOC(nband,(nkpt*nsppol))
# # # # #
# # # ]
1779 [ # # # # : 0 : ABI_MALLOC(occ,(bantot))
# # # # #
# # # ]
1780 [ # # ]: 0 : occ(1:bantot)=hdr%occ(1:bantot)
1781 [ # # ]: 0 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
1782 : 0 : nsym=hdr%nsym
1783 [ # # # # : 0 : ABI_MALLOC(symrel,(3,3,nsym))
# # # # #
# # # ]
1784 [ # # # # : 0 : symrel(:,:,:)=hdr%symrel(:,:,:)
# # ]
1785 : :
1786 : : !Get mband, as the maximum value of nband(nkpt)
1787 [ # # ]: 0 : mband=maxval(nband(:))
1788 : :
1789 : : !get ucvol etc.
1790 : 0 : iout = -1
1791 : 0 : call metric(gmet,gprimd,iout,rmet,rprimd,ucvol)
1792 : :
1793 : 0 : write(std_out,*)
1794 : 0 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1)
1795 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2)
1796 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3)
1797 : 0 : write(std_out,*)
1798 : 0 : write(std_out,'(a,3f10.5,a)' )' rprimdinv =',rprimdinv(1:3,1)
1799 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimdinv(1:3,2)
1800 : 0 : write(std_out,'(a,3f10.5,a)' )' ',rprimdinv(1:3,3)
1801 : 0 : write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband
1802 : :
1803 : : !get eigen0
1804 [ # # # # : 0 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
# # # # #
# # # ]
1805 : 0 : read(wfunt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
1806 : :
1807 : 0 : read(inpunt,*)mbpt_sciss
1808 : 0 : read(inpunt,*)dom,omin,omax,mom
1809 : 0 : close(inpunt)
1810 : :
1811 [ # # # # : 0 : ABI_MALLOC(oml1,(mom))
# # # # #
# # # ]
1812 [ # # # # : 0 : ABI_MALLOC(e1,(3,3,mom))
# # # # #
# # # ]
1813 [ # # # # : 0 : ABI_MALLOC(e2,(2,3,3,mom))
# # # # #
# # # ]
1814 [ # # # # : 0 : ABI_MALLOC(epsilon_tot,(2,3,3,mom))
# # # # ]
1815 [ # # # # : 0 : ABI_MALLOC(eps_work,(mom))
# # # # #
# # # ]
1816 : 0 : del=(omax-omin)/(mom-1)
1817 [ # # ]: 0 : do iom=1,mom
1818 : 0 : oml1(iom)=omin+dble(iom-1)*del
1819 : : end do
1820 : 0 : write(std_out,'(a,i8,4f10.5,a)')' npts,omin,omax,width,mbpt_sciss =',mom,omin,omax,dom,mbpt_sciss,' Ha'
1821 : :
1822 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,mband,mband))
# # # # #
# # # # #
# # ]
1823 : :
1824 : : !loop over spin components
1825 [ # # ]: 0 : do isppol=1,nsppol
1826 : : bdtot_index = 0
1827 : : ! loop over k points
1828 [ # # ]: 0 : do ikpt=1,nkpt
1829 : : !
1830 : : ! number of bands for this k point
1831 : 0 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1832 [ # # # # : 0 : ABI_MALLOC(eig0_k,(nband_k))
# # # # #
# # # ]
1833 [ # # # # : 0 : ABI_MALLOC(occ_k,(nband_k))
# # # # ]
1834 : : ! eigenvalues for this k-point
1835 [ # # ]: 0 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
1836 : : ! occupation numbers for this k-point
1837 [ # # ]: 0 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1838 : : ! values of -i*nabla matrix elements for this k point
1839 [ # # # # : 0 : psinablapsi=zero
# # # # ]
1840 [ # # # # : 0 : read(wfunt)((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k)
# # # # ]
1841 [ # # # # : 0 : read(wfunt)((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k)
# # # # ]
1842 [ # # # # : 0 : read(wfunt)((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k)
# # # # ]
1843 : :
1844 : : ! occupation numbers for k-point
1845 [ # # ]: 0 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1846 : : ! accumulate e2 for this k point, Eq. 17 from PRB 73, 045112 (2006) [[cite:Gajdo2006]]
1847 [ # # ]: 0 : do iband = 1, nband_k
1848 [ # # ]: 0 : do jband = 1, nband_k
1849 : 0 : fij = occ_k(iband) - occ_k(jband) !occ number difference
1850 : 0 : wij = eig0_k(iband) - eig0_k(jband) !energy difference
1851 [ # # ]: 0 : if (abs(fij) > zero) then ! only consider states of differing occupation numbers
1852 [ # # ]: 0 : do ii = 1, 3
1853 [ # # ]: 0 : do jj = 1, 3
1854 : : paijpbij(1) = psinablapsi(1,ii,iband,jband)*psinablapsi(1,jj,iband,jband) + &
1855 : 0 : & psinablapsi(2,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
1856 : : paijpbij(2) = psinablapsi(2,ii,iband,jband)*psinablapsi(1,jj,iband,jband) - &
1857 : 0 : & psinablapsi(1,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
1858 [ # # ]: 0 : do iom = 1, mom
1859 : : ! original version
1860 : : ! diffw = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
1861 : : ! gdelta = exp(-diffw*diffw/(4.0*dom*dom))/(2.0*dom*sqrt(pi)) ! delta fnc resolved as Gaussian
1862 : : ! 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))
1863 : : ! 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))
1864 : 0 : diffwm = wij - mbpt_sciss + oml1(iom) ! apply scissors term here
1865 : 0 : diffwp = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
1866 : 0 : gdelta = exp(-diffwp*diffwp/(4.0*dom*dom))/(2.0*dom*sqrt(pi))
1867 : 0 : e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(wij*wij)
1868 : 0 : e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(wij*wij)
1869 : : end do ! end loop over spectral points
1870 : : end do ! end loop over jj = 1, 3
1871 : : end do ! end loop over ii = 1, 3
1872 : : end if ! end selection on fij /= 0
1873 : : end do ! end loop over jband
1874 : : end do ! end loop over iband
1875 : :
1876 [ # # ]: 0 : ABI_FREE(eig0_k)
1877 [ # # ]: 0 : ABI_FREE(occ_k)
1878 : 0 : bdtot_index=bdtot_index+nband_k
1879 : : end do ! end loop over k points
1880 : : end do ! end loop over spin polarizations
1881 : :
1882 : : !here apply nsym symrel transformations to reconstruct full tensor from IBZ part
1883 [ # # # # : 0 : epsilon_tot(:,:,:,:) = zero
# # # # ]
1884 [ # # ]: 0 : do isym = 1, nsym
1885 [ # # # # : 0 : symd(:,:)=matmul(rprimd(:,:),matmul(symrel(:,:,isym),rprimdinv(:,:)))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1886 : 0 : symdinv(:,:)=symd(:,:)
1887 : 0 : call matrginv(symdinv,3,3)
1888 [ # # ]: 0 : do iom = 1, mom
1889 [ # # # # : 0 : e2rot(:,:)=matmul(symdinv(:,:),matmul(e2(1,:,:,iom),symd(:,:)))
# # # # #
# # # # #
# # # # #
# ]
1890 [ # # # # ]: 0 : epsilon_tot(2,:,:,iom) = epsilon_tot(2,:,:,iom)+e2rot(:,:)/nsym
1891 : : end do
1892 : : end do
1893 : :
1894 : : !generate e1 from e2 via KK transforma
1895 : 0 : method=0 ! use naive integration ( = 1 for simpson)
1896 : 0 : only_check=0 ! compute real part of eps in kk routine
1897 [ # # ]: 0 : do ii = 1, 3
1898 [ # # ]: 0 : do jj = 1, 3
1899 [ # # ]: 0 : eps_work(:) = cmplx(0.0,epsilon_tot(2,ii,jj,:), kind=dp)
1900 : 0 : call kramerskronig(mom,oml1,eps_work,method,only_check)
1901 [ # # ]: 0 : epsilon_tot(1,ii,jj,:) = real(eps_work(:))
1902 [ # # # # ]: 0 : if (ii /= jj) epsilon_tot(1,ii,jj,:) = epsilon_tot(1,ii,jj,:)- 1.0
1903 : : end do ! end loop over jj
1904 : : end do ! end loop over ii
1905 : :
1906 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_imag',msg,newunit=reunt,form='formatted') /= 0) then
1907 : 0 : ABI_ERROR(msg)
1908 : : end if
1909 : :
1910 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_real',msg,unit=imunt,form='formatted') /= 0) then
1911 : 0 : ABI_ERROR(msg)
1912 : : end if
1913 : :
1914 : 0 : write(reunt,'(a12,6a13)')' # Energy/Ha ','eps_2_xx','eps_2_yy','eps_2_zz',&
1915 : 0 : & 'eps_2_yz','eps_2_xz','eps_2_xy'
1916 : 0 : write(imunt,'(a12,6a13)')' # Energy/Ha ','eps_1_xx','eps_1_yy','eps_1_zz',&
1917 : 0 : & 'eps_1_yz','eps_1_xz','eps_1_xy'
1918 : :
1919 [ # # ]: 0 : do iom = 1, mom
1920 : 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),' ',&
1921 : 0 : & epsilon_tot(2,1,1,iom),' ',epsilon_tot(2,2,2,iom),' ',epsilon_tot(2,3,3,iom),' ',&
1922 : 0 : & epsilon_tot(2,2,3,iom),' ',epsilon_tot(2,1,3,iom),' ',epsilon_tot(2,1,2,iom)
1923 : 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),' ',&
1924 : 0 : & epsilon_tot(1,1,1,iom),' ',epsilon_tot(1,2,2,iom),' ',epsilon_tot(1,3,3,iom),' ',&
1925 : 0 : & epsilon_tot(1,2,3,iom),' ',epsilon_tot(1,1,3,iom),' ',epsilon_tot(1,1,2,iom)
1926 : : end do
1927 : :
1928 : 0 : close(reunt)
1929 : 0 : close(imunt)
1930 : :
1931 [ # # ]: 0 : ABI_FREE(nband)
1932 [ # # ]: 0 : ABI_FREE(oml1)
1933 [ # # ]: 0 : ABI_FREE(e2)
1934 [ # # ]: 0 : ABI_FREE(e1)
1935 [ # # ]: 0 : ABI_FREE(occ)
1936 [ # # ]: 0 : ABI_FREE(psinablapsi)
1937 [ # # ]: 0 : ABI_FREE(eigen0)
1938 [ # # ]: 0 : ABI_FREE(wtk)
1939 [ # # ]: 0 : ABI_FREE(kpts)
1940 : :
1941 : 0 : call hdr%free()
1942 : :
1943 : : ABI_NVTX_END_RANGE()
1944 : : DBG_EXIT("COLL")
1945 : :
1946 [ # # # # : 0 : end subroutine linear_optics_paw
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
1947 : : !!***
1948 : :
1949 : : !----------------------------------------------------------------------
1950 : :
1951 : : !!****f* m_paw_optics/pawnabla_soc_init
1952 : : !! NAME
1953 : : !! pawnabla_soc_init
1954 : : !!
1955 : : !! FUNCTION
1956 : : !! Compute the PAW SOC contribution(s) to the momentum PAW matrix elements,
1957 : : !! i.e. <Phi_i|1/4 Alpha^2 dV/dr (Sigma X vec(r)/r) |Phi_j>
1958 : : !! where:
1959 : : !! {Phi_i}= AE partial waves
1960 : : !! Alpha = inverse of fine structure constant
1961 : : !! Sigma^alpha,beta= Pauli matrices
1962 : : !! X = cross product
1963 : : !!
1964 : : !! There are 2 typical uses:
1965 : : !! - Valence-valence terms: Phi_i and Phi_j are PAW AE partial waves (unpolarized)
1966 : : !! - Core-valence terms: Phi_i is are AE partial waves and Phi_j are core AE wave-functions (spinors)
1967 : : !!
1968 : : !! In practice we compute:
1969 : : !! (Sigma^up-up X g_ij) and (Sigma^up-dn X g_ij) (X = vector cross product)
1970 : : !! (Sigma^dn-up X g_ij) and (Sigma^dn-dn X g_ij)
1971 : : !! where:
1972 : : !! g_ij= 1/4 Alpha^2 Int_[Phi_i(r)/r Phi_j(r)/r dV(r)/dr r^2 dr] . Gvec_ij
1973 : : !! and Gvec_ij= Int[S_limi S_ljmj vec(r)/r dOmega] (Gaunt coefficients)
1974 : : !!
1975 : : !! COPYRIGHT
1976 : : !! Copyright (C) 2021-2026 ABINIT group (NBrouwer,MT)
1977 : : !! This file is distributed under the terms of the
1978 : : !! GNU General Public License, see ~ABINIT/COPYING
1979 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1980 : : !!
1981 : : !! INPUTS
1982 : : !! el_temp=electronic temperature (hartree)
1983 : : !! ixc= choice of exchange-correlation scheme (see above, and below)
1984 : : !! my_natom=number of atoms treated by current processor
1985 : : !! natom=total number of atoms in cell
1986 : : !! nspden=number of spin-density components
1987 : : !! ntypat=number of types of atoms in unit cell.
1988 : : !! option_core=Type of calculation: 0=valence-valence, 1=core-valence
1989 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
1990 : : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
1991 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
1992 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
1993 : : !! pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
1994 : : !! spnorbscl=scaling factor for spin-orbit coupling
1995 : : !! typat(natom) =Type of each atoms
1996 : : !! xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
1997 : : !! xc_taupos= lowest allowed kinetic energy density (for mGGA XC functionals)
1998 : : !! znucl(ntypat)=gives the nuclear charge for all types of atoms
1999 : : !! [atm <type(paw_atomorb_type)>]=--optional-- structure containing core info
2000 : : !! [mpi_atmtab(:)]=--optional-- indexes of the atoms treated by current proc
2001 : : !! [comm_atom]=--optional-- MPI communicator over atoms
2002 : : !!
2003 : : !! OUTPUT
2004 : : ! phisocphj(dtset%natom) <type(coeff4_type)>= stores soc coefficients:
2005 : : !!! If option_core==0 or nspinor_cor==1:
2006 : : !! phisocphj(iat)%value(1,1,idir,ilmn,jlmn) is real part of (Sigma^up-up X g_ij)
2007 : : !! phisocphj(iat)%value(2,1,idir,ilmn,jlmn) is imaginary part of (Sigma^up-up X g_ij)
2008 : : !! phisocphj(iat)%value(1,2,idir,ilmn,jlmn) is real part of (Sigma^up-dn X g_ij)
2009 : : !! phisocphj(iat)%value(2,2,idir,ilmn,jlmn) is imaginary part of (Sigma^up-dn X g_ij)
2010 : : !! If option_core==1 and nspinor_cor==2 (core-valence with spinorial core WF):
2011 : : !! phisocphj(iat)%value(1,1,idir,ilmn,2*jlmn-1) is real part of (Sigma^up-up X g_ij^up)
2012 : : !! phisocphj(iat)%value(2,1,idir,ilmn,2*jlmn-1) is imaginary part of (Sigma^up-up X g_ij^up)
2013 : : !! phisocphj(iat)%value(1,2,idir,ilmn,2*jlmn-1) is real part of (Sigma^dn-up X g_ij^dn)
2014 : : !! phisocphj(iat)%value(2,2,idir,ilmn,2*jlmn-1) is imaginary part of (Sigma^dn-up X g_ij^dn)
2015 : : !! phisocphj(iat)%value(1,1,idir,ilmn,2*jlmn ) is real part of (Sigma^up-dn X g_ij^up)
2016 : : !! phisocphj(iat)%value(2,1,idir,ilmn,2*jlmn ) is imaginary part of (Sigma^up-dn X g_ij^up)
2017 : : !! phisocphj(iat)%value(1,2,idir,ilmn,2*jlmn ) is real part of (Sigma^dn-dn X g_ij^dn)
2018 : : !! phisocphj(iat)%value(2,2,idir,ilmn,2*jlmn ) is imaginary part of (Sigma^dn-dn X g_ij^dn)
2019 : : !! (idir=cartesian direction)
2020 : : !!
2021 : : !! SIDE EFFECTS
2022 : : !!
2023 : : !! NOTES
2024 : : !! If Phi_j is not polarized,
2025 : : !! (Sigma^dn-dn X g_ij)=-(Sigma^up-up X g_ij)
2026 : : !! (Sigma^dn-up X g_ij)= (Sigma^up-dn X g_ij)^*
2027 : : !! So, we store only 2 components.
2028 : : !! If Phi_j is polarized, the spin component is included in the last dimension of
2029 : : !! phisocphj(iat)%value, i.e. lmn_size_cor=2*lmn_size
2030 : : !!
2031 : : !! SOURCE
2032 : :
2033 : 3 : subroutine pawnabla_soc_init(el_temp,phisocphj,option_core,ixc,my_natom,natom,nspden,ntypat,pawang, &
2034 : 6 : & pawrad,pawrhoij,pawtab,pawxcdev,spnorbscl,typat,xc_denpos,xc_taupos,znucl, &
2035 [ + - - + : 3 : & atm,mpi_atmtab,comm_atom) ! Optional arguments
- - ]
2036 : :
2037 : : !Arguments ------------------------------------
2038 : : !scalars
2039 : : integer,intent(in) :: ixc,my_natom,natom,nspden,ntypat,option_core,pawxcdev
2040 : : integer,optional,intent(in) :: comm_atom
2041 : : real(dp),intent(in) :: el_temp,spnorbscl,xc_denpos,xc_taupos
2042 : : type(pawang_type),intent(in) :: pawang
2043 : : !arrays
2044 : : integer,intent(in) :: typat(natom)
2045 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
2046 : : real(dp),intent(in) :: znucl(ntypat)
2047 : : type(coeff5_type),allocatable,target,intent(inout) :: phisocphj(:)
2048 : : type(pawrad_type),target,intent(in) :: pawrad(ntypat)
2049 : : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
2050 : : type(pawtab_type),target,intent(in) :: pawtab(ntypat)
2051 : : type(atomorb_type), intent(in), target, optional :: atm(ntypat)
2052 : :
2053 : : !Local variables-------------------------------
2054 : : !scalars
2055 : : real(dp),parameter :: one_over_fourpi = one/sqrt(four_pi)
2056 : : real(dp),parameter :: sqr_fourpi_over_3 = sqrt(four_pi/3)
2057 : : real(dp),parameter :: QuarterFineStruct2=(half/InvFineStruct)**2
2058 : : real(dp),parameter :: hyb_mixing_ = 0.0_dp ! Fake value to be updated in the future
2059 : : integer :: iatom,iatom_tot,itypat,ii,jj,ierr,ipts,ignt,sgnkappa
2060 : : integer :: idum,option,usenhat,usekden,usecore,xclevel,nkxc,my_comm_atom
2061 : : integer :: mesh_size,lmn_size,lmn2_size,lmn_size_j,lmn_size_cor
2062 : : integer :: lm_size,ln_size,ln_size_j
2063 : : integer :: ilmn,ilm,iln,jl,jm,jm_re,jm_im,jlmn,jlm,jlm_re,jlm_im,jln,js,klm_re,klm_im
2064 : : logical :: my_atmtab_allocated,paral_atom
2065 : : real(dp) :: avg,cgc,compch_sph_dum,eexc_dum,ssxc_dum,eexcdc_dum,jmj
2066 : : real(dp) :: fact_re,fact_im,gx_re,gx_im,gy_re,gy_im,gz_re,gz_im,if3
2067 : : character(len=500) :: msg
2068 : : !arrays
2069 : 3 : integer,pointer :: my_atmtab(:),indlmn(:,:),indlmn_j(:,:)
2070 : 3 : logical,allocatable :: lmselect(:)
2071 : : real(dp) :: nhat_dum(1,1,1),trho_dum(1,1,1),kxc_dum(1,1,1),k3xc_dum(1,1,1)
2072 : 3 : real(dp),allocatable :: rho1(:,:,:),tau1(:,:,:),rhosph(:),vhartree(:),vxc(:,:,:)
2073 : 3 : real(dp),allocatable :: intf3(:,:),potsph(:),dVdr(:),func(:)
2074 : 3 : real(dp),pointer :: phi_j(:,:),soc_ij(:,:,:,:,:)
2075 : : type(pawrad_type),pointer :: pawrd
2076 : : type(pawtab_type),pointer :: pawtb
2077 : :
2078 : : ! ************************************************************************
2079 : : ABI_NVTX_START_RANGE(NVTX_PAWNABLA_SOC_INIT)
2080 : :
2081 : : !Some checks in case of core-valence (option_core=1)
2082 : 3 : if (option_core/=0.and.option_core/=1) then
2083 : 0 : msg='Wrong option_core value!'
2084 : 0 : ABI_BUG(msg)
2085 : : end if
2086 [ + + ]: 3 : if (option_core==1) then
2087 : : ! Check if we have the optional arguments
2088 [ - + ]: 2 : if (.not.present(atm)) then
2089 : 0 : msg='For core-valence calculation, need atm data!'
2090 : 0 : ABI_BUG(msg)
2091 : : end if
2092 [ + + ]: 4 : do itypat=1,ntypat
2093 : : ! Check if we have relativistic core wave functions
2094 [ - + ]: 4 : if (size(atm(itypat)%indlmn,1)<8) then
2095 : 0 : write(msg,'(a)') 'Wrong 1st dim. of indlmn_cor in pawnabla_soc_init (need spinors)!'
2096 : 0 : ABI_BUG(msg)
2097 : : end if
2098 : : enddo
2099 : : endif
2100 : :
2101 : : !Some useful variables
2102 : 3 : usekden=pawxc_get_usekden(ixc)
2103 : 3 : usecore=1 ; nkxc=0 ; usenhat=0
2104 : 3 : xclevel=pawxc_get_xclevel(ixc)
2105 [ + + ]: 3 : if (option_core==1) then
2106 : : lmn_size_cor=0
2107 [ + + ]: 4 : do itypat=1,ntypat
2108 : 4 : lmn_size_cor=max(lmn_size_cor,atm(itypat)%lmn_size) !Includes spinors
2109 : : enddo
2110 : : end if
2111 : :
2112 : : !Prepare output arrays
2113 [ - + ]: 3 : if (allocated(phisocphj)) then
2114 [ # # ]: 0 : do iatom=1,natom
2115 [ # # ]: 0 : if (allocated(phisocphj(iatom)%value)) then
2116 : 0 : ABI_FREE(phisocphj(iatom)%value)
2117 : : end if
2118 : : end do
2119 [ # # # # : 0 : ABI_FREE(phisocphj)
# # # # ]
2120 : : end if
2121 [ - + + - : 18 : ABI_MALLOC(phisocphj,(natom))
+ - + - -
+ - + - +
+ + ]
2122 [ + + ]: 12 : do iatom=1,natom
2123 : 9 : lmn_size=pawtab(typat(iatom))%lmn_size
2124 [ + + ]: 9 : if (option_core==0) then
2125 [ - + + - : 16 : ABI_MALLOC(phisocphj(iatom)%value,(2,2,3,lmn_size,lmn_size))
- + + - +
- + - - +
- + - + ]
2126 : : else
2127 : : !lmn_size_cor is double because it contains the spin component
2128 [ - + + - : 20 : ABI_MALLOC(phisocphj(iatom)%value,(2,2,3,lmn_size,lmn_size_cor))
- + + - +
- + - - +
- + - + ]
2129 : : end if
2130 [ + + + + : 27776 : phisocphj(iatom)%value=zero
+ + + + +
+ ]
2131 : : end do
2132 : :
2133 : : !Set up parallelism over atoms
2134 [ + - + - ]: 3 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2135 [ - + ]: 3 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2136 [ + - ]: 3 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2137 : : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
2138 : 3 : & my_natom_ref=my_natom)
2139 : :
2140 : : !-------------------------------------------------------------------
2141 : : !Loop over atoms
2142 : :
2143 [ + + ]: 12 : do iatom=1,my_natom
2144 : :
2145 : : ! Atom-dependent data
2146 : 9 : itypat=typat(iatom)
2147 : 9 : pawrd => pawrad(itypat)
2148 : 9 : pawtb => pawtab(itypat)
2149 : 9 : lmn_size=pawtb%lmn_size
2150 : 9 : lmn2_size=pawtb%lmn2_size
2151 : 9 : lm_size=pawtb%lcut_size**2
2152 : 9 : ln_size=pawtb%basis_size
2153 : 9 : mesh_size=pawtb%mesh_size
2154 : 9 : indlmn => pawtb%indlmn
2155 [ - + + - : 27 : ABI_MALLOC(lmselect,(lm_size))
+ - + - -
+ - + ]
2156 [ + + ]: 106 : lmselect(:)=.true.
2157 : :
2158 : : ! Distinguish valence-valence and core-valence cases
2159 [ + + ]: 9 : if (option_core==0) then
2160 : 4 : ln_size_j=ln_size
2161 : 4 : lmn_size_j=lmn_size
2162 : 4 : indlmn_j => pawtb%indlmn
2163 : 4 : phi_j => pawtb%phi
2164 : : else
2165 : 5 : ln_size_j=atm(itypat)%ln_size
2166 : 5 : lmn_size_j=atm(itypat)%lmn_size
2167 : 5 : indlmn_j => atm(itypat)%indlmn
2168 : 5 : phi_j => atm(itypat)%phi(:,:,1)
2169 [ + - ]: 5 : if (atm(itypat)%mesh_size<mesh_size) then
2170 : 0 : msg='mesh_size and mesh_size_cor not compatible!'
2171 : 0 : ABI_BUG(msg)
2172 : : end if
2173 : : endif
2174 : :
2175 : : ! Manage parallelism
2176 [ - + ]: 9 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
2177 : 9 : soc_ij => phisocphj(iatom_tot)%value(:,:,:,:,:)
2178 : :
2179 : : !-------------------------------------------------------------------
2180 : : !Compute all-electron density
2181 : :
2182 [ - + + - : 45 : ABI_MALLOC(rho1,(mesh_size,lm_size,nspden))
- + + - -
+ + - + -
+ - - + -
+ ]
2183 [ + + + + : 107505 : option=2 ; rho1=zero
+ + ]
2184 : : call pawdensities(compch_sph_dum,1,iatom_tot,lmselect,lmselect,lm_size,&
2185 : : & nhat_dum,nspden,-1,0,option,-1,0,pawang,0,pawrd,&
2186 : 9 : & pawrhoij(iatom),pawtb,rho1,trho_dum)
2187 [ - + ]: 9 : if (usekden==1) then
2188 [ # # # # : 0 : ABI_MALLOC(tau1,(mesh_size,lm_size,nspden))
# # # # #
# # # # #
# # # # ]
2189 [ # # # # : 0 : tau1=zero
# # ]
2190 : : call pawkindensities(1,lmselect,lm_size,nspden,-1,option,-1,&
2191 : 0 : & pawang,pawrd,pawrhoij(iatom),pawtb,tau1,trho_dum)
2192 : : end if
2193 : :
2194 : : !-------------------------------------------------------------------
2195 : : !Compute spherical potential and compute its first derivative dV/dr
2196 : :
2197 [ - + + - : 27 : ABI_MALLOC(potsph,(mesh_size))
+ - + - -
+ - + -
+ ]
2198 [ - + + - : 18 : ABI_MALLOC(dVdr,(mesh_size))
- + - + -
+ ]
2199 [ + + + + ]: 18789 : potsph=zero ; dVdr=zero
2200 : :
2201 : : ! Compute XC potential
2202 : 9 : option=1
2203 [ + - ]: 9 : if (pawxcdev/=0) then
2204 [ - + + - : 36 : ABI_MALLOC(vxc,(mesh_size,lm_size,nspden))
- + + - -
+ + - - +
- + ]
2205 [ + + + + : 107505 : vxc=zero
+ + ]
2206 : : call pawxcm(pawtb%coredens,eexc_dum,eexcdc_dum,ssxc_dum,idum,hyb_mixing_,ixc,kxc_dum,lm_size,&
2207 : : & lmselect,nhat_dum,nkxc,.false.,mesh_size,nspden,option,pawang,pawrd,&
2208 : 9 : & pawxcdev,rho1,usecore,usenhat,vxc,xclevel,xc_denpos,el_temp)
2209 [ + + ]: 9399 : potsph(1:mesh_size)=half*(vxc(1:mesh_size,1,1)+vxc(1:mesh_size,1,nspden))
2210 : : else
2211 [ # # # # : 0 : ABI_MALLOC(vxc,(mesh_size,pawang%angl_size,nspden))
# # # # #
# # # # #
# # # # #
# ]
2212 [ # # # # : 0 : vxc=zero
# # ]
2213 : : call pawxc(pawtb%coredens,eexc_dum,eexcdc_dum,ssxc_dum,hyb_mixing_,ixc,kxc_dum,k3xc_dum,&
2214 : : & lm_size,lmselect,nhat_dum,nkxc,nkxc,.false.,mesh_size,nspden,option,pawang,&
2215 : : & pawrd,rho1,usecore,usenhat,vxc,xclevel,xc_denpos,el_temp,&
2216 : 0 : & coretau=pawtb%coretau,taur=tau1,xc_taupos=xc_taupos)
2217 [ # # ]: 0 : potsph(1:mesh_size)=zero
2218 [ # # ]: 0 : do ipts=1,pawang%angl_size
2219 : : potsph(1:mesh_size)=potsph(1:mesh_size) &
2220 [ # # ]: 0 : & +half*(vxc(1:mesh_size,ipts,1)+vxc(1:mesh_size,ipts,nspden))*pawang%angwgth(ipts)
2221 : : end do
2222 [ # # ]: 0 : potsph(1:mesh_size)=sqrt(four_pi)*potsph(1:mesh_size)
2223 : : end if
2224 : :
2225 : : ! Compute Hartree potentialHalfFineStruct2
2226 [ - + + - : 18 : ABI_MALLOC(vhartree,(mesh_size))
- + - + -
+ ]
2227 [ - + + - : 18 : ABI_MALLOC(rhosph,(mesh_size))
- + - + -
+ ]
2228 [ + + + + ]: 18789 : vhartree=zero ; rhosph=zero
2229 [ + + ]: 9399 : rhosph(1:mesh_size)=rho1(1:mesh_size,1,1)
2230 [ + + ]: 9399 : if (usecore==1) rhosph(1:mesh_size)=rhosph(1:mesh_size)+sqrt(four_pi)*pawtb%coredens(1:mesh_size)
2231 [ + + ]: 9399 : rhosph(1:mesh_size)=rhosph(1:mesh_size)*four_pi*pawrd%rad(1:mesh_size)**2
2232 : 9 : call poisson(rhosph,0,pawrd,vhartree)
2233 [ + + ]: 9390 : vhartree(2:mesh_size)=(vhartree(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrd%rad(2:mesh_size)
2234 : 9 : call pawrad_deducer0(vhartree,mesh_size,pawrd)
2235 [ + + ]: 9399 : potsph(1:mesh_size)=potsph(1:mesh_size)+vhartree(1:mesh_size)
2236 : :
2237 : : ! Apply angular scaling factor
2238 [ + + ]: 9399 : potsph(1:mesh_size)=one_over_fourpi*potsph(1:mesh_size)
2239 : :
2240 : : ! Compute 1st derivative of potential
2241 : 9 : call nderiv_gen(dVdr,potsph,pawrd)
2242 : :
2243 : : ! Multiply by relativistic factor
2244 [ + + ]: 9399 : dVdr(1:mesh_size)=dVdr(1:mesh_size)*(one/(one-potsph(1:mesh_size)*half/InvFineStruct**2)**2)
2245 : :
2246 [ - + ]: 9 : ABI_FREE(vxc)
2247 [ - + ]: 9 : ABI_FREE(vhartree)
2248 [ - + ]: 9 : ABI_FREE(potsph)
2249 [ - + ]: 9 : ABI_FREE(rhosph)
2250 [ - + ]: 9 : ABI_FREE(lmselect)
2251 [ - + ]: 9 : ABI_FREE(rho1)
2252 [ - + ]: 9 : if (usekden==1) then
2253 [ # # ]: 0 : ABI_FREE(tau1)
2254 : : end if
2255 : :
2256 : : !-------------------------------------------------------------------
2257 : : !Compute radial and angular contributions
2258 : :
2259 [ - + + - : 36 : ABI_MALLOC(intf3,(ln_size,ln_size_j))
- + + - +
- + - - +
- + ]
2260 [ + + + + ]: 197 : intf3=zero
2261 : :
2262 : : ! >>>> Calculate f_3= alpha^2/4 int[dr ui*(r) uj(r) dV/dr]
2263 [ - + + - : 18 : ABI_MALLOC(func,(mesh_size))
- + - + -
+ ]
2264 [ + + ]: 45 : do jln=1,ln_size_j
2265 [ + + ]: 197 : do iln=1,ln_size
2266 [ + + ]: 161832 : func(1:mesh_size)=dVdr(1:mesh_size)*pawtb%phi(1:mesh_size,iln)*phi_j(1:mesh_size,jln)
2267 : 188 : call simp_gen(intf3(iln,jln),func,pawrd)
2268 : : end do
2269 : : end do
2270 [ + + + + ]: 197 : intf3(:,:)=QuarterFineStruct2*spnorbscl*intf3(:,:)
2271 [ - + ]: 9 : ABI_FREE(func)
2272 : :
2273 : : ! Loop over initial states (valence or core, according to option_core)
2274 [ + + ]: 141 : do jlmn=1,lmn_size_j
2275 : 132 : jl=indlmn_j(1,jlmn)
2276 : 132 : jm=indlmn_j(2,jlmn)
2277 : 132 : jlm=indlmn_j(4,jlmn)
2278 : 132 : jln=indlmn_j(5,jlmn)
2279 : :
2280 : : ! In case of spinorial core wave function, we have to handle imaginary
2281 : : ! spherical harmonics as linear combination of real spherical harmonics
2282 : : ! See Brouwer et al, CPC 266, 108029 (2021), equation (21)
2283 [ + + ]: 132 : if (option_core==1) then
2284 : 100 : jm_re= abs(jm)
2285 : 100 : jm_im=-abs(jm)
2286 : 100 : jlm_re=jl*(jl+1)+jm_re+1
2287 : 100 : jlm_im=jl*(jl+1)+jm_im+1
2288 : : ! Calculate spinor dependent coefficients
2289 : 100 : sgnkappa=indlmn_j(3,jlmn) !sign of kappa
2290 : 100 : jmj=half*indlmn_j(8,jlmn) !2mj is stored in indlmn_cor
2291 : 100 : js=indlmn_j(6,jlmn) !1 is up, 2 is down
2292 [ + + ]: 100 : if (sgnkappa==1) then
2293 [ + + ]: 20 : if(js==1) then
2294 : 10 : cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
2295 : : else
2296 : 10 : cgc=-sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
2297 : : endif
2298 : : else
2299 [ + + ]: 80 : if(js==1) then
2300 : 25 : cgc= sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
2301 : : else
2302 : 55 : cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
2303 : : endif
2304 : : endif
2305 : :
2306 : : ! Calculate factors to convert from complex to real sph. harm.
2307 [ + + ]: 100 : if (jm<0) then
2308 : 45 : fact_re=sqr_fourpi_over_3*half_sqrt2*cgc
2309 : 45 : fact_im=-fact_re
2310 [ + + ]: 55 : else if (jm>0) then
2311 : 15 : fact_re=sqr_fourpi_over_3*half_sqrt2*cgc*(-1)**jm
2312 : 15 : fact_im=fact_re
2313 : : else
2314 : : fact_re=sqr_fourpi_over_3
2315 : : fact_im=0
2316 : : end if
2317 : : else ! valence-valence case (real)
2318 : : js=1
2319 : : jlm_re=jlm ; jlm_im=jlm
2320 : : fact_re=sqr_fourpi_over_3
2321 : : fact_im=0
2322 : : end if
2323 : :
2324 : : ! Loop over final states
2325 [ + + ]: 1397 : do ilmn=1,lmn_size
2326 : 1256 : ilm=indlmn(4,ilmn)
2327 : 1256 : iln=indlmn(5,ilmn)
2328 : :
2329 : : ! >>>> Calculate g_ij=(g_x,g_y,g_z) = sqrt(4pi/3) int dOmega Ylm Ylm' S1-1,0,1
2330 : : ! using real Gaunt coefficients
2331 : 1256 : gx_re=zero;gy_re=zero;gz_re=zero
2332 : 1256 : gx_im=zero;gy_im=zero;gz_im=zero
2333 : 1256 : if3=zero
2334 : :
2335 : : ! jl was set as a flag for invalid combinations
2336 : : ! i.e. m=-(l+1) or m=(l+1)
2337 : : ! In these cases, cgc=0 ; so gx=gy=gz=0
2338 [ + + ]: 1256 : if (jl/=-1) then
2339 : 956 : if3=intf3(iln,jln)
2340 [ + + ]: 956 : klm_re=merge((jlm_re*(jlm_re-1))/2+ilm,(ilm*(ilm-1))/2+jlm_re,ilm<=jlm_re)
2341 : :
2342 : : ! Real parts
2343 : : ! M=-1
2344 : 956 : ignt=pawang%gntselect(2,klm_re) !get index for L=1 M =-1 ilm jlm_re
2345 [ + + ]: 956 : if (ignt/=0) gy_re=fact_re*pawang%realgnt(ignt)
2346 : : ! M=0
2347 : 956 : ignt=pawang%gntselect(3,klm_re) !get index for L=1 M = 0 ilm jlm_re
2348 [ + + ]: 956 : if (ignt/=0) gz_re=fact_re*pawang%realgnt(ignt)
2349 : : ! M=1
2350 : 956 : ignt=pawang%gntselect(4,klm_re) !get index for L=1 M = 1 ilm jlm_re
2351 [ + + ]: 956 : if (ignt/=0) gx_re=fact_re*pawang%realgnt(ignt)
2352 : :
2353 : : ! Imaginary parts
2354 [ + + ]: 956 : if (option_core==1) then
2355 [ + + ]: 700 : klm_im=merge((jlm_im*(jlm_im-1))/2+ilm,(ilm*(ilm-1))/2+jlm_im,ilm<=jlm_im)
2356 : : ! M=-1
2357 : 700 : ignt=pawang%gntselect(2,klm_im) !get index for L=1 M =-1 ilm jlm_im
2358 [ + + ]: 700 : if (ignt/=0) gy_im=fact_im*pawang%realgnt(ignt)
2359 : : ! M=0
2360 : 700 : ignt=pawang%gntselect(3,klm_im) !get index for L=1 M = 0 ilm jlm_im
2361 [ + + ]: 700 : if (ignt/=0) gz_im=fact_im*pawang%realgnt(ignt)
2362 : : ! M=1
2363 : 700 : ignt=pawang%gntselect(4,klm_im) !get index for L=1 M = 1 ilm jlm_im
2364 [ + + ]: 700 : if (ignt/=0) gx_im=fact_im*pawang%realgnt(ignt)
2365 : : end if
2366 : : end if
2367 : :
2368 : : ! >>>> Calculate Sigma X g_ij
2369 : :
2370 [ + + + + ]: 1388 : if (option_core==0.or.js==1) then
2371 : : !(Sigma^up-up X gij)_x = -gy*f_3
2372 : 606 : soc_ij(1,1,1,ilmn,jlmn)=-if3*gy_re ! real part
2373 : 606 : soc_ij(2,1,1,ilmn,jlmn)=-if3*gy_im ! imag part
2374 : : !(Sigma^up-up X gij)_y = gx*f_3
2375 : 606 : soc_ij(1,1,2,ilmn,jlmn)= if3*gx_re ! real part
2376 : 606 : soc_ij(2,1,2,ilmn,jlmn)= if3*gx_im ! imag part
2377 : : !(Sigma^up-up X gij)_z = 0
2378 : 606 : soc_ij(1,1,3,ilmn,jlmn)= zero ! real part
2379 : 606 : soc_ij(2,1,3,ilmn,jlmn)= zero ! imag part
2380 : :
2381 : : !(Sigma^dn-up X gij)_x = i.gz*f_3
2382 : 606 : soc_ij(1,2,1,ilmn,jlmn)=-if3*gz_im ! real part
2383 : 606 : soc_ij(2,2,1,ilmn,jlmn)= if3*gz_re ! imag part
2384 : : !(Sigma^dn-up X gij)_y = -gz*f_3
2385 : 606 : soc_ij(1,2,2,ilmn,jlmn)=-if3*gz_re ! real part
2386 : 606 : soc_ij(2,2,2,ilmn,jlmn)=-if3*gz_im ! imag part
2387 : : !(Sigma^dn-up X gij)_z = (gy-i.gx)*f_3
2388 : 606 : soc_ij(1,2,3,ilmn,jlmn)= if3*(gy_re+gx_im) ! real part
2389 : 606 : soc_ij(2,2,3,ilmn,jlmn)= if3*(gy_im-gx_re) ! imag part
2390 : :
2391 [ + - + + ]: 650 : else if (option_core==1.and.js==2) then
2392 : : !(Sigma^up-dn X gij^dn)_x = -i.gz^dn*f_3
2393 : 350 : soc_ij(1,1,1,ilmn,jlmn)= if3*gz_im ! real part
2394 : 350 : soc_ij(2,1,1,ilmn,jlmn)=-if3*gz_re ! imag part
2395 : : !(Sigma^up-dn X gij^dn)_y = -gz^dn*f_3
2396 : 350 : soc_ij(1,1,2,ilmn,jlmn)=-if3*gz_re ! real part
2397 : 350 : soc_ij(2,1,2,ilmn,jlmn)=-if3*gz_im ! imag part
2398 : : !(Sigma^up-dn X gij^dn)_z = (gy^dn+i.gx^dn)*f_3
2399 : 350 : soc_ij(1,1,3,ilmn,jlmn)= if3*(gy_re-gx_im) ! real part
2400 : 350 : soc_ij(2,1,3,ilmn,jlmn)= if3*(gy_im+gx_re) ! imag part
2401 : :
2402 : : !(Sigma^dn-dn X gij^dn)_x = gy^dn*f_3
2403 : 350 : soc_ij(1,2,1,ilmn,jlmn)= if3*gy_re ! real part
2404 : 350 : soc_ij(2,2,1,ilmn,jlmn)= if3*gy_im ! imag part
2405 : : !(Sigma^dn-dn X gij^dn)_y = -gx^dn*f_3
2406 : 350 : soc_ij(1,2,2,ilmn,jlmn)=-if3*gx_re ! real part
2407 : 350 : soc_ij(2,2,2,ilmn,jlmn)=-if3*gx_im ! imag part
2408 : : !(Sigma^dn-dn X gij^dn)_z = 0
2409 : 350 : soc_ij(1,2,3,ilmn,jlmn)= zero ! real part
2410 : 350 : soc_ij(2,2,3,ilmn,jlmn)= zero ! imag part
2411 : : end if
2412 : :
2413 : : end do ! ilmn
2414 : : end do ! jlmn
2415 : :
2416 : : ! Symetrization
2417 [ + + + - ]: 9 : if (option_core==0.and.lmn_size>1) then
2418 [ + + ]: 32 : do jlmn=2,lmn_size
2419 [ + + ]: 144 : do ilmn=1,jlmn-1
2420 [ + + ]: 476 : do ii=1,3
2421 [ + + ]: 1120 : do jj=1,2
2422 : 672 : avg=half*(soc_ij(1,jj,ii,ilmn,jlmn)+soc_ij(1,jj,ii,jlmn,ilmn))
2423 : 672 : soc_ij(1,jj,ii,ilmn,jlmn)=avg ; soc_ij(1,jj,ii,jlmn,ilmn)=avg
2424 : 672 : avg=half*(soc_ij(2,jj,ii,ilmn,jlmn)+soc_ij(2,jj,ii,ilmn,jlmn))
2425 : 1008 : soc_ij(2,jj,ii,ilmn,jlmn)=avg ; soc_ij(2,jj,ii,jlmn,ilmn)=avg
2426 : : end do
2427 : : end do
2428 : : end do
2429 : : end do
2430 : : end if
2431 : :
2432 [ - + ]: 9 : ABI_FREE(dVdr)
2433 [ - + ]: 12 : ABI_FREE(intf3)
2434 : :
2435 : : end do ! iatom
2436 : :
2437 : : !Reduction in case of parallelism
2438 [ - + ]: 3 : if (paral_atom) then
2439 : 0 : call xmpi_sum(phisocphj,my_comm_atom,ierr)
2440 : : end if
2441 : :
2442 : : !Destroy atom table(s) used for parallelism
2443 : 3 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2444 : : ABI_NVTX_END_RANGE()
2445 : :
2446 [ - + - + : 6 : end subroutine pawnabla_soc_init
- + - + -
+ - + -
+ ]
2447 : : !!***
2448 : :
2449 : : !----------------------------------------------------------------------
2450 : :
2451 : 0 : END MODULE m_paw_optics
2452 : : !!***
|
