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