(git:34ef472)
molsym.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Molecular symmetry routines
10 !> \par History
11 !> 2008 adapted from older routines by Matthias Krack
12 !> \author jgh
13 ! **************************************************************************************************
14 MODULE molsym
15 
16  USE kinds, ONLY: dp
17  USE mathconstants, ONLY: pi
18  USE mathlib, ONLY: angle,&
19  build_rotmat,&
20  jacobi,&
21  reflect_vector,&
22  rotate_vector,&
23  unit_matrix,&
24  vector_product
25  USE physcon, ONLY: angstrom
26  USE util, ONLY: sort
27 #include "./base/base_uses.f90"
28 
29  IMPLICIT NONE
30  PRIVATE
31 
32  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
33  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molsym'
34 
35  INTEGER, PARAMETER :: maxcn = 20, &
36  maxsec = maxcn + 1, &
37  maxses = 2*maxcn + 1, &
38  maxsig = maxcn + 1, &
39  maxsn = 2*maxcn
40 
41  PUBLIC :: molsym_type
42  PUBLIC :: release_molsym, molecular_symmetry, print_symmetry
43 
44 ! **************************************************************************************************
45 !> \brief Container for information about molecular symmetry
46 !> \param ain : Atomic identification numbers (symmetry code).
47 !> \param aw : Atomic weights for the symmetry analysis.
48 !> \param eps_geo : Accuracy requested for analysis
49 !> \param inptostd : Transformation matrix for input to standard orientation.
50 !> \param point_group_symbol: Point group symbol.
51 !> \param rotmat : Rotation matrix.
52 !> \param sec : List of C axes
53 !> (sec(:,i,j) => x,y,z of the ith j-fold C axis).
54 !> \param ses : List of S axes
55 !> (ses(:,i,j) => x,y,z of the ith j-fold S axis).
56 !> \param sig : List of mirror planes
57 !> (sig(:,i) => x,y,z of the ith mirror plane).
58 !> \param center_of_mass : Shift vector for the center of mass.
59 !> \param tenmat : Molecular tensor of inertia.
60 !> \param tenval : Eigenvalues of the molecular tensor of inertia.
61 !> \param tenvec : Eigenvectors of the molecular tensor of inertia.
62 !> \param group_of : Group of equivalent atom.
63 !> \param llequatom : Lower limit of a group in symequ_list.
64 !> \param ncn : Degree of the C axis with highest degree.
65 !> \param ndod : Number of substituted dodecahedral angles.
66 !> \param nequatom : Number of equivalent atoms.
67 !> \param ngroup : Number of groups of equivalent atoms.
68 !> \param nico : Number of substituted icosahedral angles.
69 !> \param nlin : Number of substituted angles of 180 degrees.
70 !> \param nsec : Number of C axes.
71 !> \param nses : Number of S axes.
72 !> \param nsig : Number of mirror planes.
73 !> \param nsn : Degree of the S axis with highest degree.
74 !> \param ntet : Number of substituted tetrahedral angles.
75 !> \param point_group_order : Group order.
76 !> \param symequ_list : List of all atoms ordered in groups of equivalent atoms.
77 !> \param ulequatom : Upper limit of a group in symequ_list.
78 !> \param cubic : .TRUE., if a cubic point group was found (T,Th,Td,O,Oh,I,Ih).
79 !> \param dgroup: .TRUE., if a point group of D symmetry was found (Dn,Dnh,Dnd)
80 !> \param igroup: .TRUE., if a point group of icosahedral symmetry was found (I,Ih).
81 !> \param invers: .TRUE., if the molecule has a center of inversion.
82 !> \param linear: .TRUE., if the molecule is linear.
83 !> \param maxis : .TRUE., if the molecule has a main axis.
84 !> \param ogroup: .TRUE., if a point group of octahedral symmetry was found (O,Oh)
85 !> \param sgroup: .TRUE., if a point group of S symmetry was found.
86 !> \param sigmad: .TRUE., if there is a sigma_d mirror plane.
87 !> \param sigmah: .TRUE., if there is a sigma_h mirror plane.
88 !> \param sigmav: .TRUE., if there is a sigma_v mirror plane.
89 !> \param tgroup: .TRUE., if a point group of tetrahedral symmetry was found (T,Th,Td).
90 !> \par History
91 !> 05.2008 created
92 !> \author jgh
93 ! **************************************************************************************************
94  TYPE molsym_type
95  CHARACTER(LEN=4) :: point_group_symbol
96  INTEGER :: point_group_order
97  INTEGER :: ncn, ndod, ngroup, nico, nlin, nsig, nsn, ntet
98  LOGICAL :: cubic, dgroup, igroup, invers, linear, maxis, &
99  ogroup, sgroup, sigmad, sigmah, sigmav, tgroup
100  REAL(KIND=dp) :: eps_geo
101  REAL(KIND=dp), DIMENSION(3) :: center_of_mass, tenval
102  REAL(KIND=dp), DIMENSION(3) :: x_axis, y_axis, z_axis
103  REAL(KIND=dp), DIMENSION(:), POINTER :: ain, aw
104  REAL(KIND=dp), DIMENSION(3, 3) :: inptostd, rotmat, tenmat, tenvec
105  REAL(KIND=dp), DIMENSION(3, maxsig) :: sig
106  REAL(KIND=dp), DIMENSION(3, maxsec, 2:maxcn) :: sec
107  REAL(KIND=dp), DIMENSION(3, maxses, 2:maxsn) :: ses
108  INTEGER, DIMENSION(maxcn) :: nsec
109  INTEGER, DIMENSION(maxsn) :: nses
110  INTEGER, DIMENSION(:), POINTER :: group_of, llequatom, nequatom, &
111  symequ_list, ulequatom
112  END TYPE molsym_type
113 
114 ! **************************************************************************************************
115 
116 CONTAINS
117 
118 ! **************************************************************************************************
119 !> \brief Create an object of molsym type
120 !> \param sym ...
121 !> \param natoms ...
122 !> \author jgh
123 ! **************************************************************************************************
124  SUBROUTINE create_molsym(sym, natoms)
125  TYPE(molsym_type), POINTER :: sym
126  INTEGER, INTENT(IN) :: natoms
127 
128  IF (ASSOCIATED(sym)) CALL release_molsym(sym)
129 
130  ALLOCATE (sym)
131 
132  ALLOCATE (sym%ain(natoms), sym%aw(natoms), sym%group_of(natoms), sym%llequatom(natoms), &
133  sym%nequatom(natoms), sym%symequ_list(natoms), sym%ulequatom(natoms))
134 
135  END SUBROUTINE create_molsym
136 
137 ! **************************************************************************************************
138 !> \brief release an object of molsym type
139 !> \param sym ...
140 !> \author jgh
141 ! **************************************************************************************************
142  SUBROUTINE release_molsym(sym)
143  TYPE(molsym_type), POINTER :: sym
144 
145  cpassert(ASSOCIATED(sym))
146 
147  IF (ASSOCIATED(sym%aw)) THEN
148  DEALLOCATE (sym%aw)
149  END IF
150  IF (ASSOCIATED(sym%ain)) THEN
151  DEALLOCATE (sym%ain)
152  END IF
153  IF (ASSOCIATED(sym%group_of)) THEN
154  DEALLOCATE (sym%group_of)
155  END IF
156  IF (ASSOCIATED(sym%llequatom)) THEN
157  DEALLOCATE (sym%llequatom)
158  END IF
159  IF (ASSOCIATED(sym%nequatom)) THEN
160  DEALLOCATE (sym%nequatom)
161  END IF
162  IF (ASSOCIATED(sym%symequ_list)) THEN
163  DEALLOCATE (sym%symequ_list)
164  END IF
165  IF (ASSOCIATED(sym%ulequatom)) THEN
166  DEALLOCATE (sym%ulequatom)
167  END IF
168 
169  DEALLOCATE (sym)
170 
171  END SUBROUTINE release_molsym
172 
173 ! **************************************************************************************************
174 
175 ! **************************************************************************************************
176 !> \brief Add a new C_n axis to the list sec, but first check, if the
177 !> Cn axis is already in the list.
178 !> \param n ...
179 !> \param a ...
180 !> \param sym ...
181 !> \par History
182 !> Creation (19.10.98, Matthias Krack)
183 !> \author jgh
184 ! **************************************************************************************************
185  SUBROUTINE addsec(n, a, sym)
186  INTEGER, INTENT(IN) :: n
187  REAL(dp), DIMENSION(3), INTENT(IN) :: a
188  TYPE(molsym_type), INTENT(inout) :: sym
189 
190  INTEGER :: isec
191  REAL(dp) :: length_of_a, scapro
192  REAL(dp), DIMENSION(3) :: d
193 
194  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
195  d(:) = a(:)/length_of_a
196 
197  ! Check, if the current Cn axis is already in the list
198  DO isec = 1, sym%nsec(n)
199  scapro = sym%sec(1, isec, n)*d(1) + sym%sec(2, isec, n)*d(2) + sym%sec(3, isec, n)*d(3)
200  IF (abs(abs(scapro) - 1.0_dp) < sym%eps_geo) RETURN
201  END DO
202  sym%ncn = max(sym%ncn, n)
203 
204  ! Add the new Cn axis to the list sec
205  cpassert(sym%nsec(n) < maxsec)
206  sym%nsec(1) = sym%nsec(1) + 1
207  sym%nsec(n) = sym%nsec(n) + 1
208  sym%sec(:, sym%nsec(n), n) = d(:)
209 
210  END SUBROUTINE addsec
211 
212 ! **************************************************************************************************
213 !> \brief Add a new Sn axis to the list ses, but first check, if the
214 !> Sn axis is already in the list.
215 !> \param n ...
216 !> \param a ...
217 !> \param sym ...
218 !> \par History
219 !> Creation (19.10.98, Matthias Krack)
220 !> \author jgh
221 ! **************************************************************************************************
222  SUBROUTINE addses(n, a, sym)
223  INTEGER, INTENT(IN) :: n
224  REAL(dp), DIMENSION(3), INTENT(IN) :: a
225  TYPE(molsym_type), INTENT(inout) :: sym
226 
227  INTEGER :: ises
228  REAL(dp) :: length_of_a, scapro
229  REAL(dp), DIMENSION(3) :: d
230 
231  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
232  d(:) = a(:)/length_of_a
233 
234  ! Check, if the current Sn axis is already in the list
235  DO ises = 1, sym%nses(n)
236  scapro = sym%ses(1, ises, n)*d(1) + sym%ses(2, ises, n)*d(2) + sym%ses(3, ises, n)*d(3)
237  IF (abs(abs(scapro) - 1.0_dp) < sym%eps_geo) RETURN
238  END DO
239  sym%nsn = max(sym%nsn, n)
240 
241  ! Add the new Sn axis to the list ses
242  cpassert(sym%nses(n) < maxses)
243  sym%nses(1) = sym%nses(1) + 1
244  sym%nses(n) = sym%nses(n) + 1
245  sym%ses(:, sym%nses(n), n) = d(:)
246 
247  END SUBROUTINE addses
248 
249 ! **************************************************************************************************
250 !> \brief Add a new mirror plane to the list sig, but first check, if the
251 !> normal vector of the mirror plane is already in the list.
252 !> \param a ...
253 !> \param sym ...
254 !> \par History
255 !> Creation (19.10.98, Matthias Krack)
256 !> \author jgh
257 ! **************************************************************************************************
258  SUBROUTINE addsig(a, sym)
259  REAL(dp), DIMENSION(3), INTENT(IN) :: a
260  TYPE(molsym_type), INTENT(inout) :: sym
261 
262  INTEGER :: isig
263  REAL(dp) :: length_of_a, scapro
264  REAL(dp), DIMENSION(3) :: d
265 
266  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
267  d(:) = a(:)/length_of_a
268 
269  ! Check, if the normal vector of the current mirror plane is already in the list
270  DO isig = 1, sym%nsig
271  scapro = sym%sig(1, isig)*d(1) + sym%sig(2, isig)*d(2) + sym%sig(3, isig)*d(3)
272  IF (abs(abs(scapro) - 1.0_dp) < sym%eps_geo) RETURN
273  END DO
274 
275  ! Add the normal vector of the new mirror plane to the list sig
276  cpassert(sym%nsig < maxsig)
277  sym%nsig = sym%nsig + 1
278  sym%sig(:, sym%nsig) = d(:)
279 
280  END SUBROUTINE addsig
281 
282 ! **************************************************************************************************
283 !> \brief Symmetry handling for only one atom.
284 !> \param sym ...
285 !> \par History
286 !> Creation (19.10.98, Matthias Krack)
287 !> \author jgh
288 ! **************************************************************************************************
289  SUBROUTINE atomsym(sym)
290  TYPE(molsym_type), INTENT(inout) :: sym
291 
292 ! Set point group symbol
293 
294  sym%point_group_symbol = "R(3)"
295 
296  ! Set variables like D*h
297  sym%linear = .true.
298  sym%invers = .true.
299  sym%dgroup = .true.
300  sym%sigmah = .true.
301 
302  END SUBROUTINE atomsym
303 
304 ! **************************************************************************************************
305 !> \brief Search for point groups with AXial SYMmetry (Cn,Cnh,Cnv,Dn,Dnh,Dnd,S2n).
306 !> \param coord ...
307 !> \param sym ...
308 !> \par History
309 !> Creation (19.10.98, Matthias Krack)
310 !> \author jgh
311 ! **************************************************************************************************
312  SUBROUTINE axsym(coord, sym)
313  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
314  TYPE(molsym_type), INTENT(inout) :: sym
315 
316  INTEGER :: iatom, isig, jatom, m, n, natoms, nx, &
317  ny, nz
318  REAL(dp) :: phi
319  REAL(dp), DIMENSION(3) :: a, b, d
320 
321 ! Find the correct main axis and rotate main axis to z axis
322 
323  IF ((sym%ncn == 2) .AND. (sym%nses(2) == 3)) THEN
324  IF (sym%nsn == 4) THEN
325  ! Special case: D2d
326  phi = angle(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
327  d(:) = vector_product(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
328  CALL rotate_molecule(phi, d(:), sym, coord)
329  ELSE
330  ! Special cases: D2 and D2h
331  phi = 0.5_dp*pi
332  nx = naxis(sym%x_axis(:), coord, sym)
333  ny = naxis(sym%y_axis(:), coord, sym)
334  nz = naxis(sym%z_axis(:), coord, sym)
335  IF ((nx > ny) .AND. (nx > nz)) THEN
336  CALL rotate_molecule(-phi, sym%y_axis(:), sym, coord)
337  ELSE IF ((ny > nz) .AND. (ny > nx)) THEN
338  CALL rotate_molecule(phi, sym%x_axis(:), sym, coord)
339  END IF
340  END IF
341  ELSE
342  phi = angle(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
343  d(:) = vector_product(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
344  CALL rotate_molecule(phi, d(:), sym, coord)
345  END IF
346 
347  ! Search for C2 axes perpendicular to the main axis
348  ! Loop over all pairs of atoms (except on the z axis)
349  natoms = SIZE(coord, 2)
350  DO iatom = 1, natoms
351  a(:) = coord(:, iatom)
352  IF ((abs(a(1)) > sym%eps_geo) .OR. (abs(a(2)) > sym%eps_geo)) THEN
353  a(3) = 0.0_dp
354  IF (caxis(2, a(:), sym, coord)) CALL addsec(2, a(:), sym)
355  d(:) = vector_product(a(:), sym%z_axis(:))
356  IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
357  DO jatom = iatom + 1, natoms
358  b(:) = coord(:, jatom)
359  IF ((abs(b(1)) > sym%eps_geo) .OR. (abs(b(2)) > sym%eps_geo)) THEN
360  b(3) = 0.0_dp
361  phi = 0.5_dp*angle(a(:), b(:))
362  d(:) = vector_product(a(:), b(:))
363  b(:) = rotate_vector(a(:), phi, d(:))
364  IF (caxis(2, b(:), sym, coord)) CALL addsec(2, b(:), sym)
365  d(:) = vector_product(b(:), sym%z_axis)
366  IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
367  END IF
368  END DO
369  END IF
370  END DO
371 
372  ! Check the xy plane for mirror plane
373  IF (sigma(sym%z_axis(:), sym, coord)) THEN
374  CALL addsig(sym%z_axis(:), sym)
375  sym%sigmah = .true.
376  END IF
377 
378  ! Set logical variables for point group determination ***
379  IF (sym%nsec(2) > 1) THEN
380  sym%dgroup = .true.
381  IF ((.NOT. sym%sigmah) .AND. (sym%nsig > 0)) sym%sigmad = .true.
382  ELSE
383  IF ((.NOT. sym%sigmah) .AND. (sym%nsig > 0)) sym%sigmav = .true.
384  ! Cnh groups with n>3 were wrong
385  ! IF (sym%nsn > 3) sym%sgroup = .TRUE.
386  IF ((.NOT. sym%sigmah) .AND. (sym%nsn > 3)) sym%sgroup = .true.
387  END IF
388 
389  ! Rotate to standard orientation
390  n = 0
391  m = 0
392  IF ((sym%dgroup .AND. sym%sigmah) .OR. (sym%dgroup .AND. sym%sigmad) .OR. sym%sigmav) THEN
393  ! Dnh, Dnd or Cnv
394  DO isig = 1, sym%nsig
395  IF (abs(sym%sig(3, isig)) < sym%eps_geo) THEN
396  n = nsigma(sym%sig(:, isig), sym, coord)
397  IF (n > m) THEN
398  m = n
399  a(:) = sym%sig(:, isig)
400  END IF
401  END IF
402  END DO
403  IF (m > 0) THEN
404  ! Check for special case: C2v with all atoms in a plane
405  IF (sym%sigmav .AND. (m == natoms)) THEN
406  phi = angle(a(:), sym%x_axis(:))
407  d(:) = vector_product(a(:), sym%x_axis(:))
408  ELSE
409  phi = angle(a(:), sym%y_axis(:))
410  d(:) = vector_product(a(:), sym%y_axis(:))
411  END IF
412  CALL rotate_molecule(phi, d(:), sym, coord)
413  END IF
414  ELSE IF (sym%sigmah) THEN
415  ! Cnh
416  DO iatom = 1, natoms
417  IF (abs(coord(3, iatom)) < sym%eps_geo) THEN
418  n = naxis(coord(:, iatom), coord, sym)
419  IF (n > m) THEN
420  m = n
421  a(:) = coord(:, iatom)
422  END IF
423  END IF
424  END DO
425  IF (m > 0) THEN
426  phi = angle(a(:), sym%x_axis(:))
427  d(:) = vector_product(a(:), sym%x_axis(:))
428  CALL rotate_molecule(phi, d(:), sym, coord)
429  END IF
430  ! No action for Dn, Cn or S2n ***
431  END IF
432 
433  END SUBROUTINE axsym
434 
435 ! **************************************************************************************************
436 !> \brief Generate a symmetry list to identify equivalent atoms.
437 !> \param sym ...
438 !> \param coord ...
439 !> \par History
440 !> Creation (19.10.98, Matthias Krack)
441 !> \author jgh
442 ! **************************************************************************************************
443  SUBROUTINE build_symequ_list(sym, coord)
444  TYPE(molsym_type), INTENT(inout) :: sym
445  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
446 
447  INTEGER :: i, iatom, icn, iequatom, incr, isec, &
448  ises, isig, isn, jatom, jcn, jsn, &
449  natoms
450  REAL(kind=dp) :: phi
451  REAL(kind=dp), DIMENSION(3) :: a
452 
453  natoms = SIZE(coord, 2)
454 
455  IF (sym%sigmah) THEN
456  incr = 1
457  ELSE
458  incr = 2
459  END IF
460 
461  ! Initialize the atom and the group counter
462  iatom = 1
463  sym%ngroup = 1
464 
465  loop: DO
466 
467  ! Loop over all mirror planes
468  DO isig = 1, sym%nsig
469  a(:) = reflect_vector(coord(:, iatom), sym%sig(:, isig))
470  iequatom = equatom(iatom, a(:), sym, coord)
471  IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
472  sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
473  sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
474  sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
475  END IF
476  END DO
477 
478  ! Loop over all Cn axes
479  DO icn = 2, sym%ncn
480  DO isec = 1, sym%nsec(icn)
481  DO jcn = 1, icn - 1
482  IF (newse(jcn, icn)) THEN
483  phi = 2.0_dp*pi*real(jcn, kind=dp)/real(icn, kind=dp)
484  a(:) = rotate_vector(coord(:, iatom), phi, sym%sec(:, isec, icn))
485  iequatom = equatom(iatom, a(:), sym, coord)
486  IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
487  sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
488  sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
489  sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
490  END IF
491  END IF
492  END DO
493  END DO
494  END DO
495 
496  ! Loop over all Sn axes
497  DO isn = 2, sym%nsn
498  DO ises = 1, sym%nses(isn)
499  DO jsn = 1, isn - 1, incr
500  IF (newse(jsn, isn)) THEN
501  phi = 2.0_dp*pi*real(jsn, kind=dp)/real(isn, kind=dp)
502  a(:) = rotate_vector(coord(:, iatom), phi, sym%ses(:, ises, isn))
503  a(:) = reflect_vector(a(:), sym%ses(:, ises, isn))
504  iequatom = equatom(iatom, a(:), sym, coord)
505  IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
506  sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
507  sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
508  sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
509  END IF
510  END IF
511  END DO
512  END DO
513  END DO
514 
515  ! Exit loop, if all atoms are in the list
516  IF (sym%ulequatom(sym%ngroup) == natoms) EXIT
517 
518  ! Search for the next atom which is not in the list yet
519  DO jatom = 2, natoms
520  IF (.NOT. in_symequ_list(jatom, sym)) THEN
521  iatom = jatom
522  sym%ngroup = sym%ngroup + 1
523  sym%llequatom(sym%ngroup) = sym%ulequatom(sym%ngroup - 1) + 1
524  sym%ulequatom(sym%ngroup) = sym%llequatom(sym%ngroup)
525  sym%symequ_list(sym%llequatom(sym%ngroup)) = iatom
526  cycle loop
527  END IF
528  END DO
529 
530  END DO loop
531 
532  ! Generate list vector group_of
533  DO i = 1, sym%ngroup
534  DO iequatom = sym%llequatom(i), sym%ulequatom(i)
535  sym%group_of(sym%symequ_list(iequatom)) = i
536  END DO
537  END DO
538 
539  END SUBROUTINE build_symequ_list
540 
541 ! **************************************************************************************************
542 !> \brief Rotate the molecule about a n-fold axis defined by the vector a
543 !> using a rotation angle of 2*pi/n. Check, if the axis is a Cn axis.
544 !> \param n ...
545 !> \param a ...
546 !> \param sym ...
547 !> \param coord ...
548 !> \return ...
549 !> \par History
550 !> Creation (19.10.98, Matthias Krack)
551 !> \author jgh
552 ! **************************************************************************************************
553  FUNCTION caxis(n, a, sym, coord)
554  INTEGER, INTENT(IN) :: n
555  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
556  TYPE(molsym_type), INTENT(inout) :: sym
557  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
558  LOGICAL :: caxis
559 
560  INTEGER :: iatom, natoms
561  REAL(kind=dp) :: length_of_a, phi
562  REAL(kind=dp), DIMENSION(3) :: b
563 
564  caxis = .false.
565  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
566 
567  ! Check the length of the axis vector a
568  natoms = SIZE(coord, 2)
569  IF (length_of_a > sym%eps_geo) THEN
570  ! Calculate the rotation angle phi and build up the rotation matrix rotmat
571  phi = 2.0_dp*pi/real(n, dp)
572  CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
573  ! Check all atoms
574  DO iatom = 1, natoms
575  ! Rotate the actual atom by 2*pi/n about a
576  b(:) = matmul(sym%rotmat(:, :), coord(:, iatom))
577  ! Check if the coordinates of the rotated atom
578  ! are in the coordinate set of the molecule
579  IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
580  END DO
581  ! The axis defined by vector a is a Cn axis, thus return with caxis = .TRUE.
582  caxis = .true.
583  END IF
584 
585  END FUNCTION caxis
586 
587 ! **************************************************************************************************
588 !> \brief Check for a point group of cubic symmetry (I,Ih,O,Oh,T,Th,Td).
589 !> \param sym ...
590 !> \param coord ...
591 !> \param failed ...
592 !> \par History
593 !> Creation (19.10.98, Matthias Krack)
594 !> \author jgh
595 ! **************************************************************************************************
596  SUBROUTINE cubsym(sym, coord, failed)
597  TYPE(molsym_type), INTENT(inout) :: sym
598  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
599  LOGICAL, INTENT(INOUT) :: failed
600 
601  INTEGER :: i, iatom, iax, ic3, isec, jatom, jc3, &
602  jsec, katom, natoms, nc3
603  REAL(kind=dp) :: phi1, phi2, phidd, phidi, phiii
604  REAL(kind=dp), DIMENSION(3) :: a, b, c, d
605 
606  ! Angle between two adjacent atoms of the dodecahedron => <(C3,C3)
607  phidd = atan(0.4_dp*sqrt(5.0_dp))
608  ! Angle between two adjacent atoms of the icosahedron and the dodecahedron => <(C5,C3)
609  phidi = atan(3.0_dp - sqrt(5.0_dp))
610  ! Angle between two adjacent atoms of the icosahedron <(C5,C5)
611  phiii = atan(2.0_dp)
612 
613  ! Find a pair of C3 axes
614  natoms = SIZE(coord, 2)
615  DO iatom = 1, natoms
616  ! Check all atomic vectors for C3 axis
617  IF (caxis(3, coord(:, iatom), sym, coord)) THEN
618  CALL addsec(3, coord(:, iatom), sym)
619  IF (sym%nsec(3) > 1) EXIT
620  END IF
621  END DO
622 
623  ! Check the center of mass vector of a triangle
624  ! defined by three atomic vectors for C3 axis
625  IF (sym%nsec(3) < 2) THEN
626 
627  loop: DO iatom = 1, natoms
628  a(:) = coord(:, iatom)
629  DO jatom = iatom + 1, natoms
630  DO katom = jatom + 1, natoms
631  IF ((abs(dist(coord(:, iatom), coord(:, jatom)) &
632  - dist(coord(:, iatom), coord(:, katom))) < sym%eps_geo) .AND. &
633  (abs(dist(coord(:, iatom), coord(:, jatom)) &
634  - dist(coord(:, jatom), coord(:, katom))) < sym%eps_geo)) THEN
635  b(:) = a(:) + coord(:, jatom) + coord(:, katom)
636  IF (caxis(3, b(:), sym, coord)) THEN
637  CALL addsec(3, b(:), sym)
638  IF (sym%nsec(3) > 1) EXIT loop
639  END IF
640  END IF
641  END DO
642  END DO
643  END DO loop
644 
645  END IF
646 
647  ! Return after unsuccessful search for a pair of C3 axes
648  IF (sym%nsec(3) < 2) THEN
649  failed = .true.
650  RETURN
651  END IF
652 
653  ! Generate the remaining C3 axes
654  nc3 = 0
655  DO
656  nc3 = sym%nsec(3)
657  DO ic3 = 1, nc3
658  a(:) = sym%sec(:, ic3, 3)
659  DO jc3 = 1, nc3
660  IF (ic3 /= jc3) THEN
661  d(:) = sym%sec(:, jc3, 3)
662  DO i = 1, 2
663  phi1 = 2.0_dp*real(i, kind=dp)*pi/3.0_dp
664  b(:) = rotate_vector(a(:), phi1, d(:))
665  CALL addsec(3, b(:), sym)
666  END DO
667  END IF
668  END DO
669  END DO
670 
671  IF (sym%nsec(3) > nc3) THEN
672  ! New C3 axes were found
673  cycle
674  ELSE
675  ! In the case of I or Ih there have to be more C3 axes
676  IF (sym%nsec(3) == 4) THEN
677  a(:) = sym%sec(:, 1, 3)
678  b(:) = sym%sec(:, 2, 3)
679  phi1 = 0.5_dp*angle(a(:), b(:))
680  IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
681  d(:) = vector_product(a(:), b(:))
682  b(:) = rotate_vector(a(:), phi1, d(:))
683  c(:) = sym%sec(:, 3, 3)
684  phi1 = 0.5_dp*angle(a(:), c(:))
685  IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
686  d(:) = vector_product(a(:), c(:))
687  c(:) = rotate_vector(a(:), phi1, d(:))
688  d(:) = vector_product(b(:), c(:))
689  phi1 = 0.5_dp*phidd
690  a(:) = rotate_vector(b(:), phi1, d(:))
691  IF (caxis(3, a(:), sym, coord)) THEN
692  CALL addsec(3, a(:), sym)
693  ELSE
694  phi2 = 0.5_dp*pi - phi1
695  a(:) = rotate_vector(b(:), phi2, d(:))
696  IF (caxis(3, a(:), sym, coord)) CALL addsec(3, a(:), sym)
697  END IF
698 
699  IF (sym%nsec(3) > 4) cycle
700 
701  ELSE IF (sym%nsec(3) /= 10) THEN
702 
703  ! C3 axes found, but only 4 or 10 are possible
704  failed = .true.
705  RETURN
706 
707  END IF
708 
709  ! Exit loop, if 4 or 10 C3 axes were found
710  EXIT
711 
712  END IF
713 
714  END DO
715 
716  ! Loop over all pairs of C3 axes to find all C5, C4 and C2 axes
717  DO isec = 1, sym%nsec(3)
718 
719  a(:) = sym%sec(:, isec, 3)
720  DO jsec = isec + 1, sym%nsec(3)
721  phi1 = 0.5_dp*angle(a(:), sym%sec(:, jsec, 3))
722  d(:) = vector_product(a(:), sym%sec(:, jsec, 3))
723 
724  ! Check for C5 axis (I,Ih)
725  IF (sym%nsec(3) == 10) THEN
726  b(:) = rotate_vector(a(:), phidi, d(:))
727  IF (caxis(5, b(:), sym, coord)) THEN
728  CALL addsec(5, b(:), sym)
729  phi1 = phidi + phiii
730  b(:) = rotate_vector(a(:), phi1, d(:))
731  IF (caxis(5, b(:), sym, coord)) CALL addsec(5, b(:), sym)
732  END IF
733  END IF
734 
735  ! Check for C4 (O,Oh), C2 and S2 (center of inversion) axis
736  DO i = 0, 1
737  phi2 = phi1 - 0.5_dp*real(i, kind=dp)*pi
738  b(:) = rotate_vector(a(:), phi2, d(:))
739  IF (caxis(2, b(:), sym, coord)) THEN
740  CALL addsec(2, b(:), sym)
741  IF (sym%nsec(3) == 4) THEN
742  IF (caxis(4, b(:), sym, coord)) CALL addsec(4, b(:), sym)
743  END IF
744  IF (saxis(2, b(:), sym, coord)) THEN
745  CALL addses(2, b(:), sym)
746  sym%invers = .true.
747  END IF
748  END IF
749  IF (sigma(b(:), sym, coord)) CALL addsig(b(:), sym)
750  END DO
751 
752  ! Check for mirror plane
753  IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
754 
755  END DO
756 
757  END DO
758 
759  ! Set the logical variables for point group determination
760  iax = sym%ncn
761 
762  IF ((sym%ncn == 5) .AND. (sym%nsec(5) > 1)) THEN
763  sym%igroup = .true.
764  ELSE IF ((sym%ncn == 4) .AND. (sym%nsec(4) > 1)) THEN
765  sym%ogroup = .true.
766  ELSE IF ((sym%ncn == 3) .AND. (sym%nsec(3) > 1)) THEN
767  sym%tgroup = .true.
768  IF ((.NOT. sym%invers) .AND. (sym%nsig > 0)) sym%sigmad = .true.
769  iax = 2
770  ELSE
771  ! No main axis found
772  failed = .true.
773  RETURN
774  END IF
775 
776  ! Rotate molecule to standard orientation
777  phi1 = angle(sym%sec(:, 1, iax), sym%z_axis(:))
778  d(:) = vector_product(sym%sec(:, 1, iax), sym%z_axis(:))
779  CALL rotate_molecule(phi1, d(:), sym, coord)
780 
781  a(:) = sym%sec(:, 2, iax)
782  a(3) = 0.0_dp
783  phi2 = angle(a(:), sym%x_axis(:))
784  d(:) = vector_product(a(:), sym%x_axis(:))
785  CALL rotate_molecule(phi2, d(:), sym, coord)
786 
787  END SUBROUTINE cubsym
788 
789 ! **************************************************************************************************
790 !> \brief The number of a equivalent atom to atoma is returned. If there
791 !> is no equivalent atom, then zero is returned. The cartesian
792 !> coordinates of the equivalent atom are defined by the vector a.
793 !> \param atoma ...
794 !> \param a ...
795 !> \param sym ...
796 !> \param coord ...
797 !> \return ...
798 !> \par History
799 !> Creation (19.10.98, Matthias Krack)
800 !> \author jgh
801 ! **************************************************************************************************
802  FUNCTION equatom(atoma, a, sym, coord)
803  INTEGER, INTENT(IN) :: atoma
804  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
805  TYPE(molsym_type), INTENT(inout) :: sym
806  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
807  INTEGER :: equatom
808 
809  INTEGER :: iatom, natoms
810 
811  equatom = 0
812  natoms = SIZE(coord, 2)
813  DO iatom = 1, natoms
814  IF ((abs(sym%ain(iatom) - sym%ain(atoma)) < tiny(0.0_dp)) .AND. &
815  (abs(a(1) - coord(1, iatom)) < sym%eps_geo) .AND. &
816  (abs(a(2) - coord(2, iatom)) < sym%eps_geo) .AND. &
817  (abs(a(3) - coord(3, iatom)) < sym%eps_geo)) THEN
818  equatom = iatom
819  RETURN
820  END IF
821  END DO
822 
823  END FUNCTION equatom
824 
825 ! **************************************************************************************************
826 !> \brief Calculate the order of the point group.
827 !> \param sym ...
828 !> \par History
829 !> Creation (19.10.98, Matthias Krack)
830 !> \author jgh
831 ! **************************************************************************************************
832  SUBROUTINE get_point_group_order(sym)
833  TYPE(molsym_type), INTENT(inout) :: sym
834 
835  INTEGER :: icn, incr, isec, ises, isn, jcn, jsn
836 
837  ! Count all symmetry elements of the symmetry group
838  ! First E and all mirror planes
839  sym%point_group_order = 1 + sym%nsig
840 
841  ! Loop over all C axes
842  DO icn = 2, sym%ncn
843  DO isec = 1, sym%nsec(icn)
844  DO jcn = 1, icn - 1
845  IF (newse(jcn, icn)) sym%point_group_order = sym%point_group_order + 1
846  END DO
847  END DO
848  END DO
849 
850  ! Loop over all S axes
851  IF (sym%sigmah) THEN
852  incr = 1
853  ELSE
854  incr = 2
855  END IF
856 
857  DO isn = 2, sym%nsn
858  DO ises = 1, sym%nses(isn)
859  DO jsn = 1, isn - 1, incr
860  IF (newse(jsn, isn)) sym%point_group_order = sym%point_group_order + 1
861  END DO
862  END DO
863  END DO
864 
865  END SUBROUTINE get_point_group_order
866 
867 ! **************************************************************************************************
868 !> \brief Get the point group symbol.
869 !> \param sym ...
870 !> \par History
871 !> Creation (19.10.98, Matthias Krack)
872 !> \author jgh
873 ! **************************************************************************************************
874  SUBROUTINE get_point_group_symbol(sym)
875  TYPE(molsym_type), INTENT(inout) :: sym
876 
877  CHARACTER(LEN=4) :: degree
878 
879  IF (sym%linear) THEN
880  IF (sym%invers) THEN
881  sym%point_group_symbol = "D*h"
882  ELSE
883  sym%point_group_symbol = "C*v"
884  END IF
885  ELSE IF (sym%cubic) THEN
886  IF (sym%igroup) THEN
887  sym%point_group_symbol = "I"
888  ELSE IF (sym%ogroup) THEN
889  sym%point_group_symbol = "O"
890  ELSE IF (sym%tgroup) THEN
891  sym%point_group_symbol = "T"
892  END IF
893  IF (sym%invers) THEN
894  sym%point_group_symbol = trim(sym%point_group_symbol)//"h"
895  ELSE IF (sym%sigmad) THEN
896  sym%point_group_symbol = trim(sym%point_group_symbol)//"d"
897  END IF
898  ELSE IF (sym%maxis) THEN
899  IF (sym%dgroup) THEN
900  WRITE (degree, "(I3)") sym%ncn
901  sym%point_group_symbol = "D"//trim(adjustl(degree))
902  IF (sym%sigmah) THEN
903  sym%point_group_symbol = trim(sym%point_group_symbol)//"h"
904  ELSE IF (sym%sigmad) THEN
905  sym%point_group_symbol = trim(sym%point_group_symbol)//"d"
906  END IF
907  ELSE IF (sym%sgroup) THEN
908  WRITE (degree, "(I3)") sym%nsn
909  sym%point_group_symbol = "S"//trim(adjustl(degree))
910  ELSE
911  WRITE (degree, "(I3)") sym%ncn
912  sym%point_group_symbol = "C"//trim(adjustl(degree))
913  IF (sym%sigmah) THEN
914  sym%point_group_symbol = trim(sym%point_group_symbol)//"h"
915  ELSE IF (sym%sigmav) THEN
916  sym%point_group_symbol = trim(sym%point_group_symbol)//"v"
917  END IF
918  END IF
919  ELSE
920  IF (sym%invers) THEN
921  sym%point_group_symbol = "Ci"
922  ELSE IF (sym%sigmah) THEN
923  sym%point_group_symbol = "Cs"
924  ELSE
925  sym%point_group_symbol = "C1"
926  END IF
927  END IF
928 
929  END SUBROUTINE get_point_group_symbol
930 
931 ! **************************************************************************************************
932 !> \brief Initialization of the global variables of module symmetry.
933 !> \param sym ...
934 !> \param atype ...
935 !> \param weight ...
936 !> \par History
937 !> Creation (19.10.98, Matthias Krack)
938 !> \author jgh
939 ! **************************************************************************************************
940  SUBROUTINE init_symmetry(sym, atype, weight)
941  TYPE(molsym_type), INTENT(inout) :: sym
942  INTEGER, DIMENSION(:), INTENT(IN) :: atype
943  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: weight
944 
945  INTEGER :: i, iatom, natoms
946 
947  ! Define the Cartesian axis vectors
948  sym%x_axis = (/1.0_dp, 0.0_dp, 0.0_dp/)
949  sym%y_axis = (/0.0_dp, 1.0_dp, 0.0_dp/)
950  sym%z_axis = (/0.0_dp, 0.0_dp, 1.0_dp/)
951 
952  ! Initialize lists for symmetry elements
953  sym%sec(:, :, :) = 0.0_dp
954  sym%ses(:, :, :) = 0.0_dp
955  sym%sig(:, :) = 0.0_dp
956 
957  sym%center_of_mass(:) = 0.0_dp
958 
959  ! Initialize symmetry element counters
960  sym%ncn = 1
961  sym%nsn = 1
962  sym%nsec(:) = 0
963  sym%nses(:) = 0
964  sym%nsig = 0
965 
966  ! Initialize logical variables
967  sym%cubic = .false.
968  sym%dgroup = .false.
969  sym%igroup = .false.
970  sym%invers = .false.
971  sym%linear = .false.
972  sym%maxis = .false.
973  sym%ogroup = .false.
974  sym%sgroup = .false.
975  sym%sigmad = .false.
976  sym%sigmah = .false.
977  sym%sigmav = .false.
978  sym%tgroup = .false.
979 
980  ! Initialize list of equivalent atoms (C1 symmetry)
981  natoms = SIZE(sym%nequatom)
982  sym%ngroup = natoms
983  sym%nequatom(:) = 1
984  DO i = 1, sym%ngroup
985  sym%group_of(i) = i
986  sym%llequatom(i) = i
987  sym%symequ_list(i) = i
988  sym%ulequatom(i) = i
989  END DO
990 
991  sym%point_group_order = -1
992 
993  DO iatom = 1, natoms
994  sym%aw(iatom) = weight(iatom)
995  END DO
996 
997  ! Generate atomic identification numbers (symmetry code) ***
998  sym%ain(:) = real(atype(:), kind=dp) + sym%aw(:)
999 
1000  ! Initialize the transformation matrix for input orientation -> standard orientation
1001  CALL unit_matrix(sym%inptostd(:, :))
1002 
1003  END SUBROUTINE init_symmetry
1004 
1005 ! **************************************************************************************************
1006 !> \brief Return .TRUE., if the atom atoma is in the list symequ_list.
1007 !> \param atoma ...
1008 !> \param sym ...
1009 !> \return ...
1010 !> \par History
1011 !> Creation (21.04.95, Matthias Krack)
1012 !> \author jgh
1013 ! **************************************************************************************************
1014  FUNCTION in_symequ_list(atoma, sym)
1015  INTEGER, INTENT(IN) :: atoma
1016  TYPE(molsym_type), INTENT(inout) :: sym
1017  LOGICAL :: in_symequ_list
1018 
1019  INTEGER :: iatom
1020 
1021  in_symequ_list = .false.
1022 
1023  DO iatom = 1, sym%ulequatom(sym%ngroup)
1024  IF (sym%symequ_list(iatom) == atoma) THEN
1025  in_symequ_list = .true.
1026  RETURN
1027  END IF
1028  END DO
1029 
1030  END FUNCTION in_symequ_list
1031 
1032 ! **************************************************************************************************
1033 !> \brief Search for point groups of LOW SYMmetry (Ci,Cs).
1034 !> \param sym ...
1035 !> \param coord ...
1036 !> \par History
1037 !> Creation (21.04.95, Matthias Krack)
1038 !> \author jgh
1039 ! **************************************************************************************************
1040  SUBROUTINE lowsym(sym, coord)
1041  TYPE(molsym_type), INTENT(inout) :: sym
1042  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1043 
1044  REAL(kind=dp) :: phi
1045  REAL(kind=dp), DIMENSION(3) :: d
1046 
1047  IF (sym%nsn == 2) THEN
1048  ! Ci
1049  sym%invers = .true.
1050  phi = angle(sym%ses(:, 1, 2), sym%z_axis(:))
1051  d(:) = vector_product(sym%ses(:, 1, 2), sym%z_axis(:))
1052  CALL rotate_molecule(phi, d(:), sym, coord)
1053  ELSE IF (sym%nsig == 1) THEN
1054  ! Cs
1055  sym%sigmah = .true.
1056  phi = angle(sym%sig(:, 1), sym%z_axis(:))
1057  d(:) = vector_product(sym%sig(:, 1), sym%z_axis(:))
1058  CALL rotate_molecule(phi, d(:), sym, coord)
1059  END IF
1060 
1061  END SUBROUTINE lowsym
1062 
1063 ! **************************************************************************************************
1064 !> \brief Main program for symmetry analysis.
1065 !> \param sym ...
1066 !> \param eps_geo ...
1067 !> \param coord ...
1068 !> \param atype ...
1069 !> \param weight ...
1070 !> \par History
1071 !> Creation (14.11.98, Matthias Krack)
1072 !> \author jgh
1073 ! **************************************************************************************************
1074  SUBROUTINE molecular_symmetry(sym, eps_geo, coord, atype, weight)
1075  TYPE(molsym_type), POINTER :: sym
1076  REAL(kind=dp), INTENT(IN) :: eps_geo
1077  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1078  INTEGER, DIMENSION(:), INTENT(IN) :: atype
1079  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: weight
1080 
1081  CHARACTER(LEN=*), PARAMETER :: routinen = 'molecular_symmetry'
1082 
1083  INTEGER :: handle, natoms
1084 
1085  CALL timeset(routinen, handle)
1086 
1087  ! Perform memory allocation for the symmetry analysis
1088  natoms = SIZE(coord, 2)
1089  CALL create_molsym(sym, natoms)
1090  sym%eps_geo = eps_geo
1091 
1092  ! Initialization of symmetry analysis
1093  CALL init_symmetry(sym, atype, weight)
1094 
1095  IF (natoms == 1) THEN
1096  ! Special case: only one atom
1097  CALL atomsym(sym)
1098  ELSE
1099  ! Find molecular symmetry
1100  CALL moleculesym(sym, coord)
1101  ! Get point group and load point group table
1102  CALL get_point_group_symbol(sym)
1103  END IF
1104 
1105  ! Calculate group order
1106  IF (.NOT. sym%linear) CALL get_point_group_order(sym)
1107 
1108  ! Generate a list of equivalent atoms
1109  CALL build_symequ_list(sym, coord)
1110 
1111  CALL timestop(handle)
1112 
1113  END SUBROUTINE molecular_symmetry
1114 
1115 ! **************************************************************************************************
1116 !> \brief Find the molecular symmetry.
1117 !> \param sym ...
1118 !> \param coord ...
1119 !> \par History
1120 !> Creation (14.11.98, Matthias Krack)
1121 !> \author jgh
1122 ! **************************************************************************************************
1123  SUBROUTINE moleculesym(sym, coord)
1124  TYPE(molsym_type), INTENT(inout) :: sym
1125  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1126 
1127  INTEGER :: icn, isec, isn
1128  LOGICAL :: failed
1129  REAL(kind=dp) :: eps_tenval
1130 
1131 ! Tolerance of degenerate eigenvalues for the molecular tensor of inertia
1132 
1133  eps_tenval = 0.01_dp*sym%eps_geo
1134 
1135  ! Calculate the molecular tensor of inertia
1136  CALL tensor(sym, coord)
1137  ! Use symmetry information from the eigenvalues of the molecular tensor of inertia
1138  IF ((sym%tenval(3) - sym%tenval(1)) < eps_tenval) THEN
1139  ! 0 < tenval(1) = tenval(2) = tenval(3)
1140  sym%cubic = .true.
1141  ELSE IF ((sym%tenval(3) - sym%tenval(2)) < eps_tenval) THEN
1142  ! 0 < tenval(1) < tenval(2) = tenval(3)
1143  ! special case: 0 = tenval(1) < tenval(2) = tenval(3)
1144  IF (sym%tenval(1) < eps_tenval) sym%linear = .true.
1145  CALL tracar(sym, coord, (/3, 1, 2/))
1146  ELSE
1147  ! 0 < tenval(1) = tenval(2) < tenval(3) or 0 < tenval(1) < tenval(2) < tenval(3)
1148  CALL tracar(sym, coord, (/1, 2, 3/))
1149  END IF
1150 
1151  ! Continue with the new coordinate axes
1152  DO
1153  failed = .false.
1154  IF (sym%cubic) THEN
1155  CALL cubsym(sym, coord, failed)
1156  IF (failed) THEN
1157  sym%cubic = .false.
1158  cycle
1159  END IF
1160 
1161  ELSE IF (sym%linear) THEN
1162 
1163  ! Linear molecule
1164  IF (saxis(2, sym%z_axis(:), sym, coord)) THEN
1165  sym%invers = .true.
1166  sym%dgroup = .true.
1167  sym%sigmah = .true.
1168  END IF
1169 
1170  ELSE
1171 
1172  ! Check the new coordinate axes for Cn axes
1173  DO icn = 2, maxcn
1174  IF (caxis(icn, sym%z_axis(:), sym, coord)) CALL addsec(icn, sym%z_axis(:), sym)
1175  IF (caxis(icn, sym%x_axis(:), sym, coord)) CALL addsec(icn, sym%x_axis(:), sym)
1176  IF (caxis(icn, sym%y_axis(:), sym, coord)) CALL addsec(icn, sym%y_axis(:), sym)
1177  END DO
1178 
1179  ! Check the new coordinate axes for Sn axes
1180  DO isn = 2, maxsn
1181  IF (saxis(isn, sym%z_axis(:), sym, coord)) CALL addses(isn, sym%z_axis(:), sym)
1182  IF (saxis(isn, sym%x_axis(:), sym, coord)) CALL addses(isn, sym%x_axis(:), sym)
1183  IF (saxis(isn, sym%y_axis(:), sym, coord)) CALL addses(isn, sym%y_axis(:), sym)
1184  END DO
1185 
1186  ! Check the new coordinate planes for mirror planes
1187  IF (sigma(sym%z_axis(:), sym, coord)) CALL addsig(sym%z_axis(:), sym)
1188  IF (sigma(sym%x_axis(:), sym, coord)) CALL addsig(sym%x_axis(:), sym)
1189  IF (sigma(sym%y_axis(:), sym, coord)) CALL addsig(sym%y_axis(:), sym)
1190 
1191  ! There is a main axis (MAXIS = .TRUE.)
1192  IF ((sym%ncn > 1) .OR. (sym%nsn > 3)) THEN
1193  sym%maxis = .true.
1194  CALL axsym(coord, sym)
1195  ELSE
1196  ! Only low or no symmetry (Ci, Cs or C1)
1197  CALL lowsym(sym, coord)
1198  END IF
1199 
1200  END IF
1201 
1202  IF (.NOT. failed) EXIT
1203 
1204  END DO
1205 
1206  ! Find the remaining S axes
1207  IF (.NOT. sym%linear) THEN
1208  DO icn = 2, sym%ncn
1209  DO isec = 1, sym%nsec(icn)
1210  IF (saxis(icn, sym%sec(:, isec, icn), sym, coord)) &
1211  CALL addses(icn, sym%sec(:, isec, icn), sym)
1212  IF (saxis(2*icn, sym%sec(:, isec, icn), sym, coord)) &
1213  CALL addses(2*icn, sym%sec(:, isec, icn), sym)
1214  END DO
1215  END DO
1216  END IF
1217 
1218  ! Remove redundant S2 axes
1219  IF (sym%nses(2) > 0) THEN
1220  sym%nses(1) = sym%nses(1) - sym%nses(2)
1221  sym%nses(2) = 0
1222  CALL addses(2, sym%z_axis(:), sym)
1223  END IF
1224 
1225  END SUBROUTINE moleculesym
1226 
1227 ! **************************************************************************************************
1228 !> \brief The number of atoms which are placed on an axis defined by the vector a is returned.
1229 !> \param a ...
1230 !> \param coord ...
1231 !> \param sym ...
1232 !> \return ...
1233 !> \par History
1234 !> Creation (16.11.98, Matthias Krack)
1235 !> \author jgh
1236 ! **************************************************************************************************
1237  FUNCTION naxis(a, coord, sym)
1238  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1239  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1240  TYPE(molsym_type), INTENT(inout) :: sym
1241  INTEGER :: naxis
1242 
1243  INTEGER :: iatom, natoms
1244  REAL(kind=dp) :: length_of_a, length_of_b, scapro
1245  REAL(kind=dp), DIMENSION(3) :: a_norm, b, b_norm
1246 
1247  naxis = 0
1248  natoms = SIZE(coord, 2)
1249 
1250  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1251 
1252  ! Check the length of vector a
1253  IF (length_of_a > sym%eps_geo) THEN
1254 
1255  DO iatom = 1, natoms
1256  b(:) = coord(:, iatom)
1257  length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1258  ! An atom in the origin counts for each axis
1259  IF (length_of_b < sym%eps_geo) THEN
1260  naxis = naxis + 1
1261  ELSE
1262  a_norm = a(:)/length_of_a
1263  b_norm = b(:)/length_of_b
1264  scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1265  IF (abs(abs(scapro) - 1.0_dp) < sym%eps_geo) naxis = naxis + 1
1266  END IF
1267  END DO
1268 
1269  END IF
1270 
1271  END FUNCTION naxis
1272 
1273 ! **************************************************************************************************
1274 !> \brief Return .FALSE., if the symmetry elements se1 and se2 are a pair of
1275 !> redundant symmetry elements.
1276 !> \param se1 ...
1277 !> \param se2 ...
1278 !> \return ...
1279 !> \par History
1280 !> Creation (16.11.98, Matthias Krack)
1281 !> \author jgh
1282 ! **************************************************************************************************
1283  FUNCTION newse(se1, se2)
1284  INTEGER, INTENT(IN) :: se1, se2
1285  LOGICAL :: newse
1286 
1287  INTEGER :: ise
1288 
1289  newse = .true.
1290 
1291  DO ise = 2, min(se1, se2)
1292  IF ((modulo(se1, ise) == 0) .AND. (modulo(se2, ise) == 0)) THEN
1293  newse = .false.
1294  RETURN
1295  END IF
1296  END DO
1297 
1298  END FUNCTION newse
1299 
1300 ! **************************************************************************************************
1301 !> \brief The number of atoms which are placed in a plane defined by the
1302 !> the normal vector a is returned.
1303 !> \param a ...
1304 !> \param sym ...
1305 !> \param coord ...
1306 !> \return ...
1307 !> \par History
1308 !> Creation (16.11.98, Matthias Krack)
1309 !> \author jgh
1310 ! **************************************************************************************************
1311  FUNCTION nsigma(a, sym, coord)
1312  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1313  TYPE(molsym_type), INTENT(inout) :: sym
1314  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1315  INTEGER :: nsigma
1316 
1317  INTEGER :: iatom, natoms
1318  REAL(kind=dp) :: length_of_a, length_of_b, scapro
1319  REAL(kind=dp), DIMENSION(3) :: a_norm, b, b_norm
1320 
1321  nsigma = 0
1322 
1323  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1324 
1325  ! Check the length of vector a
1326  IF (length_of_a > sym%eps_geo) THEN
1327  natoms = SIZE(coord, 2)
1328  DO iatom = 1, natoms
1329  b(:) = coord(:, iatom)
1330  length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1331  ! An atom in the origin counts for each mirror plane
1332  IF (length_of_b < sym%eps_geo) THEN
1333  nsigma = nsigma + 1
1334  ELSE
1335  a_norm = a(:)/length_of_a
1336  b_norm = b(:)/length_of_b
1337  scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1338  IF (abs(scapro) < sym%eps_geo) nsigma = nsigma + 1
1339  END IF
1340  END DO
1341  END IF
1342 
1343  END FUNCTION nsigma
1344 
1345 ! **************************************************************************************************
1346 !> \brief Style the output of the symmetry elements.
1347 !> \param se ...
1348 !> \param eps ...
1349 !> \par History
1350 !> Creation (16.11.98, Matthias Krack)
1351 !> \author jgh
1352 ! **************************************************************************************************
1353  SUBROUTINE outse(se, eps)
1354  REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: se
1355  REAL(kind=dp), INTENT(IN) :: eps
1356 
1357 ! Positive z component for all vectors
1358 
1359  IF (abs(se(3)) < eps) THEN
1360  IF (abs(se(1)) < eps) THEN
1361  se(2) = -se(2)
1362  ELSE
1363  IF (se(1) < 0.0_dp) se(1:2) = -se(1:2)
1364  END IF
1365  ELSE
1366  IF (se(3) < 0.0_dp) se(:) = -se(:)
1367  END IF
1368 
1369  END SUBROUTINE outse
1370 
1371 ! **************************************************************************************************
1372 !> \brief Print the output of the symmetry analysis.
1373 !> \param sym ...
1374 !> \param coord ...
1375 !> \param atype ...
1376 !> \param element ...
1377 !> \param z ...
1378 !> \param weight ...
1379 !> \param iw ...
1380 !> \param plevel ...
1381 !> \par History
1382 !> Creation (16.11.98, Matthias Krack)
1383 !> \author jgh
1384 ! **************************************************************************************************
1385  SUBROUTINE print_symmetry(sym, coord, atype, element, z, weight, iw, plevel)
1386  TYPE(molsym_type), INTENT(inout) :: sym
1387  REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: coord
1388  INTEGER, DIMENSION(:), INTENT(in) :: atype
1389  CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: element
1390  INTEGER, DIMENSION(:), INTENT(in) :: z
1391  REAL(kind=dp), DIMENSION(:), INTENT(in) :: weight
1392  INTEGER, INTENT(IN) :: iw, plevel
1393 
1394  REAL(kind=dp), PARAMETER :: formatmaxval = 1.0e5_dp
1395 
1396  CHARACTER(LEN=3) :: string
1397  INTEGER :: i, icn, iequatom, isec, ises, isig, isn, &
1398  j, nequatom, secount
1399  INTEGER, ALLOCATABLE, DIMENSION(:) :: equatomlist, idxlist
1400 
1401  ! Print point group symbol and point group order
1402  WRITE (iw, "(/,(T2,A))") &
1403  "MOLSYM| Molecular symmetry information", &
1404  "MOLSYM|"
1405  WRITE (iw, "(T2,A,T77,A4)") &
1406  "MOLSYM| Point group", adjustr(sym%point_group_symbol)
1407  IF (sym%point_group_order > -1) THEN
1408  WRITE (iw, "(T2,A,T77,I4)") &
1409  "MOLSYM| Group order ", sym%point_group_order
1410  END IF
1411 
1412  IF (mod(plevel, 10) == 1) THEN
1413  WRITE (iw, "(T2,A)") &
1414  "MOLSYM|", &
1415  "MOLSYM| Groups of atoms equivalent by symmetry"
1416  WRITE (iw, "(T2,A,T31,A)") &
1417  "MOLSYM| Group number", "Group members (atomic indices)"
1418  DO i = 1, sym%ngroup
1419  nequatom = sym%ulequatom(i) - sym%llequatom(i) + 1
1420  cpassert(nequatom > 0)
1421  ALLOCATE (equatomlist(nequatom), idxlist(nequatom))
1422  equatomlist(1:nequatom) = sym%symequ_list(sym%llequatom(i):sym%ulequatom(i))
1423  idxlist = 0
1424  CALL sort(equatomlist, nequatom, idxlist)
1425  WRITE (iw, "(T2,A,T10,I5,T16,I6,T31,10(1X,I4),/,(T2,'MOLSYM|',T31,10(1X,I4)))") &
1426  "MOLSYM|", i, nequatom, (equatomlist(iequatom), iequatom=1, nequatom)
1427  DEALLOCATE (equatomlist, idxlist)
1428  END DO
1429  END IF
1430 
1431  IF (mod(plevel/100, 10) == 1) THEN
1432  ! Print all symmetry elements
1433  WRITE (iw, "(T2,A,/,T2,A,/,T2,A,T31,A,T45,A,T60,A,T75,A)") &
1434  "MOLSYM|", &
1435  "MOLSYM| Symmetry elements", &
1436  "MOLSYM| Element number", "Type", "X", "Y", "Z"
1437  secount = 1
1438  WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A1,T36,3(1X,F14.8))") &
1439  "MOLSYM|", secount, secount, "E", 0.0_dp, 0.0_dp, 0.0_dp
1440  ! Mirror planes
1441  string = "@ "
1442  DO isig = 1, sym%nsig
1443  secount = secount + 1
1444  CALL outse(sym%sig(:, isig), sym%eps_geo)
1445  WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
1446  "MOLSYM|", secount, isig, string, (sym%sig(i, isig), i=1, 3)
1447  END DO
1448  ! C axes
1449  DO icn = 2, sym%ncn
1450  IF (icn < 10) THEN
1451  WRITE (string, "(A1,I1)") "C", icn
1452  ELSE
1453  WRITE (string, "(A1,I2)") "C", icn
1454  END IF
1455  DO isec = 1, sym%nsec(icn)
1456  secount = secount + 1
1457  CALL outse(sym%sec(:, isec, icn), sym%eps_geo)
1458  WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
1459  "MOLSYM|", secount, isec, string, (sym%sec(i, isec, icn), i=1, 3)
1460  END DO
1461  END DO
1462  ! S axes
1463  DO isn = 2, sym%nsn
1464  IF (isn == 2) THEN
1465  WRITE (string, "(A3)") "i "
1466  ELSE IF (isn < 10) THEN
1467  WRITE (string, "(A1,I1)") "S", isn
1468  ELSE
1469  WRITE (string, "(A1,I2)") "S", isn
1470  END IF
1471  DO ises = 1, sym%nses(isn)
1472  secount = secount + 1
1473  CALL outse(sym%ses(:, ises, isn), sym%eps_geo)
1474  WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
1475  "MOLSYM|", secount, ises, string, (sym%ses(i, ises, icn), i=1, 3)
1476  END DO
1477  END DO
1478  END IF
1479 
1480  IF (mod(plevel/10, 10) == 1 .AND. SIZE(coord, 2) > 1) THEN
1481  ! Print the moments of the molecular inertia tensor
1482  WRITE (iw, "(T2,A,/,T2,A,/,T2,A,T43,A,T58,A,T73,A)") &
1483  "MOLSYM|", &
1484  "MOLSYM| Moments of the molecular inertia tensor [a.u.]", &
1485  "MOLSYM|", "I(x)", "I(y)", "I(z)"
1486  string = "xyz"
1487  IF (maxval(sym%tenmat) >= formatmaxval) THEN
1488  DO i = 1, 3
1489  WRITE (iw, "(T2,A,T32,A,T36,3(1X,ES14.4))") &
1490  "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1491  END DO
1492  ELSE
1493  DO i = 1, 3
1494  WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1495  "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1496  END DO
1497  END IF
1498  WRITE (iw, "(T2,A)") &
1499  "MOLSYM|", &
1500  "MOLSYM| Principal moments and axes of inertia [a.u.]"
1501  IF (maxval(sym%tenmat) >= formatmaxval) THEN
1502  WRITE (iw, "(T2,A,T36,3(1X,ES14.4))") &
1503  "MOLSYM|", (sym%tenval(i), i=1, 3)
1504  ELSE
1505  WRITE (iw, "(T2,A,T36,3(1X,F14.8))") &
1506  "MOLSYM|", (sym%tenval(i), i=1, 3)
1507  END IF
1508  DO i = 1, 3
1509  WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1510  "MOLSYM|", string(i:i), (sym%tenvec(i, j), j=1, 3)
1511  END DO
1512  END IF
1513 
1514  IF (mod(plevel, 10) == 1) THEN
1515  ! Print the Cartesian coordinates of the standard orientation
1516  CALL write_coordinates(coord, atype, element, z, weight, iw)
1517  END IF
1518 
1519  END SUBROUTINE print_symmetry
1520 
1521 ! **************************************************************************************************
1522 !> \brief Rotate the molecule about an axis defined by the vector a. The
1523 !> rotation angle is phi (radians).
1524 !> \param phi ...
1525 !> \param a ...
1526 !> \param sym ...
1527 !> \param coord ...
1528 !> \par History
1529 !> Creation (16.11.98, Matthias Krack)
1530 !> \author jgh
1531 ! **************************************************************************************************
1532  SUBROUTINE rotate_molecule(phi, a, sym, coord)
1533  REAL(kind=dp), INTENT(IN) :: phi
1534  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1535  TYPE(molsym_type), INTENT(inout) :: sym
1536  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1537 
1538  INTEGER :: i
1539  REAL(kind=dp) :: length_of_a
1540 
1541  ! Check the length of vector a
1542  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1543  IF (length_of_a > sym%eps_geo) THEN
1544 
1545  ! Build up the rotation matrix
1546  CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1547 
1548  ! Rotate the molecule by phi around a
1549  coord(:, :) = matmul(sym%rotmat(:, :), coord(:, :))
1550 
1551  ! Rotate the normal vectors of the mirror planes which are already found
1552  sym%sig(:, 1:sym%nsig) = matmul(sym%rotmat(:, :), sym%sig(:, 1:sym%nsig))
1553 
1554  ! Rotate the Cn axes which are already found
1555  DO i = 2, sym%ncn
1556  sym%sec(:, 1:sym%nsec(i), i) = matmul(sym%rotmat(:, :), sym%sec(:, 1:sym%nsec(i), i))
1557  END DO
1558 
1559  ! Rotate the Sn axes which are already found
1560  DO i = 2, sym%nsn
1561  sym%ses(:, 1:sym%nses(i), i) = matmul(sym%rotmat(:, :), sym%ses(:, 1:sym%nses(i), i))
1562  END DO
1563 
1564  ! Store current rotation
1565  sym%inptostd(:, :) = matmul(sym%rotmat(:, :), sym%inptostd(:, :))
1566 
1567  END IF
1568 
1569  END SUBROUTINE rotate_molecule
1570 
1571 ! **************************************************************************************************
1572 !> \brief Rotate the molecule about a n-fold axis defined by the vector a
1573 !> using a rotation angle of 2*pi/n. Check, if the axis is a Cn axis.
1574 !> \param n ...
1575 !> \param a ...
1576 !> \param sym ...
1577 !> \param coord ...
1578 !> \return ...
1579 !> \par History
1580 !> Creation (16.11.98, Matthias Krack)
1581 !> \author jgh
1582 ! **************************************************************************************************
1583  FUNCTION saxis(n, a, sym, coord)
1584  INTEGER, INTENT(IN) :: n
1585  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1586  TYPE(molsym_type), INTENT(inout) :: sym
1587  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1588  LOGICAL :: saxis
1589 
1590  INTEGER :: iatom, natoms
1591  REAL(kind=dp) :: length_of_a, phi
1592  REAL(kind=dp), DIMENSION(3) :: b
1593 
1594  saxis = .false.
1595 
1596  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1597 
1598  natoms = SIZE(coord, 2)
1599 
1600  ! Check the length of the axis vector a
1601  IF (length_of_a > sym%eps_geo) THEN
1602  ! Calculate the rotation angle phi and build up the rotation matrix rotmat
1603  phi = 2.0_dp*pi/real(n, kind=dp)
1604  CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1605  ! Check all atoms
1606  DO iatom = 1, natoms
1607  ! Rotate the actual atom by 2*pi/n about a
1608  b(:) = matmul(sym%rotmat(:, :), coord(:, iatom))
1609  ! Reflect the coordinates of the rotated atom
1610  b(:) = reflect_vector(b(:), a(:))
1611  ! Check if the coordinates of the rotated atom are in the coordinate set of the molecule
1612  IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
1613  END DO
1614  ! The axis defined by a is a Sn axis, thus return with saxis = .TRUE.
1615  saxis = .true.
1616  END IF
1617 
1618  END FUNCTION saxis
1619 
1620 ! **************************************************************************************************
1621 !> \brief Reflect all atoms of the molecule through a mirror plane defined
1622 !> by the normal vector a.
1623 !> \param a ...
1624 !> \param sym ...
1625 !> \param coord ...
1626 !> \return ...
1627 !> \par History
1628 !> Creation (16.11.98, Matthias Krack)
1629 !> \author jgh
1630 ! **************************************************************************************************
1631  FUNCTION sigma(a, sym, coord)
1632  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1633  TYPE(molsym_type), INTENT(inout) :: sym
1634  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1635  LOGICAL :: sigma
1636 
1637  INTEGER :: iatom, natoms
1638  REAL(kind=dp) :: length_of_a
1639  REAL(kind=dp), DIMENSION(3) :: b
1640 
1641  sigma = .false.
1642 
1643  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1644 
1645  ! Check the length of vector a
1646  IF (length_of_a > sym%eps_geo) THEN
1647  natoms = SIZE(coord, 2)
1648  DO iatom = 1, natoms
1649  ! Reflect the actual atom
1650  b(:) = reflect_vector(coord(:, iatom), a(:))
1651  ! Check if the coordinates of the reflected atom are in the coordinate set of the molecule
1652  IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
1653  END DO
1654  ! The plane defined by a is a mirror plane, thus return with sigma = .TRUE.
1655  sigma = .true.
1656  END IF
1657 
1658  END FUNCTION sigma
1659 
1660 ! **************************************************************************************************
1661 !> \brief Calculate the molecular tensor of inertia and the vector to the
1662 !> center of mass of the molecule.
1663 !> \param sym ...
1664 !> \param coord ...
1665 !> \par History
1666 !> Creation (16.11.98, Matthias Krack)
1667 !> \author jgh
1668 ! **************************************************************************************************
1669  SUBROUTINE tensor(sym, coord)
1670  TYPE(molsym_type), INTENT(inout) :: sym
1671  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1672 
1673  INTEGER :: natoms
1674  REAL(kind=dp) :: tt
1675 
1676  ! Find the vector center_of_mass to molecular center of mass
1677  natoms = SIZE(coord, 2)
1678  sym%center_of_mass(:) = matmul(coord(:, 1:natoms), sym%aw(1:natoms))/sum(sym%aw(1:natoms))
1679 
1680  ! Translate the center of mass of the molecule to the origin
1681  coord(:, 1:natoms) = coord(:, 1:natoms) - spread(sym%center_of_mass(:), 2, natoms)
1682 
1683  ! Build up the molecular tensor of inertia
1684 
1685  sym%tenmat(1, 1) = dot_product(sym%aw(1:natoms), (coord(2, 1:natoms)**2 + coord(3, 1:natoms)**2))
1686  sym%tenmat(2, 2) = dot_product(sym%aw(1:natoms), (coord(3, 1:natoms)**2 + coord(1, 1:natoms)**2))
1687  sym%tenmat(3, 3) = dot_product(sym%aw(1:natoms), (coord(1, 1:natoms)**2 + coord(2, 1:natoms)**2))
1688 
1689  sym%tenmat(1, 2) = -dot_product(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(2, 1:natoms)))
1690  sym%tenmat(1, 3) = -dot_product(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(3, 1:natoms)))
1691  sym%tenmat(2, 3) = -dot_product(sym%aw(1:natoms), (coord(2, 1:natoms)*coord(3, 1:natoms)))
1692 
1693  ! Symmetrize tensor matrix
1694  sym%tenmat(2, 1) = sym%tenmat(1, 2)
1695  sym%tenmat(3, 1) = sym%tenmat(1, 3)
1696  sym%tenmat(3, 2) = sym%tenmat(2, 3)
1697 
1698  ! Diagonalize the molecular tensor of inertia
1699  CALL jacobi(sym%tenmat(:, :), sym%tenval(:), sym%tenvec(:, :))
1700 
1701  ! Secure that the principal axes are right-handed
1702  sym%tenvec(:, 3) = vector_product(sym%tenvec(:, 1), sym%tenvec(:, 2))
1703 
1704  tt = sqrt(sym%tenval(1)**2 + sym%tenval(2)**2 + sym%tenval(3)**2)
1705  cpassert(tt /= 0)
1706 
1707  END SUBROUTINE tensor
1708 
1709 ! **************************************************************************************************
1710 !> \brief Transformation the Cartesian coodinates with the matrix tenvec.
1711 !> The coordinate axes can be switched according to the index
1712 !> vector idx. If idx(1) is equal to 3 instead to 1, then the first
1713 !> eigenvector (smallest eigenvalue) of the molecular tensor of
1714 !> inertia is connected to the z axis instead of the x axis. In
1715 !> addition the local atomic coordinate systems are initialized,
1716 !> if the symmetry is turned off.
1717 !> \param sym ...
1718 !> \param coord ...
1719 !> \param idx ...
1720 !> \par History
1721 !> Creation (16.11.98, Matthias Krack)
1722 !> \author jgh
1723 ! **************************************************************************************************
1724  SUBROUTINE tracar(sym, coord, idx)
1725  TYPE(molsym_type), INTENT(inout) :: sym
1726  REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1727  INTEGER, DIMENSION(3), INTENT(IN) :: idx
1728 
1729  INTEGER :: iatom, natom
1730  REAL(kind=dp), DIMENSION(3, 3) :: tenvec
1731 
1732  ! Transformation of the atomic coordinates ***
1733  natom = SIZE(coord, 2)
1734  tenvec = transpose(sym%tenvec)
1735  DO iatom = 1, natom
1736  coord(idx, iatom) = matmul(tenvec, coord(:, iatom))
1737  END DO
1738 
1739  ! Modify the transformation matrix for input orientation -> standard orientation
1740  sym%inptostd(idx, :) = tenvec
1741 
1742  END SUBROUTINE tracar
1743 
1744 ! **************************************************************************************************
1745 !> \brief Distance between two points
1746 !> \param a ...
1747 !> \param b ...
1748 !> \return ...
1749 !> \author jgh
1750 ! **************************************************************************************************
1751  FUNCTION dist(a, b) RESULT(d)
1752  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a, b
1753  REAL(kind=dp) :: d
1754 
1755  d = sqrt(sum((a - b)**2))
1756 
1757  END FUNCTION
1758 ! **************************************************************************************************
1759 !> \brief Write atomic coordinates to output unit iw.
1760 !> \param coord ...
1761 !> \param atype ...
1762 !> \param element ...
1763 !> \param z ...
1764 !> \param weight ...
1765 !> \param iw ...
1766 !> \date 08.05.2008
1767 !> \author JGH
1768 ! **************************************************************************************************
1769  SUBROUTINE write_coordinates(coord, atype, element, z, weight, iw)
1770  REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: coord
1771  INTEGER, DIMENSION(:), INTENT(in) :: atype
1772  CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: element
1773  INTEGER, DIMENSION(:), INTENT(in) :: z
1774  REAL(kind=dp), DIMENSION(:), INTENT(in) :: weight
1775  INTEGER, INTENT(IN) :: iw
1776 
1777  INTEGER :: iatom, natom
1778 
1779  IF (iw > 0) THEN
1780  natom = SIZE(coord, 2)
1781  WRITE (unit=iw, fmt="(T2,A)") &
1782  "MOLSYM|", &
1783  "MOLSYM| Atomic coordinates of the standard orientation in BOHR"
1784  WRITE (unit=iw, fmt="(T2,A,T37,A,T51,A,T65,A,T77,A)") &
1785  "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1786  DO iatom = 1, natom
1787  WRITE (unit=iw, fmt="(T2,A,I7,1X,I4,1X,A2,1X,I4,3(1X,F13.8),1X,F9.4)") &
1788  "MOLSYM|", iatom, atype(iatom), adjustl(element(iatom)), z(iatom), &
1789  coord(1:3, iatom), weight(iatom)
1790  END DO
1791  WRITE (unit=iw, fmt="(T2,A)") &
1792  "MOLSYM|", &
1793  "MOLSYM| Atomic coordinates of the standard orientation in ANGSTROM"
1794  WRITE (unit=iw, fmt="(T2,A,T37,A,T51,A,T65,A,T77,A)") &
1795  "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1796  DO iatom = 1, natom
1797  WRITE (unit=iw, fmt="(T2,A,I7,1X,I4,1X,A2,1X,I4,3(1X,F13.8),1X,F9.4)") &
1798  "MOLSYM|", iatom, atype(iatom), adjustl(element(iatom)), z(iatom), &
1799  coord(1:3, iatom)*angstrom, weight(iatom)
1800  END DO
1801  END IF
1802 
1803  END SUBROUTINE write_coordinates
1804 
1805 END MODULE molsym
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition: grid_common.h:153
tensor
struct tensor_ tensor
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::jacobi
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition: mathlib.F:1479
mathlib::reflect_vector
pure real(kind=dp) function, dimension(3), public reflect_vector(a, b)
Reflection of the vector a through a mirror plane defined by the normal vector b. The reflected vecto...
Definition: mathlib.F:1097
mathlib::angle
pure real(kind=dp) function, public angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition: mathlib.F:175
mathlib::build_rotmat
pure subroutine, public build_rotmat(phi, a, rotmat)
The rotation matrix rotmat which rotates a vector about a rotation axis defined by the vector a is bu...
Definition: mathlib.F:285
mathlib::rotate_vector
pure real(kind=dp) function, dimension(3), public rotate_vector(a, phi, b)
Rotation of the vector a about an rotation axis defined by the vector b. The rotation angle is phi (r...
Definition: mathlib.F:1135
mathlib::vector_product
pure real(kind=dp) function, dimension(3), public vector_product(a, b)
Calculation of the vector product c = a x b.
Definition: mathlib.F:1274
molsym
Molecular symmetry routines.
Definition: molsym.F:14
molsym::molecular_symmetry
subroutine, public molecular_symmetry(sym, eps_geo, coord, atype, weight)
Main program for symmetry analysis.
Definition: molsym.F:1075
molsym::release_molsym
subroutine, public release_molsym(sym)
release an object of molsym type
Definition: molsym.F:143
molsym::print_symmetry
subroutine, public print_symmetry(sym, coord, atype, element, z, weight, iw, plevel)
Print the output of the symmetry analysis.
Definition: molsym.F:1386
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
util
All kind of helpful little routines.
Definition: util.F:14