1! 2! PCMSolver, an API for the Polarizable Continuum Model 3! Copyright (C) 2020 Roberto Di Remigio, Luca Frediani and contributors. 4! 5! This file is part of PCMSolver. 6! 7! PCMSolver is free software: you can redistribute it and/or modify 8! it under the terms of the GNU Lesser General Public License as published by 9! the Free Software Foundation, either version 3 of the License, or 10! (at your option) any later version. 11! 12! PCMSolver is distributed in the hope that it will be useful, 13! but WITHOUT ANY WARRANTY; without even the implied warranty of 14! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! GNU Lesser General Public License for more details. 16! 17! You should have received a copy of the GNU Lesser General Public License 18! along with PCMSolver. If not, see <http://www.gnu.org/licenses/>. 19! 20! For information on the complete list of contributors to the 21! PCMSolver API, see: <http://pcmsolver.readthedocs.io/> 22! 23 24 module pedra_symmetry 25! NOTE: the ibtfun module is gone, as we can safely use Fortran 26! standard intrinsic functions. 27! Mapping between intrinsics and ibtfun: 28! ibtand(i, j) <-> iand(i, j) 29! ibtor(i, j) <-> ior(i, j) 30! ibtshl(i, j) <-> ishft(i, j) 31! ibtshr(i, j) <-> ishft(i, -j) !WARNING! 32! ibtxor(i, j) <-> ieor(i, j) 33 34 use, intrinsic :: iso_c_binding 35 use pedra_precision 36 37 implicit none 38 39 public get_pt 40 public get_point_group 41! MINI-MANUAL 42! A. Indexing of symmetry operations and their mapping to a bitstring: 43! zyx Parity 44! 0 000 E 1.0 45! 1 001 Oyz -1.0 46! 2 010 Oxz -1.0 47! 3 011 C2z 1.0 48! 4 100 Oxy -1.0 49! 5 101 C2y 1.0 50! 6 110 C2x 1.0 51! 7 111 i -1.0 52! B. Indexing of isymax matrix 53! The isymax array contains the irrep to which the linear functions 54! (first column) and the rotations (second column) belong. The indexing 55! is given above. 56! C. Indexing of the jsop array 57! The jsop array contains the position at which the operation 58! i (index given in point A) appears 59! 60 61 type, public :: point_group 62 ! A type containing all you need to know about the point group 63 64 ! String with the group name: 65 ! C1, C2, Cs, Ci, D2, C2v, C2h, D2h 66 character(len=3) :: group_name 67 ! Integer representing the group 68 ! 0, 1, 2, 3, 4, 5, 6, 7 69 integer(kind=regint_k) :: group_int 70 ! Number of generators 71 integer(kind=regint_k) :: nr_generators 72 ! Number of not-trivial symmetry operations (2**nr_generators - 1) 73 integer(kind=regint_k) :: maxrep 74 ! group%isymax(i, 1): behaviour of principal axes under basic operations 75 ! (x-y-z) 76 ! group%isymax(i, 2): behaviour of principal rotations under basic operations 77 ! (Rx-Ry-Rz) 78 integer(kind=regint_k) :: isymax(3, 2) 79 ! Symmetry operations in the Abelian groups. 80 ! Bitstring: 1 coordinate changes sign under operation; 81 ! 0 coordinate does not change sign. 82 ! Of course, that's also the binary representation of 83 ! numbers from 0 to 7! 84 integer(kind=regint_k) :: jsop(0:7) 85 ! 86 integer(kind=regint_k) :: nr_rotations 87 ! 88 integer(kind=regint_k) :: nr_reflections 89 ! 90 integer(kind=regint_k) :: nr_inversion 91 end type point_group 92 93 private 94 95 integer(kind=regint_k) :: global_print_unit 96 97 contains 98 99 real(kind=dp) function get_pt(bit_rep) 100 101 integer(kind=regint_k), intent(in) :: bit_rep 102 103 real(kind=dp) :: pt(0:7) 104 105! 106! PT is the parity of a bitstring: 107! 1 for an even number of ones: 000,011,110,101 108! -1 for an odd number of ones: 001,010,100,111 109! 110 pt(0) = 1.0d0 111 pt(1) = -1.0d0 112 pt(2) = -1.0d0 113 pt(3) = 1.0d0 114 pt(4) = -1.0d0 115 pt(5) = 1.0d0 116 pt(6) = 1.0d0 117 pt(7) = -1.0d0 118 119 get_pt = pt(bit_rep) 120 121 end function get_pt 122 123 subroutine get_point_group(print_unit, pgroup, nr_gen, gen1, gen2, gen3) 124 125 type(point_group), intent(inout) :: pgroup 126 integer(kind=regint_k), intent(in) :: print_unit 127 integer(kind=regint_k), intent(in) :: nr_gen 128 integer(kind=regint_k), intent(in) :: gen1, gen2, gen3 129 130 global_print_unit = print_unit 131 pgroup = build_point_group(nr_gen, gen1, gen2, gen3) 132 133 end subroutine get_point_group 134 135 type(point_group) function build_point_group(nr_gen, gen1, gen2, gen3) 136! 137! Builds point group given the generators 138! Originally written by Trond Saue for DALTON/DIRAC 139! Copied and adapted by Roberto Di Remigio 140! 141 142 integer(kind=regint_k), intent(in) :: nr_gen 143 integer(kind=regint_k), intent(in) :: gen1, gen2, gen3 144! Local variables 145 integer(kind=regint_k) :: maxrep 146 integer(kind=regint_k) :: isymax(3, 2) 147 integer(kind=regint_k) :: igen(3) 148! Integer representation of the rotations bitmaps 149 integer(kind=regint_k), parameter :: irots(3) = [3, 5, 6] 150 integer(kind=regint_k), parameter :: rots(3) = [6, 5, 3] 151! Integer representation of the reflections bitmaps 152 integer(kind=regint_k), parameter :: irefl(3) = [4, 2, 1] 153! Parity of the symmetry operations bitmaps 154 integer(kind=regint_k), parameter :: jpar(0:7) = [1, -1, -1, 1, -1, 1, 1, -1] 155 integer(kind=regint_k) :: i, j, k, l, i0, i1, i2, ind, ipos, bitmap 156 integer(kind=regint_k) :: nrots, nrefl, ninvc, igroup 157 integer(kind=regint_k) :: char_tab(0:7, 0:7) 158 logical :: lsymop(0:7) 159 integer(kind=regint_k) :: jsop(0:7), ipar(0:7) 160 character(3) :: group, rep(0:7) 161 character(3), parameter :: groups(0:7) = ['C1 ', 'C2 ', 'Cs ', & 162 'Ci ', 'D2 ', 'C2v', & 163 'C2h', 'D2h'] 164 character(3), parameter :: symop(0:7) = [' E ', 'Oyz', 'Oxz', & 165 'C2z', 'Oxy', 'C2y', & 166 'C2x', ' i '] 167 168 isymax = 0 169 igen = 0 170 maxrep = 2**nr_gen - 1 171! igen contains the bitmap for the generators 172 igen = [gen1, gen2, gen3] 173! Build isymax(:, 1) 174! determine to which irrep the translations belong to 175! Loop over Cartesian axes 176 do i = 1, 3 177 bitmap = 0 178! Loop over generators 179 do j = 1, nr_gen 180! Apply generators on Cartesian axes rots(i) and check the character 181 if (nint(get_pt(ior(igen(j), rots(i)))) == -1) then 182! Set the bitmap 183 bitmap = ibset(bitmap, j) 184 end if 185 end do 186! Right-shift the bitmap and assign to isymax 187 isymax(i, 1) = ishft(bitmap, -1) 188 end do 189 190! Build isymax(:, 2) 191! determine to which irrep the rotations belong to 192! R_x = (y XOR z) and cyclic permutations 193 isymax(1, 2) = ieor(isymax(2, 1), isymax(3, 1)) 194 isymax(2, 2) = ieor(isymax(3, 1), isymax(1, 1)) 195 isymax(3, 2) = ieor(isymax(1, 1), isymax(2, 1)) 196 197! Build the character table 198 lsymop = .false. 199! Activate all symmetry operations of the group 200 lsymop(0) = .true. 201 jsop(0) = 0 202 ipar(0) = 1 203 do i = 1, maxrep 204 i0 = iand(1_regint_k, i) * igen(1) 205 i1 = iand(1_regint_k, ishft(i, -1)) * igen(2) 206 i2 = iand(1_regint_k, ishft(i, -2)) * igen(3) 207 ind = ieor(ieor(i0, i1),i2) 208 lsymop(ind) = .true. 209 ipar(i) = jpar(ind) 210 end do 211! List group operations in preferred order 212! Identity, E 213 ind = 0 214 jsop(ind) = 0 215! Rotations 216 nrots = 0 217 do i = 1, 3 218 if (lsymop(irots(i))) then 219 ind = ind + 1 220 jsop(ind) = irots(i) 221 nrots = nrots + 1 222 end if 223 end do 224! Inversion 225 ninvc = 0 226 if (lsymop(7)) then 227 ind = ind + 1 228 jsop(ind) = 7 229 ninvc = 1 230 end if 231! Reflections 232 nrefl = 0 233 do i = 1, 3 234 if (lsymop(irefl(i))) then 235 ind = ind + 1 236 jsop(ind) = irefl(i) 237 nrefl = nrefl + 1 238 end if 239 end do 240! Classify group 241! ============== 242! tsaue - Here I have devised a highly empirical formula, but it works !!! 243 igroup = min(7, nint((4 * nrots + 8 * ninvc + 6 * nrefl) / 3.0)) 244 group = groups(igroup) 245 char_tab = 0 246! Now generate the character table 247 do i = 0, maxrep 248 ! The character of the identity is always +1 249 char_tab(0, i) = 1 250 do j = 1, nr_gen 251 char_tab(igen(j), i) = nint(get_pt(iand(ishft(i,-(j-1)), 1_regint_k))) 252 do k = 1, (j-1) 253 ind = ieor(igen(j), igen(k)) 254 char_tab(ind, i) = char_tab(igen(j), i) * char_tab(igen(k), i) 255 do l = 1, (k-1) 256 char_tab(ieor(ind, igen(l)), i) = char_tab(ind, i) * char_tab(igen(l), i) 257 end do 258 end do 259 end do 260 end do 261! Classify irrep 262 do i = 0, maxrep 263 rep(i) = 'A ' 264 ipos = 2 265! 1. Rotational symmetry 266 if (nrots == 3) then 267 ind = (1 - char_tab(jsop(1), i)) + (1 - char_tab(jsop(2), i))/2 268 if (ind /= 0) then 269 rep(i)(1:1) = 'B' 270 rep(i)(2:2) = char(ichar('0') + ind) 271 ipos = 3 272 end if 273 else if (nrots == 1) then 274 if (char_tab(jsop(1), i) == -1) then 275 rep(i)(1:1) = 'B' 276 end if 277 if (nrefl == 2) then 278 if (iand(ishft(jsop(1), -1), 1_regint_k) == 1) then 279 ind = 2 280 else 281 ind = 3 282 end if 283 if (char_tab(jsop(ind), i) == 1) then 284 rep(i)(2:2) = '1' 285 else 286 rep(i)(2:2) = '2' 287 end if 288 end if 289 else if (nrefl == 1) then 290! 2. Mirror symmetry 291 if (char_tab(jsop(1), i) == 1) then 292 rep(i)(2:2) = '''' 293 else if (char_tab(jsop(1), i) == -1) then 294 rep(i)(2:2) = '"' 295 end if 296 end if 297! 3. Inversion symmetry 298 if (ninvc == 1) then 299 ind = nrots + 1 300 if (char_tab(jsop(ind), i) == 1) then 301 rep(i)(ipos:ipos) = 'g' 302 else 303 rep(i)(ipos:ipos) = 'u' 304 end if 305 end if 306 end do 307! Output 308! 1. Group name and generators 309 write(global_print_unit, '(a, a3)') 'Point group: ', group 310 if (nr_gen > 0) then 311 write(global_print_unit, '(/3x, a/)') '* The point group was generated by:' 312 do i = 1, nr_gen 313 if (symop(igen(i))(1:1) == 'C') then 314 write(global_print_unit, '(6x, 3a)') 'Rotation about the ', symop(igen(i))(3:3),'-axis' 315 else if (symop(igen(i))(1:1) == 'O') then 316 write(global_print_unit, '(6x, 3a)') 'Reflection in the ', symop(igen(i))(2:3),'-plane' 317 else 318 write(global_print_unit, '(6x, a)') 'Inversion center' 319 end if 320 end do 321! 2. Group multiplication table 322 write(global_print_unit,'(/3x, a/)') '* Group multiplication table' 323 write(global_print_unit,'(8x, a1, 8(1x, a3, 1x))') '|', (symop(jsop(i)), i = 0, maxrep) 324 write(global_print_unit,'(3x,a6,8a5)') '-----+', ('-----', i = 0, maxrep) 325 do i = 0, maxrep 326 write(global_print_unit,'(4x, a3, 1x, a1, 8(1x, a3, 1x))') symop(jsop(i)), '|', & 327 (symop(ieor(jsop(i), jsop(j))), j = 0, maxrep) 328 end do 329! 3. Character table 330 write(global_print_unit,'(/3x, a/)') '* Character table' 331 write(global_print_unit,'(8x, a1, 8(1x, a3, 1x))') '|', (symop(jsop(j)), j = 0, maxrep) 332 write(global_print_unit,'(3x, a6, 8a5)') '-----+', ('-----', i = 0, maxrep) 333 do i = 0, maxrep 334 write(global_print_unit,'(4x, a3, 1x, a1, 8(1x, i3, 1x))') rep(i), '|', (char_tab(jsop(j), i), j=0, maxrep) 335 end do 336! 4. Direct product table 337 write(global_print_unit,'(/3x, a/)') '* Direct product table' 338 write(global_print_unit,'(8x, a1, 8(1x, a3, 1x))') '|', (rep(i), i = 0, maxrep) 339 write(global_print_unit,'(3x, a6, 8a5)') '-----+', ('-----', i = 0, maxrep) 340 do i = 0, maxrep 341 write(global_print_unit,'(3x, 1x, a3, 1x, a1, 8(1x, a3, 1x))') rep(i), '|', (rep(ieor(i, j)), j = 0, maxrep) 342 end do 343 end if 344! Fields: group name, group integer(kind=regint_k), number of generators, 345! number of nontrivial operations, isymax, jsop, 346! number of rotations, number of reflections, 347! number of inversions. 348 build_point_group = point_group(group, igroup, nr_gen, maxrep, isymax, jsop, nrots, nrefl, ninvc) 349 350 end function build_point_group 351 352 end module pedra_symmetry 353