(git:374b731)
Loading...
Searching...
No Matches
rmsd.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
27CONTAINS
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 SUBROUTINE rmsd3(particle_set, r, r0, output_unit, weights, my_val, &
52 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w
70 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: r0p, rp
71
72 cpassert(SIZE(r) == SIZE(r0))
73 natom = SIZE(particle_set)
74 my_rotate = .false.
75 IF (PRESENT(rotate)) my_rotate = rotate
76 ALLOCATE (w(natom))
77 IF (PRESENT(weights)) THEN
78 w(:) = weights
79 ELSE
80 ! All atoms have a weight proportional to their mass in the RMSD unless
81 ! differently requested
82 DO i = 1, natom
83 w(i) = particle_set(i)%atomic_kind%mass
84 END DO
85 mtot = minval(w)
86 IF (mtot /= 0.0_dp) w(:) = w(:)/mtot
87 END IF
88 ALLOCATE (rp(3, natom))
89 ALLOCATE (r0p(3, natom))
90 ! Molecule given by coordinates R
91 ! Find COM and center molecule in COM
92 xx = 0.0_dp
93 yy = 0.0_dp
94 zz = 0.0_dp
95 mtot = 0.0_dp
96 DO i = 1, natom
97 mtot = mtot + particle_set(i)%atomic_kind%mass
98 xx = xx + r((i - 1)*3 + 1)*particle_set(i)%atomic_kind%mass
99 yy = yy + r((i - 1)*3 + 2)*particle_set(i)%atomic_kind%mass
100 zz = zz + r((i - 1)*3 + 3)*particle_set(i)%atomic_kind%mass
101 END DO
102 xx = xx/mtot
103 yy = yy/mtot
104 zz = zz/mtot
105 DO i = 1, natom
106 rp(1, i) = r((i - 1)*3 + 1) - xx
107 rp(2, i) = r((i - 1)*3 + 2) - yy
108 rp(3, i) = r((i - 1)*3 + 3) - zz
109 END DO
110 IF (PRESENT(transl)) THEN
111 transl(1) = xx
112 transl(2) = yy
113 transl(3) = zz
114 END IF
115 ! Molecule given by coordinates R0
116 ! Find COM and center molecule in COM
117 xx = 0.0_dp
118 yy = 0.0_dp
119 zz = 0.0_dp
120 DO i = 1, natom
121 xx = xx + r0((i - 1)*3 + 1)*particle_set(i)%atomic_kind%mass
122 yy = yy + r0((i - 1)*3 + 2)*particle_set(i)%atomic_kind%mass
123 zz = zz + r0((i - 1)*3 + 3)*particle_set(i)%atomic_kind%mass
124 END DO
125 xx = xx/mtot
126 yy = yy/mtot
127 zz = zz/mtot
128 DO i = 1, natom
129 r0p(1, i) = r0((i - 1)*3 + 1) - xx
130 r0p(2, i) = r0((i - 1)*3 + 2) - yy
131 r0p(3, i) = r0((i - 1)*3 + 3) - zz
132 END DO
133 loc_tr(1) = xx
134 loc_tr(2) = yy
135 loc_tr(3) = zz
136 ! Give back the translational vector
137 IF (PRESENT(transl)) THEN
138 transl(1) = transl(1) - xx
139 transl(2) = transl(2) - yy
140 transl(3) = transl(3) - zz
141 END IF
142 m = 0.0_dp
143 !
144 DO i = 1, natom
145 IF (w(i) == 0.0_dp) cycle
146 rr(1) = rp(1, i)
147 rr(2) = rp(2, i)
148 rr(3) = rp(3, i)
149 rr0(1) = r0p(1, i)
150 rr0(2) = r0p(2, i)
151 rr0(3) = r0p(3, i)
152 rrsq = w(i)*(rr0(1)**2 + rr0(2)**2 + rr0(3)**2 + rr(1)**2 + rr(2)**2 + rr(3)**2)
153 rr0(1) = w(i)*rr0(1)
154 rr0(2) = w(i)*rr0(2)
155 rr0(3) = w(i)*rr0(3)
156 m(1, 1) = m(1, 1) + rrsq + 2.0_dp*(-rr0(1)*rr(1) - rr0(2)*rr(2) - rr0(3)*rr(3))
157 m(2, 2) = m(2, 2) + rrsq + 2.0_dp*(-rr0(1)*rr(1) + rr0(2)*rr(2) + rr0(3)*rr(3))
158 m(3, 3) = m(3, 3) + rrsq + 2.0_dp*(rr0(1)*rr(1) - rr0(2)*rr(2) + rr0(3)*rr(3))
159 m(4, 4) = m(4, 4) + rrsq + 2.0_dp*(rr0(1)*rr(1) + rr0(2)*rr(2) - rr0(3)*rr(3))
160 m(1, 2) = m(1, 2) + 2.0_dp*(-rr0(2)*rr(3) + rr0(3)*rr(2))
161 m(1, 3) = m(1, 3) + 2.0_dp*(rr0(1)*rr(3) - rr0(3)*rr(1))
162 m(1, 4) = m(1, 4) + 2.0_dp*(-rr0(1)*rr(2) + rr0(2)*rr(1))
163 m(2, 3) = m(2, 3) - 2.0_dp*(rr0(1)*rr(2) + rr0(2)*rr(1))
164 m(2, 4) = m(2, 4) - 2.0_dp*(rr0(1)*rr(3) + rr0(3)*rr(1))
165 m(3, 4) = m(3, 4) - 2.0_dp*(rr0(2)*rr(3) + rr0(3)*rr(2))
166 END DO
167 ! Symmetrize
168 m(2, 1) = m(1, 2)
169 m(3, 1) = m(1, 3)
170 m(3, 2) = m(2, 3)
171 m(4, 1) = m(1, 4)
172 m(4, 2) = m(2, 4)
173 m(4, 3) = m(3, 4)
174 ! Solve the eigenvalue problem for M
175 z = m
176 CALL diamat_all(z, lambda)
177 ! Pick the correct eigenvectors
178 s = 1.0_dp
179 IF (z(1, 1) .LT. 0.0_dp) s = -1.0_dp
180 q(0) = s*z(1, 1)
181 q(1) = s*z(2, 1)
182 q(2) = s*z(3, 1)
183 q(3) = s*z(4, 1)
184 IF (PRESENT(my_val)) THEN
185 IF (abs(lambda(1)) < epsi) THEN
186 my_val = 0.0_dp
187 ELSE
188 my_val = lambda(1)/real(natom, kind=dp)
189 END IF
190 END IF
191 IF (abs(lambda(1) - lambda(2)) < epsi) THEN
192 IF (output_unit > 0) WRITE (output_unit, fmt='(T2,"RMSD3|",A)') &
193 'NORMAL EXECUTION, NON-UNIQUE SOLUTION'
194 END IF
195 ! Computes derivatives w.r.t. the positions if requested
196 IF (PRESENT(drmsd3)) THEN
197 DO i = 1, natom
198 IF (w(i) == 0.0_dp) cycle
199 rr(1) = w(i)*2.0_dp*rp(1, i)
200 rr(2) = w(i)*2.0_dp*rp(2, i)
201 rr(3) = w(i)*2.0_dp*rp(3, i)
202 rr0(1) = w(i)*2.0_dp*r0p(1, i)
203 rr0(2) = w(i)*2.0_dp*r0p(2, i)
204 rr0(3) = w(i)*2.0_dp*r0p(3, i)
205
206 dm_r(1, 1, 1) = (rr(1) - rr0(1))
207 dm_r(1, 1, 2) = (rr(2) - rr0(2))
208 dm_r(1, 1, 3) = (rr(3) - rr0(3))
209
210 dm_r(1, 2, 1) = 0.0_dp
211 dm_r(1, 2, 2) = rr0(3)
212 dm_r(1, 2, 3) = -rr0(2)
213
214 dm_r(1, 3, 1) = -rr0(3)
215 dm_r(1, 3, 2) = 0.0_dp
216 dm_r(1, 3, 3) = rr0(1)
217
218 dm_r(1, 4, 1) = rr0(2)
219 dm_r(1, 4, 2) = -rr0(1)
220 dm_r(1, 4, 3) = 0.0_dp
221
222 dm_r(2, 2, 1) = (rr(1) - rr0(1))
223 dm_r(2, 2, 2) = (rr(2) + rr0(2))
224 dm_r(2, 2, 3) = (rr(3) + rr0(3))
225
226 dm_r(2, 3, 1) = -rr0(2)
227 dm_r(2, 3, 2) = -rr0(1)
228 dm_r(2, 3, 3) = 0.0_dp
229
230 dm_r(2, 4, 1) = -rr0(3)
231 dm_r(2, 4, 2) = 0.0_dp
232 dm_r(2, 4, 3) = -rr0(1)
233
234 dm_r(3, 3, 1) = (rr(1) + rr0(1))
235 dm_r(3, 3, 2) = (rr(2) - rr0(2))
236 dm_r(3, 3, 3) = (rr(3) + rr0(3))
237
238 dm_r(3, 4, 1) = 0.0_dp
239 dm_r(3, 4, 2) = -rr0(3)
240 dm_r(3, 4, 3) = -rr0(2)
241
242 dm_r(4, 4, 1) = (rr(1) + rr0(1))
243 dm_r(4, 4, 2) = (rr(2) + rr0(2))
244 dm_r(4, 4, 3) = (rr(3) - rr0(3))
245
246 DO ix = 1, 3
247 dm_r(2, 1, ix) = dm_r(1, 2, ix)
248 dm_r(3, 1, ix) = dm_r(1, 3, ix)
249 dm_r(4, 1, ix) = dm_r(1, 4, ix)
250 dm_r(3, 2, ix) = dm_r(2, 3, ix)
251 dm_r(4, 2, ix) = dm_r(2, 4, ix)
252 dm_r(4, 3, ix) = dm_r(3, 4, ix)
253 END DO
254 !
255 DO ix = 1, 3
256 drmsd3(ix, i) = 0.0_dp
257 DO k = 1, 4
258 DO j = 1, 4
259 drmsd3(ix, i) = drmsd3(ix, i) + q(k - 1)*q(j - 1)*dm_r(j, k, ix)
260 END DO
261 END DO
262 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 my_rot(1, 1) = -2.0_dp*q(2)**2 - 2.0_dp*q(3)**2 + 1.0_dp
268 my_rot(1, 2) = 2.0_dp*(-q(0)*q(3) + q(1)*q(2))
269 my_rot(1, 3) = 2.0_dp*(q(0)*q(2) + q(1)*q(3))
270 my_rot(2, 1) = 2.0_dp*(q(0)*q(3) + q(1)*q(2))
271 my_rot(2, 2) = -2.0_dp*q(1)**2 - 2.0_dp*q(3)**2 + 1.0_dp
272 my_rot(2, 3) = 2.0_dp*(-q(0)*q(1) + q(2)*q(3))
273 my_rot(3, 1) = 2.0_dp*(-q(0)*q(2) + q(1)*q(3))
274 my_rot(3, 2) = 2.0_dp*(q(0)*q(1) + q(2)*q(3))
275 my_rot(3, 3) = -2.0_dp*q(1)**2 - 2.0_dp*q(2)**2 + 1.0_dp
276 IF (PRESENT(rot)) rot = my_rot
277 ! Give back coordinates rotated in order to minimize the RMSD
278 IF (my_rotate) THEN
279 DO i = 1, natom
280 r((i - 1)*3 + 1:(i - 1)*3 + 3) = matmul(transpose(my_rot), rp(:, i)) + loc_tr
281 END DO
282 END IF
283 DEALLOCATE (w)
284 DEALLOCATE (rp)
285 DEALLOCATE (r0p)
286 END SUBROUTINE rmsd3
287
288END MODULE rmsd
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:372
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
rmsd
Defines functions to perform rmsd in 3D.
Definition rmsd.F:12
rmsd::rmsd3
subroutine, public rmsd3(particle_set, r, r0, output_unit, weights, my_val, rotate, transl, rot, drmsd3)
Computes the RMSD in 3D. Provides also derivatives.
Definition rmsd.F:53
particle_types::particle_type
Definition particle_types.F:35