Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Molecular symmetry routines
10 : !> \par History
11 : !> 2008 adapted from older routines by Matthias Krack
12 : !> \author jgh
13 : ! **************************************************************************************************
14 : MODULE molsym
15 :
16 : USE kinds, ONLY: dp
17 : USE mathconstants, ONLY: pi
18 : USE mathlib, ONLY: angle,&
19 : build_rotmat,&
20 : jacobi,&
21 : reflect_vector,&
22 : rotate_vector,&
23 : unit_matrix,&
24 : vector_product
25 : USE physcon, ONLY: angstrom
26 : USE util, ONLY: sort
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 : PRIVATE
31 :
32 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molsym'
34 :
35 : INTEGER, PARAMETER :: maxcn = 20, &
36 : maxsec = maxcn + 1, &
37 : maxses = 2*maxcn + 1, &
38 : maxsig = maxcn + 1, &
39 : maxsn = 2*maxcn
40 :
41 : PUBLIC :: molsym_type
42 : PUBLIC :: release_molsym, molecular_symmetry, print_symmetry
43 :
44 : ! **************************************************************************************************
45 : !> \brief Container for information about molecular symmetry
46 : !> \param ain : Atomic identification numbers (symmetry code).
47 : !> \param aw : Atomic weights for the symmetry analysis.
48 : !> \param eps_geo : Accuracy requested for analysis
49 : !> \param inptostd : Transformation matrix for input to standard orientation.
50 : !> \param point_group_symbol: Point group symbol.
51 : !> \param rotmat : Rotation matrix.
52 : !> \param sec : List of C axes
53 : !> (sec(:,i,j) => x,y,z of the ith j-fold C axis).
54 : !> \param ses : List of S axes
55 : !> (ses(:,i,j) => x,y,z of the ith j-fold S axis).
56 : !> \param sig : List of mirror planes
57 : !> (sig(:,i) => x,y,z of the ith mirror plane).
58 : !> \param center_of_mass : Shift vector for the center of mass.
59 : !> \param tenmat : Molecular tensor of inertia.
60 : !> \param tenval : Eigenvalues of the molecular tensor of inertia.
61 : !> \param tenvec : Eigenvectors of the molecular tensor of inertia.
62 : !> \param group_of : Group of equivalent atom.
63 : !> \param llequatom : Lower limit of a group in symequ_list.
64 : !> \param ncn : Degree of the C axis with highest degree.
65 : !> \param ndod : Number of substituted dodecahedral angles.
66 : !> \param nequatom : Number of equivalent atoms.
67 : !> \param ngroup : Number of groups of equivalent atoms.
68 : !> \param nico : Number of substituted icosahedral angles.
69 : !> \param nlin : Number of substituted angles of 180 degrees.
70 : !> \param nsec : Number of C axes.
71 : !> \param nses : Number of S axes.
72 : !> \param nsig : Number of mirror planes.
73 : !> \param nsn : Degree of the S axis with highest degree.
74 : !> \param ntet : Number of substituted tetrahedral angles.
75 : !> \param point_group_order : Group order.
76 : !> \param symequ_list : List of all atoms ordered in groups of equivalent atoms.
77 : !> \param ulequatom : Upper limit of a group in symequ_list.
78 : !> \param cubic : .TRUE., if a cubic point group was found (T,Th,Td,O,Oh,I,Ih).
79 : !> \param dgroup: .TRUE., if a point group of D symmetry was found (Dn,Dnh,Dnd)
80 : !> \param igroup: .TRUE., if a point group of icosahedral symmetry was found (I,Ih).
81 : !> \param invers: .TRUE., if the molecule has a center of inversion.
82 : !> \param linear: .TRUE., if the molecule is linear.
83 : !> \param maxis : .TRUE., if the molecule has a main axis.
84 : !> \param ogroup: .TRUE., if a point group of octahedral symmetry was found (O,Oh)
85 : !> \param sgroup: .TRUE., if a point group of S symmetry was found.
86 : !> \param sigmad: .TRUE., if there is a sigma_d mirror plane.
87 : !> \param sigmah: .TRUE., if there is a sigma_h mirror plane.
88 : !> \param sigmav: .TRUE., if there is a sigma_v mirror plane.
89 : !> \param tgroup: .TRUE., if a point group of tetrahedral symmetry was found (T,Th,Td).
90 : !> \par History
91 : !> 05.2008 created
92 : !> \author jgh
93 : ! **************************************************************************************************
94 : TYPE molsym_type
95 : CHARACTER(LEN=4) :: point_group_symbol
96 : INTEGER :: point_group_order
97 : INTEGER :: ncn, ndod, ngroup, nico, nlin, nsig, nsn, ntet
98 : LOGICAL :: cubic, dgroup, igroup, invers, linear, maxis, &
99 : ogroup, sgroup, sigmad, sigmah, sigmav, tgroup
100 : REAL(KIND=dp) :: eps_geo
101 : REAL(KIND=dp), DIMENSION(3) :: center_of_mass, tenval
102 : REAL(KIND=dp), DIMENSION(3) :: x_axis, y_axis, z_axis
103 : REAL(KIND=dp), DIMENSION(:), POINTER :: ain, aw
104 : REAL(KIND=dp), DIMENSION(3, 3) :: inptostd, rotmat, tenmat, tenvec
105 : REAL(KIND=dp), DIMENSION(3, maxsig) :: sig
106 : REAL(KIND=dp), DIMENSION(3, maxsec, 2:maxcn) :: sec
107 : REAL(KIND=dp), DIMENSION(3, maxses, 2:maxsn) :: ses
108 : INTEGER, DIMENSION(maxcn) :: nsec
109 : INTEGER, DIMENSION(maxsn) :: nses
110 : INTEGER, DIMENSION(:), POINTER :: group_of, llequatom, nequatom, &
111 : symequ_list, ulequatom
112 : END TYPE molsym_type
113 :
114 : ! **************************************************************************************************
115 :
116 : CONTAINS
117 :
118 : ! **************************************************************************************************
119 : !> \brief Create an object of molsym type
120 : !> \param sym ...
121 : !> \param natoms ...
122 : !> \author jgh
123 : ! **************************************************************************************************
124 58 : SUBROUTINE create_molsym(sym, natoms)
125 : TYPE(molsym_type), POINTER :: sym
126 : INTEGER, INTENT(IN) :: natoms
127 :
128 58 : IF (ASSOCIATED(sym)) CALL release_molsym(sym)
129 :
130 58 : ALLOCATE (sym)
131 :
132 : ALLOCATE (sym%ain(natoms), sym%aw(natoms), sym%group_of(natoms), sym%llequatom(natoms), &
133 870 : sym%nequatom(natoms), sym%symequ_list(natoms), sym%ulequatom(natoms))
134 :
135 58 : END SUBROUTINE create_molsym
136 :
137 : ! **************************************************************************************************
138 : !> \brief release an object of molsym type
139 : !> \param sym ...
140 : !> \author jgh
141 : ! **************************************************************************************************
142 58 : SUBROUTINE release_molsym(sym)
143 : TYPE(molsym_type), POINTER :: sym
144 :
145 58 : CPASSERT(ASSOCIATED(sym))
146 :
147 58 : IF (ASSOCIATED(sym%aw)) THEN
148 58 : DEALLOCATE (sym%aw)
149 : END IF
150 58 : IF (ASSOCIATED(sym%ain)) THEN
151 58 : DEALLOCATE (sym%ain)
152 : END IF
153 58 : IF (ASSOCIATED(sym%group_of)) THEN
154 58 : DEALLOCATE (sym%group_of)
155 : END IF
156 58 : IF (ASSOCIATED(sym%llequatom)) THEN
157 58 : DEALLOCATE (sym%llequatom)
158 : END IF
159 58 : IF (ASSOCIATED(sym%nequatom)) THEN
160 58 : DEALLOCATE (sym%nequatom)
161 : END IF
162 58 : IF (ASSOCIATED(sym%symequ_list)) THEN
163 58 : DEALLOCATE (sym%symequ_list)
164 : END IF
165 58 : IF (ASSOCIATED(sym%ulequatom)) THEN
166 58 : DEALLOCATE (sym%ulequatom)
167 : END IF
168 :
169 58 : DEALLOCATE (sym)
170 :
171 58 : 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 1419 : 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 1419 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
195 5676 : d(:) = a(:)/length_of_a
196 :
197 : ! Check, if the current Cn axis is already in the list
198 6384 : DO isec = 1, sym%nsec(n)
199 6088 : scapro = sym%sec(1, isec, n)*d(1) + sym%sec(2, isec, n)*d(2) + sym%sec(3, isec, n)*d(3)
200 6384 : IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) RETURN
201 : END DO
202 296 : sym%ncn = MAX(sym%ncn, n)
203 :
204 : ! Add the new Cn axis to the list sec
205 296 : CPASSERT(sym%nsec(n) < maxsec)
206 296 : sym%nsec(1) = sym%nsec(1) + 1
207 296 : sym%nsec(n) = sym%nsec(n) + 1
208 1184 : 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 284 : 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 284 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
232 1136 : d(:) = a(:)/length_of_a
233 :
234 : ! Check, if the current Sn axis is already in the list
235 1051 : DO ises = 1, sym%nses(n)
236 888 : scapro = sym%ses(1, ises, n)*d(1) + sym%ses(2, ises, n)*d(2) + sym%ses(3, ises, n)*d(3)
237 1051 : IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) RETURN
238 : END DO
239 163 : sym%nsn = MAX(sym%nsn, n)
240 :
241 : ! Add the new Sn axis to the list ses
242 163 : CPASSERT(sym%nses(n) < maxses)
243 163 : sym%nses(1) = sym%nses(1) + 1
244 163 : sym%nses(n) = sym%nses(n) + 1
245 652 : 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 597 : 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 597 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
267 2388 : d(:) = a(:)/length_of_a
268 :
269 : ! Check, if the normal vector of the current mirror plane is already in the list
270 2426 : DO isig = 1, sym%nsig
271 2273 : scapro = sym%sig(1, isig)*d(1) + sym%sig(2, isig)*d(2) + sym%sig(3, isig)*d(3)
272 2426 : 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 153 : CPASSERT(sym%nsig < maxsig)
277 153 : sym%nsig = sym%nsig + 1
278 612 : 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 1 : SUBROUTINE atomsym(sym)
290 : TYPE(molsym_type), INTENT(inout) :: sym
291 :
292 : ! Set point group symbol
293 :
294 1 : sym%point_group_symbol = "R(3)"
295 :
296 : ! Set variables like D*h
297 1 : sym%linear = .TRUE.
298 1 : sym%invers = .TRUE.
299 1 : sym%dgroup = .TRUE.
300 1 : sym%sigmah = .TRUE.
301 :
302 1 : 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 45 : 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 45 : IF ((sym%ncn == 2) .AND. (sym%nses(2) == 3)) THEN
324 2 : IF (sym%nsn == 4) THEN
325 : ! Special case: D2d
326 0 : phi = angle(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
327 0 : d(:) = vector_product(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
328 0 : CALL rotate_molecule(phi, d(:), sym, coord)
329 : ELSE
330 : ! Special cases: D2 and D2h
331 2 : phi = 0.5_dp*pi
332 2 : nx = naxis(sym%x_axis(:), coord, sym)
333 2 : ny = naxis(sym%y_axis(:), coord, sym)
334 2 : nz = naxis(sym%z_axis(:), coord, sym)
335 2 : IF ((nx > ny) .AND. (nx > nz)) THEN
336 0 : CALL rotate_molecule(-phi, sym%y_axis(:), sym, coord)
337 2 : ELSE IF ((ny > nz) .AND. (ny > nx)) THEN
338 0 : CALL rotate_molecule(phi, sym%x_axis(:), sym, coord)
339 : END IF
340 : END IF
341 : ELSE
342 43 : phi = angle(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
343 43 : d(:) = vector_product(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
344 43 : 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 45 : natoms = SIZE(coord, 2)
350 459 : DO iatom = 1, natoms
351 1656 : a(:) = coord(:, iatom)
352 459 : IF ((ABS(a(1)) > sym%eps_geo) .OR. (ABS(a(2)) > sym%eps_geo)) THEN
353 388 : a(3) = 0.0_dp
354 388 : IF (caxis(2, a(:), sym, coord)) CALL addsec(2, a(:), sym)
355 388 : d(:) = vector_product(a(:), sym%z_axis(:))
356 388 : IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
357 2283 : DO jatom = iatom + 1, natoms
358 7580 : b(:) = coord(:, jatom)
359 2309 : IF ((ABS(b(1)) > sym%eps_geo) .OR. (ABS(b(2)) > sym%eps_geo)) THEN
360 1871 : b(3) = 0.0_dp
361 1871 : phi = 0.5_dp*angle(a(:), b(:))
362 1871 : d(:) = vector_product(a(:), b(:))
363 1871 : b(:) = rotate_vector(a(:), phi, d(:))
364 1871 : IF (caxis(2, b(:), sym, coord)) CALL addsec(2, b(:), sym)
365 1871 : d(:) = vector_product(b(:), sym%z_axis)
366 1895 : 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 45 : IF (sigma(sym%z_axis(:), sym, coord)) THEN
374 14 : CALL addsig(sym%z_axis(:), sym)
375 14 : sym%sigmah = .TRUE.
376 : END IF
377 :
378 : ! Set logical variables for point group determination ***
379 45 : IF (sym%nsec(2) > 1) THEN
380 21 : sym%dgroup = .TRUE.
381 21 : IF ((.NOT. sym%sigmah) .AND. (sym%nsig > 0)) sym%sigmad = .TRUE.
382 : ELSE
383 24 : 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 24 : IF ((.NOT. sym%sigmah) .AND. (sym%nsn > 3)) sym%sgroup = .TRUE.
387 : END IF
388 :
389 : ! Rotate to standard orientation
390 45 : n = 0
391 45 : m = 0
392 45 : IF ((sym%dgroup .AND. sym%sigmah) .OR. (sym%dgroup .AND. sym%sigmad) .OR. sym%sigmav) THEN
393 : ! Dnh, Dnd or Cnv
394 133 : DO isig = 1, sym%nsig
395 133 : IF (ABS(sym%sig(3, isig)) < sym%eps_geo) THEN
396 105 : n = nsigma(sym%sig(:, isig), sym, coord)
397 105 : IF (n > m) THEN
398 21 : m = n
399 84 : a(:) = sym%sig(:, isig)
400 : END IF
401 : END IF
402 : END DO
403 21 : IF (m > 0) THEN
404 : ! Check for special case: C2v with all atoms in a plane
405 21 : IF (sym%sigmav .AND. (m == natoms)) THEN
406 1 : phi = angle(a(:), sym%x_axis(:))
407 1 : d(:) = vector_product(a(:), sym%x_axis(:))
408 : ELSE
409 20 : phi = angle(a(:), sym%y_axis(:))
410 20 : d(:) = vector_product(a(:), sym%y_axis(:))
411 : END IF
412 21 : CALL rotate_molecule(phi, d(:), sym, coord)
413 : END IF
414 24 : ELSE IF (sym%sigmah) THEN
415 : ! Cnh
416 81 : DO iatom = 1, natoms
417 81 : IF (ABS(coord(3, iatom)) < sym%eps_geo) THEN
418 70 : n = naxis(coord(:, iatom), coord, sym)
419 70 : IF (n > m) THEN
420 28 : m = n
421 28 : a(:) = coord(:, iatom)
422 : END IF
423 : END IF
424 : END DO
425 7 : IF (m > 0) THEN
426 7 : phi = angle(a(:), sym%x_axis(:))
427 7 : d(:) = vector_product(a(:), sym%x_axis(:))
428 7 : CALL rotate_molecule(phi, d(:), sym, coord)
429 : END IF
430 : ! No action for Dn, Cn or S2n ***
431 : END IF
432 :
433 45 : 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 58 : 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 58 : natoms = SIZE(coord, 2)
454 :
455 58 : IF (sym%sigmah) THEN
456 : incr = 1
457 : ELSE
458 41 : incr = 2
459 : END IF
460 :
461 : ! Initialize the atom and the group counter
462 58 : iatom = 1
463 58 : sym%ngroup = 1
464 :
465 : loop: DO
466 :
467 : ! Loop over all mirror planes
468 332 : DO isig = 1, sym%nsig
469 226 : a(:) = reflect_vector(coord(:, iatom), sym%sig(:, isig))
470 226 : iequatom = equatom(iatom, a(:), sym, coord)
471 332 : IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
472 111 : sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
473 111 : sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
474 111 : sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
475 : END IF
476 : END DO
477 :
478 : ! Loop over all Cn axes
479 445 : DO icn = 2, sym%ncn
480 829 : DO isec = 1, sym%nsec(icn)
481 1465 : DO jcn = 1, icn - 1
482 1126 : IF (newse(jcn, icn)) THEN
483 644 : phi = 2.0_dp*pi*REAL(jcn, KIND=dp)/REAL(icn, KIND=dp)
484 644 : a(:) = rotate_vector(coord(:, iatom), phi, sym%sec(:, isec, icn))
485 644 : iequatom = equatom(iatom, a(:), sym, coord)
486 644 : IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
487 328 : sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
488 328 : sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
489 328 : 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 339 : DO isn = 2, sym%nsn
498 442 : DO ises = 1, sym%nses(isn)
499 643 : DO jsn = 1, isn - 1, incr
500 410 : IF (newse(jsn, isn)) THEN
501 243 : phi = 2.0_dp*pi*REAL(jsn, KIND=dp)/REAL(isn, KIND=dp)
502 243 : a(:) = rotate_vector(coord(:, iatom), phi, sym%ses(:, ises, isn))
503 972 : a(:) = reflect_vector(a(:), sym%ses(:, ises, isn))
504 243 : iequatom = equatom(iatom, a(:), sym, coord)
505 243 : IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
506 30 : sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
507 30 : sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
508 30 : 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 106 : IF (sym%ulequatom(sym%ngroup) == natoms) EXIT
517 :
518 : ! Search for the next atom which is not in the list yet
519 156 : DO jatom = 2, natoms
520 98 : IF (.NOT. in_symequ_list(jatom, sym)) THEN
521 48 : iatom = jatom
522 48 : sym%ngroup = sym%ngroup + 1
523 48 : sym%llequatom(sym%ngroup) = sym%ulequatom(sym%ngroup - 1) + 1
524 48 : sym%ulequatom(sym%ngroup) = sym%llequatom(sym%ngroup)
525 48 : sym%symequ_list(sym%llequatom(sym%ngroup)) = iatom
526 48 : CYCLE loop
527 : END IF
528 : END DO
529 :
530 : END DO loop
531 :
532 : ! Generate list vector group_of
533 164 : DO i = 1, sym%ngroup
534 739 : DO iequatom = sym%llequatom(i), sym%ulequatom(i)
535 681 : sym%group_of(sym%symequ_list(iequatom)) = i
536 : END DO
537 : END DO
538 :
539 58 : 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 5543 : 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 5543 : caxis = .FALSE.
565 5543 : 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 5543 : natoms = SIZE(coord, 2)
569 5543 : IF (length_of_a > sym%eps_geo) THEN
570 : ! Calculate the rotation angle phi and build up the rotation matrix rotmat
571 5409 : phi = 2.0_dp*pi/REAL(n, dp)
572 5409 : CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
573 : ! Check all atoms
574 18830 : DO iatom = 1, natoms
575 : ! Rotate the actual atom by 2*pi/n about a
576 234507 : 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 18830 : 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 7 : 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 7 : phidd = ATAN(0.4_dp*SQRT(5.0_dp))
608 : ! Angle between two adjacent atoms of the icosahedron and the dodecahedron => <(C5,C3)
609 7 : phidi = ATAN(3.0_dp - SQRT(5.0_dp))
610 : ! Angle between two adjacent atoms of the icosahedron <(C5,C5)
611 7 : phiii = ATAN(2.0_dp)
612 :
613 : ! Find a pair of C3 axes
614 7 : natoms = SIZE(coord, 2)
615 146 : DO iatom = 1, natoms
616 : ! Check all atomic vectors for C3 axis
617 146 : IF (caxis(3, coord(:, iatom), sym, coord)) THEN
618 2 : CALL addsec(3, coord(:, iatom), sym)
619 2 : 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 7 : IF (sym%nsec(3) < 2) THEN
626 :
627 6 : loop: DO iatom = 1, natoms
628 24 : a(:) = coord(:, iatom)
629 24 : DO jatom = iatom + 1, natoms
630 734 : DO katom = jatom + 1, natoms
631 : IF ((ABS(dist(coord(:, iatom), coord(:, jatom)) &
632 716 : - dist(coord(:, iatom), coord(:, katom))) < sym%eps_geo) .AND. &
633 : (ABS(dist(coord(:, iatom), coord(:, jatom)) &
634 18 : - dist(coord(:, jatom), coord(:, katom))) < sym%eps_geo)) THEN
635 48 : b(:) = a(:) + coord(:, jatom) + coord(:, katom)
636 12 : IF (caxis(3, b(:), sym, coord)) THEN
637 12 : CALL addsec(3, b(:), sym)
638 12 : 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 7 : IF (sym%nsec(3) < 2) THEN
649 0 : failed = .TRUE.
650 0 : RETURN
651 : END IF
652 :
653 : ! Generate the remaining C3 axes
654 16 : nc3 = 0
655 : DO
656 16 : nc3 = sym%nsec(3)
657 82 : DO ic3 = 1, nc3
658 264 : a(:) = sym%sec(:, ic3, 3)
659 462 : DO jc3 = 1, nc3
660 446 : IF (ic3 /= jc3) THEN
661 1256 : d(:) = sym%sec(:, jc3, 3)
662 942 : DO i = 1, 2
663 628 : phi1 = 2.0_dp*REAL(i, KIND=dp)*pi/3.0_dp
664 628 : b(:) = rotate_vector(a(:), phi1, d(:))
665 942 : CALL addsec(3, b(:), sym)
666 : END DO
667 : END IF
668 : END DO
669 : END DO
670 :
671 16 : 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 7 : IF (sym%nsec(3) == 4) THEN
677 20 : a(:) = sym%sec(:, 1, 3)
678 20 : b(:) = sym%sec(:, 2, 3)
679 5 : phi1 = 0.5_dp*angle(a(:), b(:))
680 5 : IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
681 5 : d(:) = vector_product(a(:), b(:))
682 5 : b(:) = rotate_vector(a(:), phi1, d(:))
683 20 : c(:) = sym%sec(:, 3, 3)
684 5 : phi1 = 0.5_dp*angle(a(:), c(:))
685 5 : IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
686 5 : d(:) = vector_product(a(:), c(:))
687 5 : c(:) = rotate_vector(a(:), phi1, d(:))
688 5 : d(:) = vector_product(b(:), c(:))
689 5 : phi1 = 0.5_dp*phidd
690 5 : a(:) = rotate_vector(b(:), phi1, d(:))
691 5 : IF (caxis(3, a(:), sym, coord)) THEN
692 0 : CALL addsec(3, a(:), sym)
693 : ELSE
694 5 : phi2 = 0.5_dp*pi - phi1
695 5 : a(:) = rotate_vector(b(:), phi2, d(:))
696 5 : IF (caxis(3, a(:), sym, coord)) CALL addsec(3, a(:), sym)
697 : END IF
698 :
699 5 : IF (sym%nsec(3) > 4) CYCLE
700 :
701 2 : ELSE IF (sym%nsec(3) /= 10) THEN
702 :
703 : ! C3 axes found, but only 4 or 10 are possible
704 0 : failed = .TRUE.
705 0 : RETURN
706 :
707 : END IF
708 :
709 : ! Exit loop, if 4 or 10 C3 axes were found
710 9 : 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 47 : DO isec = 1, sym%nsec(3)
718 :
719 160 : a(:) = sym%sec(:, isec, 3)
720 167 : DO jsec = isec + 1, sym%nsec(3)
721 120 : phi1 = 0.5_dp*angle(a(:), sym%sec(:, jsec, 3))
722 120 : d(:) = vector_product(a(:), sym%sec(:, jsec, 3))
723 :
724 : ! Check for C5 axis (I,Ih)
725 120 : IF (sym%nsec(3) == 10) THEN
726 90 : b(:) = rotate_vector(a(:), phidi, d(:))
727 90 : IF (caxis(5, b(:), sym, coord)) THEN
728 14 : CALL addsec(5, b(:), sym)
729 14 : phi1 = phidi + phiii
730 14 : b(:) = rotate_vector(a(:), phi1, d(:))
731 14 : 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 360 : DO i = 0, 1
737 240 : phi2 = phi1 - 0.5_dp*REAL(i, KIND=dp)*pi
738 240 : b(:) = rotate_vector(a(:), phi2, d(:))
739 240 : IF (caxis(2, b(:), sym, coord)) THEN
740 134 : CALL addsec(2, b(:), sym)
741 134 : IF (sym%nsec(3) == 4) THEN
742 42 : IF (caxis(4, b(:), sym, coord)) CALL addsec(4, b(:), sym)
743 : END IF
744 134 : IF (saxis(2, b(:), sym, coord)) THEN
745 64 : CALL addses(2, b(:), sym)
746 64 : sym%invers = .TRUE.
747 : END IF
748 : END IF
749 360 : IF (sigma(b(:), sym, coord)) CALL addsig(b(:), sym)
750 : END DO
751 :
752 : ! Check for mirror plane
753 160 : 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 7 : iax = sym%ncn
761 :
762 7 : IF ((sym%ncn == 5) .AND. (sym%nsec(5) > 1)) THEN
763 2 : sym%igroup = .TRUE.
764 5 : ELSE IF ((sym%ncn == 4) .AND. (sym%nsec(4) > 1)) THEN
765 2 : sym%ogroup = .TRUE.
766 3 : ELSE IF ((sym%ncn == 3) .AND. (sym%nsec(3) > 1)) THEN
767 3 : sym%tgroup = .TRUE.
768 3 : IF ((.NOT. sym%invers) .AND. (sym%nsig > 0)) sym%sigmad = .TRUE.
769 : iax = 2
770 : ELSE
771 : ! No main axis found
772 0 : failed = .TRUE.
773 0 : RETURN
774 : END IF
775 :
776 : ! Rotate molecule to standard orientation
777 7 : phi1 = angle(sym%sec(:, 1, iax), sym%z_axis(:))
778 7 : d(:) = vector_product(sym%sec(:, 1, iax), sym%z_axis(:))
779 7 : CALL rotate_molecule(phi1, d(:), sym, coord)
780 :
781 28 : a(:) = sym%sec(:, 2, iax)
782 7 : a(3) = 0.0_dp
783 7 : phi2 = angle(a(:), sym%x_axis(:))
784 7 : d(:) = vector_product(a(:), sym%x_axis(:))
785 7 : 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 38170 : 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 38170 : equatom = 0
812 38170 : natoms = SIZE(coord, 2)
813 413408 : 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 400638 : (ABS(a(2) - coord(2, iatom)) < sym%eps_geo) .AND. &
817 12770 : (ABS(a(3) - coord(3, iatom)) < sym%eps_geo)) THEN
818 38170 : 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 55 : 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 55 : sym%point_group_order = 1 + sym%nsig
840 :
841 : ! Loop over all C axes
842 249 : DO icn = 2, sym%ncn
843 545 : DO isec = 1, sym%nsec(icn)
844 1021 : DO jcn = 1, icn - 1
845 827 : 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 55 : IF (sym%sigmah) THEN
852 : incr = 1
853 : ELSE
854 40 : incr = 2
855 : END IF
856 :
857 212 : DO isn = 2, sym%nsn
858 285 : DO ises = 1, sym%nses(isn)
859 452 : DO jsn = 1, isn - 1, incr
860 295 : IF (newse(jsn, isn)) sym%point_group_order = sym%point_group_order + 1
861 : END DO
862 : END DO
863 : END DO
864 :
865 55 : 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 57 : SUBROUTINE get_point_group_symbol(sym)
875 : TYPE(molsym_type), INTENT(inout) :: sym
876 :
877 : CHARACTER(LEN=4) :: degree
878 :
879 57 : IF (sym%linear) THEN
880 2 : IF (sym%invers) THEN
881 1 : sym%point_group_symbol = "D*h"
882 : ELSE
883 1 : sym%point_group_symbol = "C*v"
884 : END IF
885 55 : ELSE IF (sym%cubic) THEN
886 7 : IF (sym%igroup) THEN
887 2 : sym%point_group_symbol = "I"
888 5 : ELSE IF (sym%ogroup) THEN
889 2 : sym%point_group_symbol = "O"
890 3 : ELSE IF (sym%tgroup) THEN
891 3 : sym%point_group_symbol = "T"
892 : END IF
893 7 : IF (sym%invers) THEN
894 3 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
895 4 : ELSE IF (sym%sigmad) THEN
896 1 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"d"
897 : END IF
898 48 : ELSE IF (sym%maxis) THEN
899 45 : IF (sym%dgroup) THEN
900 21 : WRITE (degree, "(I3)") sym%ncn
901 21 : sym%point_group_symbol = "D"//TRIM(ADJUSTL(degree))
902 21 : IF (sym%sigmah) THEN
903 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
904 14 : ELSE IF (sym%sigmad) THEN
905 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"d"
906 : END IF
907 24 : ELSE IF (sym%sgroup) THEN
908 3 : WRITE (degree, "(I3)") sym%nsn
909 3 : sym%point_group_symbol = "S"//TRIM(ADJUSTL(degree))
910 : ELSE
911 21 : WRITE (degree, "(I3)") sym%ncn
912 21 : sym%point_group_symbol = "C"//TRIM(ADJUSTL(degree))
913 21 : IF (sym%sigmah) THEN
914 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
915 14 : ELSE IF (sym%sigmav) THEN
916 7 : sym%point_group_symbol = TRIM(sym%point_group_symbol)//"v"
917 : END IF
918 : END IF
919 : ELSE
920 3 : IF (sym%invers) THEN
921 1 : sym%point_group_symbol = "Ci"
922 2 : ELSE IF (sym%sigmah) THEN
923 1 : sym%point_group_symbol = "Cs"
924 : ELSE
925 1 : sym%point_group_symbol = "C1"
926 : END IF
927 : END IF
928 :
929 57 : 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 58 : 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 232 : sym%x_axis = (/1.0_dp, 0.0_dp, 0.0_dp/)
949 232 : sym%y_axis = (/0.0_dp, 1.0_dp, 0.0_dp/)
950 232 : sym%z_axis = (/0.0_dp, 0.0_dp, 1.0_dp/)
951 :
952 : ! Initialize lists for symmetry elements
953 93728 : sym%sec(:, :, :) = 0.0_dp
954 373288 : sym%ses(:, :, :) = 0.0_dp
955 4930 : sym%sig(:, :) = 0.0_dp
956 :
957 232 : sym%center_of_mass(:) = 0.0_dp
958 :
959 : ! Initialize symmetry element counters
960 58 : sym%ncn = 1
961 58 : sym%nsn = 1
962 1218 : sym%nsec(:) = 0
963 2378 : sym%nses(:) = 0
964 58 : sym%nsig = 0
965 :
966 : ! Initialize logical variables
967 58 : sym%cubic = .FALSE.
968 58 : sym%dgroup = .FALSE.
969 58 : sym%igroup = .FALSE.
970 58 : sym%invers = .FALSE.
971 58 : sym%linear = .FALSE.
972 58 : sym%maxis = .FALSE.
973 58 : sym%ogroup = .FALSE.
974 58 : sym%sgroup = .FALSE.
975 58 : sym%sigmad = .FALSE.
976 58 : sym%sigmah = .FALSE.
977 58 : sym%sigmav = .FALSE.
978 58 : sym%tgroup = .FALSE.
979 :
980 : ! Initialize list of equivalent atoms (C1 symmetry)
981 58 : natoms = SIZE(sym%nequatom)
982 58 : sym%ngroup = natoms
983 633 : sym%nequatom(:) = 1
984 633 : DO i = 1, sym%ngroup
985 575 : sym%group_of(i) = i
986 575 : sym%llequatom(i) = i
987 575 : sym%symequ_list(i) = i
988 633 : sym%ulequatom(i) = i
989 : END DO
990 :
991 58 : sym%point_group_order = -1
992 :
993 633 : DO iatom = 1, natoms
994 633 : sym%aw(iatom) = weight(iatom)
995 : END DO
996 :
997 : ! Generate atomic identification numbers (symmetry code) ***
998 633 : sym%ain(:) = REAL(atype(:), KIND=dp) + sym%aw(:)
999 :
1000 : ! Initialize the transformation matrix for input orientation -> standard orientation
1001 58 : CALL unit_matrix(sym%inptostd(:, :))
1002 :
1003 58 : 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 1211 : 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 1211 : in_symequ_list = .FALSE.
1022 :
1023 7969 : DO iatom = 1, sym%ulequatom(sym%ngroup)
1024 7969 : IF (sym%symequ_list(iatom) == atoma) THEN
1025 1211 : 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 3 : 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 3 : IF (sym%nsn == 2) THEN
1048 : ! Ci
1049 1 : sym%invers = .TRUE.
1050 1 : phi = angle(sym%ses(:, 1, 2), sym%z_axis(:))
1051 1 : d(:) = vector_product(sym%ses(:, 1, 2), sym%z_axis(:))
1052 1 : CALL rotate_molecule(phi, d(:), sym, coord)
1053 2 : ELSE IF (sym%nsig == 1) THEN
1054 : ! Cs
1055 1 : sym%sigmah = .TRUE.
1056 1 : phi = angle(sym%sig(:, 1), sym%z_axis(:))
1057 1 : d(:) = vector_product(sym%sig(:, 1), sym%z_axis(:))
1058 1 : CALL rotate_molecule(phi, d(:), sym, coord)
1059 : END IF
1060 :
1061 3 : 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 58 : 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 58 : CALL timeset(routineN, handle)
1086 :
1087 : ! Perform memory allocation for the symmetry analysis
1088 58 : natoms = SIZE(coord, 2)
1089 58 : CALL create_molsym(sym, natoms)
1090 58 : sym%eps_geo = eps_geo
1091 :
1092 : ! Initialization of symmetry analysis
1093 58 : CALL init_symmetry(sym, atype, weight)
1094 :
1095 58 : IF (natoms == 1) THEN
1096 : ! Special case: only one atom
1097 1 : CALL atomsym(sym)
1098 : ELSE
1099 : ! Find molecular symmetry
1100 57 : CALL moleculesym(sym, coord)
1101 : ! Get point group and load point group table
1102 57 : CALL get_point_group_symbol(sym)
1103 : END IF
1104 :
1105 : ! Calculate group order
1106 58 : IF (.NOT. sym%linear) CALL get_point_group_order(sym)
1107 :
1108 : ! Generate a list of equivalent atoms
1109 58 : CALL build_symequ_list(sym, coord)
1110 :
1111 58 : CALL timestop(handle)
1112 :
1113 58 : 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 57 : 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 57 : eps_tenval = 0.01_dp*sym%eps_geo
1134 :
1135 : ! Calculate the molecular tensor of inertia
1136 57 : CALL tensor(sym, coord)
1137 : ! Use symmetry information from the eigenvalues of the molecular tensor of inertia
1138 57 : IF ((sym%tenval(3) - sym%tenval(1)) < eps_tenval) THEN
1139 : ! 0 < tenval(1) = tenval(2) = tenval(3)
1140 7 : sym%cubic = .TRUE.
1141 50 : 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 14 : IF (sym%tenval(1) < eps_tenval) sym%linear = .TRUE.
1145 14 : 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 36 : CALL tracar(sym, coord, (/1, 2, 3/))
1149 : END IF
1150 :
1151 : ! Continue with the new coordinate axes
1152 : DO
1153 57 : failed = .FALSE.
1154 57 : IF (sym%cubic) THEN
1155 7 : CALL cubsym(sym, coord, failed)
1156 7 : IF (failed) THEN
1157 0 : sym%cubic = .FALSE.
1158 0 : CYCLE
1159 : END IF
1160 :
1161 50 : ELSE IF (sym%linear) THEN
1162 :
1163 : ! Linear molecule
1164 2 : IF (saxis(2, sym%z_axis(:), sym, coord)) THEN
1165 1 : sym%invers = .TRUE.
1166 1 : sym%dgroup = .TRUE.
1167 1 : sym%sigmah = .TRUE.
1168 : END IF
1169 :
1170 : ELSE
1171 :
1172 : ! Check the new coordinate axes for Cn axes
1173 960 : DO icn = 2, maxcn
1174 912 : IF (caxis(icn, sym%z_axis(:), sym, coord)) CALL addsec(icn, sym%z_axis(:), sym)
1175 912 : IF (caxis(icn, sym%x_axis(:), sym, coord)) CALL addsec(icn, sym%x_axis(:), sym)
1176 960 : 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 1920 : DO isn = 2, maxsn
1181 1872 : IF (saxis(isn, sym%z_axis(:), sym, coord)) CALL addses(isn, sym%z_axis(:), sym)
1182 1872 : IF (saxis(isn, sym%x_axis(:), sym, coord)) CALL addses(isn, sym%x_axis(:), sym)
1183 1920 : 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 48 : IF (sigma(sym%z_axis(:), sym, coord)) CALL addsig(sym%z_axis(:), sym)
1188 48 : IF (sigma(sym%x_axis(:), sym, coord)) CALL addsig(sym%x_axis(:), sym)
1189 48 : IF (sigma(sym%y_axis(:), sym, coord)) CALL addsig(sym%y_axis(:), sym)
1190 :
1191 : ! There is a main axis (MAXIS = .TRUE.)
1192 48 : IF ((sym%ncn > 1) .OR. (sym%nsn > 3)) THEN
1193 45 : sym%maxis = .TRUE.
1194 45 : CALL axsym(coord, sym)
1195 : ELSE
1196 : ! Only low or no symmetry (Ci, Cs or C1)
1197 3 : CALL lowsym(sym, coord)
1198 : END IF
1199 :
1200 : END IF
1201 :
1202 57 : IF (.NOT. failed) EXIT
1203 :
1204 : END DO
1205 :
1206 : ! Find the remaining S axes
1207 57 : IF (.NOT. sym%linear) THEN
1208 249 : DO icn = 2, sym%ncn
1209 545 : DO isec = 1, sym%nsec(icn)
1210 296 : IF (saxis(icn, sym%sec(:, isec, icn), sym, coord)) &
1211 89 : CALL addses(icn, sym%sec(:, isec, icn), sym)
1212 296 : IF (saxis(2*icn, sym%sec(:, isec, icn), sym, coord)) &
1213 243 : 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 57 : IF (sym%nses(2) > 0) THEN
1220 16 : sym%nses(1) = sym%nses(1) - sym%nses(2)
1221 16 : sym%nses(2) = 0
1222 16 : CALL addses(2, sym%z_axis(:), sym)
1223 : END IF
1224 :
1225 57 : 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 76 : 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 76 : naxis = 0
1248 76 : natoms = SIZE(coord, 2)
1249 :
1250 76 : 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 76 : IF (length_of_a > sym%eps_geo) THEN
1254 :
1255 956 : DO iatom = 1, natoms
1256 3520 : b(:) = coord(:, iatom)
1257 880 : 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 956 : IF (length_of_b < sym%eps_geo) THEN
1260 0 : naxis = naxis + 1
1261 : ELSE
1262 3520 : a_norm = a(:)/length_of_a
1263 3520 : b_norm = b(:)/length_of_b
1264 880 : scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1265 880 : 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 76 : 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 1802 : FUNCTION newse(se1, se2)
1284 : INTEGER, INTENT(IN) :: se1, se2
1285 : LOGICAL :: newse
1286 :
1287 : INTEGER :: ise
1288 :
1289 1802 : newse = .TRUE.
1290 :
1291 3878 : DO ise = 2, MIN(se1, se2)
1292 3878 : IF ((MODULO(se1, ise) == 0) .AND. (MODULO(se2, ise) == 0)) THEN
1293 1802 : 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 105 : 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 105 : nsigma = 0
1322 :
1323 105 : 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 105 : IF (length_of_a > sym%eps_geo) THEN
1327 105 : natoms = SIZE(coord, 2)
1328 998 : DO iatom = 1, natoms
1329 3572 : b(:) = coord(:, iatom)
1330 893 : 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 998 : IF (length_of_b < sym%eps_geo) THEN
1333 0 : nsigma = nsigma + 1
1334 : ELSE
1335 3572 : a_norm = a(:)/length_of_a
1336 3572 : b_norm = b(:)/length_of_b
1337 893 : scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
1338 893 : IF (ABS(scapro) < sym%eps_geo) nsigma = nsigma + 1
1339 : END IF
1340 : END DO
1341 : END IF
1342 :
1343 105 : 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 522 : 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 522 : IF (ABS(se(3)) < eps) THEN
1360 255 : IF (ABS(se(1)) < eps) THEN
1361 47 : se(2) = -se(2)
1362 : ELSE
1363 448 : IF (se(1) < 0.0_dp) se(1:2) = -se(1:2)
1364 : END IF
1365 : ELSE
1366 471 : IF (se(3) < 0.0_dp) se(:) = -se(:)
1367 : END IF
1368 :
1369 522 : 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 58 : 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 58 : INTEGER, ALLOCATABLE, DIMENSION(:) :: equatomlist, idxlist
1400 :
1401 : ! Print point group symbol and point group order
1402 : WRITE (iw, "(/,(T2,A))") &
1403 58 : "MOLSYM| Molecular symmetry information", &
1404 116 : "MOLSYM|"
1405 : WRITE (iw, "(T2,A,T77,A4)") &
1406 58 : "MOLSYM| Point group", ADJUSTR(sym%point_group_symbol)
1407 58 : IF (sym%point_group_order > -1) THEN
1408 : WRITE (iw, "(T2,A,T77,I4)") &
1409 55 : "MOLSYM| Group order ", sym%point_group_order
1410 : END IF
1411 :
1412 58 : IF (MOD(plevel, 10) == 1) THEN
1413 : WRITE (iw, "(T2,A)") &
1414 58 : "MOLSYM|", &
1415 116 : "MOLSYM| Groups of atoms equivalent by symmetry"
1416 : WRITE (iw, "(T2,A,T31,A)") &
1417 58 : "MOLSYM| Group number", "Group members (atomic indices)"
1418 164 : DO i = 1, sym%ngroup
1419 106 : nequatom = sym%ulequatom(i) - sym%llequatom(i) + 1
1420 106 : CPASSERT(nequatom > 0)
1421 530 : ALLOCATE (equatomlist(nequatom), idxlist(nequatom))
1422 681 : equatomlist(1:nequatom) = sym%symequ_list(sym%llequatom(i):sym%ulequatom(i))
1423 681 : idxlist = 0
1424 106 : 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 106 : "MOLSYM|", i, nequatom, (equatomlist(iequatom), iequatom=1, nequatom)
1427 164 : DEALLOCATE (equatomlist, idxlist)
1428 : END DO
1429 : END IF
1430 :
1431 58 : 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 58 : "MOLSYM|", &
1435 58 : "MOLSYM| Symmetry elements", &
1436 116 : "MOLSYM| Element number", "Type", "X", "Y", "Z"
1437 58 : secount = 1
1438 : WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A1,T36,3(1X,F14.8))") &
1439 58 : "MOLSYM|", secount, secount, "E", 0.0_dp, 0.0_dp, 0.0_dp
1440 : ! Mirror planes
1441 58 : string = "@ "
1442 211 : DO isig = 1, sym%nsig
1443 153 : secount = secount + 1
1444 153 : 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 670 : "MOLSYM|", secount, isig, string, (sym%sig(i, isig), i=1, 3)
1447 : END DO
1448 : ! C axes
1449 252 : DO icn = 2, sym%ncn
1450 194 : IF (icn < 10) THEN
1451 194 : WRITE (string, "(A1,I1)") "C", icn
1452 : ELSE
1453 0 : WRITE (string, "(A1,I2)") "C", icn
1454 : END IF
1455 548 : DO isec = 1, sym%nsec(icn)
1456 296 : secount = secount + 1
1457 296 : 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 1378 : "MOLSYM|", secount, isec, string, (sym%sec(i, isec, icn), i=1, 3)
1460 : END DO
1461 : END DO
1462 : ! S axes
1463 215 : DO isn = 2, sym%nsn
1464 157 : IF (isn == 2) THEN
1465 29 : WRITE (string, "(A3)") "i "
1466 128 : ELSE IF (isn < 10) THEN
1467 111 : WRITE (string, "(A1,I1)") "S", isn
1468 : ELSE
1469 17 : WRITE (string, "(A1,I2)") "S", isn
1470 : END IF
1471 288 : DO ises = 1, sym%nses(isn)
1472 73 : secount = secount + 1
1473 73 : 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 449 : "MOLSYM|", secount, ises, string, (sym%ses(i, ises, icn), i=1, 3)
1476 : END DO
1477 : END DO
1478 : END IF
1479 :
1480 58 : 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 57 : "MOLSYM|", &
1484 57 : "MOLSYM| Moments of the molecular inertia tensor [a.u.]", &
1485 114 : "MOLSYM|", "I(x)", "I(y)", "I(z)"
1486 57 : string = "xyz"
1487 741 : IF (MAXVAL(sym%tenmat) >= formatmaxval) THEN
1488 0 : DO i = 1, 3
1489 : WRITE (iw, "(T2,A,T32,A,T36,3(1X,ES14.4))") &
1490 0 : "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1491 : END DO
1492 : ELSE
1493 228 : DO i = 1, 3
1494 : WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1495 741 : "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
1496 : END DO
1497 : END IF
1498 : WRITE (iw, "(T2,A)") &
1499 57 : "MOLSYM|", &
1500 114 : "MOLSYM| Principal moments and axes of inertia [a.u.]"
1501 741 : IF (MAXVAL(sym%tenmat) >= formatmaxval) THEN
1502 : WRITE (iw, "(T2,A,T36,3(1X,ES14.4))") &
1503 0 : "MOLSYM|", (sym%tenval(i), i=1, 3)
1504 : ELSE
1505 : WRITE (iw, "(T2,A,T36,3(1X,F14.8))") &
1506 228 : "MOLSYM|", (sym%tenval(i), i=1, 3)
1507 : END IF
1508 228 : DO i = 1, 3
1509 : WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
1510 741 : "MOLSYM|", string(i:i), (sym%tenvec(i, j), j=1, 3)
1511 : END DO
1512 : END IF
1513 :
1514 58 : IF (MOD(plevel, 10) == 1) THEN
1515 : ! Print the Cartesian coordinates of the standard orientation
1516 58 : CALL write_coordinates(coord, atype, element, z, weight, iw)
1517 : END IF
1518 :
1519 58 : 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 87 : 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 87 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1543 87 : IF (length_of_a > sym%eps_geo) THEN
1544 :
1545 : ! Build up the rotation matrix
1546 42 : CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1547 :
1548 : ! Rotate the molecule by phi around a
1549 10689 : coord(:, :) = MATMUL(sym%rotmat(:, :), coord(:, :))
1550 :
1551 : ! Rotate the normal vectors of the mirror planes which are already found
1552 3339 : sym%sig(:, 1:sym%nsig) = MATMUL(sym%rotmat(:, :), sym%sig(:, 1:sym%nsig))
1553 :
1554 : ! Rotate the Cn axes which are already found
1555 186 : DO i = 2, sym%ncn
1556 9696 : 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 134 : DO i = 2, sym%nsn
1561 3014 : 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 2688 : sym%inptostd(:, :) = MATMUL(sym%rotmat(:, :), sym%inptostd(:, :))
1566 :
1567 : END IF
1568 :
1569 87 : 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 6344 : 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 6344 : saxis = .FALSE.
1595 :
1596 6344 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
1597 :
1598 6344 : natoms = SIZE(coord, 2)
1599 :
1600 : ! Check the length of the axis vector a
1601 6344 : IF (length_of_a > sym%eps_geo) THEN
1602 : ! Calculate the rotation angle phi and build up the rotation matrix rotmat
1603 6344 : phi = 2.0_dp*pi/REAL(n, KIND=dp)
1604 6344 : CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
1605 : ! Check all atoms
1606 9816 : DO iatom = 1, natoms
1607 : ! Rotate the actual atom by 2*pi/n about a
1608 124111 : b(:) = MATMUL(sym%rotmat(:, :), coord(:, iatom))
1609 : ! Reflect the coordinates of the rotated atom
1610 38188 : b(:) = reflect_vector(b(:), a(:))
1611 : ! Check if the coordinates of the rotated atom are in the coordinate set of the molecule
1612 9816 : 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 2808 : 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 2808 : sigma = .FALSE.
1642 :
1643 2808 : 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 2808 : IF (length_of_a > sym%eps_geo) THEN
1647 2674 : natoms = SIZE(coord, 2)
1648 10068 : DO iatom = 1, natoms
1649 : ! Reflect the actual atom
1650 9471 : b(:) = reflect_vector(coord(:, iatom), a(:))
1651 : ! Check if the coordinates of the reflected atom are in the coordinate set of the molecule
1652 10068 : 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 57 : 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 57 : natoms = SIZE(coord, 2)
1678 3269 : 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 2353 : coord(:, 1:natoms) = coord(:, 1:natoms) - SPREAD(sym%center_of_mass(:), 2, natoms)
1682 :
1683 : ! Build up the molecular tensor of inertia
1684 :
1685 631 : sym%tenmat(1, 1) = DOT_PRODUCT(sym%aw(1:natoms), (coord(2, 1:natoms)**2 + coord(3, 1:natoms)**2))
1686 631 : sym%tenmat(2, 2) = DOT_PRODUCT(sym%aw(1:natoms), (coord(3, 1:natoms)**2 + coord(1, 1:natoms)**2))
1687 631 : sym%tenmat(3, 3) = DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)**2 + coord(2, 1:natoms)**2))
1688 :
1689 631 : sym%tenmat(1, 2) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(2, 1:natoms)))
1690 631 : sym%tenmat(1, 3) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(3, 1:natoms)))
1691 631 : sym%tenmat(2, 3) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(2, 1:natoms)*coord(3, 1:natoms)))
1692 :
1693 : ! Symmetrize tensor matrix
1694 57 : sym%tenmat(2, 1) = sym%tenmat(1, 2)
1695 57 : sym%tenmat(3, 1) = sym%tenmat(1, 3)
1696 57 : sym%tenmat(3, 2) = sym%tenmat(2, 3)
1697 :
1698 : ! Diagonalize the molecular tensor of inertia
1699 57 : CALL jacobi(sym%tenmat(:, :), sym%tenval(:), sym%tenvec(:, :))
1700 :
1701 : ! Secure that the principal axes are right-handed
1702 228 : sym%tenvec(:, 3) = vector_product(sym%tenvec(:, 1), sym%tenvec(:, 2))
1703 :
1704 57 : tt = SQRT(sym%tenval(1)**2 + sym%tenval(2)**2 + sym%tenval(3)**2)
1705 57 : CPASSERT(tt /= 0)
1706 :
1707 57 : 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 50 : 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 50 : natom = SIZE(coord, 2)
1734 650 : tenvec = TRANSPOSE(sym%tenvec)
1735 482 : DO iatom = 1, natom
1736 3074 : coord(idx, iatom) = MATMUL(tenvec, coord(:, iatom))
1737 : END DO
1738 :
1739 : ! Modify the transformation matrix for input orientation -> standard orientation
1740 650 : sym%inptostd(idx, :) = tenvec
1741 :
1742 50 : END SUBROUTINE tracar
1743 :
1744 : ! **************************************************************************************************
1745 : !> \brief Distance between two points
1746 : !> \param a ...
1747 : !> \param b ...
1748 : !> \return ...
1749 : !> \author jgh
1750 : ! **************************************************************************************************
1751 2864 : FUNCTION dist(a, b) RESULT(d)
1752 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b
1753 : REAL(KIND=dp) :: d
1754 :
1755 11456 : d = SQRT(SUM((a - b)**2))
1756 :
1757 2864 : 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 58 : 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 58 : IF (iw > 0) THEN
1780 58 : natom = SIZE(coord, 2)
1781 : WRITE (UNIT=iw, FMT="(T2,A)") &
1782 58 : "MOLSYM|", &
1783 116 : "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 58 : "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1786 633 : 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 575 : "MOLSYM|", iatom, atype(iatom), ADJUSTL(element(iatom)), z(iatom), &
1789 1208 : coord(1:3, iatom), weight(iatom)
1790 : END DO
1791 : WRITE (UNIT=iw, FMT="(T2,A)") &
1792 58 : "MOLSYM|", &
1793 116 : "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 58 : "MOLSYM| Atom Kind Element", "X", "Y", "Z", "Mass"
1796 633 : 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 575 : "MOLSYM|", iatom, atype(iatom), ADJUSTL(element(iatom)), z(iatom), &
1799 2933 : coord(1:3, iatom)*angstrom, weight(iatom)
1800 : END DO
1801 : END IF
1802 :
1803 58 : END SUBROUTINE write_coordinates
1804 :
1805 0 : END MODULE molsym
|
