1!>==========================================================================
2!!
3!! Module sort_m:
4!!
5!! (1) gcutoff
6!!
7!!     Given G-vectors sorted by kinetic energy and an energy cutoff,
8!!     find the corresponding G-vector cutoff.
9!!
10!! (2,3) sortrx, sortix
11!!
12!!     Sorts an array by the quicksort method. real(DP) and integer versions.
13!!     See included sort_inc.f90.
14!!
15!! (4) make_identity_symmetry_first
16!!
17!!     The identity must always be the first symmetry, as assumed in various places
18!!     in the code. We enforce this by swapping it with op #1 if it is not first.
19!!
20!! (5) sort_symmetries
21!!
22!!     Bring symmetries into a standardized order. We are not currently using this.
23!!
24!!==========================================================================
25
26#include "f_defs.h"
27
28module sort_m
29  use global_m
30
31  implicit none
32
33  private
34
35  public ::   &
36    gcutoff,  &
37    sortrx, &
38    sortix, &
39    make_identity_symmetry_first, &
40    sort_symmetries, &
41    sort_sequential, &
42    sort_threaded, &
43    get_GK_array_from_gvecs
44
45  interface sortrx
46    module procedure sortrx_gvec, sortrx_no_gvec
47  end interface
48  interface sortix
49    module procedure sortix_gvec, sortix_no_gvec
50  end interface
51  interface sort_sequential
52    module procedure &
53      sortrx_no_gvec, sortrx_gvec_GK, sortrx_gvec, &
54      sortix_no_gvec, sortix_gvec_GK, sortix_gvec
55  end interface sort_sequential
56  interface sort_threaded
57    module procedure &
58      sortrx_no_gvec_threaded, &
59      sortrx_gvec_threaded_GK, &
60      sortrx_gvec_threaded,    &
61      sortix_no_gvec_threaded, &
62      sortix_gvec_threaded_GK, &
63      sortix_gvec_threaded
64  end interface sort_threaded
65
66
67contains
68
69  !> Given G-vectors sorted by kinetic energy and an energy cutoff, find the corresponding G-vector cutoff
70  !! such that all(ekin(isrtrq(ig)) <= ecutoff) for ig <= gcutoff.
71  integer function gcutoff(ng, ekin, isrtrq, ecutoff)
72    integer, intent(in) :: ng !< number of G-vectors
73    real(DP), intent(in) :: ekin(:) !< (ng) kinetic energies, should be sorted already
74    integer, intent(in) :: isrtrq(:) !< (ng) this is the index array returned by sorting ekin
75    real(DP), intent(in) :: ecutoff  !< energy cutoff, in same units as ekin (Ry typically)
76
77    integer :: gup, gdn, gmid, ig
78
79    PUSH_SUB(gcutoff)
80
81    ! perhaps all G-vectors fall within the cutoff
82    if(ekin(isrtrq(ng)) < ecutoff) then
83      gcutoff = ng
84      POP_SUB(gcutoff)
85      return
86    endif
87
88    ! otherwise, use bisection
89    gup = ng
90    gdn = 1
91
92    do ig = 1, ng
93      gmid = (gup + gdn) / 2
94      if(gmid == gdn) exit
95      if(ekin(isrtrq(gmid)) > ecutoff) then
96        gup = gmid
97      else
98        gdn = gmid
99      endif
100    enddo
101    gcutoff = gdn
102
103    POP_SUB(gcutoff)
104    return
105  end function gcutoff
106
107!=====================================================================
108!> The identity must always be the first symmetry, as assumed in various places in the code.
109  subroutine make_identity_symmetry_first(nsyms, mtrx, tnp)
110    integer,  intent(in)    :: nsyms
111    integer,  intent(inout) :: mtrx(3, 3, 48)
112    real(DP), intent(inout) :: tnp(3, 48)
113
114    integer :: isym, mtrx_temp(3, 3), identity(3,3)
115    real(DP) :: tnp_temp(3)
116    logical :: found
117
118    PUSH_SUB(make_identity_symmetry_first)
119
120    identity = reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), shape(identity))
121    found = all(mtrx(1:3, 1:3, 1) == identity(1:3, 1:3))
122
123    do isym = 2, nsyms
124      if(all(mtrx(1:3, 1:3, isym) == identity(1:3, 1:3))) then
125        if(.not. found) then
126          ! if identity is not first, swap
127          mtrx_temp(1:3, 1:3)  = mtrx(1:3, 1:3, 1)
128          mtrx(1:3, 1:3, 1)    = mtrx(1:3, 1:3, isym)
129          mtrx(1:3, 1:3, isym) = mtrx_temp(1:3, 1:3)
130
131          tnp_temp(1:3)  = tnp(1:3, 1)
132          tnp(1:3, 1)    = tnp(1:3, isym)
133          tnp(1:3, isym) = tnp_temp(1:3)
134
135          found = .true.
136          write(0,'(a,i2)') 'WARNING: making identity first by swapping with symmetry op #', isym
137        else
138          call die("There is a duplicate identity in the symmetry operations.")
139        endif
140      endif
141    enddo
142
143    if(.not. found) then
144      call die("Identity is not present in the list of symmetries.")
145    endif
146
147    POP_SUB(make_identity_symmetry_first)
148    return
149  end subroutine make_identity_symmetry_first
150
151!=====================================================================
152!> Bring symmetries into a standardized order.
153!! The identity is always the first one.
154  subroutine sort_symmetries(nsyms, mtrx, tnp)
155    integer,  intent(in)    :: nsyms
156    integer,  intent(inout) :: mtrx(3, 3, 48)
157    real(DP), intent(inout) :: tnp(3, 48)
158
159    integer :: isym, ii, jj, factor, hash(48), order(48), mtrx_temp(3, 3, 48), identity(3,3)
160    real(DP) :: tnp_temp(3, 48)
161
162    PUSH_SUB(sort_symmetries)
163
164    identity = reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), shape(identity))
165
166    do isym = 1, nsyms
167
168      ! make sure the identity comes first
169      if(all(mtrx(1:3, 1:3, isym) == identity(1:3, 1:3))) then
170        hash(isym) = -1d9
171        cycle
172      endif
173
174      hash(isym) = 0
175      factor = 1
176      do jj = 1, 3
177        if(jj > 1) factor = factor * 3
178        do ii = 1, 3
179          if(ii > 1) factor = factor * 3
180          hash(isym) = hash(isym) + mtrx(4 - ii, 4 - jj, isym) * factor
181        enddo
182      enddo
183    enddo
184
185    call sortix(nsyms, hash, order)
186
187    do isym = 1, nsyms
188      mtrx_temp(1:3, 1:3, isym) = mtrx(1:3, 1:3, order(isym))
189      tnp_temp(1:3, isym) = tnp(1:3, order(isym))
190    enddo
191
192    mtrx(1:3, 1:3, 1:nsyms) = mtrx_temp(1:3, 1:3, 1:nsyms)
193    tnp(1:3, 1:nsyms) = tnp_temp(1:3, 1:nsyms)
194
195    POP_SUB(sort_symmetries)
196    return
197  end subroutine sort_symmetries
198
199!=====================================================================
200
201!> Internal subroutine to get the ranking array GK from the G-vectors
202!! array gvec. This is used internally by the sorting routines, and should
203!! be transparent to most developers.
204subroutine get_GK_array_from_gvecs(NVAL, gvec, GK)
205  integer, intent(in) :: NVAL
206  integer, intent(in) :: gvec(3,NVAL)
207  integer, intent(out) :: GK(NVAL)
208
209  integer :: eff_grid(3), ii
210
211  PUSH_SUB(get_GK_array_from_gvecs)
212
213  eff_grid(1:3) = maxval(gvec(1:3,1:NVAL), 2) - minval(gvec(1:3,1:NVAL), 2) + 1
214  !!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ii)
215  do ii = 1, NVAL
216    GK(ii) = gvec(3,ii) + eff_grid(3)*(gvec(2,ii) + eff_grid(2)*gvec(1,ii))
217  enddo
218  !!$OMP END PARALLEL DO
219
220  POP_SUB(get_GK_array_from_gvecs)
221
222end subroutine get_GK_array_from_gvecs
223
224! FHJ: Use the preprocessor to create the following routines:
225! sortix_gvec, sortix_no_gvec, sortrx_gvec, sortrx_no_gvec
226
227#define DTYPE integer
228#define LABEL sortix_gvec
229#define HAS_GVEC
230#include "sort_inc.f90"
231#undef LABEL
232#undef HAS_GVEC
233#define LABEL sortix_no_gvec
234#include "sort_inc.f90"
235#undef LABEL
236#undef DTYPE
237
238#define DTYPE real(DP)
239#define LABEL sortrx_gvec
240#define HAS_GVEC
241#include "sort_inc.f90"
242#undef LABEL
243#undef HAS_GVEC
244#define LABEL sortrx_no_gvec
245#include "sort_inc.f90"
246#undef LABEL
247#undef DTYPE
248
249end module sort_m
250
251!! Local Variables:
252!! mode: f90
253!! coding: utf-8
254!! End:
255