Line data Source code
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 : ! **************************************************************************************************
14 : MODULE molsym
15 :
16 : USE kinds, ONLY: dp
17 : USE mathconstants, ONLY: pi
18 : USE mathlib, ONLY: angle,&
19 : build_rotmat,&
20 : jacobi,&
21 : reflect_vector,&
22 : rotate_vector,&
23 : unit_matrix,&
24 : vector_product
25 : USE physcon, ONLY: angstrom
26 : USE util, ONLY: sort
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 : PRIVATE
31 :
32 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molsym'
34 :
35 : INTEGER, PARAMETER :: maxcn = 20, &
36 : maxsec = maxcn + 1, &
37 : maxses = 2*maxcn + 1, &
38 : maxsig = maxcn + 1, &
39 : maxsn = 2*maxcn
40 :
41 : PUBLIC :: molsym_type
42 : PUBLIC :: release_molsym, molecular_symmetry, print_symmetry
43 :
44 : ! **************************************************************************************************
45 : !> \brief Container for information about molecular symmetry
46 : !> \param ain : Atomic identification numbers (symmetry code).
47 : !> \param aw : Atomic weights for the symmetry analysis.
48 : !> \param eps_geo : Accuracy requested for analysis
49 : !> \param inptostd : Transformation matrix for input to standard orientation.
50 : !> \param point_group_symbol: Point group symbol.
51 : !> \param rotmat : Rotation matrix.
52 : !> \param sec : List of C axes
53 : !> (sec(:,i,j) => x,y,z of the ith j-fold C axis).
54 : !> \param ses : List of S axes
55 : !> (ses(:,i,j) => x,y,z of the ith j-fold S axis).
56 : !> \param sig : List of mirror planes
57 : !> (sig(:,i) => x,y,z of the ith mirror plane).
58 : !> \param center_of_mass : Shift vector for the center of mass.
59 : !> \param tenmat : Molecular tensor of inertia.
60 : !> \param tenval : Eigenvalues of the molecular tensor of inertia.
61 : !> \param tenvec : Eigenvectors of the molecular tensor of inertia.
62 : !> \param group_of : Group of equivalent atom.
63 : !> \param llequatom : Lower limit of a group in symequ_list.
64 : !> \param ncn : Degree of the C axis with highest degree.
65 : !> \param ndod : Number of substituted dodecahedral angles.
66 : !> \param nequatom : Number of equivalent atoms.
67 : !> \param ngroup : Number of groups of equivalent atoms.
68 : !> \param nico : Number of substituted icosahedral angles.
69 : !> \param nlin : Number of substituted angles of 180 degrees.
70 : !> \param nsec : Number of C axes.
71 : !> \param nses : Number of S axes.
72 : !> \param nsig : Number of mirror planes.
73 : !> \param nsn : Degree of the S axis with highest degree.
74 : !> \param ntet : Number of substituted tetrahedral angles.
75 : !> \param point_group_order : Group order.
76 : !> \param symequ_list : List of all atoms ordered in groups of equivalent atoms.
77 : !> \param ulequatom : Upper limit of a group in symequ_list.
78 : !> \param cubic : .TRUE., if a cubic point group was found (T,Th,Td,O,Oh,I,Ih).
79 : !> \param dgroup: .TRUE., if a point group of D symmetry was found (Dn,Dnh,Dnd)
80 : !> \param igroup: .TRUE., if a point group of icosahedral symmetry was found (I,Ih).
81 : !> \param invers: .TRUE., if the molecule has a center of inversion.
82 : !> \param linear: .TRUE., if the molecule is linear.
83 : !> \param maxis : .TRUE., if the molecule has a main axis.
84 : !> \param ogroup: .TRUE., if a point group of octahedral symmetry was found (O,Oh)
85 : !> \param sgroup: .TRUE., if a point group of S symmetry was found.
86 : !> \param sigmad: .TRUE., if there is a sigma_d mirror plane.
87 : !> \param sigmah: .TRUE., if there is a sigma_h mirror plane.
88 : !> \param sigmav: .TRUE., if there is a sigma_v mirror plane.
89 : !> \param tgroup: .TRUE., if a point group of tetrahedral symmetry was found (T,Th,Td).
90 : !> \par History
91 : !> 05.2008 created
92 : !> \author jgh
93 : ! **************************************************************************************************
94 : TYPE molsym_type
95 : CHARACTER(LEN=4) :: point_group_symbol = ""
96 : INTEGER :: point_group_order = -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 :
119 : CONTAINS
120 :
121 : ! **************************************************************************************************
122 : !> \brief Create an object of molsym type
123 : !> \param sym ...
124 : !> \param natoms ...
125 : !> \author jgh
126 : ! **************************************************************************************************
127 58 : SUBROUTINE create_molsym(sym, natoms)
128 : TYPE(molsym_type), POINTER :: sym
129 : INTEGER, INTENT(IN) :: natoms
130 :
131 58 : IF (ASSOCIATED(sym)) CALL release_molsym(sym)
132 :
133 479718 : ALLOCATE (sym)
134 :
135 : ALLOCATE (sym%ain(natoms), sym%aw(natoms), sym%group_of(natoms), sym%llequatom(natoms), &
136 870 : sym%nequatom(natoms), sym%symequ_list(natoms), sym%ulequatom(natoms))
137 :
138 58 : END SUBROUTINE create_molsym
139 :
140 : ! **************************************************************************************************
141 : !> \brief release an object of molsym type
142 : !> \param sym ...
143 : !> \author jgh
144 : ! **************************************************************************************************
145 58 : SUBROUTINE release_molsym(sym)
146 : TYPE(molsym_type), POINTER :: sym
147 :
148 58 : CPASSERT(ASSOCIATED(sym))
149 :
150 58 : IF (ASSOCIATED(sym%aw)) THEN
151 58 : DEALLOCATE (sym%aw)
152 : END IF
153 58 : IF (ASSOCIATED(sym%ain)) THEN
154 58 : DEALLOCATE (sym%ain)
155 : END IF
156 58 : IF (ASSOCIATED(sym%group_of)) THEN
157 58 : DEALLOCATE (sym%group_of)
158 : END IF
159 58 : IF (ASSOCIATED(sym%llequatom)) THEN
160 58 : DEALLOCATE (sym%llequatom)
161 : END IF
162 58 : IF (ASSOCIATED(sym%nequatom)) THEN
163 58 : DEALLOCATE (sym%nequatom)
164 : END IF
165 58 : IF (ASSOCIATED(sym%symequ_list)) THEN
166 58 : DEALLOCATE (sym%symequ_list)
167 : END IF
168 58 : IF (ASSOCIATED(sym%ulequatom)) THEN
169 58 : DEALLOCATE (sym%ulequatom)
170 : END IF
171 :
172 58 : DEALLOCATE (sym)
173 :
174 58 : 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 1419 : 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 1419 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
198 5676 : d(:) = a(:)/length_of_a
199 :
200 : ! Check, if the current Cn axis is already in the list
201 6384 : DO isec = 1, sym%nsec(n)
202 6088 : scapro = sym%sec(1, isec, n)*d(1) + sym%sec(2, isec, n)*d(2) + sym%sec(3, isec, n)*d(3)
203 6384 : IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) RETURN
204 : END DO
205 296 : sym%ncn = MAX(sym%ncn, n)
206 :
207 : ! Add the new Cn axis to the list sec
208 296 : CPASSERT(sym%nsec(n) < maxsec)
209 296 : sym%nsec(1) = sym%nsec(1) + 1
210 296 : sym%nsec(n) = sym%nsec(n) + 1
211 1184 : 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 284 : 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 284 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
235 1136 : d(:) = a(:)/length_of_a
236 :
237 : ! Check, if the current Sn axis is already in the list
238 1051 : DO ises = 1, sym%nses(n)
239 888 : scapro = sym%ses(1, ises, n)*d(1) + sym%ses(2, ises, n)*d(2) + sym%ses(3, ises, n)*d(3)
240 1051 : IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) RETURN
241 : END DO
242 163 : sym%nsn = MAX(sym%nsn, n)
243 :
244 : ! Add the new Sn axis to the list ses
245 163 : CPASSERT(sym%nses(n) < maxses)
246 163 : sym%nses(1) = sym%nses(1) + 1
247 163 : sym%nses(n) = sym%nses(n) + 1
248 652 : 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 597 : 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 597 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
270 2388 : d(:) = a(:)/length_of_a
271 :
272 : ! Check, if the normal vector of the current mirror plane is already in the list
273 2426 : DO isig = 1, sym%nsig
274 2273 : scapro = sym%sig(1, isig)*d(1) + sym%sig(2, isig)*d(2) + sym%sig(3, isig)*d(3)
275 2426 : 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 153 : CPASSERT(sym%nsig < maxsig)
280 153 : sym%nsig = sym%nsig + 1
281 612 : 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 1 : SUBROUTINE atomsym(sym)
293 : TYPE(molsym_type), INTENT(inout) :: sym
294 :
295 : ! Set point group symbol
296 :
297 1 : sym%point_group_symbol = "R(3)"
298 :
299 : ! Set variables like D*h
300 1 : sym%linear = .TRUE.
301 1 : sym%invers = .TRUE.
302 1 : sym%dgroup = .TRUE.
303 1 : sym%sigmah = .TRUE.
304 :
305 1 : 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 45 : 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 45 : IF ((sym%ncn == 2) .AND. (sym%nses(2) == 3)) THEN
327 2 : IF (sym%nsn == 4) THEN
328 : ! Special case: D2d
329 0 : phi = angle(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
330 0 : d(:) = vector_product(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
331 0 : CALL rotate_molecule(phi, d(:), sym, coord)
332 : ELSE
333 : ! Special cases: D2 and D2h
334 2 : phi = 0.5_dp*pi
335 2 : nx = naxis(sym%x_axis(:), coord, sym)
336 2 : ny = naxis(sym%y_axis(:), coord, sym)
337 2 : nz = naxis(sym%z_axis(:), coord, sym)
338 2 : IF ((nx > ny) .AND. (nx > nz)) THEN
339 0 : CALL rotate_molecule(-phi, sym%y_axis(:), sym, coord)
340 2 : ELSE IF ((ny > nz) .AND. (ny > nx)) THEN
341 0 : CALL rotate_molecule(phi, sym%x_axis(:), sym, coord)
342 : END IF
343 : END IF
344 : ELSE
345 43 : phi = angle(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
346 43 : d(:) = vector_product(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
347 43 : 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 45 : natoms = SIZE(coord, 2)
353 459 : DO iatom = 1, natoms
354 1656 : a(:) = coord(:, iatom)
355 459 : IF ((ABS(a(1)) > sym%eps_geo) .OR. (ABS(a(2)) > sym%eps_geo)) THEN
356 388 : a(3) = 0.0_dp
357 388 : IF (caxis(2, a(:), sym, coord)) CALL addsec(2, a(:), sym)
358 388 : d(:) = vector_product(a(:), sym%z_axis(:))
359 388 : IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
360 2283 : DO jatom = iatom + 1, natoms
361 7580 : b(:) = coord(:, jatom)
362 2309 : IF ((ABS(b(1)) > sym%eps_geo) .OR. (ABS(b(2)) > sym%eps_geo)) THEN
363 1871 : b(3) = 0.0_dp
364 1871 : phi = 0.5_dp*angle(a(:), b(:))
365 1871 : d(:) = vector_product(a(:), b(:))
366 1871 : b(:) = rotate_vector(a(:), phi, d(:))
367 1871 : IF (caxis(2, b(:), sym, coord)) CALL addsec(2, b(:), sym)
368 1871 : d(:) = vector_product(b(:), sym%z_axis)
369 1895 : 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 45 : IF (sigma(sym%z_axis(:), sym, coord)) THEN
377 14 : CALL addsig(sym%z_axis(:), sym)
378 14 : sym%sigmah = .TRUE.
379 : END IF
380 :
381 : ! Set logical variables for point group determination ***
382 45 : IF (sym%nsec(2) > 1) THEN
383 21 : sym%dgroup = .TRUE.
384 21 : IF ((.NOT. sym%sigmah) .AND. (sym%nsig > 0)) sym%sigmad = .TRUE.
385 : ELSE
386 24 : 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 24 : IF ((.NOT. sym%sigmah) .AND. (sym%nsn > 3)) sym%sgroup = .TRUE.
390 : END IF
391 :
392 : ! Rotate to standard orientation
393 45 : n = 0
394 45 : m = 0
395 45 : IF ((sym%dgroup .AND. sym%sigmah) .OR. (sym%dgroup .AND. sym%sigmad) .OR. sym%sigmav) THEN
396 : ! Dnh, Dnd or Cnv
397 133 : DO isig = 1, sym%nsig
398 133 : IF (ABS(sym%sig(3, isig)) < sym%eps_geo) THEN
399 105 : n = nsigma(sym%sig(:, isig), sym, coord)
400 105 : IF (n > m) THEN
401 21 : m = n
402 84 : a(:) = sym%sig(:, isig)
403 : END IF
404 : END IF
405 : END DO
406 21 : IF (m > 0) THEN
407 : ! Check for special case: C2v with all atoms in a plane
408 21 : IF (sym%sigmav .AND. (m == natoms)) THEN
409 1 : phi = angle(a(:), sym%x_axis(:))
410 1 : d(:) = vector_product(a(:), sym%x_axis(:))
411 : ELSE
412 20 : phi = angle(a(:), sym%y_axis(:))
413 20 : d(:) = vector_product(a(:), sym%y_axis(:))
414 : END IF
415 21 : CALL rotate_molecule(phi, d(:), sym, coord)
416 : END IF
417 24 : ELSE IF (sym%sigmah) THEN
418 : ! Cnh
419 81 : DO iatom = 1, natoms
420 81 : IF (ABS(coord(3, iatom)) < sym%eps_geo) THEN
421 70 : n = naxis(coord(:, iatom), coord, sym)
422 70 : IF (n > m) THEN
423 28 : m = n
424 28 : a(:) = coord(:, iatom)
425 : END IF
426 : END IF
427 : END DO
428 7 : IF (m > 0) THEN
429 7 : phi = angle(a(:), sym%x_axis(:))
430 7 : d(:) = vector_product(a(:), sym%x_axis(:))
431 7 : CALL rotate_molecule(phi, d(:), sym, coord)
432 : END IF
433 : ! No action for Dn, Cn or S2n ***
434 : END IF
435 :
436 45 : 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 58 : 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 58 : natoms = SIZE(coord, 2)
457 :
458 58 : IF (sym%sigmah) THEN
459 : incr = 1
460 : ELSE
461 41 : incr = 2
462 : END IF
463 :
464 : ! Initialize the atom and the group counter
465 58 : iatom = 1
466 58 : sym%ngroup = 1
467 :
468 : loop: DO
469 :
470 : ! Loop over all mirror planes
471 332 : DO isig = 1, sym%nsig
472 226 : a(:) = reflect_vector(coord(:, iatom), sym%sig(:, isig))
473 226 : iequatom = equatom(iatom, a(:), sym, coord)
474 332 : IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
475 111 : sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
476 111 : sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
477 111 : sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
478 : END IF
479 : END DO
480 :
481 : ! Loop over all Cn axes
482 445 : DO icn = 2, sym%ncn
483 829 : DO isec = 1, sym%nsec(icn)
484 1465 : DO jcn = 1, icn - 1
485 1126 : IF (newse(jcn, icn)) THEN
486 644 : phi = 2.0_dp*pi*REAL(jcn, KIND=dp)/REAL(icn, KIND=dp)
487 644 : a(:) = rotate_vector(coord(:, iatom), phi, sym%sec(:, isec, icn))
488 644 : iequatom = equatom(iatom, a(:), sym, coord)
489 644 : IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
490 328 : sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
491 328 : sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
492 328 : 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 339 : DO isn = 2, sym%nsn
501 442 : DO ises = 1, sym%nses(isn)
502 643 : DO jsn = 1, isn - 1, incr
503 410 : IF (newse(jsn, isn)) THEN
504 243 : phi = 2.0_dp*pi*REAL(jsn, KIND=dp)/REAL(isn, KIND=dp)
505 243 : a(:) = rotate_vector(coord(:, iatom), phi, sym%ses(:, ises, isn))
506 972 : a(:) = reflect_vector(a(:), sym%ses(:, ises, isn))
507 243 : iequatom = equatom(iatom, a(:), sym, coord)
508 243 : IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
509 30 : sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
510 30 : sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
511 30 : 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 106 : IF (sym%ulequatom(sym%ngroup) == natoms) EXIT
520 :
521 : ! Search for the next atom which is not in the list yet
522 156 : DO jatom = 2, natoms
523 98 : IF (.NOT. in_symequ_list(jatom, sym)) THEN
524 48 : iatom = jatom
525 48 : sym%ngroup = sym%ngroup + 1
526 48 : sym%llequatom(sym%ngroup) = sym%ulequatom(sym%ngroup - 1) + 1
527 48 : sym%ulequatom(sym%ngroup) = sym%llequatom(sym%ngroup)
528 48 : sym%symequ_list(sym%llequatom(sym%ngroup)) = iatom
529 48 : CYCLE loop
530 : END IF
531 : END DO
532 :
533 : END DO loop
534 :
535 : ! Generate list vector group_of
536 164 : DO i = 1, sym%ngroup
537 739 : DO iequatom = sym%llequatom(i), sym%ulequatom(i)
538 681 : sym%group_of(sym%symequ_list(iequatom)) = i
539 : END DO
540 : END DO
541 :
542 58 : 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 5543 : 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 5543 : caxis = .FALSE.
568 5543 : 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 5543 : natoms = SIZE(coord, 2)
572 5543 : IF (length_of_a > sym%eps_geo) THEN
573 : ! Calculate the rotation angle phi and build up the rotation matrix rotmat
574 5409 : phi = 2.0_dp*pi/REAL(n, dp)
575 5409 : CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
576 : ! Check all atoms
577 18830 : DO iatom = 1, natoms
578 : ! Rotate the actual atom by 2*pi/n about a
579 234507 : 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 18830 : 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 7 : 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 7 : phidd = ATAN(0.4_dp*SQRT(5.0_dp))
611 : ! Angle between two adjacent atoms of the icosahedron and the dodecahedron => <(C5,C3)
612 7 : phidi = ATAN(3.0_dp - SQRT(5.0_dp))
613 : ! Angle between two adjacent atoms of the icosahedron <(C5,C5)
614 7 : phiii = ATAN(2.0_dp)
615 :
616 : ! Find a pair of C3 axes
617 7 : natoms = SIZE(coord, 2)
618 146 : DO iatom = 1, natoms
619 : ! Check all atomic vectors for C3 axis
620 146 : IF (caxis(3, coord(:, iatom), sym, coord)) THEN
621 2 : CALL addsec(3, coord(:, iatom), sym)
622 2 : 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 7 : IF (sym%nsec(3) < 2) THEN
629 :
630 6 : loop: DO iatom = 1, natoms
631 24 : a(:) = coord(:, iatom)
632 24 : DO jatom = iatom + 1, natoms
633 734 : DO katom = jatom + 1, natoms
634 : IF ((ABS(dist(coord(:, iatom), coord(:, jatom)) &
635 716 : - dist(coord(:, iatom), coord(:, katom))) < sym%eps_geo) .AND. &
636 : (ABS(dist(coord(:, iatom), coord(:, jatom)) &
637 18 : - dist(coord(:, jatom), coord(:, katom))) < sym%eps_geo)) THEN
638 48 : b(:) = a(:) + coord(:, jatom) + coord(:, katom)
639 12 : IF (caxis(3, b(:), sym, coord)) THEN
640 12 : CALL addsec(3, b(:), sym)
641 12 : 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 7 : IF (sym%nsec(3) < 2) THEN
652 0 : failed = .TRUE.
653 0 : RETURN
654 : END IF
655 :
656 : ! Generate the remaining C3 axes
657 16 : nc3 = 0
658 : DO
659 16 : nc3 = sym%nsec(3)
660 82 : DO ic3 = 1, nc3
661 264 : a(:) = sym%sec(:, ic3, 3)
662 462 : DO jc3 = 1, nc3
663 446 : IF (ic3 /= jc3) THEN
664 1256 : d(:) = sym%sec(:, jc3, 3)
665 942 : DO i = 1, 2
666 628 : phi1 = 2.0_dp*REAL(i, KIND=dp)*pi/3.0_dp
667 628 : b(:) = rotate_vector(a(:), phi1, d(:))
668 942 : CALL addsec(3, b(:), sym)
669 : END DO
670 : END IF
671 : END DO
672 : END DO
673 :
674 16 : 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 7 : IF (sym%nsec(3) == 4) THEN
680 20 : a(:) = sym%sec(:, 1, 3)
681 20 : b(:) = sym%sec(:, 2, 3)
682 5 : phi1 = 0.5_dp*angle(a(:), b(:))
683 5 : IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
684 5 : d(:) = vector_product(a(:), b(:))
685 5 : b(:) = rotate_vector(a(:), phi1, d(:))
686 20 : c(:) = sym%sec(:, 3, 3)
687 5 : phi1 = 0.5_dp*angle(a(:), c(:))
688 5 : IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
689 5 : d(:) = vector_product(a(:), c(:))
690 5 : c(:) = rotate_vector(a(:), phi1, d(:))
691 5 : d(:) = vector_product(b(:), c(:))
692 5 : phi1 = 0.5_dp*phidd
693 5 : a(:) = rotate_vector(b(:), phi1, d(:))
694 5 : IF (caxis(3, a(:), sym, coord)) THEN
695 0 : CALL addsec(3, a(:), sym)
696 : ELSE
697 5 : phi2 = 0.5_dp*pi - phi1
698 5 : a(:) = rotate_vector(b(:), phi2, d(:))
699 5 : IF (caxis(3, a(:), sym, coord)) CALL addsec(3, a(:), sym)
700 : END IF
701 :
702 5 : IF (sym%nsec(3) > 4) CYCLE
703 :
704 2 : ELSE IF (sym%nsec(3) /= 10) THEN
705 :
706 : ! C3 axes found, but only 4 or 10 are possible
707 0 : failed = .TRUE.
708 0 : RETURN
709 :
710 : END IF
711 :
712 : ! Exit loop, if 4 or 10 C3 axes were found
713 9 : 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 47 : DO isec = 1, sym%nsec(3)
721 :
722 160 : a(:) = sym%sec(:, isec, 3)
723 167 : DO jsec = isec + 1, sym%nsec(3)
724 120 : phi1 = 0.5_dp*angle(a(:), sym%sec(:, jsec, 3))
725 120 : d(:) = vector_product(a(:), sym%sec(:, jsec, 3))
726 :
727 : ! Check for C5 axis (I,Ih)
728 120 : IF (sym%nsec(3) == 10) THEN
729 90 : b(:) = rotate_vector(a(:), phidi, d(:))
730 90 : IF (caxis(5, b(:), sym, coord)) THEN
731 14 : CALL addsec(5, b(:), sym)
732 14 : phi1 = phidi + phiii
733 14 : b(:) = rotate_vector(a(:), phi1, d(:))
734 14 : 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 360 : DO i = 0, 1
740 240 : phi2 = phi1 - 0.5_dp*REAL(i, KIND=dp)*pi
741 240 : b(:) = rotate_vector(a(:), phi2, d(:))
742 240 : IF (caxis(2, b(:), sym, coord)) THEN
743 134 : CALL addsec(2, b(:), sym)
744 134 : IF (sym%nsec(3) == 4) THEN
745 42 : IF (caxis(4, b(:), sym, coord)) CALL addsec(4, b(:), sym)
746 : END IF
747 134 : IF (saxis(2, b(:), sym, coord)) THEN
748 64 : CALL addses(2, b(:), sym)
749 64 : sym%invers = .TRUE.
750 : END IF
751 : END IF
752 360 : IF (sigma(b(:), sym, coord)) CALL addsig(b(:), sym)
753 : END DO
754 :
755 : ! Check for mirror plane
756 160 : 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 7 : iax = sym%ncn
764 :
765 7 : IF ((sym%ncn == 5) .AND. (sym%nsec(5) > 1)) THEN
766 2 : sym%igroup = .TRUE.
767 5 : ELSE IF ((sym%ncn == 4) .AND. (sym%nsec(4) > 1)) THEN
768 2 : sym%ogroup = .TRUE.
769 3 : ELSE IF ((sym%ncn == 3) .AND. (sym%nsec(3) > 1)) THEN
770 3 : sym%tgroup = .TRUE.
771 3 : IF ((.NOT. sym%invers) .AND. (sym%nsig > 0)) sym%sigmad = .TRUE.
772 : iax = 2
773 : ELSE
774 : ! No main axis found
775 0 : failed = .TRUE.
776 0 : RETURN
777 : END IF
778 :
779 : ! Rotate molecule to standard orientation
780 7 : phi1 = angle(sym%sec(:, 1, iax), sym%z_axis(:))
781 7 : d(:) = vector_product(sym%sec(:, 1, iax), sym%z_axis(:))
782 7 : CALL rotate_molecule(phi1, d(:), sym, coord)
783 :
784 28 : a(:) = sym%sec(:, 2, iax)
785 7 : a(3) = 0.0_dp
786 7 : phi2 = angle(a(:), sym%x_axis(:))
787 7 : d(:) = vector_product(a(:), sym%x_axis(:))
788 7 : 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 38170 : 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 38170 : equatom = 0
815 38170 : natoms = SIZE(coord, 2)
816 413408 : 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 400638 : (ABS(a(2) - coord(2, iatom)) < sym%eps_geo) .AND. &
820 12770 : (ABS(a(3) - coord(3, iatom)) < sym%eps_geo)) THEN
821 38170 : 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 55 : 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 55 : sym%point_group_order = 1 + sym%nsig
843 :
844 : ! Loop over all C axes
845 249 : DO icn = 2, sym%ncn
846 545 : DO isec = 1, sym%nsec(icn)
847 1021 : DO jcn = 1, icn - 1
848 827 : 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 55 : IF (sym%sigmah) THEN
855 : incr = 1
856 : ELSE
857 40 : incr = 2
858 : END IF
859 :
860 212 : DO isn = 2, sym%nsn
861 285 : DO ises = 1, sym%nses(isn)
862 452 : DO jsn = 1, isn - 1, incr
863 295 : IF (newse(jsn, isn)) sym%point_group_order = sym%point_group_order + 1
864 : END DO
865 : END DO
866 : END DO
867 :
868 55 : 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 57 : SUBROUTINE get_point_group_symbol(sym)
878 : TYPE(molsym_type), INTENT(inout) :: sym
879 :
880 : CHARACTER(LEN=4) :: degree
881 :
882 57 : IF (sym%linear) THEN
883 2 : IF (sym%invers) THEN
884 1 : sym%point_group_symbol = "D*h"
885 : ELSE
886 1 : sym%point_group_symbol = "C*v"
887 : END IF
888 55 : ELSE IF (sym%cubic) THEN
889 7 : IF (sym%igroup) THEN
890 2 : sym%point_group_symbol = "I"
891 5 : ELSE IF (sym%ogroup) THEN
892 2 : sym%point_group_symbol = "O"
893 3 : ELSE IF (sym%tgroup) THEN
894 3 : sym%point_group_symbol = "T"
895 : END IF
896 7 : IF (sym%invers) THEN
897 3 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
898 4 : ELSE IF (sym%sigmad) THEN
899 1 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"d"
900 : END IF
901 48 : ELSE IF (sym%maxis) THEN
902 45 : IF (sym%dgroup) THEN
903 21 : WRITE (degree, "(I3)") sym%ncn
904 21 : sym%point_group_symbol = "D"//TRIM(ADJUSTL(degree))
905 21 : IF (sym%sigmah) THEN
906 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
907 14 : ELSE IF (sym%sigmad) THEN
908 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"d"
909 : END IF
910 24 : ELSE IF (sym%sgroup) THEN
911 3 : WRITE (degree, "(I3)") sym%nsn
912 3 : sym%point_group_symbol = "S"//TRIM(ADJUSTL(degree))
913 : ELSE
914 21 : WRITE (degree, "(I3)") sym%ncn
915 21 : sym%point_group_symbol = "C"//TRIM(ADJUSTL(degree))
916 21 : IF (sym%sigmah) THEN
917 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
918 14 : ELSE IF (sym%sigmav) THEN
919 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"v"
920 : END IF
921 : END IF
922 : ELSE
923 3 : IF (sym%invers) THEN
924 1 : sym%point_group_symbol = "Ci"
925 2 : ELSE IF (sym%sigmah) THEN
926 1 : sym%point_group_symbol = "Cs"
927 : ELSE
928 1 : sym%point_group_symbol = "C1"
929 : END IF
930 : END IF
931 :
932 57 : 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 58 : 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 232 : sym%x_axis = (/1.0_dp, 0.0_dp, 0.0_dp/)
952 232 : sym%y_axis = (/0.0_dp, 1.0_dp, 0.0_dp/)
953 232 : sym%z_axis = (/0.0_dp, 0.0_dp, 1.0_dp/)
954 :
955 : ! Initialize lists for symmetry elements
956 93728 : sym%sec(:, :, :) = 0.0_dp
957 373288 : sym%ses(:, :, :) = 0.0_dp
958 4930 : sym%sig(:, :) = 0.0_dp
959 :
960 232 : sym%center_of_mass(:) = 0.0_dp
961 :
962 : ! Initialize symmetry element counters
963 58 : sym%ncn = 1
964 58 : sym%nsn = 1
965 1218 : sym%nsec(:) = 0
966 2378 : sym%nses(:) = 0
967 58 : sym%nsig = 0
968 :
969 : ! Initialize logical variables
970 58 : sym%cubic = .FALSE.
971 58 : sym%dgroup = .FALSE.
972 58 : sym%igroup = .FALSE.
973 58 : sym%invers = .FALSE.
974 58 : sym%linear = .FALSE.
975 58 : sym%maxis = .FALSE.
976 58 : sym%ogroup = .FALSE.
977 58 : sym%sgroup = .FALSE.
978 58 : sym%sigmad = .FALSE.
979 58 : sym%sigmah = .FALSE.
980 58 : sym%sigmav = .FALSE.
981 58 : sym%tgroup = .FALSE.
982 :
983 : ! Initialize list of equivalent atoms (C1 symmetry)
984 58 : natoms = SIZE(sym%nequatom)
985 58 : sym%ngroup = natoms
986 633 : sym%nequatom(:) = 1
987 633 : DO i = 1, sym%ngroup
988 575 : sym%group_of(i) = i
989 575 : sym%llequatom(i) = i
990 575 : sym%symequ_list(i) = i
991 633 : sym%ulequatom(i) = i
992 : END DO
993 :
994 58 : sym%point_group_order = -1
995 :
996 633 : DO iatom = 1, natoms
997 633 : sym%aw(iatom) = weight(iatom)
998 : END DO
999 :
1000 : ! Generate atomic identification numbers (symmetry code) ***
1001 633 : sym%ain(:) = REAL(atype(:), KIND=dp) + sym%aw(:)
1002 :
1003 : ! Initialize the transformation matrix for input orientation -> standard orientation
1004 58 : CALL unit_matrix(sym%inptostd(:, :))
1005 :
1006 58 : 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 1211 : 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 1211 : in_symequ_list = .FALSE.
1025 :
1026 7969 : DO iatom = 1, sym%ulequatom(sym%ngroup)
1027 7969 : IF (sym%symequ_list(iatom) == atoma) THEN
1028 1211 : 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 3 : 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 3 : IF (sym%nsn == 2) THEN
1051 : ! Ci
1052 1 : sym%invers = .TRUE.
1053 1 : phi = angle(sym%ses(:, 1, 2), sym%z_axis(:))
1054 1 : d(:) = vector_product(sym%ses(:, 1, 2), sym%z_axis(:))
1055 1 : CALL rotate_molecule(phi, d(:), sym, coord)
1056 2 : ELSE IF (sym%nsig == 1) THEN
1057 : ! Cs
1058 1 : sym%sigmah = .TRUE.
1059 1 : phi = angle(sym%sig(:, 1), sym%z_axis(:))
1060 1 : d(:) = vector_product(sym%sig(:, 1), sym%z_axis(:))
1061 1 : CALL rotate_molecule(phi, d(:), sym, coord)
1062 : END IF
1063 :
1064 3 : 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 58 : 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 58 : CALL timeset(routineN, handle)
1089 :
1090 : ! Perform memory allocation for the symmetry analysis
1091 58 : natoms = SIZE(coord, 2)
1092 58 : CALL create_molsym(sym, natoms)
1093 58 : sym%eps_geo = eps_geo
1094 :
1095 : ! Initialization of symmetry analysis
1096 58 : CALL init_symmetry(sym, atype, weight)
1097 :
1098 58 : IF (natoms == 1) THEN
1099 : ! Special case: only one atom
1100 1 : CALL atomsym(sym)
1101 : ELSE
1102 : ! Find molecular symmetry
1103 57 : CALL moleculesym(sym, coord)
1104 : ! Get point group and load point group table
1105 57 : CALL get_point_group_symbol(sym)
1106 : END IF
1107 :
1108 : ! Calculate group order
1109 58 : IF (.NOT. sym%linear) CALL get_point_group_order(sym)
1110 :
1111 : ! Generate a list of equivalent atoms
1112 58 : CALL build_symequ_list(sym, coord)
1113 :
1114 58 : CALL timestop(handle)
1115 :
1116 58 : 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 57 : 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 57 : eps_tenval = 0.01_dp*sym%eps_geo
1137 :
1138 : ! Calculate the molecular tensor of inertia
1139 57 : CALL tensor(sym, coord)
1140 : ! Use symmetry information from the eigenvalues of the molecular tensor of inertia
1141 57 : IF ((sym%tenval(3) - sym%tenval(1)) < eps_tenval) THEN
1142 : ! 0 < tenval(1) = tenval(2) = tenval(3)
1143 7 : sym%cubic = .TRUE.
1144 50 : 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 14 : IF (sym%tenval(1) < eps_tenval) sym%linear = .TRUE.
1148 14 : 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 36 : CALL tracar(sym, coord, (/1, 2, 3/))
1152 : END IF
1153 :
1154 : ! Continue with the new coordinate axes
1155 : DO
1156 57 : failed = .FALSE.
1157 57 : IF (sym%cubic) THEN
1158 7 : CALL cubsym(sym, coord, failed)
1159 7 : IF (failed) THEN
1160 0 : sym%cubic = .FALSE.
1161 0 : CYCLE
1162 : END IF
1163 :
1164 50 : ELSE IF (sym%linear) THEN
1165 :
1166 : ! Linear molecule
1167 2 : IF (saxis(2, sym%z_axis(:), sym, coord)) THEN
1168 1 : sym%invers = .TRUE.
1169 1 : sym%dgroup = .TRUE.
1170 1 : sym%sigmah = .TRUE.
1171 : END IF
1172 :
1173 : ELSE
1174 :
1175 : ! Check the new coordinate axes for Cn axes
1176 960 : DO icn = 2, maxcn
1177 912 : IF (caxis(icn, sym%z_axis(:), sym, coord)) CALL addsec(icn, sym%z_axis(:), sym)
1178 912 : IF (caxis(icn, sym%x_axis(:), sym, coord)) CALL addsec(icn, sym%x_axis(:), sym)
1179 960 : 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 1920 : DO isn = 2, maxsn
1184 1872 : IF (saxis(isn, sym%z_axis(:), sym, coord)) CALL addses(isn, sym%z_axis(:), sym)
1185 1872 : IF (saxis(isn, sym%x_axis(:), sym, coord)) CALL addses(isn, sym%x_axis(:), sym)
1186 1920 : 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 48 : IF (sigma(sym%z_axis(:), sym, coord)) CALL addsig(sym%z_axis(:), sym)
1191 48 : IF (sigma(sym%x_axis(:), sym, coord)) CALL addsig(sym%x_axis(:), sym)
1192 48 : IF (sigma(sym%y_axis(:), sym, coord)) CALL addsig(sym%y_axis(:), sym)
1193 :
1194 : ! There is a main axis (MAXIS = .TRUE.)
1195 48 : IF ((sym%ncn > 1) .OR. (sym%nsn > 3)) THEN
1196 45 : sym%maxis = .TRUE.
1197 45 : CALL axsym(coord, sym)
1198 : ELSE
1199 : ! Only low or no symmetry (Ci, Cs or C1)
1200 3 : CALL lowsym(sym, coord)
1201 : END IF
1202 :
1203 : END IF
1204 :
1205 57 : IF (.NOT. failed) EXIT
1206 :
1207 : END DO
1208 :
1209 : ! Find the remaining S axes
1210 57 : IF (.NOT. sym%linear) THEN
1211 249 : DO icn = 2, sym%ncn
1212 545 : DO isec = 1, sym%nsec(icn)
1213 296 : IF (saxis(icn, sym%sec(:, isec, icn), sym, coord)) &
1214 89 : CALL addses(icn, sym%sec(:, isec, icn), sym)
1215 296 : IF (saxis(2*icn, sym%sec(:, isec, icn), sym, coord)) &
1216 243 : 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 57 : IF (sym%nses(2) > 0) THEN
1223 16 : sym%nses(1) = sym%nses(1) - sym%nses(2)
1224 16 : sym%nses(2) = 0
1225 16 : CALL addses(2, sym%z_axis(:), sym)
1226 : END IF
1227 :
1228 57 : 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 76 : 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 76 : naxis = 0
1251 76 : natoms = SIZE(coord, 2)
1252 :
1253 76 : 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 76 : IF (length_of_a > sym%eps_geo) THEN
1257 :
1258 956 : DO iatom = 1, natoms
1259 3520 : b(:) = coord(:, iatom)
1260 880 : 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 956 : IF (length_of_b < sym%eps_geo) THEN
1263 0 : naxis = naxis + 1
1264 : ELSE
1265 3520 : a_norm = a(:)/length_of_a
1266 3520 : b_norm = b(:)/length_of_b
1267 880 : scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1268 880 : 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 76 : 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 1802 : FUNCTION newse(se1, se2)
1287 : INTEGER, INTENT(IN) :: se1, se2
1288 : LOGICAL :: newse
1289 :
1290 : INTEGER :: ise
1291 :
1292 1802 : newse = .TRUE.
1293 :
1294 3878 : DO ise = 2, MIN(se1, se2)
1295 3878 : IF ((MODULO(se1, ise) == 0) .AND. (MODULO(se2, ise) == 0)) THEN
1296 1802 : 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 105 : 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 105 : nsigma = 0
1325 :
1326 105 : 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 105 : IF (length_of_a > sym%eps_geo) THEN
1330 105 : natoms = SIZE(coord, 2)
1331 998 : DO iatom = 1, natoms
1332 3572 : b(:) = coord(:, iatom)
1333 893 : 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 998 : IF (length_of_b < sym%eps_geo) THEN
1336 0 : nsigma = nsigma + 1
1337 : ELSE
1338 3572 : a_norm = a(:)/length_of_a
1339 3572 : b_norm = b(:)/length_of_b
1340 893 : scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1341 893 : IF (ABS(scapro) < sym%eps_geo) nsigma = nsigma + 1
1342 : END IF
1343 : END DO
1344 : END IF
1345 :
1346 105 : 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 522 : 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 522 : IF (ABS(se(3)) < eps) THEN
1363 255 : IF (ABS(se(1)) < eps) THEN
1364 47 : se(2) = -se(2)
1365 : ELSE
1366 448 : IF (se(1) < 0.0_dp) se(1:2) = -se(1:2)
1367 : END IF
1368 : ELSE
1369 471 : IF (se(3) < 0.0_dp) se(:) = -se(:)
1370 : END IF
1371 :
1372 522 : 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 58 : 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 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: equatomlist, idxlist
1403 :
1404 : ! Print point group symbol and point group order
1405 : WRITE (iw, "(/,(T2,A))") &
1406 58 : "MOLSYM| Molecular symmetry information", &
1407 116 : "MOLSYM|"
1408 : WRITE (iw, "(T2,A,T77,A4)") &
1409 58 : "MOLSYM| Point group", ADJUSTR(sym%point_group_symbol)
1410 58 : IF (sym%point_group_order > -1) THEN
1411 : WRITE (iw, "(T2,A,T77,I4)") &
1412 55 : "MOLSYM| Group order ", sym%point_group_order
1413 : END IF
1414 :
1415 58 : IF (MOD(plevel, 10) == 1) THEN
1416 : WRITE (iw, "(T2,A)") &
1417 58 : "MOLSYM|", &
1418 116 : "MOLSYM| Groups of atoms equivalent by symmetry"
1419 : WRITE (iw, "(T2,A,T31,A)") &
1420 58 : "MOLSYM| Group number", "Group members (atomic indices)"
1421 164 : DO i = 1, sym%ngroup
1422 106 : nequatom = sym%ulequatom(i) - sym%llequatom(i) + 1
1423 106 : CPASSERT(nequatom > 0)
1424 424 : ALLOCATE (equatomlist(nequatom), idxlist(nequatom))
1425 681 : equatomlist(1:nequatom) = sym%symequ_list(sym%llequatom(i):sym%ulequatom(i))
1426 681 : idxlist = 0
1427 106 : 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 106 : "MOLSYM|", i, nequatom, (equatomlist(iequatom), iequatom=1, nequatom)
1430 164 : DEALLOCATE (equatomlist, idxlist)
1431 : END DO
1432 : END IF
1433 :
1434 58 : 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 58 : "MOLSYM|", &
1438 58 : "MOLSYM| Symmetry elements", &
1439 116 : "MOLSYM| Element number", "Type", "X", "Y", "Z"
1440 58 : secount = 1
1441 : WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A1,T36,3(1X,F14.8))") &
1442 58 : "MOLSYM|", secount, secount, "E", 0.0_dp, 0.0_dp, 0.0_dp
1443 : ! Mirror planes
1444 58 : string = "@ "
1445 211 : DO isig = 1, sym%nsig
1446 153 : secount = secount + 1
1447 153 : 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 670 : "MOLSYM|", secount, isig, string, (sym%sig(i, isig), i=1, 3)
1450 : END DO
1451 : ! C axes
1452 252 : DO icn = 2, sym%ncn
1453 194 : IF (icn < 10) THEN
1454 194 : WRITE (string, "(A1,I1)") "C", icn
1455 : ELSE
1456 0 : WRITE (string, "(A1,I2)") "C", icn
1457 : END IF
1458 548 : DO isec = 1, sym%nsec(icn)
1459 296 : secount = secount + 1
1460 296 : 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 1378 : "MOLSYM|", secount, isec, string, (sym%sec(i, isec, icn), i=1, 3)
1463 : END DO
1464 : END DO
1465 : ! S axes
1466 215 : DO isn = 2, sym%nsn
1467 157 : IF (isn == 2) THEN
1468 29 : WRITE (string, "(A3)") "i "
1469 128 : ELSE IF (isn < 10) THEN
1470 111 : WRITE (string, "(A1,I1)") "S", isn
1471 : ELSE
1472 17 : WRITE (string, "(A1,I2)") "S", isn
1473 : END IF
1474 288 : DO ises = 1, sym%nses(isn)
1475 73 : secount = secount + 1
1476 73 : 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 449 : "MOLSYM|", secount, ises, string, (sym%ses(i, ises, icn), i=1, 3)
1479 : END DO
1480 : END DO
1481 : END IF
1482 :
1483 58 : 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 57 : "MOLSYM|", &
1487 57 : "MOLSYM| Moments of the molecular inertia tensor [a.u.]", &
1488 114 : "MOLSYM|", "I(x)", "I(y)", "I(z)"
1489 57 : string = "xyz"
1490 741 : IF (MAXVAL(sym%tenmat) >= formatmaxval) THEN
1491 0 : DO i = 1, 3
1492 : WRITE (iw, "(T2,A,T32,A,T36,3(1X,ES14.4))") &
1493 0 : "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1494 : END DO
1495 : ELSE
1496 228 : DO i = 1, 3
1497 : WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1498 741 : "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1499 : END DO
1500 : END IF
1501 : WRITE (iw, "(T2,A)") &
1502 57 : "MOLSYM|", &
1503 114 : "MOLSYM| Principal moments and axes of inertia [a.u.]"
1504 741 : IF (MAXVAL(sym%tenmat) >= formatmaxval) THEN
1505 : WRITE (iw, "(T2,A,T36,3(1X,ES14.4))") &
1506 0 : "MOLSYM|", (sym%tenval(i), i=1, 3)
1507 : ELSE
1508 : WRITE (iw, "(T2,A,T36,3(1X,F14.8))") &
1509 228 : "MOLSYM|", (sym%tenval(i), i=1, 3)
1510 : END IF
1511 228 : DO i = 1, 3
1512 : WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1513 741 : "MOLSYM|", string(i:i), (sym%tenvec(i, j), j=1, 3)
1514 : END DO
1515 : END IF
1516 :
1517 58 : IF (MOD(plevel, 10) == 1) THEN
1518 : ! Print the Cartesian coordinates of the standard orientation
1519 58 : CALL write_coordinates(coord, atype, element, z, weight, iw)
1520 : END IF
1521 :
1522 58 : 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 87 : 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 87 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1546 87 : IF (length_of_a > sym%eps_geo) THEN
1547 :
1548 : ! Build up the rotation matrix
1549 42 : CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1550 :
1551 : ! Rotate the molecule by phi around a
1552 10689 : coord(:, :) = MATMUL(sym%rotmat(:, :), coord(:, :))
1553 :
1554 : ! Rotate the normal vectors of the mirror planes which are already found
1555 3339 : sym%sig(:, 1:sym%nsig) = MATMUL(sym%rotmat(:, :), sym%sig(:, 1:sym%nsig))
1556 :
1557 : ! Rotate the Cn axes which are already found
1558 186 : DO i = 2, sym%ncn
1559 6843 : 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 134 : DO i = 2, sym%nsn
1564 2150 : 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 2688 : sym%inptostd(:, :) = MATMUL(sym%rotmat(:, :), sym%inptostd(:, :))
1569 :
1570 : END IF
1571 :
1572 87 : 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 6344 : 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 6344 : saxis = .FALSE.
1598 :
1599 6344 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1600 :
1601 6344 : natoms = SIZE(coord, 2)
1602 :
1603 : ! Check the length of the axis vector a
1604 6344 : IF (length_of_a > sym%eps_geo) THEN
1605 : ! Calculate the rotation angle phi and build up the rotation matrix rotmat
1606 6344 : phi = 2.0_dp*pi/REAL(n, KIND=dp)
1607 6344 : CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1608 : ! Check all atoms
1609 9816 : DO iatom = 1, natoms
1610 : ! Rotate the actual atom by 2*pi/n about a
1611 124111 : b(:) = MATMUL(sym%rotmat(:, :), coord(:, iatom))
1612 : ! Reflect the coordinates of the rotated atom
1613 38188 : b(:) = reflect_vector(b(:), a(:))
1614 : ! Check if the coordinates of the rotated atom are in the coordinate set of the molecule
1615 9816 : 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 2808 : 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 2808 : sigma = .FALSE.
1645 :
1646 2808 : 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 2808 : IF (length_of_a > sym%eps_geo) THEN
1650 2674 : natoms = SIZE(coord, 2)
1651 10068 : DO iatom = 1, natoms
1652 : ! Reflect the actual atom
1653 9471 : b(:) = reflect_vector(coord(:, iatom), a(:))
1654 : ! Check if the coordinates of the reflected atom are in the coordinate set of the molecule
1655 10068 : 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 57 : 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 57 : natoms = SIZE(coord, 2)
1681 3269 : 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 2353 : coord(:, 1:natoms) = coord(:, 1:natoms) - SPREAD(sym%center_of_mass(:), 2, natoms)
1685 :
1686 : ! Build up the molecular tensor of inertia
1687 :
1688 631 : sym%tenmat(1, 1) = DOT_PRODUCT(sym%aw(1:natoms), (coord(2, 1:natoms)**2 + coord(3, 1:natoms)**2))
1689 631 : sym%tenmat(2, 2) = DOT_PRODUCT(sym%aw(1:natoms), (coord(3, 1:natoms)**2 + coord(1, 1:natoms)**2))
1690 631 : sym%tenmat(3, 3) = DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)**2 + coord(2, 1:natoms)**2))
1691 :
1692 631 : sym%tenmat(1, 2) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(2, 1:natoms)))
1693 631 : sym%tenmat(1, 3) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(3, 1:natoms)))
1694 631 : sym%tenmat(2, 3) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(2, 1:natoms)*coord(3, 1:natoms)))
1695 :
1696 : ! Symmetrize tensor matrix
1697 57 : sym%tenmat(2, 1) = sym%tenmat(1, 2)
1698 57 : sym%tenmat(3, 1) = sym%tenmat(1, 3)
1699 57 : sym%tenmat(3, 2) = sym%tenmat(2, 3)
1700 :
1701 : ! Diagonalize the molecular tensor of inertia
1702 57 : CALL jacobi(sym%tenmat(:, :), sym%tenval(:), sym%tenvec(:, :))
1703 :
1704 : ! Secure that the principal axes are right-handed
1705 228 : sym%tenvec(:, 3) = vector_product(sym%tenvec(:, 1), sym%tenvec(:, 2))
1706 :
1707 57 : tt = SQRT(sym%tenval(1)**2 + sym%tenval(2)**2 + sym%tenval(3)**2)
1708 57 : CPASSERT(tt /= 0)
1709 :
1710 57 : 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 50 : 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 50 : natom = SIZE(coord, 2)
1737 650 : tenvec = TRANSPOSE(sym%tenvec)
1738 482 : DO iatom = 1, natom
1739 3074 : coord(idx, iatom) = MATMUL(tenvec, coord(:, iatom))
1740 : END DO
1741 :
1742 : ! Modify the transformation matrix for input orientation -> standard orientation
1743 650 : sym%inptostd(idx, :) = tenvec
1744 :
1745 50 : END SUBROUTINE tracar
1746 :
1747 : ! **************************************************************************************************
1748 : !> \brief Distance between two points
1749 : !> \param a ...
1750 : !> \param b ...
1751 : !> \return ...
1752 : !> \author jgh
1753 : ! **************************************************************************************************
1754 2864 : FUNCTION dist(a, b) RESULT(d)
1755 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b
1756 : REAL(KIND=dp) :: d
1757 :
1758 11456 : d = SQRT(SUM((a - b)**2))
1759 :
1760 2864 : 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 58 : 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 58 : IF (iw > 0) THEN
1783 58 : natom = SIZE(coord, 2)
1784 : WRITE (UNIT=iw, FMT="(T2,A)") &
1785 58 : "MOLSYM|", &
1786 116 : "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 58 : "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1789 633 : 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 575 : "MOLSYM|", iatom, atype(iatom), ADJUSTL(element(iatom)), z(iatom), &
1792 1208 : coord(1:3, iatom), weight(iatom)
1793 : END DO
1794 : WRITE (UNIT=iw, FMT="(T2,A)") &
1795 58 : "MOLSYM|", &
1796 116 : "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 58 : "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1799 633 : 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 575 : "MOLSYM|", iatom, atype(iatom), ADJUSTL(element(iatom)), z(iatom), &
1802 2933 : coord(1:3, iatom)*angstrom, weight(iatom)
1803 : END DO
1804 : END IF
1805 :
1806 58 : END SUBROUTINE write_coordinates
1807 :
1808 0 : END MODULE molsym
|
