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