TABLE OF CONTENTS


ABINIT/m_optics_vloc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_optics_vloc

FUNCTION

COPYRIGHT

  Copyright (C) 2010-2018 ABINIT group (SM,VR,FJ,MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_optics_vloc
28 
29  implicit none
30 
31  private

ABINIT/optics_vloc [ Functions ]

[ Top ] [ Functions ]

NAME

 optics_vloc

FUNCTION

 Compute matrix elements need for optical conductivity in a LOCAL potential
 and store them in a file.
 Matrix elements = <Phi_i|Nabla|Phi_j>

INPUTS

  cg(2,mcg)=planewave coefficients of wavefunctions.
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  fildata= name of the output file
  gprimd(3,3)=dimensional reciprocal space primitive translations
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw.
  nkpt=number of k points.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nsppol=1 for unpolarized, 2 for spin-polarized

OUTPUT

  (only writing in a file)

SIDE EFFECTS

NOTES

PARENTS

      outscfcv

CHILDREN

      hdr_io,timab,wffclose,wffopen,xmpi_exch,xmpi_sum_master

SOURCE

 80  subroutine optics_vloc(cg,dtfil,dtset,eigen0,gprimd,hdr,kg,mband,mcg,mkmem,mpi_enreg,mpw,&
 81 &                       nkpt,npwarr,nsppol)
 82 
 83  use defs_basis
 84  use defs_abitypes
 85  use m_abicore
 86  use m_errors
 87  use m_wffile
 88  use m_xmpi
 89  use m_hdr
 90 
 91  use m_time,         only : timab
 92  use m_io_tools,     only : get_unit
 93  use m_mpinfo,       only : proc_distrb_cycle
 94 
 95 !This section has been created automatically by the script Abilint (TD).
 96 !Do not modify the following lines by hand.
 97 #undef ABI_FUNC
 98 #define ABI_FUNC 'optics_vloc'
 99 !End of the abilint section
100 
101  implicit none
102 
103 !Arguments ------------------------------------
104 !scalars
105  integer,intent(in) :: mband,mcg,mkmem,mpw,nkpt,nsppol
106  type(datafiles_type),intent(in) :: dtfil
107  type(dataset_type),intent(in) :: dtset
108  type(MPI_type),intent(in) :: mpi_enreg
109  type(hdr_type),intent(inout) :: hdr
110 !arrays
111  integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
112  real(dp),intent(in) :: gprimd(3,3)
113  real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
114  real(dp),intent(inout) :: cg(2,mcg)
115 
116 !Local variables-------------------------------
117 !scalars
118  integer :: iomode,bdtot_index,cplex,etiq,fformopt,ib,icg,ierr,ikg,ikpt
119  integer :: ipw,isppol,istwf_k,iwavef,jb,jwavef
120  integer :: me,me_kpt,my_nspinor,nband_k,npw_k,sender,ount
121  integer :: spaceComm_band,spaceComm_bandfftspin,spaceComm_fft,spaceComm_k
122  logical :: mykpt
123  real(dp) :: cgnm1,cgnm2
124  character(len=500) :: msg
125 !arrays
126  integer :: tmp_shape(3)
127  integer,allocatable :: kg_k(:,:)
128  real(dp) :: kpoint(3),tsec(2)
129  real(dp),allocatable :: eig0_k(:),kpg_k(:,:)
130  real(dp),allocatable :: psinablapsi(:,:,:,:),tnm(:,:,:,:)
131  type(wffile_type) :: wff1
132 ! ************************************************************************
133 
134  DBG_ENTER("COLL")
135 
136 !Compatibility tests
137 
138 !----------------------------------------------------------------------------------
139 !2- Computation of <psi_n|-i.nabla|psi_m> for each k
140 !----------------------------------------------------------------------------------
141 
142 !Init parallelism
143  if (mpi_enreg%paral_kgb==1) then
144    spaceComm_k=mpi_enreg%comm_kpt
145    spaceComm_fft=mpi_enreg%comm_fft
146    spaceComm_band=mpi_enreg%comm_band
147    spaceComm_bandfftspin=mpi_enreg%comm_bandspinorfft
148    me_kpt=mpi_enreg%me_kpt
149  else
150    spaceComm_band=0;spaceComm_fft=0;spaceComm_bandfftspin=0
151    spaceComm_k=mpi_enreg%comm_cell
152    me=xmpi_comm_rank(spaceComm_k)
153    me_kpt=me
154  end if
155  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
156 
157 !Initialize main variables
158  ABI_ALLOCATE(psinablapsi,(2,3,mband,mband))
159  psinablapsi=zero
160 
161  iomode= IO_MODE_FORTRAN_MASTER
162  fformopt=612
163  ount = get_unit()
164  call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt,ierr,wff1,0,me,ount)
165  call hdr_io(fformopt,hdr,2,wff1)
166  !if (me == 0) then
167  !call hdr_fort_write(hdr, wff1%unwff, fformopt,ierr)
168  !ABI_CHECK(ierr /= 0, "hdr_fort_write returned ierr = 0")
169  !end if
170 
171 !LOOP OVER SPINS
172  icg=0
173  do isppol=1,nsppol
174 
175 !  LOOP OVER k POINTS
176    ikg=0
177    do ikpt=1,nkpt
178      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
179      etiq=ikpt+(isppol-1)*nkpt
180      if (me==0) then
181        ABI_ALLOCATE(eig0_k,(nband_k))
182        eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
183      end if
184 
185 !    Select my k-points
186      mykpt=.true.
187      mykpt=(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)))
188      if (mykpt) then
189 
190 !      Allocations depending on k-point
191        kpoint(:)=dtset%kptns(:,ikpt)
192        istwf_k=dtset%istwfk(ikpt)
193        npw_k=npwarr(ikpt)
194        cplex=2;if (istwf_k>1) cplex=1
195        ABI_ALLOCATE(kg_k,(3,npw_k))
196        ABI_ALLOCATE(kpg_k,(npw_k*dtset%nspinor,3))
197 
198 !      Get G-vectors for this k-point
199        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
200        ikg=ikg+npw_k
201 
202 !      Calculation of k+G in cartesian coordinates
203        do ipw=1,npw_k
204          kpg_k(ipw,1)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1)&
205 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2)&
206 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3)
207          kpg_k(ipw,2)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1)&
208 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2)&
209 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3)
210          kpg_k(ipw,3)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1)&
211 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2)&
212 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3)
213        end do !ipw
214        kpg_k=two_pi*kpg_k
215        if (dtset%nspinor==2) kpg_k(npw_k+1:2*npw_k,1:3)=kpg_k(1:npw_k,1:3)
216        ABI_DEALLOCATE(kg_k)
217 
218 !      2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m>
219 !      ----------------------------------------------------------------------------------
220 !      Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
221 
222        ABI_ALLOCATE(tnm,(2,3,nband_k,nband_k))
223        tnm=zero
224 
225 !      Loops on bands
226        do jb=1,nband_k
227          jwavef=(jb-1)*npw_k*my_nspinor+icg
228          if (mpi_enreg%paral_kgb/=1) then
229            tmp_shape = shape(mpi_enreg%proc_distrb)
230            if (ikpt > tmp_shape(1)) then
231              msg='  ikpt out of bounds '
232              MSG_BUG(msg)
233            end if
234            if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle
235          end if
236          do ib=1,jb
237            iwavef=(ib-1)*npw_k*my_nspinor+icg
238 
239 !          Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
240            if (cplex==1) then
241              do ipw=1,npw_k*my_nspinor
242                cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)
243                tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
244              end do
245            else
246              do ipw=1,npw_k*my_nspinor
247                cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef)
248                cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef)
249                tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
250                tnm(2,1:3,ib,jb)=tnm(2,1:3,ib,jb)+cgnm2*kpg_k(ipw,1:3)
251              end do
252            end if
253 
254 !          Second half of the (n,m) matrix
255            if (ib/=jb) then
256              tnm(1,1:3,jb,ib)= tnm(1,1:3,ib,jb)
257              tnm(2,1:3,jb,ib)=-tnm(2,1:3,ib,jb)
258            end if
259 
260          end do ! ib
261        end do ! jb
262 
263 !      Reduction in case of parallelism
264        if (mpi_enreg%paral_kgb == 1) then
265          call timab(48,1,tsec)
266          call xmpi_sum_master(tnm,0,spaceComm_bandfftspin,ierr)
267          call timab(48,2,tsec)
268        end if
269 
270        psinablapsi(:,:,:,:)=tnm(:,:,:,:)
271 
272        ABI_DEALLOCATE(tnm)
273 
274        if (mkmem/=0) then
275          icg = icg + npw_k*my_nspinor*nband_k
276        end if
277 
278        ABI_DEALLOCATE(kpg_k)
279 
280        if (me==0) then
281          write(ount)(eig0_k(ib),ib=1,nband_k)
282          write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
283          write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
284          write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
285        elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then
286          call xmpi_exch(psinablapsi,etiq,me_kpt,psinablapsi,0,spaceComm_k,ierr)
287        end if
288 
289      elseif (me==0) then
290        sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
291        call xmpi_exch(psinablapsi,etiq,sender,psinablapsi,0,spaceComm_k,ierr)
292        write(ount)(eig0_k(ib),ib=1,nband_k)
293        write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
294        write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
295        write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
296      end if ! mykpt
297 
298      bdtot_index=bdtot_index+nband_k
299      if (me==0)  then
300        ABI_DEALLOCATE(eig0_k)
301      end if
302 !    End loop on spin,kpt
303    end do ! ikpt
304  end do !isppol
305 
306 !Close file
307  call WffClose(wff1,ierr)
308 
309 !Datastructures deallocations
310  ABI_DEALLOCATE(psinablapsi)
311 
312  DBG_EXIT("COLL")
313 
314 end subroutine optics_vloc