Branch data 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 : : !!***
|
