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