2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "pbcmethods.h"
42 #include "gromacs/pbcutil/pbc.h"
43 #include "gromacs/topology/topology.h"
44 #include "gromacs/utility/fatalerror.h"
45 #include "gromacs/utility/smalloc.h"
47 void calc_pbc_cluster(int ecenter
, int nrefat
, t_topology
* top
, PbcType pbcType
, rvec x
[], const int index
[], matrix box
)
49 int m
, i
, j
, j0
, j1
, jj
, ai
, aj
;
52 rvec dx
, xtest
, box_center
;
53 int nmol
, imol_center
;
55 gmx_bool
*bMol
, *bTmp
;
56 rvec
* m_com
, *m_shift
;
63 calc_box_center(ecenter
, box
, box_center
);
65 /* Initiate the pbc structure */
66 std::memset(&pbc
, 0, sizeof(pbc
));
67 set_pbc(&pbc
, pbcType
, box
);
69 /* Convert atom index to molecular */
71 molind
= top
->mols
.index
;
77 snew(bTmp
, top
->atoms
.nr
);
79 for (i
= 0; (i
< nrefat
); i
++)
81 /* Mark all molecules in the index */
84 /* Binary search assuming the molecules are sorted */
89 if (ai
< molind
[j0
+ 1])
93 else if (ai
>= molind
[j1
])
100 if (ai
< molind
[jj
+ 1])
112 /* Double check whether all atoms in all molecules that are marked are part
113 * of the cluster. Simultaneously compute the center of geometry.
115 min_dist2
= 10 * gmx::square(trace(box
));
118 for (i
= 0; i
< nmol
; i
++)
120 for (j
= molind
[i
]; j
< molind
[i
+ 1]; j
++)
122 if (bMol
[i
] && !bTmp
[j
])
125 "Molecule %d marked for clustering but not atom %d in it - check your "
129 else if (!bMol
[i
] && bTmp
[j
])
132 "Atom %d marked for clustering but not molecule %d - this is an internal "
138 /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
141 pbc_dx(&pbc
, x
[j
], x
[j
- 1], dx
);
142 rvec_add(x
[j
- 1], dx
, x
[j
]);
144 /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
145 rvec_inc(m_com
[i
], x
[j
]);
150 /* Normalize center of geometry */
151 fac
= 1.0 / (molind
[i
+ 1] - molind
[i
]);
152 for (m
= 0; (m
< DIM
); m
++)
156 /* Determine which molecule is closest to the center of the box */
157 pbc_dx(&pbc
, box_center
, m_com
[i
], dx
);
158 tmp_r2
= iprod(dx
, dx
);
160 if (tmp_r2
< min_dist2
)
165 cluster
[ncluster
++] = i
;
172 fprintf(stderr
, "No molecules selected in the cluster\n");
175 else if (imol_center
== -1)
177 fprintf(stderr
, "No central molecules could be found\n");
182 added
[nadded
++] = imol_center
;
183 bMol
[imol_center
] = FALSE
;
185 while (nadded
< ncluster
)
187 /* Find min distance between cluster molecules and those remaining to be added */
188 min_dist2
= 10 * gmx::square(trace(box
));
191 /* Loop over added mols */
192 for (i
= 0; i
< nadded
; i
++)
195 /* Loop over all mols */
196 for (j
= 0; j
< ncluster
; j
++)
199 /* check those remaining to be added */
202 pbc_dx(&pbc
, m_com
[aj
], m_com
[ai
], dx
);
203 tmp_r2
= iprod(dx
, dx
);
204 if (tmp_r2
< min_dist2
)
214 /* Add the best molecule */
215 added
[nadded
++] = jmin
;
217 /* Calculate the shift from the ai molecule */
218 pbc_dx(&pbc
, m_com
[jmin
], m_com
[imin
], dx
);
219 rvec_add(m_com
[imin
], dx
, xtest
);
220 rvec_sub(xtest
, m_com
[jmin
], m_shift
[jmin
]);
221 rvec_inc(m_com
[jmin
], m_shift
[jmin
]);
223 for (j
= molind
[jmin
]; j
< molind
[jmin
+ 1]; j
++)
225 rvec_inc(x
[j
], m_shift
[jmin
]);
227 fprintf(stdout
, "\rClustering iteration %d of %d...", nadded
, ncluster
);
237 fprintf(stdout
, "\n");
240 void put_molecule_com_in_box(int unitcell_enum
,
251 rvec com
, shift
, box_center
;
256 calc_box_center(ecenter
, box
, box_center
);
257 set_pbc(&pbc
, pbcType
, box
);
261 "There are no molecule descriptions. I need a .tpr file for this pbc option.");
263 for (i
= 0; (i
< mols
->nr
); i
++)
268 for (j
= mols
->index
[i
]; (j
< mols
->index
[i
+ 1] && j
< natoms
); j
++)
271 for (d
= 0; d
< DIM
; d
++)
273 com
[d
] += m
* x
[j
][d
];
277 /* calculate final COM */
278 svmul(1.0 / mtot
, com
, com
);
280 /* check if COM is outside box */
282 copy_rvec(com
, newCom
);
283 auto newComArrayRef
= gmx::arrayRefFromArray(&newCom
, 1);
284 switch (unitcell_enum
)
286 case euRect
: put_atoms_in_box(pbcType
, box
, newComArrayRef
); break;
287 case euTric
: put_atoms_in_triclinic_unitcell(ecenter
, box
, newComArrayRef
); break;
289 put_atoms_in_compact_unitcell(pbcType
, ecenter
, box
, newComArrayRef
);
292 rvec_sub(newCom
, com
, shift
);
293 if (norm2(shift
) > 0)
298 "\nShifting position of molecule %d "
299 "by %8.3f %8.3f %8.3f\n",
300 i
+ 1, shift
[XX
], shift
[YY
], shift
[ZZ
]);
302 for (j
= mols
->index
[i
]; (j
< mols
->index
[i
+ 1] && j
< natoms
); j
++)
304 rvec_inc(x
[j
], shift
);
310 void put_residue_com_in_box(int unitcell_enum
,
318 int i
, j
, res_start
, res_end
;
322 rvec box_center
, com
, shift
;
323 static const int NOTSET
= -12347;
324 calc_box_center(ecenter
, box
, box_center
);
330 for (i
= 0; i
< natoms
+ 1; i
++)
332 if (i
== natoms
|| (presnr
!= atom
[i
].resind
&& presnr
!= NOTSET
))
334 /* calculate final COM */
336 svmul(1.0 / mtot
, com
, com
);
338 /* check if COM is outside box */
340 copy_rvec(com
, newCom
);
341 auto newComArrayRef
= gmx::arrayRefFromArray(&newCom
, 1);
342 switch (unitcell_enum
)
344 case euRect
: put_atoms_in_box(pbcType
, box
, newComArrayRef
); break;
345 case euTric
: put_atoms_in_triclinic_unitcell(ecenter
, box
, newComArrayRef
); break;
347 put_atoms_in_compact_unitcell(pbcType
, ecenter
, box
, newComArrayRef
);
350 rvec_sub(newCom
, com
, shift
);
351 if (norm2(shift
) != 0.0F
)
356 "\nShifting position of residue %d (atoms %d-%d) "
358 atom
[res_start
].resind
+ 1, res_start
+ 1, res_end
+ 1, shift
[XX
],
359 shift
[YY
], shift
[ZZ
]);
361 for (j
= res_start
; j
< res_end
; j
++)
363 rvec_inc(x
[j
], shift
);
369 /* remember start of new residue */
376 for (d
= 0; d
< DIM
; d
++)
378 com
[d
] += m
* x
[i
][d
];
382 presnr
= atom
[i
].resind
;
387 void center_x(int ecenter
, rvec x
[], matrix box
, int n
, int nc
, const int ci
[])
390 rvec cmin
, cmax
, box_center
, dx
;
394 copy_rvec(x
[ci
[0]], cmin
);
395 copy_rvec(x
[ci
[0]], cmax
);
396 for (i
= 0; i
< nc
; i
++)
399 for (m
= 0; m
< DIM
; m
++)
401 if (x
[ai
][m
] < cmin
[m
])
405 else if (x
[ai
][m
] > cmax
[m
])
411 calc_box_center(ecenter
, box
, box_center
);
412 for (m
= 0; m
< DIM
; m
++)
414 dx
[m
] = box_center
[m
] - (cmin
[m
] + cmax
[m
]) * 0.5;
417 for (i
= 0; i
< n
; i
++)