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, 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.
37 /* This file is completely threadsafe - keep it that way! */
44 #include "gromacs/linearalgebra/nrjac.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/topology/topology.h"
47 #include "gromacs/utility/smalloc.h"
52 static void m_op(matrix mat
, rvec x
)
57 for (m
= 0; (m
< DIM
); m
++)
59 xp
[m
] = mat
[m
][XX
]*x
[XX
]+mat
[m
][YY
]*x
[YY
]+mat
[m
][ZZ
]*x
[ZZ
];
61 fprintf(stderr
, "x %8.3f %8.3f %8.3f\n", x
[XX
], x
[YY
], x
[ZZ
]);
62 fprintf(stderr
, "xp %8.3f %8.3f %8.3f\n", xp
[XX
], xp
[YY
], xp
[ZZ
]);
63 fprintf(stderr
, "fac %8.3f %8.3f %8.3f\n", xp
[XX
]/x
[XX
], xp
[YY
]/x
[YY
],
67 static void ptrans(char *s
, real
**inten
, real d
[], real e
[])
71 for (m
= 1; (m
< NDIM
); m
++)
77 fprintf(stderr
, "%8s %8.3f %8.3f %8.3f, norm:%8.3f, d:%8.3f, e:%8.3f\n",
78 s
, x
, y
, z
, std::sqrt(n
), d
[m
], e
[m
]);
80 fprintf(stderr
, "\n");
83 void t_trans(matrix trans
, real d
[], real
**ev
)
87 for (j
= 0; (j
< DIM
); j
++)
93 fprintf(stderr
, "d[%d]=%g\n", j
, d
[j
+1]);
98 void principal_comp(int n
, int index
[], t_atom atom
[], rvec x
[],
101 int i
, j
, ai
, m
, nrot
;
103 double **inten
, dd
[NDIM
], tvec
[NDIM
], **ev
;
111 for (i
= 0; (i
< NDIM
); i
++)
113 snew(inten
[i
], NDIM
);
121 for (i
= 0; (i
< NDIM
); i
++)
123 for (m
= 0; (m
< NDIM
); m
++)
128 for (i
= 0; (i
< n
); i
++)
135 inten
[0][0] += mm
*(sqr(ry
)+sqr(rz
));
136 inten
[1][1] += mm
*(sqr(rx
)+sqr(rz
));
137 inten
[2][2] += mm
*(sqr(rx
)+sqr(ry
));
138 inten
[1][0] -= mm
*(ry
*rx
);
139 inten
[2][0] -= mm
*(rx
*rz
);
140 inten
[2][1] -= mm
*(rz
*ry
);
142 inten
[0][1] = inten
[1][0];
143 inten
[0][2] = inten
[2][0];
144 inten
[1][2] = inten
[2][1];
146 ptrans("initial", inten
, dd
, e
);
149 for (i
= 0; (i
< DIM
); i
++)
151 for (m
= 0; (m
< DIM
); m
++)
153 trans
[i
][m
] = inten
[i
][m
];
157 /* Call numerical recipe routines */
158 jacobi(inten
, 3, dd
, ev
, &nrot
);
160 ptrans("jacobi", ev
, dd
, e
);
163 /* Sort eigenvalues in ascending order */
165 if (std::abs(dd[i+1]) < std::abs(dd[i])) { \
167 for (j = 0; (j < NDIM); j++) { tvec[j] = ev[j][i]; } \
169 for (j = 0; (j < NDIM); j++) { ev[j][i] = ev[j][i+1]; } \
171 for (j = 0; (j < NDIM); j++) { ev[j][i+1] = tvec[j]; } \
177 ptrans("swap", ev
, dd
, e
);
178 t_trans(trans
, dd
, ev
);
181 for (i
= 0; (i
< DIM
); i
++)
184 for (m
= 0; (m
< DIM
); m
++)
186 trans
[i
][m
] = ev
[m
][i
];
190 for (i
= 0; (i
< NDIM
); i
++)
199 void rotate_atoms(int gnx
, int *index
, rvec x
[], matrix trans
)
204 for (i
= 0; (i
< gnx
); i
++)
206 ii
= index
? index
[i
] : i
;
210 x
[ii
][XX
] = trans
[XX
][XX
]*xt
+trans
[XX
][YY
]*yt
+trans
[XX
][ZZ
]*zt
;
211 x
[ii
][YY
] = trans
[YY
][XX
]*xt
+trans
[YY
][YY
]*yt
+trans
[YY
][ZZ
]*zt
;
212 x
[ii
][ZZ
] = trans
[ZZ
][XX
]*xt
+trans
[ZZ
][YY
]*yt
+trans
[ZZ
][ZZ
]*zt
;
216 real
calc_xcm(rvec x
[], int gnx
, int *index
, t_atom
*atom
, rvec xcm
,
224 for (i
= 0; (i
< gnx
); i
++)
226 ii
= index
? index
[i
] : i
;
231 m0
= std::abs(atom
[ii
].q
);
243 for (m
= 0; (m
< DIM
); m
++)
245 xcm
[m
] += m0
*x
[ii
][m
];
248 for (m
= 0; (m
< DIM
); m
++)
256 real
sub_xcm(rvec x
[], int gnx
, int *index
, t_atom atom
[], rvec xcm
,
262 tm
= calc_xcm(x
, gnx
, index
, atom
, xcm
, bQ
);
263 for (i
= 0; (i
< gnx
); i
++)
265 ii
= index
? index
[i
] : i
;
266 rvec_dec(x
[ii
], xcm
);
271 void add_xcm(rvec x
[], int gnx
, int *index
, rvec xcm
)
275 for (i
= 0; (i
< gnx
); i
++)
277 ii
= index
? index
[i
] : i
;
278 rvec_inc(x
[ii
], xcm
);
282 void orient_princ(t_atoms
*atoms
, int isize
, int *index
,
283 int natoms
, rvec x
[], rvec
*v
, rvec d
)
289 calc_xcm(x
, isize
, index
, atoms
->atom
, xcm
, FALSE
);
290 for (i
= 0; i
< natoms
; i
++)
294 principal_comp(isize
, index
, atoms
->atom
, x
, trans
, prcomp
);
297 copy_rvec(prcomp
, d
);
300 /* Check whether this trans matrix mirrors the molecule */
303 for (m
= 0; (m
< DIM
); m
++)
305 trans
[ZZ
][m
] = -trans
[ZZ
][m
];
308 rotate_atoms(natoms
, NULL
, x
, trans
);
311 rotate_atoms(natoms
, NULL
, v
, trans
);
314 for (i
= 0; i
< natoms
; i
++)