Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 6184 : 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 2 : my_val = 0.0_dp
187 : ELSE
188 1054 : 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
|
