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