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