LCOV - code coverage report
Current view: top level - src - rmsd.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 162 171 94.7 %
Date: 2024-03-28 07:31:50 Functions: 1 1 100.0 %

          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 Defines functions to perform rmsd in 3D
      10             : !> \author Teodoro Laino 09.2006
      11             : ! **************************************************************************************************
      12             : MODULE rmsd
      13             : 
      14             :    USE kinds,                           ONLY: dp
      15             :    USE mathlib,                         ONLY: diamat_all
      16             :    USE particle_types,                  ONLY: particle_type
      17             : #include "./base/base_uses.f90"
      18             : 
      19             :    IMPLICIT NONE
      20             :    PRIVATE
      21             : 
      22             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rmsd'
      23             :    REAL(KIND=dp), PARAMETER, PRIVATE    :: Epsi = EPSILON(0.0_dp)*1.0E4_dp
      24             : 
      25             :    PUBLIC :: rmsd3
      26             : 
      27             : CONTAINS
      28             : 
      29             : ! **************************************************************************************************
      30             : !> \brief Computes the RMSD in 3D. Provides also derivatives.
      31             : !> \param particle_set ...
      32             : !> \param r ...
      33             : !> \param r0 ...
      34             : !> \param output_unit ...
      35             : !> \param weights ...
      36             : !> \param my_val ...
      37             : !> \param rotate ...
      38             : !> \param transl ...
      39             : !> \param rot ...
      40             : !> \param drmsd3 ...
      41             : !> \author Teodoro Laino 08.2006
      42             : !> \note
      43             : !>      Optional arguments:
      44             : !>          my_val -> gives back the value of the RMSD
      45             : !>          transl -> provides the translational vector
      46             : !>          rotate -> if .true. gives back in r the coordinates rotated
      47             : !>                    in order to minimize the RMSD3
      48             : !>          rot    -> provides the rotational matrix
      49             : !>          drmsd3 -> derivatives of RMSD3 w.r.t. atomic positions
      50             : ! **************************************************************************************************
      51        3092 :    SUBROUTINE rmsd3(particle_set, r, r0, output_unit, weights, my_val, &
      52        3092 :                     rotate, transl, rot, drmsd3)
      53             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      54             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: r, r0
      55             :       INTEGER, INTENT(IN)                                :: output_unit
      56             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: weights
      57             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: my_val
      58             :       LOGICAL, INTENT(IN), OPTIONAL                      :: rotate
      59             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: transl
      60             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
      61             :          OPTIONAL                                        :: rot, drmsd3
      62             : 
      63             :       INTEGER                                            :: i, ix, j, k, natom
      64             :       LOGICAL                                            :: my_rotate
      65             :       REAL(KIND=dp)                                      :: dm_r(4, 4, 3), lambda(4), loc_tr(3), &
      66             :                                                             M(4, 4), mtot, my_rot(3, 3), Q(0:3), &
      67             :                                                             rr(3), rr0(3), rrsq, s, xx, yy, &
      68             :                                                             Z(4, 4), zz
      69        3092 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: w
      70        3092 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r0p, rp
      71             : 
      72        3092 :       CPASSERT(SIZE(r) == SIZE(r0))
      73        3092 :       natom = SIZE(particle_set)
      74        3092 :       my_rotate = .FALSE.
      75        3092 :       IF (PRESENT(rotate)) my_rotate = rotate
      76        9276 :       ALLOCATE (w(natom))
      77        3092 :       IF (PRESENT(weights)) THEN
      78         468 :          w(:) = weights
      79             :       ELSE
      80             :          ! All atoms have a weight proportional to their mass in the RMSD unless
      81             :          ! differently requested
      82       61920 :          DO i = 1, natom
      83       61920 :             w(i) = particle_set(i)%atomic_kind%mass
      84             :          END DO
      85       64976 :          mtot = MINVAL(w)
      86       64976 :          IF (mtot /= 0.0_dp) w(:) = w(:)/mtot
      87             :       END IF
      88        9276 :       ALLOCATE (rp(3, natom))
      89        9276 :       ALLOCATE (r0p(3, natom))
      90             :       ! Molecule given by coordinates R
      91             :       ! Find COM and center molecule in COM
      92       62388 :       xx = 0.0_dp
      93       62388 :       yy = 0.0_dp
      94       62388 :       zz = 0.0_dp
      95       62388 :       mtot = 0.0_dp
      96       62388 :       DO i = 1, natom
      97       59296 :          mtot = mtot + particle_set(i)%atomic_kind%mass
      98       59296 :          xx = xx + r((i - 1)*3 + 1)*particle_set(i)%atomic_kind%mass
      99       59296 :          yy = yy + r((i - 1)*3 + 2)*particle_set(i)%atomic_kind%mass
     100       62388 :          zz = zz + r((i - 1)*3 + 3)*particle_set(i)%atomic_kind%mass
     101             :       END DO
     102        3092 :       xx = xx/mtot
     103        3092 :       yy = yy/mtot
     104        3092 :       zz = zz/mtot
     105       62388 :       DO i = 1, natom
     106       59296 :          rp(1, i) = r((i - 1)*3 + 1) - xx
     107       59296 :          rp(2, i) = r((i - 1)*3 + 2) - yy
     108       62388 :          rp(3, i) = r((i - 1)*3 + 3) - zz
     109             :       END DO
     110        3092 :       IF (PRESENT(transl)) THEN
     111           0 :          transl(1) = xx
     112           0 :          transl(2) = yy
     113           0 :          transl(3) = zz
     114             :       END IF
     115             :       ! Molecule given by coordinates R0
     116             :       ! Find COM and center molecule in COM
     117        3092 :       xx = 0.0_dp
     118        3092 :       yy = 0.0_dp
     119        3092 :       zz = 0.0_dp
     120       62388 :       DO i = 1, natom
     121       59296 :          xx = xx + r0((i - 1)*3 + 1)*particle_set(i)%atomic_kind%mass
     122       59296 :          yy = yy + r0((i - 1)*3 + 2)*particle_set(i)%atomic_kind%mass
     123       62388 :          zz = zz + r0((i - 1)*3 + 3)*particle_set(i)%atomic_kind%mass
     124             :       END DO
     125        3092 :       xx = xx/mtot
     126        3092 :       yy = yy/mtot
     127        3092 :       zz = zz/mtot
     128       62388 :       DO i = 1, natom
     129       59296 :          r0p(1, i) = r0((i - 1)*3 + 1) - xx
     130       59296 :          r0p(2, i) = r0((i - 1)*3 + 2) - yy
     131       62388 :          r0p(3, i) = r0((i - 1)*3 + 3) - zz
     132             :       END DO
     133        3092 :       loc_tr(1) = xx
     134        3092 :       loc_tr(2) = yy
     135        3092 :       loc_tr(3) = zz
     136             :       ! Give back the translational vector
     137        3092 :       IF (PRESENT(transl)) THEN
     138           0 :          transl(1) = transl(1) - xx
     139           0 :          transl(2) = transl(2) - yy
     140           0 :          transl(3) = transl(3) - zz
     141             :       END IF
     142       62388 :       M = 0.0_dp
     143             :       !
     144       62388 :       DO i = 1, natom
     145       59296 :          IF (w(i) == 0.0_dp) CYCLE
     146       59080 :          rr(1) = rp(1, I)
     147       59080 :          rr(2) = rp(2, I)
     148       59080 :          rr(3) = rp(3, I)
     149       59080 :          rr0(1) = r0p(1, I)
     150       59080 :          rr0(2) = r0p(2, I)
     151       59080 :          rr0(3) = r0p(3, I)
     152       59080 :          rrsq = w(I)*(rr0(1)**2 + rr0(2)**2 + rr0(3)**2 + rr(1)**2 + rr(2)**2 + rr(3)**2)
     153       59080 :          rr0(1) = w(I)*rr0(1)
     154       59080 :          rr0(2) = w(I)*rr0(2)
     155       59080 :          rr0(3) = w(I)*rr0(3)
     156       59080 :          M(1, 1) = M(1, 1) + rrsq + 2.0_dp*(-rr0(1)*rr(1) - rr0(2)*rr(2) - rr0(3)*rr(3))
     157       59080 :          M(2, 2) = M(2, 2) + rrsq + 2.0_dp*(-rr0(1)*rr(1) + rr0(2)*rr(2) + rr0(3)*rr(3))
     158       59080 :          M(3, 3) = M(3, 3) + rrsq + 2.0_dp*(rr0(1)*rr(1) - rr0(2)*rr(2) + rr0(3)*rr(3))
     159       59080 :          M(4, 4) = M(4, 4) + rrsq + 2.0_dp*(rr0(1)*rr(1) + rr0(2)*rr(2) - rr0(3)*rr(3))
     160       59080 :          M(1, 2) = M(1, 2) + 2.0_dp*(-rr0(2)*rr(3) + rr0(3)*rr(2))
     161       59080 :          M(1, 3) = M(1, 3) + 2.0_dp*(rr0(1)*rr(3) - rr0(3)*rr(1))
     162       59080 :          M(1, 4) = M(1, 4) + 2.0_dp*(-rr0(1)*rr(2) + rr0(2)*rr(1))
     163       59080 :          M(2, 3) = M(2, 3) - 2.0_dp*(rr0(1)*rr(2) + rr0(2)*rr(1))
     164       59080 :          M(2, 4) = M(2, 4) - 2.0_dp*(rr0(1)*rr(3) + rr0(3)*rr(1))
     165       62388 :          M(3, 4) = M(3, 4) - 2.0_dp*(rr0(2)*rr(3) + rr0(3)*rr(2))
     166             :       END DO
     167             :       ! Symmetrize
     168        3092 :       M(2, 1) = M(1, 2)
     169        3092 :       M(3, 1) = M(1, 3)
     170        3092 :       M(3, 2) = M(2, 3)
     171        3092 :       M(4, 1) = M(1, 4)
     172        3092 :       M(4, 2) = M(2, 4)
     173        3092 :       M(4, 3) = M(3, 4)
     174             :       ! Solve the eigenvalue problem for M
     175        3092 :       Z = M
     176        3092 :       CALL diamat_all(Z, lambda)
     177             :       ! Pick the correct eigenvectors
     178        3092 :       S = 1.0_dp
     179        3092 :       IF (Z(1, 1) .LT. 0.0_dp) S = -1.0_dp
     180        3092 :       Q(0) = S*Z(1, 1)
     181        3092 :       Q(1) = S*Z(2, 1)
     182        3092 :       Q(2) = S*Z(3, 1)
     183        3092 :       Q(3) = S*Z(4, 1)
     184        3092 :       IF (PRESENT(my_val)) THEN
     185        1056 :          IF (ABS(lambda(1)) < epsi) THEN
     186           0 :             my_val = 0.0_dp
     187             :          ELSE
     188        1056 :             my_val = lambda(1)/REAL(natom, KIND=dp)
     189             :          END IF
     190             :       END IF
     191        3092 :       IF (ABS(lambda(1) - lambda(2)) < epsi) THEN
     192           0 :          IF (output_unit > 0) WRITE (output_unit, FMT='(T2,"RMSD3|",A)') &
     193           0 :             'NORMAL EXECUTION, NON-UNIQUE SOLUTION'
     194             :       END IF
     195             :       ! Computes derivatives w.r.t. the positions if requested
     196        3092 :       IF (PRESENT(drmsd3)) THEN
     197         468 :          DO I = 1, natom
     198         432 :             IF (W(I) == 0.0_dp) CYCLE
     199         216 :             rr(1) = W(I)*2.0_dp*rp(1, I)
     200         216 :             rr(2) = W(I)*2.0_dp*rp(2, I)
     201         216 :             rr(3) = W(I)*2.0_dp*rp(3, I)
     202         216 :             rr0(1) = W(I)*2.0_dp*r0p(1, I)
     203         216 :             rr0(2) = W(I)*2.0_dp*r0p(2, I)
     204         216 :             rr0(3) = W(I)*2.0_dp*r0p(3, I)
     205             : 
     206         216 :             dm_r(1, 1, 1) = (rr(1) - rr0(1))
     207         216 :             dm_r(1, 1, 2) = (rr(2) - rr0(2))
     208         216 :             dm_r(1, 1, 3) = (rr(3) - rr0(3))
     209             : 
     210         216 :             dm_r(1, 2, 1) = 0.0_dp
     211         216 :             dm_r(1, 2, 2) = rr0(3)
     212         216 :             dm_r(1, 2, 3) = -rr0(2)
     213             : 
     214         216 :             dm_r(1, 3, 1) = -rr0(3)
     215         216 :             dm_r(1, 3, 2) = 0.0_dp
     216         216 :             dm_r(1, 3, 3) = rr0(1)
     217             : 
     218         216 :             dm_r(1, 4, 1) = rr0(2)
     219         216 :             dm_r(1, 4, 2) = -rr0(1)
     220         216 :             dm_r(1, 4, 3) = 0.0_dp
     221             : 
     222         216 :             dm_r(2, 2, 1) = (rr(1) - rr0(1))
     223         216 :             dm_r(2, 2, 2) = (rr(2) + rr0(2))
     224         216 :             dm_r(2, 2, 3) = (rr(3) + rr0(3))
     225             : 
     226         216 :             dm_r(2, 3, 1) = -rr0(2)
     227         216 :             dm_r(2, 3, 2) = -rr0(1)
     228         216 :             dm_r(2, 3, 3) = 0.0_dp
     229             : 
     230         216 :             dm_r(2, 4, 1) = -rr0(3)
     231         216 :             dm_r(2, 4, 2) = 0.0_dp
     232         216 :             dm_r(2, 4, 3) = -rr0(1)
     233             : 
     234         216 :             dm_r(3, 3, 1) = (rr(1) + rr0(1))
     235         216 :             dm_r(3, 3, 2) = (rr(2) - rr0(2))
     236         216 :             dm_r(3, 3, 3) = (rr(3) + rr0(3))
     237             : 
     238         216 :             dm_r(3, 4, 1) = 0.0_dp
     239         216 :             dm_r(3, 4, 2) = -rr0(3)
     240         216 :             dm_r(3, 4, 3) = -rr0(2)
     241             : 
     242         216 :             dm_r(4, 4, 1) = (rr(1) + rr0(1))
     243         216 :             dm_r(4, 4, 2) = (rr(2) + rr0(2))
     244         216 :             dm_r(4, 4, 3) = (rr(3) - rr0(3))
     245             : 
     246         864 :             DO ix = 1, 3
     247         648 :                dm_r(2, 1, ix) = dm_r(1, 2, ix)
     248         648 :                dm_r(3, 1, ix) = dm_r(1, 3, ix)
     249         648 :                dm_r(4, 1, ix) = dm_r(1, 4, ix)
     250         648 :                dm_r(3, 2, ix) = dm_r(2, 3, ix)
     251         648 :                dm_r(4, 2, ix) = dm_r(2, 4, ix)
     252         864 :                dm_r(4, 3, ix) = dm_r(3, 4, ix)
     253             :             END DO
     254             :             !
     255         900 :             DO ix = 1, 3
     256         648 :                drmsd3(ix, I) = 0.0_dp
     257        3240 :                DO k = 1, 4
     258       13608 :                   DO j = 1, 4
     259       12960 :                      drmsd3(ix, i) = drmsd3(ix, i) + Q(K - 1)*Q(j - 1)*dm_r(j, k, ix)
     260             :                   END DO
     261             :                END DO
     262        1080 :                drmsd3(ix, I) = drmsd3(ix, I)/REAL(natom, KIND=dp)
     263             :             END DO
     264             :          END DO
     265             :       END IF
     266             :       ! Computes the rotation matrix in terms of quaternions
     267        3092 :       my_rot(1, 1) = -2.0_dp*Q(2)**2 - 2.0_dp*Q(3)**2 + 1.0_dp
     268        3092 :       my_rot(1, 2) = 2.0_dp*(-Q(0)*Q(3) + Q(1)*Q(2))
     269        3092 :       my_rot(1, 3) = 2.0_dp*(Q(0)*Q(2) + Q(1)*Q(3))
     270        3092 :       my_rot(2, 1) = 2.0_dp*(Q(0)*Q(3) + Q(1)*Q(2))
     271        3092 :       my_rot(2, 2) = -2.0_dp*Q(1)**2 - 2.0_dp*Q(3)**2 + 1.0_dp
     272        3092 :       my_rot(2, 3) = 2.0_dp*(-Q(0)*Q(1) + Q(2)*Q(3))
     273        3092 :       my_rot(3, 1) = 2.0_dp*(-Q(0)*Q(2) + Q(1)*Q(3))
     274        3092 :       my_rot(3, 2) = 2.0_dp*(Q(0)*Q(1) + Q(2)*Q(3))
     275        3092 :       my_rot(3, 3) = -2.0_dp*Q(1)**2 - 2.0_dp*Q(2)**2 + 1.0_dp
     276       10112 :       IF (PRESENT(rot)) rot = my_rot
     277             :       ! Give back coordinates rotated in order to minimize the RMSD
     278        3092 :       IF (my_rotate) THEN
     279       61920 :          DO i = 1, natom
     280      297376 :             r((i - 1)*3 + 1:(i - 1)*3 + 3) = MATMUL(TRANSPOSE(my_rot), rp(:, i)) + loc_tr
     281             :          END DO
     282             :       END IF
     283        3092 :       DEALLOCATE (w)
     284        3092 :       DEALLOCATE (rp)
     285        3092 :       DEALLOCATE (r0p)
     286        3092 :    END SUBROUTINE rmsd3
     287             : 
     288             : END MODULE rmsd

Generated by: LCOV version 1.15