(git:b195825)
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 ! **************************************************************************************************
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  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 
288 END 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:376
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