1! ---
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt .
6! See Docs/Contributors.txt for a list of contributors.
7! ---
8module m_ts_kpoints
9!
10! Routines that are related to TS kpoint sampling
11!
12!==============================================================================
13! CONTAINS:
14!          1) setup_ts_scf_kscell
15!          2) setup_ts_kpoint_grid
16!          3) ts_write_k_points
17!          4) ts_iokp
18
19  use precision, only : dp
20
21  implicit none
22
23  private
24  save
25
26!===================== K-POINT RELATED VARIABLES ==============================
27!==============================================================================
28! Contains data structures and routines to deal with the kpoint-grid
29! for the self-consistent calculation with Green Functions.
30! Original verison modified so that the kpoint sampling is done for the
31! direction perpendicular to the transport directins
32! Version created to integrate TRANSIESTA in siesta2.3
33!==============================================================================
34! The kpoint mesh parameters that are public may be accessed in other
35! parts of the code. With respesct to the original names of the variables used
36! in the m_kpoint_grid module, a "ts_" was added. Also in the SIESTA module
37! kscell and kdispl are not public, but are made public (ts_kscell and
38! ts_kdispl) here.
39
40!==============================================================================
41! DETAILS: To obtain the kpoints for the GFs calculations it uses the same
42! scheme as for SIESTA but puts ts_kscell(3,3)=1 and ts_kdispl(3)=0.0
43!
44!==============================================================================
45
46  logical, public  :: ts_scf_kgrid_first_time = .true.
47  logical, public  :: ts_Gamma
48  integer, public  :: ts_nkpnt             ! Total number of k-points
49  real(dp), public :: ts_eff_kgrid_cutoff  ! Effective kgrid_cutoff
50
51  real(dp), pointer, public :: ts_kweight(:)
52  real(dp), pointer, public :: ts_kpoint(:,:)
53
54  integer,  public :: ts_kscell(3,3) = 0
55  real(dp), public :: ts_kdispl(3) = 0.0_dp
56
57  logical, public :: ts_firm_displ = .false.
58  logical, public :: ts_user_requested_mp = .false.
59  logical, public :: ts_user_requested_cutoff = .false.
60
61  public :: setup_ts_scf_kscell, setup_ts_kpoint_grid
62  public :: ts_write_k_points
63
64contains
65
66  subroutine setup_ts_scf_kscell( cell )
67
68! ***************** INPUT **********************************************
69! real*8  cell(3,3)  : Unit cell vectors in real space cell(ixyz,ivec)
70! type(Elec) Elecs(:) : Electrodes
71
72!   The relevant fdf labels are kgrid_cutoff and kgrid_Monkhorst_Pack.
73!   If both are present, kgrid_Monkhorst_Pack has priority. If none is
74!   present, the cutoff default is zero, producing only the gamma point.
75!   Examples of fdf data specifications:
76!     kgrid_cutoff  50. Bohr
77!     %block kgrid_Monkhorst_Pack  # Defines kscell and kdispl
78!     4  0  0   0.50               # (kscell(i,1),i=1,3), kdispl(1)
79!     0  4  0   0.50               # (kscell(i,2),i=1,3), kdispl(2)
80!     0  0  4   0.50               # (kscell(i,3),i=1,3), kdispl(3)
81!     %endblock kgrid_Monkhorst_Pack
82! **********************************************************************
83
84!  Modules
85
86    use parallel,   only : IONode
87    use m_minvec,   only : minvec
88    use fdf
89    use sys,        only : die
90#ifdef MPI
91    use mpi_siesta
92#endif
93
94    use m_ts_tdir, only: ts_tidx
95    use m_ts_global_vars, only : TSmode
96
97    implicit          none
98
99    ! Passed variables
100    real(dp), intent(in)   :: cell(3,3)
101
102    ! Internal variables
103    integer :: i, j, factor(3,3), expansion_factor
104    real(dp) :: scmin(3,3), vmod, cutoff
105    real(dp) :: ctransf(3,3)
106
107    type(block_fdf)            :: bfdf
108    type(parsed_line), pointer :: pline
109
110    logical :: bool
111    real(dp), parameter :: defcut = 0.0_dp
112
113    bool = fdf_block('TS.kgrid_Monkhorst_Pack',bfdf)
114    if ( .not. bool ) &
115         bool = fdf_block('kgrid_Monkhorst_Pack',bfdf)
116
117    if ( bool ) then
118       ts_user_requested_mp = .true.
119       do i = 1,3
120          if (fdf_bline(bfdf,pline)) then
121            ts_kscell(1,i) = fdf_bintegers(pline,1)
122            ts_kscell(2,i) = fdf_bintegers(pline,2)
123            ts_kscell(3,i) = fdf_bintegers(pline,3)
124            if ( fdf_bnvalues(pline) > 3 ) then
125              ts_kdispl(i) = mod(fdf_bvalues(pline,4), 1._dp)
126            else
127              ts_kdispl(i) = 0._dp
128            end if
129          else
130             call die( 'setup_ts_scf_kscell: ERROR no data in' // &
131                  'kgrid_Monkhorst_Pack block' )
132          endif
133       enddo
134       call fdf_bclose(bfdf)
135       ts_firm_displ = .true.
136
137    else
138
139       if ( IONode .and. TSmode ) then
140          write(*,*) 'WARNING !!!'
141          write(*,*) 'TS kgrid determined first with 3D cell !!!'
142          write(*,*) 'Specifying only cutoff in Electrode AND Scattering calculations might lead to problems !!'
143       end if
144
145       cutoff = fdf_get('kgrid_cutoff',defcut,'Bohr')
146       ts_user_requested_cutoff = (cutoff /= defcut)
147
148       ts_kdispl(1:3) = 0.0_dp  ! Might be changed later
149       ts_firm_displ = .false.
150
151       ! Find equivalent rounded unit-cell
152       call minvec( cell, scmin, ctransf )
153
154       expansion_factor = 1
155       do j = 1,3
156          factor(j,1:3) = 0
157          vmod = sqrt(dot_product(scmin(1:3,j),scmin(1:3,j)))
158          factor(j,j) = int(2.0_dp*cutoff/vmod) + 1
159          expansion_factor = expansion_factor * factor(j,j)
160       enddo
161       ! Generate actual supercell skeleton
162       ts_kscell = matmul(ctransf, factor)
163
164       if ( expansion_factor == 1 ) then
165          ts_kscell = 0
166          ts_kscell(1,1) = 1
167          ts_kscell(2,2) = 1
168          ts_kscell(3,3) = 1
169       end if
170
171    end if
172
173    ! In case of TSmode we have a transiesta run.
174    ! This means that we truncate the k-points in the transport direction.
175    ! However, if we are dealing with an electrode calculation we simply allow it
176    ! to have the same cell (the check made will disregard the transport directions k-points)
177    if ( TSmode .and. ts_tidx > 0 ) then
178       i = ts_tidx
179       ts_kscell(:,i) = 0
180       ts_kscell(i,:) = 0
181       ts_kscell(i,i) = 1
182       ts_kdispl(i)   = 0._dp
183    end if
184
185  end subroutine setup_ts_scf_kscell
186
187  subroutine setup_ts_kpoint_grid( ucell )
188
189    ! SIESTA Modules
190    USE precision, only : dp
191    USE fdf, only       : fdf_get
192    USE m_find_kgrid, only : find_kgrid
193    USE parallel, only  : IONode
194    use m_ts_global_vars, only : TSmode
195    use kpoint_grid
196
197    ! Local Variables
198    real(dp), intent(in) :: ucell(3,3)
199
200    if ( .not. TSMode ) then
201
202       ! Same as SIESTA (mainly to be able to write the TSHS files).
203       ts_Gamma = gamma_SCF
204       ts_nkpnt = nkpnt
205       ts_eff_kgrid_cutoff = eff_kgrid_cutoff
206       ts_kweight => kweight
207       ts_kpoint => kpoint
208       ts_kscell = kscell
209       ts_kdispl = kdispl
210
211       ! Since this is not transiesta we immediately return
212       ! and then we also skip printing out transiesta information.
213       return
214
215    end if
216
217    if ( ts_scf_kgrid_first_time ) then
218
219       nullify(ts_kweight,ts_kpoint)
220       if ( fdf_get('SpinSpiral', .false.) ) then
221          call die('transiesta: Does not work with spin-spiral')
222       end if
223
224       call setup_ts_scf_kscell(ucell)
225
226       ts_scf_kgrid_first_time = .false.
227
228    else
229       if ( ts_user_requested_mp    ) then
230          ! no need to set up the kscell again
231       else
232          call setup_ts_scf_kscell(ucell)
233       endif
234    endif
235
236    call find_kgrid(ucell,ts_kscell,ts_kdispl,ts_firm_displ, &
237         .true., &
238         ts_nkpnt,ts_kpoint,ts_kweight, ts_eff_kgrid_cutoff)
239
240    ts_Gamma = ts_nkpnt == 1 .and. &
241         dot_product(ts_kpoint(:,1),ts_kpoint(:,1)) < 1.0e-20_dp
242
243    if (IONode) call ts_write_k_points()
244
245  end subroutine setup_ts_kpoint_grid
246
247  subroutine ts_write_k_points()
248    USE siesta_options, only: writek
249    USE siesta_cml
250
251    implicit none
252
253    integer  :: ik, ix, i
254
255    if ( writek ) then
256       write(*,'(/,a)') 'transiesta: ts_k-point coordinates (Bohr**-1) and weights:'
257       write(*,'(a,i4,3f12.6,3x,f12.6)')                          &
258            ('transiesta: ', ik, (ts_kpoint(ix,ik),ix=1,3), ts_kweight(ik), &
259            ik=1,ts_nkpnt)
260    endif
261
262    ! Always write the TranSIESTA k-points
263    call ts_iokp( ts_nkpnt, ts_kpoint, ts_kweight )
264
265    write(*,'(/,a,i6)')  'transiesta: k-grid: Number of Green function k-points =', ts_nkpnt
266    write(*,'(a)') 'transiesta: k-grid: Supercell and displacements'
267    write(*,'(a,3i4,3x,f8.3)') 'transiesta: k-grid: ',        &
268         (ts_kscell(i,1),i=1,3), ts_kdispl(1)
269    write(*,'(a,3i4,3x,f8.3)') 'transiesta: k-grid: ',        &
270         (ts_kscell(i,2),i=1,3), ts_kdispl(2)
271    write(*,'(a,3i4,3x,f8.3)') 'transiesta: k-grid: ',        &
272         (ts_kscell(i,3),i=1,3), ts_kdispl(3)
273    if (cml_p) then
274       call cmlStartPropertyList(xf=mainXML, title="Transiesta k-points", &
275            dictRef="siesta:ts_kpoints")
276       call cmlAddProperty(xf=mainXML, value=ts_nkpnt, dictref='siesta:ts_nkpnt', &
277            units="cmlUnits:countable")
278       do ik=1, ts_nkpnt
279          call cmlAddKPoint(xf=mainXML, coords=ts_kpoint(:,ik), weight=ts_kweight(ik))
280       enddo
281       call cmlEndPropertyList(mainXML)
282       call cmlAddProperty(xf=mainXML, value=ts_kscell, dictref='siesta:ts_kscell', &
283            units="siestaUnits:Ang")
284       call cmlAddProperty(xf=mainXML, value=ts_kdispl, dictref='siesta:ts_kdispl', &
285            units="siestaUnits:Ang")
286    endif
287  end subroutine ts_write_k_points
288
289  subroutine ts_iokp( nk, points, weight )
290! *******************************************************************
291! Saves TranSIESTA k-points (only writing) Bohr^-1
292! Emilio Artacho, Feb. 1999
293! Modified by Nick Papior Andersen to not overwrite the SIESTA k-points
294! ********** INPUT **************************************************
295! integer nk           : Number of TS k-points
296! real*8  points(3,nk) : TS k-point coordinates
297! real*8  weight(3,nk) : TS k-point weight
298! *******************************************************************
299    use fdf
300    use files,     only : slabel, label_length
301
302    integer  :: nk
303    real(dp) :: points(3,nk), weight(nk)
304    external :: io_assign, io_close
305
306! Internal
307    character(len=label_length+5) :: fname
308    integer                       :: iu, ik, ix
309! -------------------------------------------------------------------
310
311    fname = trim( slabel ) // '.TSKP'
312
313    call io_assign( iu )
314    open( iu, file=fname, form='formatted', status='unknown' )
315
316    write(iu,'(i6)') nk
317    write(iu,'(i6,3f12.6,3x,f12.6)') &
318         (ik, (points(ix,ik),ix=1,3), weight(ik), ik=1,nk)
319
320    call io_close( iu )
321
322  end subroutine ts_iokp
323
324end module m_ts_kpoints
325