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