(git:b77b4be)
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-2025 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
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)) &
1214 CALL addses(icn, sym%sec(:, isec, icn), sym)
1215 IF (saxis(2*icn, sym%sec(:, isec, icn), sym, coord)) &
1216 CALL addses(2*icn, sym%sec(:, isec, icn), sym)
1217 END DO
1218 END DO
1219 END IF
1220
1221 ! Remove redundant S2 axes
1222 IF (sym%nses(2) > 0) THEN
1223 sym%nses(1) = sym%nses(1) - sym%nses(2)
1224 sym%nses(2) = 0
1225 CALL addses(2, sym%z_axis(:), sym)
1226 END IF
1227
1228 END SUBROUTINE moleculesym
1229
1230! **************************************************************************************************
1231!> \brief The number of atoms which are placed on an axis defined by the vector a is returned.
1232!> \param a ...
1233!> \param coord ...
1234!> \param sym ...
1235!> \return ...
1236!> \par History
1237!> Creation (16.11.98, Matthias Krack)
1238!> \author jgh
1239! **************************************************************************************************
1240 FUNCTION naxis(a, coord, sym)
1241 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1242 REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1243 TYPE(molsym_type), INTENT(inout) :: sym
1244 INTEGER :: naxis
1245
1246 INTEGER :: iatom, natoms
1247 REAL(kind=dp) :: length_of_a, length_of_b, scapro
1248 REAL(kind=dp), DIMENSION(3) :: a_norm, b, b_norm
1249
1250 naxis = 0
1251 natoms = SIZE(coord, 2)
1252
1253 length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1254
1255 ! Check the length of vector a
1256 IF (length_of_a > sym%eps_geo) THEN
1257
1258 DO iatom = 1, natoms
1259 b(:) = coord(:, iatom)
1260 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1261 ! An atom in the origin counts for each axis
1262 IF (length_of_b < sym%eps_geo) THEN
1263 naxis = naxis + 1
1264 ELSE
1265 a_norm = a(:)/length_of_a
1266 b_norm = b(:)/length_of_b
1267 scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1268 IF (abs(abs(scapro) - 1.0_dp) < sym%eps_geo) naxis = naxis + 1
1269 END IF
1270 END DO
1271
1272 END IF
1273
1274 END FUNCTION naxis
1275
1276! **************************************************************************************************
1277!> \brief Return .FALSE., if the symmetry elements se1 and se2 are a pair of
1278!> redundant symmetry elements.
1279!> \param se1 ...
1280!> \param se2 ...
1281!> \return ...
1282!> \par History
1283!> Creation (16.11.98, Matthias Krack)
1284!> \author jgh
1285! **************************************************************************************************
1286 FUNCTION newse(se1, se2)
1287 INTEGER, INTENT(IN) :: se1, se2
1288 LOGICAL :: newse
1289
1290 INTEGER :: ise
1291
1292 newse = .true.
1293
1294 DO ise = 2, min(se1, se2)
1295 IF ((modulo(se1, ise) == 0) .AND. (modulo(se2, ise) == 0)) THEN
1296 newse = .false.
1297 RETURN
1298 END IF
1299 END DO
1300
1301 END FUNCTION newse
1302
1303! **************************************************************************************************
1304!> \brief The number of atoms which are placed in a plane defined by the
1305!> the normal vector a is returned.
1306!> \param a ...
1307!> \param sym ...
1308!> \param coord ...
1309!> \return ...
1310!> \par History
1311!> Creation (16.11.98, Matthias Krack)
1312!> \author jgh
1313! **************************************************************************************************
1314 FUNCTION nsigma(a, sym, coord)
1315 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1316 TYPE(molsym_type), INTENT(inout) :: sym
1317 REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1318 INTEGER :: nsigma
1319
1320 INTEGER :: iatom, natoms
1321 REAL(kind=dp) :: length_of_a, length_of_b, scapro
1322 REAL(kind=dp), DIMENSION(3) :: a_norm, b, b_norm
1323
1324 nsigma = 0
1325
1326 length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1327
1328 ! Check the length of vector a
1329 IF (length_of_a > sym%eps_geo) THEN
1330 natoms = SIZE(coord, 2)
1331 DO iatom = 1, natoms
1332 b(:) = coord(:, iatom)
1333 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1334 ! An atom in the origin counts for each mirror plane
1335 IF (length_of_b < sym%eps_geo) THEN
1336 nsigma = nsigma + 1
1337 ELSE
1338 a_norm = a(:)/length_of_a
1339 b_norm = b(:)/length_of_b
1340 scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1341 IF (abs(scapro) < sym%eps_geo) nsigma = nsigma + 1
1342 END IF
1343 END DO
1344 END IF
1345
1346 END FUNCTION nsigma
1347
1348! **************************************************************************************************
1349!> \brief Style the output of the symmetry elements.
1350!> \param se ...
1351!> \param eps ...
1352!> \par History
1353!> Creation (16.11.98, Matthias Krack)
1354!> \author jgh
1355! **************************************************************************************************
1356 SUBROUTINE outse(se, eps)
1357 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: se
1358 REAL(kind=dp), INTENT(IN) :: eps
1359
1360! Positive z component for all vectors
1361
1362 IF (abs(se(3)) < eps) THEN
1363 IF (abs(se(1)) < eps) THEN
1364 se(2) = -se(2)
1365 ELSE
1366 IF (se(1) < 0.0_dp) se(1:2) = -se(1:2)
1367 END IF
1368 ELSE
1369 IF (se(3) < 0.0_dp) se(:) = -se(:)
1370 END IF
1371
1372 END SUBROUTINE outse
1373
1374! **************************************************************************************************
1375!> \brief Print the output of the symmetry analysis.
1376!> \param sym ...
1377!> \param coord ...
1378!> \param atype ...
1379!> \param element ...
1380!> \param z ...
1381!> \param weight ...
1382!> \param iw ...
1383!> \param plevel ...
1384!> \par History
1385!> Creation (16.11.98, Matthias Krack)
1386!> \author jgh
1387! **************************************************************************************************
1388 SUBROUTINE print_symmetry(sym, coord, atype, element, z, weight, iw, plevel)
1389 TYPE(molsym_type), INTENT(inout) :: sym
1390 REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: coord
1391 INTEGER, DIMENSION(:), INTENT(in) :: atype
1392 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: element
1393 INTEGER, DIMENSION(:), INTENT(in) :: z
1394 REAL(kind=dp), DIMENSION(:), INTENT(in) :: weight
1395 INTEGER, INTENT(IN) :: iw, plevel
1396
1397 REAL(kind=dp), PARAMETER :: formatmaxval = 1.0e5_dp
1398
1399 CHARACTER(LEN=3) :: string
1400 INTEGER :: i, icn, iequatom, isec, ises, isig, isn, &
1401 j, nequatom, secount
1402 INTEGER, ALLOCATABLE, DIMENSION(:) :: equatomlist, idxlist
1403
1404 ! Print point group symbol and point group order
1405 WRITE (iw, "(/,(T2,A))") &
1406 "MOLSYM| Molecular symmetry information", &
1407 "MOLSYM|"
1408 WRITE (iw, "(T2,A,T77,A4)") &
1409 "MOLSYM| Point group", adjustr(sym%point_group_symbol)
1410 IF (sym%point_group_order > -1) THEN
1411 WRITE (iw, "(T2,A,T77,I4)") &
1412 "MOLSYM| Group order ", sym%point_group_order
1413 END IF
1414
1415 IF (mod(plevel, 10) == 1) THEN
1416 WRITE (iw, "(T2,A)") &
1417 "MOLSYM|", &
1418 "MOLSYM| Groups of atoms equivalent by symmetry"
1419 WRITE (iw, "(T2,A,T31,A)") &
1420 "MOLSYM| Group number", "Group members (atomic indices)"
1421 DO i = 1, sym%ngroup
1422 nequatom = sym%ulequatom(i) - sym%llequatom(i) + 1
1423 cpassert(nequatom > 0)
1424 ALLOCATE (equatomlist(nequatom), idxlist(nequatom))
1425 equatomlist(1:nequatom) = sym%symequ_list(sym%llequatom(i):sym%ulequatom(i))
1426 idxlist = 0
1427 CALL sort(equatomlist, nequatom, idxlist)
1428 WRITE (iw, "(T2,A,T10,I5,T16,I6,T31,10(1X,I4),/,(T2,'MOLSYM|',T31,10(1X,I4)))") &
1429 "MOLSYM|", i, nequatom, (equatomlist(iequatom), iequatom=1, nequatom)
1430 DEALLOCATE (equatomlist, idxlist)
1431 END DO
1432 END IF
1433
1434 IF (mod(plevel/100, 10) == 1) THEN
1435 ! Print all symmetry elements
1436 WRITE (iw, "(T2,A,/,T2,A,/,T2,A,T31,A,T45,A,T60,A,T75,A)") &
1437 "MOLSYM|", &
1438 "MOLSYM| Symmetry elements", &
1439 "MOLSYM| Element number", "Type", "X", "Y", "Z"
1440 secount = 1
1441 WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A1,T36,3(1X,F14.8))") &
1442 "MOLSYM|", secount, secount, "E", 0.0_dp, 0.0_dp, 0.0_dp
1443 ! Mirror planes
1444 string = "@ "
1445 DO isig = 1, sym%nsig
1446 secount = secount + 1
1447 CALL outse(sym%sig(:, isig), sym%eps_geo)
1448 WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
1449 "MOLSYM|", secount, isig, string, (sym%sig(i, isig), i=1, 3)
1450 END DO
1451 ! C axes
1452 DO icn = 2, sym%ncn
1453 IF (icn < 10) THEN
1454 WRITE (string, "(A1,I1)") "C", icn
1455 ELSE
1456 WRITE (string, "(A1,I2)") "C", icn
1457 END IF
1458 DO isec = 1, sym%nsec(icn)
1459 secount = secount + 1
1460 CALL outse(sym%sec(:, isec, icn), sym%eps_geo)
1461 WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
1462 "MOLSYM|", secount, isec, string, (sym%sec(i, isec, icn), i=1, 3)
1463 END DO
1464 END DO
1465 ! S axes
1466 DO isn = 2, sym%nsn
1467 IF (isn == 2) THEN
1468 WRITE (string, "(A3)") "i "
1469 ELSE IF (isn < 10) THEN
1470 WRITE (string, "(A1,I1)") "S", isn
1471 ELSE
1472 WRITE (string, "(A1,I2)") "S", isn
1473 END IF
1474 DO ises = 1, sym%nses(isn)
1475 secount = secount + 1
1476 CALL outse(sym%ses(:, ises, isn), sym%eps_geo)
1477 WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
1478 "MOLSYM|", secount, ises, string, (sym%ses(i, ises, icn), i=1, 3)
1479 END DO
1480 END DO
1481 END IF
1482
1483 IF (mod(plevel/10, 10) == 1 .AND. SIZE(coord, 2) > 1) THEN
1484 ! Print the moments of the molecular inertia tensor
1485 WRITE (iw, "(T2,A,/,T2,A,/,T2,A,T43,A,T58,A,T73,A)") &
1486 "MOLSYM|", &
1487 "MOLSYM| Moments of the molecular inertia tensor [a.u.]", &
1488 "MOLSYM|", "I(x)", "I(y)", "I(z)"
1489 string = "xyz"
1490 IF (maxval(sym%tenmat) >= formatmaxval) THEN
1491 DO i = 1, 3
1492 WRITE (iw, "(T2,A,T32,A,T36,3(1X,ES14.4))") &
1493 "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1494 END DO
1495 ELSE
1496 DO i = 1, 3
1497 WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1498 "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1499 END DO
1500 END IF
1501 WRITE (iw, "(T2,A)") &
1502 "MOLSYM|", &
1503 "MOLSYM| Principal moments and axes of inertia [a.u.]"
1504 IF (maxval(sym%tenmat) >= formatmaxval) THEN
1505 WRITE (iw, "(T2,A,T36,3(1X,ES14.4))") &
1506 "MOLSYM|", (sym%tenval(i), i=1, 3)
1507 ELSE
1508 WRITE (iw, "(T2,A,T36,3(1X,F14.8))") &
1509 "MOLSYM|", (sym%tenval(i), i=1, 3)
1510 END IF
1511 DO i = 1, 3
1512 WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1513 "MOLSYM|", string(i:i), (sym%tenvec(i, j), j=1, 3)
1514 END DO
1515 END IF
1516
1517 IF (mod(plevel, 10) == 1) THEN
1518 ! Print the Cartesian coordinates of the standard orientation
1519 CALL write_coordinates(coord, atype, element, z, weight, iw)
1520 END IF
1521
1522 END SUBROUTINE print_symmetry
1523
1524! **************************************************************************************************
1525!> \brief Rotate the molecule about an axis defined by the vector a. The
1526!> rotation angle is phi (radians).
1527!> \param phi ...
1528!> \param a ...
1529!> \param sym ...
1530!> \param coord ...
1531!> \par History
1532!> Creation (16.11.98, Matthias Krack)
1533!> \author jgh
1534! **************************************************************************************************
1535 SUBROUTINE rotate_molecule(phi, a, sym, coord)
1536 REAL(kind=dp), INTENT(IN) :: phi
1537 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1538 TYPE(molsym_type), INTENT(inout) :: sym
1539 REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1540
1541 INTEGER :: i
1542 REAL(kind=dp) :: length_of_a
1543
1544 ! Check the length of vector a
1545 length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1546 IF (length_of_a > sym%eps_geo) THEN
1547
1548 ! Build up the rotation matrix
1549 CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1550
1551 ! Rotate the molecule by phi around a
1552 coord(:, :) = matmul(sym%rotmat(:, :), coord(:, :))
1553
1554 ! Rotate the normal vectors of the mirror planes which are already found
1555 sym%sig(:, 1:sym%nsig) = matmul(sym%rotmat(:, :), sym%sig(:, 1:sym%nsig))
1556
1557 ! Rotate the Cn axes which are already found
1558 DO i = 2, sym%ncn
1559 sym%sec(:, 1:sym%nsec(i), i) = matmul(sym%rotmat(:, :), sym%sec(:, 1:sym%nsec(i), i))
1560 END DO
1561
1562 ! Rotate the Sn axes which are already found
1563 DO i = 2, sym%nsn
1564 sym%ses(:, 1:sym%nses(i), i) = matmul(sym%rotmat(:, :), sym%ses(:, 1:sym%nses(i), i))
1565 END DO
1566
1567 ! Store current rotation
1568 sym%inptostd(:, :) = matmul(sym%rotmat(:, :), sym%inptostd(:, :))
1569
1570 END IF
1571
1572 END SUBROUTINE rotate_molecule
1573
1574! **************************************************************************************************
1575!> \brief Rotate the molecule about a n-fold axis defined by the vector a
1576!> using a rotation angle of 2*pi/n. Check, if the axis is a Cn axis.
1577!> \param n ...
1578!> \param a ...
1579!> \param sym ...
1580!> \param coord ...
1581!> \return ...
1582!> \par History
1583!> Creation (16.11.98, Matthias Krack)
1584!> \author jgh
1585! **************************************************************************************************
1586 FUNCTION saxis(n, a, sym, coord)
1587 INTEGER, INTENT(IN) :: n
1588 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1589 TYPE(molsym_type), INTENT(inout) :: sym
1590 REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1591 LOGICAL :: saxis
1592
1593 INTEGER :: iatom, natoms
1594 REAL(kind=dp) :: length_of_a, phi
1595 REAL(kind=dp), DIMENSION(3) :: b
1596
1597 saxis = .false.
1598
1599 length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1600
1601 natoms = SIZE(coord, 2)
1602
1603 ! Check the length of the axis vector a
1604 IF (length_of_a > sym%eps_geo) THEN
1605 ! Calculate the rotation angle phi and build up the rotation matrix rotmat
1606 phi = 2.0_dp*pi/real(n, kind=dp)
1607 CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1608 ! Check all atoms
1609 DO iatom = 1, natoms
1610 ! Rotate the actual atom by 2*pi/n about a
1611 b(:) = matmul(sym%rotmat(:, :), coord(:, iatom))
1612 ! Reflect the coordinates of the rotated atom
1613 b(:) = reflect_vector(b(:), a(:))
1614 ! Check if the coordinates of the rotated atom are in the coordinate set of the molecule
1615 IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
1616 END DO
1617 ! The axis defined by a is a Sn axis, thus return with saxis = .TRUE.
1618 saxis = .true.
1619 END IF
1620
1621 END FUNCTION saxis
1622
1623! **************************************************************************************************
1624!> \brief Reflect all atoms of the molecule through a mirror plane defined
1625!> by the normal vector a.
1626!> \param a ...
1627!> \param sym ...
1628!> \param coord ...
1629!> \return ...
1630!> \par History
1631!> Creation (16.11.98, Matthias Krack)
1632!> \author jgh
1633! **************************************************************************************************
1634 FUNCTION sigma(a, sym, coord)
1635 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1636 TYPE(molsym_type), INTENT(inout) :: sym
1637 REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1638 LOGICAL :: sigma
1639
1640 INTEGER :: iatom, natoms
1641 REAL(kind=dp) :: length_of_a
1642 REAL(kind=dp), DIMENSION(3) :: b
1643
1644 sigma = .false.
1645
1646 length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1647
1648 ! Check the length of vector a
1649 IF (length_of_a > sym%eps_geo) THEN
1650 natoms = SIZE(coord, 2)
1651 DO iatom = 1, natoms
1652 ! Reflect the actual atom
1653 b(:) = reflect_vector(coord(:, iatom), a(:))
1654 ! Check if the coordinates of the reflected atom are in the coordinate set of the molecule
1655 IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
1656 END DO
1657 ! The plane defined by a is a mirror plane, thus return with sigma = .TRUE.
1658 sigma = .true.
1659 END IF
1660
1661 END FUNCTION sigma
1662
1663! **************************************************************************************************
1664!> \brief Calculate the molecular tensor of inertia and the vector to the
1665!> center of mass of the molecule.
1666!> \param sym ...
1667!> \param coord ...
1668!> \par History
1669!> Creation (16.11.98, Matthias Krack)
1670!> \author jgh
1671! **************************************************************************************************
1672 SUBROUTINE tensor(sym, coord)
1673 TYPE(molsym_type), INTENT(inout) :: sym
1674 REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1675
1676 INTEGER :: natoms
1677 REAL(kind=dp) :: tt
1678
1679 ! Find the vector center_of_mass to molecular center of mass
1680 natoms = SIZE(coord, 2)
1681 sym%center_of_mass(:) = matmul(coord(:, 1:natoms), sym%aw(1:natoms))/sum(sym%aw(1:natoms))
1682
1683 ! Translate the center of mass of the molecule to the origin
1684 coord(:, 1:natoms) = coord(:, 1:natoms) - spread(sym%center_of_mass(:), 2, natoms)
1685
1686 ! Build up the molecular tensor of inertia
1687
1688 sym%tenmat(1, 1) = dot_product(sym%aw(1:natoms), (coord(2, 1:natoms)**2 + coord(3, 1:natoms)**2))
1689 sym%tenmat(2, 2) = dot_product(sym%aw(1:natoms), (coord(3, 1:natoms)**2 + coord(1, 1:natoms)**2))
1690 sym%tenmat(3, 3) = dot_product(sym%aw(1:natoms), (coord(1, 1:natoms)**2 + coord(2, 1:natoms)**2))
1691
1692 sym%tenmat(1, 2) = -dot_product(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(2, 1:natoms)))
1693 sym%tenmat(1, 3) = -dot_product(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(3, 1:natoms)))
1694 sym%tenmat(2, 3) = -dot_product(sym%aw(1:natoms), (coord(2, 1:natoms)*coord(3, 1:natoms)))
1695
1696 ! Symmetrize tensor matrix
1697 sym%tenmat(2, 1) = sym%tenmat(1, 2)
1698 sym%tenmat(3, 1) = sym%tenmat(1, 3)
1699 sym%tenmat(3, 2) = sym%tenmat(2, 3)
1700
1701 ! Diagonalize the molecular tensor of inertia
1702 CALL jacobi(sym%tenmat(:, :), sym%tenval(:), sym%tenvec(:, :))
1703
1704 ! Secure that the principal axes are right-handed
1705 sym%tenvec(:, 3) = vector_product(sym%tenvec(:, 1), sym%tenvec(:, 2))
1706
1707 tt = sqrt(sym%tenval(1)**2 + sym%tenval(2)**2 + sym%tenval(3)**2)
1708 cpassert(tt /= 0)
1709
1710 END SUBROUTINE tensor
1711
1712! **************************************************************************************************
1713!> \brief Transformation the Cartesian coodinates with the matrix tenvec.
1714!> The coordinate axes can be switched according to the index
1715!> vector idx. If idx(1) is equal to 3 instead to 1, then the first
1716!> eigenvector (smallest eigenvalue) of the molecular tensor of
1717!> inertia is connected to the z axis instead of the x axis. In
1718!> addition the local atomic coordinate systems are initialized,
1719!> if the symmetry is turned off.
1720!> \param sym ...
1721!> \param coord ...
1722!> \param idx ...
1723!> \par History
1724!> Creation (16.11.98, Matthias Krack)
1725!> \author jgh
1726! **************************************************************************************************
1727 SUBROUTINE tracar(sym, coord, idx)
1728 TYPE(molsym_type), INTENT(inout) :: sym
1729 REAL(kind=dp), DIMENSION(:, :), INTENT(inout) :: coord
1730 INTEGER, DIMENSION(3), INTENT(IN) :: idx
1731
1732 INTEGER :: iatom, natom
1733 REAL(kind=dp), DIMENSION(3, 3) :: tenvec
1734
1735 ! Transformation of the atomic coordinates ***
1736 natom = SIZE(coord, 2)
1737 tenvec = transpose(sym%tenvec)
1738 DO iatom = 1, natom
1739 coord(idx, iatom) = matmul(tenvec, coord(:, iatom))
1740 END DO
1741
1742 ! Modify the transformation matrix for input orientation -> standard orientation
1743 sym%inptostd(idx, :) = tenvec
1744
1745 END SUBROUTINE tracar
1746
1747! **************************************************************************************************
1748!> \brief Distance between two points
1749!> \param a ...
1750!> \param b ...
1751!> \return ...
1752!> \author jgh
1753! **************************************************************************************************
1754 FUNCTION dist(a, b) RESULT(d)
1755 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a, b
1756 REAL(kind=dp) :: d
1757
1758 d = sqrt(sum((a - b)**2))
1759
1760 END FUNCTION
1761! **************************************************************************************************
1762!> \brief Write atomic coordinates to output unit iw.
1763!> \param coord ...
1764!> \param atype ...
1765!> \param element ...
1766!> \param z ...
1767!> \param weight ...
1768!> \param iw ...
1769!> \date 08.05.2008
1770!> \author JGH
1771! **************************************************************************************************
1772 SUBROUTINE write_coordinates(coord, atype, element, z, weight, iw)
1773 REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: coord
1774 INTEGER, DIMENSION(:), INTENT(in) :: atype
1775 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: element
1776 INTEGER, DIMENSION(:), INTENT(in) :: z
1777 REAL(kind=dp), DIMENSION(:), INTENT(in) :: weight
1778 INTEGER, INTENT(IN) :: iw
1779
1780 INTEGER :: iatom, natom
1781
1782 IF (iw > 0) THEN
1783 natom = SIZE(coord, 2)
1784 WRITE (unit=iw, fmt="(T2,A)") &
1785 "MOLSYM|", &
1786 "MOLSYM| Atomic coordinates of the standard orientation in BOHR"
1787 WRITE (unit=iw, fmt="(T2,A,T37,A,T51,A,T65,A,T77,A)") &
1788 "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1789 DO iatom = 1, natom
1790 WRITE (unit=iw, fmt="(T2,A,I7,1X,I4,1X,A2,1X,I4,3(1X,F13.8),1X,F9.4)") &
1791 "MOLSYM|", iatom, atype(iatom), adjustl(element(iatom)), z(iatom), &
1792 coord(1:3, iatom), weight(iatom)
1793 END DO
1794 WRITE (unit=iw, fmt="(T2,A)") &
1795 "MOLSYM|", &
1796 "MOLSYM| Atomic coordinates of the standard orientation in ANGSTROM"
1797 WRITE (unit=iw, fmt="(T2,A,T37,A,T51,A,T65,A,T77,A)") &
1798 "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1799 DO iatom = 1, natom
1800 WRITE (unit=iw, fmt="(T2,A,I7,1X,I4,1X,A2,1X,I4,3(1X,F13.8),1X,F9.4)") &
1801 "MOLSYM|", iatom, atype(iatom), adjustl(element(iatom)), z(iatom), &
1802 coord(1:3, iatom)*angstrom, weight(iatom)
1803 END DO
1804 END IF
1805
1806 END SUBROUTINE write_coordinates
1807
1808END 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:120
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:156
tensor
struct tensor_ tensor
mathlib::unit_matrix
Definition mathlib.F:82
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:1468
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:1086
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:176
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:286
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:1124
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:1263
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:1389
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