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