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