Line data Source code
1 : !!****m* ABINIT/m_optics_vloc
2 : !! NAME
3 : !! m_optics_vloc
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2010-2026 ABINIT group (SM,VR,FJ,MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_optics_vloc
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_wffile
28 : use m_xmpi
29 : use m_hdr
30 : use m_dtset
31 : use m_dtfil
32 :
33 : use defs_abitypes, only : MPI_type
34 : use m_time, only : timab
35 : use m_io_tools, only : get_unit
36 : use m_mpinfo, only : proc_distrb_cycle
37 :
38 : implicit none
39 :
40 : private
41 : !!***
42 :
43 : public :: optics_vloc
44 : !!***
45 :
46 : contains
47 : !!***
48 :
49 : !!****f* ABINIT/optics_vloc
50 : !! NAME
51 : !! optics_vloc
52 : !!
53 : !! FUNCTION
54 : !! Compute matrix elements need for optical conductivity in a LOCAL potential
55 : !! and store them in a file.
56 : !! Matrix elements = <Phi_i|Nabla|Phi_j>
57 : !!
58 : !! INPUTS
59 : !! cg(2,mcg)=planewave coefficients of wavefunctions.
60 : !! dtfil <type(datafiles_type)>=variables related to files
61 : !! dtset <type(dataset_type)>=all input variables for this dataset
62 : !! fildata= name of the output file
63 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
64 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
65 : !! mband=maximum number of bands
66 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
67 : !! mkmem =number of k points treated by this node.
68 : !! mpi_enreg=information about MPI parallelization
69 : !! mpw=maximum dimensioned size of npw.
70 : !! nkpt=number of k points.
71 : !! npwarr(nkpt)=number of planewaves in basis at this k point
72 : !! nsppol=1 for unpolarized, 2 for spin-polarized
73 : !!
74 : !! OUTPUT
75 : !! (only writing in a file)
76 : !!
77 : !! SIDE EFFECTS
78 : !!
79 : !! NOTES
80 : !!
81 : !! SOURCE
82 :
83 0 : subroutine optics_vloc(cg,dtfil,dtset,eigen0,gprimd,hdr,kg,mband,mcg,mkmem,mpi_enreg,mpw,&
84 0 : & nkpt,npwarr,nsppol)
85 : !Arguments ------------------------------------
86 : !scalars
87 : integer,intent(in) :: mband,mcg,mkmem,mpw,nkpt,nsppol
88 : type(datafiles_type),intent(in) :: dtfil
89 : type(dataset_type),intent(in) :: dtset
90 : type(MPI_type),intent(in) :: mpi_enreg
91 : type(hdr_type),intent(inout) :: hdr
92 : !arrays
93 : integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
94 : real(dp),intent(in) :: gprimd(3,3)
95 : real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
96 : real(dp),intent(inout) :: cg(2,mcg)
97 :
98 : !Local variables-------------------------------
99 : !scalars
100 : integer :: iomode,bdtot_index,cplex,etiq,fformopt,ib,icg,ierr,ikg,ikpt
101 : integer :: ipw,isppol,istwf_k,iwavef,jb,jwavef
102 : integer :: me,me_kpt,my_nspinor,nband_k,npw_k,sender,ount,pnp_size
103 : integer :: spaceComm_band,spaceComm_bandfftspin,spaceComm_fft,spaceComm_k
104 : logical :: mykpt
105 : real(dp) :: cgnm1,cgnm2
106 : character(len=500) :: msg
107 : !arrays
108 : integer :: tmp_shape(3)
109 0 : integer,allocatable :: kg_k(:,:)
110 : real(dp) :: kpoint(3),tsec(2)
111 0 : real(dp),allocatable :: eig0_k(:),kpg_k(:,:)
112 0 : real(dp),allocatable :: psinablapsi(:,:,:,:),tnm(:,:,:,:)
113 : type(wffile_type) :: wff1
114 : ! ************************************************************************
115 :
116 : DBG_ENTER("COLL")
117 :
118 : !Compatibility tests
119 :
120 : !----------------------------------------------------------------------------------
121 : !2- Computation of <psi_n|-i.nabla|psi_m> for each k
122 : !----------------------------------------------------------------------------------
123 :
124 : !Init parallelism
125 0 : if (mpi_enreg%paral_kgb==1) then
126 0 : spaceComm_k=mpi_enreg%comm_kpt
127 0 : spaceComm_fft=mpi_enreg%comm_fft
128 0 : spaceComm_band=mpi_enreg%comm_band
129 0 : spaceComm_bandfftspin=mpi_enreg%comm_bandspinorfft
130 0 : me_kpt=mpi_enreg%me_kpt
131 : else
132 0 : spaceComm_band=0;spaceComm_fft=0;spaceComm_bandfftspin=0
133 0 : spaceComm_k=mpi_enreg%comm_cell
134 0 : me=xmpi_comm_rank(spaceComm_k)
135 0 : me_kpt=me
136 : end if
137 0 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
138 :
139 : !Initialize main variables
140 0 : ABI_MALLOC(psinablapsi,(2,3,mband,mband))
141 0 : pnp_size=size(psinablapsi)
142 0 : psinablapsi=zero
143 :
144 0 : iomode= IO_MODE_FORTRAN_MASTER
145 0 : fformopt=612
146 0 : ount = get_unit()
147 0 : call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt,ierr,wff1,0,me,ount)
148 0 : call hdr_io(fformopt,hdr,2,wff1)
149 : !if (me == 0) then
150 : !call hdr_fort_write(hdr, wff1%unwff, fformopt,ierr)
151 : !ABI_CHECK(ierr /= 0, "hdr_fort_write returned ierr = 0")
152 : !end if
153 :
154 : !LOOP OVER SPINS
155 0 : icg=0
156 0 : do isppol=1,nsppol
157 :
158 : ! LOOP OVER k POINTS
159 0 : ikg=0
160 0 : do ikpt=1,nkpt
161 0 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
162 0 : etiq=ikpt+(isppol-1)*nkpt
163 0 : if (me==0) then
164 0 : ABI_MALLOC(eig0_k,(nband_k))
165 0 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
166 : end if
167 :
168 : ! Select my k-points
169 0 : mykpt=.true.
170 0 : mykpt=(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)))
171 0 : if (mykpt) then
172 :
173 : ! Allocations depending on k-point
174 0 : kpoint(:)=dtset%kptns(:,ikpt)
175 0 : istwf_k=dtset%istwfk(ikpt)
176 0 : npw_k=npwarr(ikpt)
177 0 : cplex=2;if (istwf_k>1) cplex=1
178 0 : ABI_MALLOC(kg_k,(3,npw_k))
179 0 : ABI_MALLOC(kpg_k,(npw_k*dtset%nspinor,3))
180 :
181 : ! Get G-vectors for this k-point
182 0 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
183 0 : ikg=ikg+npw_k
184 :
185 : ! Calculation of k+G in cartesian coordinates
186 0 : do ipw=1,npw_k
187 : kpg_k(ipw,1)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1)&
188 : & +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2)&
189 0 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3)
190 : kpg_k(ipw,2)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1)&
191 : & +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2)&
192 0 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3)
193 : kpg_k(ipw,3)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1)&
194 : & +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2)&
195 0 : & +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3)
196 : end do !ipw
197 0 : kpg_k=two_pi*kpg_k
198 0 : if (dtset%nspinor==2) kpg_k(npw_k+1:2*npw_k,1:3)=kpg_k(1:npw_k,1:3)
199 0 : ABI_FREE(kg_k)
200 :
201 : ! 2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m>
202 : ! ----------------------------------------------------------------------------------
203 : ! Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
204 :
205 0 : ABI_MALLOC(tnm,(2,3,nband_k,nband_k))
206 0 : tnm=zero
207 :
208 : ! Loops on bands
209 0 : do jb=1,nband_k
210 0 : jwavef=(jb-1)*npw_k*my_nspinor+icg
211 0 : if (mpi_enreg%paral_kgb/=1) then
212 0 : tmp_shape = shape(mpi_enreg%proc_distrb)
213 0 : if (ikpt > tmp_shape(1)) then
214 0 : msg=' ikpt out of bounds '
215 0 : ABI_BUG(msg)
216 : end if
217 0 : if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle
218 : end if
219 0 : do ib=1,jb
220 0 : iwavef=(ib-1)*npw_k*my_nspinor+icg
221 :
222 : ! Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
223 0 : if (cplex==1) then
224 0 : do ipw=1,npw_k*my_nspinor
225 0 : cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)
226 0 : tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
227 : end do
228 : else
229 0 : do ipw=1,npw_k*my_nspinor
230 0 : cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef)
231 0 : cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef)
232 0 : tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
233 0 : tnm(2,1:3,ib,jb)=tnm(2,1:3,ib,jb)+cgnm2*kpg_k(ipw,1:3)
234 : end do
235 : end if
236 :
237 : ! Second half of the (n,m) matrix
238 0 : if (ib/=jb) then
239 0 : tnm(1,1:3,jb,ib)= tnm(1,1:3,ib,jb)
240 0 : tnm(2,1:3,jb,ib)=-tnm(2,1:3,ib,jb)
241 : end if
242 :
243 : end do ! ib
244 : end do ! jb
245 :
246 : ! Reduction in case of parallelism
247 0 : if (mpi_enreg%paral_kgb == 1) then
248 0 : call timab(48,1,tsec)
249 0 : call xmpi_sum_master(tnm,0,spaceComm_bandfftspin,ierr)
250 0 : call timab(48,2,tsec)
251 : end if
252 :
253 0 : psinablapsi(:,:,:,:)=tnm(:,:,:,:)
254 :
255 0 : ABI_FREE(tnm)
256 :
257 0 : if (mkmem/=0) then
258 0 : icg = icg + npw_k*my_nspinor*nband_k
259 : end if
260 :
261 0 : ABI_FREE(kpg_k)
262 :
263 0 : if (me==0) then
264 0 : write(ount)(eig0_k(ib),ib=1,nband_k)
265 0 : write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
266 0 : write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
267 0 : write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
268 0 : elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then
269 0 : call xmpi_exch(psinablapsi,pnp_size,me_kpt,psinablapsi,0,spaceComm_k,etiq,ierr)
270 : end if
271 :
272 0 : elseif (me==0) then
273 0 : sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
274 0 : call xmpi_exch(psinablapsi,pnp_size,sender,psinablapsi,0,spaceComm_k,etiq,ierr)
275 0 : write(ount)(eig0_k(ib),ib=1,nband_k)
276 0 : write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
277 0 : write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
278 0 : write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
279 : end if ! mykpt
280 :
281 0 : bdtot_index=bdtot_index+nband_k
282 0 : if (me==0) then
283 0 : ABI_FREE(eig0_k)
284 : end if
285 : ! End loop on spin,kpt
286 : end do ! ikpt
287 : end do !isppol
288 :
289 : !Close file
290 0 : call WffClose(wff1,ierr)
291 :
292 : !Datastructures deallocations
293 0 : ABI_FREE(psinablapsi)
294 :
295 : DBG_EXIT("COLL")
296 :
297 0 : end subroutine optics_vloc
298 : !!***
299 :
300 : end module m_optics_vloc
301 : !!***
|
