2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2012,2014,2015,2016,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/linearalgebra/nrjac.h"
46 #include "gromacs/math/functions.h"
47 #include "gromacs/math/utilities.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/utility/smalloc.h"
52 real
calc_similar_ind(gmx_bool bRho
, int nind
, int *index
, real mass
[],
56 real m
, tm
, xs
, xd
, rs
, rd
;
61 for (j
= 0; j
< nind
; j
++)
73 for (d
= 0; d
< DIM
; d
++)
75 xd
= x
[i
][d
] - xp
[i
][d
];
76 rd
+= m
* gmx::square(xd
);
79 xs
= x
[i
][d
] + xp
[i
][d
];
80 rs
+= m
* gmx::square(xs
);
86 return 2*std::sqrt(rd
/rs
);
90 return std::sqrt(rd
/tm
);
94 real
rmsdev_ind(int nind
, int index
[], real mass
[], rvec x
[], rvec xp
[])
96 return calc_similar_ind(FALSE
, nind
, index
, mass
, x
, xp
);
99 real
rmsdev(int natoms
, real mass
[], rvec x
[], rvec xp
[])
101 return calc_similar_ind(FALSE
, natoms
, nullptr, mass
, x
, xp
);
104 real
rhodev_ind(int nind
, int index
[], real mass
[], rvec x
[], rvec xp
[])
106 return calc_similar_ind(TRUE
, nind
, index
, mass
, x
, xp
);
109 real
rhodev(int natoms
, real mass
[], rvec x
[], rvec xp
[])
111 return calc_similar_ind(TRUE
, natoms
, nullptr, mass
, x
, xp
);
114 void calc_fit_R(int ndim
, int natoms
, real
*w_rls
, const rvec
*xp
, rvec
*x
, matrix R
)
116 int c
, r
, n
, j
, i
, irot
, s
;
117 double **omega
, **om
;
118 double d
[2*DIM
], xnr
, xpc
;
124 if (ndim
!= 3 && ndim
!= 2)
126 gmx_fatal(FARGS
, "calc_fit_R called with ndim=%d instead of 3 or 2", ndim
);
131 for (i
= 0; i
< 2*ndim
; i
++)
133 snew(omega
[i
], 2*ndim
);
137 for (i
= 0; i
< 2*ndim
; i
++)
140 for (j
= 0; j
< 2*ndim
; j
++)
147 /*calculate the matrix U*/
149 for (n
= 0; (n
< natoms
); n
++)
151 if ((mn
= w_rls
[n
]) != 0.0)
153 for (c
= 0; (c
< ndim
); c
++)
156 for (r
= 0; (r
< ndim
); r
++)
159 u
[c
][r
] += mn
*xnr
*xpc
;
166 /*omega is symmetric -> omega==omega' */
167 for (r
= 0; r
< 2*ndim
; r
++)
169 for (c
= 0; c
<= r
; c
++)
171 if (r
>= ndim
&& c
< ndim
)
173 omega
[r
][c
] = u
[r
-ndim
][c
];
174 omega
[c
][r
] = u
[r
-ndim
][c
];
184 /*determine h and k*/
185 jacobi(omega
, 2*ndim
, d
, om
, &irot
);
186 /*real **omega = input matrix a[0..n-1][0..n-1] must be symmetric
187 * int natoms = number of rows and columns
188 * real NULL = d[0]..d[n-1] are the eigenvalues of a[][]
189 * real **v = v[0..n-1][0..n-1] contains the vectors in columns
190 * int *irot = number of jacobi rotations
193 if (debug
&& irot
== 0)
195 fprintf(debug
, "IROT=0\n");
198 index
= 0; /* For the compiler only */
200 /* Copy only the first ndim-1 eigenvectors */
201 for (j
= 0; j
< ndim
-1; j
++)
204 for (i
= 0; i
< 2*ndim
; i
++)
213 for (i
= 0; i
< ndim
; i
++)
215 vh
[j
][i
] = M_SQRT2
*om
[i
][index
];
216 vk
[j
][i
] = M_SQRT2
*om
[i
+ndim
][index
];
221 /* Calculate the last eigenvector as the outer-product of the first two.
222 * This insures that the conformation is not mirrored and
223 * prevents problems with completely flat reference structures.
225 cprod(vh
[0], vh
[1], vh
[2]);
226 cprod(vk
[0], vk
[1], vk
[2]);
230 /* Calculate the last eigenvector from the first one */
231 vh
[1][XX
] = -vh
[0][YY
];
232 vh
[1][YY
] = vh
[0][XX
];
233 vk
[1][XX
] = -vk
[0][YY
];
234 vk
[1][YY
] = vk
[0][XX
];
239 for (r
= 0; r
< ndim
; r
++)
241 for (c
= 0; c
< ndim
; c
++)
243 for (s
= 0; s
< ndim
; s
++)
245 R
[r
][c
] += vk
[s
][r
]*vh
[s
][c
];
249 for (r
= ndim
; r
< DIM
; r
++)
254 for (i
= 0; i
< 2*ndim
; i
++)
263 void do_fit_ndim(int ndim
, int natoms
, real
*w_rls
, const rvec
*xp
, rvec
*x
)
269 /* Calculate the rotation matrix R */
270 calc_fit_R(ndim
, natoms
, w_rls
, xp
, x
, R
);
273 for (j
= 0; j
< natoms
; j
++)
275 for (m
= 0; m
< DIM
; m
++)
279 for (r
= 0; r
< DIM
; r
++)
282 for (c
= 0; c
< DIM
; c
++)
284 x
[j
][r
] += R
[r
][c
]*x_old
[c
];
290 void do_fit(int natoms
, real
*w_rls
, const rvec
*xp
, rvec
*x
)
292 do_fit_ndim(3, natoms
, w_rls
, xp
, x
);
295 void reset_x_ndim(int ndim
, int ncm
, const int *ind_cm
,
296 int nreset
, const int *ind_reset
,
297 rvec x
[], const real mass
[])
305 gmx_incons("More than 3 dimensions not supported.");
309 if (ind_cm
!= nullptr)
311 for (i
= 0; i
< ncm
; i
++)
315 for (m
= 0; m
< ndim
; m
++)
317 xcm
[m
] += mm
*x
[ai
][m
];
324 for (i
= 0; i
< ncm
; i
++)
327 for (m
= 0; m
< ndim
; m
++)
329 xcm
[m
] += mm
*x
[i
][m
];
334 for (m
= 0; m
< ndim
; m
++)
339 if (ind_reset
!= nullptr)
341 for (i
= 0; i
< nreset
; i
++)
343 rvec_dec(x
[ind_reset
[i
]], xcm
);
348 for (i
= 0; i
< nreset
; i
++)
355 void reset_x(int ncm
, const int *ind_cm
,
356 int nreset
, const int *ind_reset
,
357 rvec x
[], const real mass
[])
359 reset_x_ndim(3, ncm
, ind_cm
, nreset
, ind_reset
, x
, mass
);