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) 2013,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.
39 #include "chargegroup.h"
43 #include "gromacs/math/vec.h"
44 #include "gromacs/pbcutil/pbc.h"
45 #include "gromacs/topology/topology.h"
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/utility/smalloc.h"
50 void calc_chargegroup_radii(const gmx_mtop_t
*mtop
, rvec
*x
,
51 real
*rvdw1
, real
*rvdw2
,
52 real
*rcoul1
, real
*rcoul2
)
54 real r2v1
, r2v2
, r2c1
, r2c2
, r2
;
55 int ntype
, i
, j
, mb
, m
, cg
, a_mol
, a0
, a1
, a
;
68 ntype
= mtop
->ffparams
.atnr
;
70 for (i
= 0; i
< ntype
; i
++)
73 for (j
= 0; j
< ntype
; j
++)
75 if (mtop
->ffparams
.iparams
[i
*ntype
+j
].lj
.c6
!= 0 ||
76 mtop
->ffparams
.iparams
[i
*ntype
+j
].lj
.c12
!= 0)
84 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
86 molb
= &mtop
->molblock
[mb
];
87 molt
= &mtop
->moltype
[molb
->type
];
89 atom
= molt
->atoms
.atom
;
90 for (m
= 0; m
< molb
->nmol
; m
++)
92 for (cg
= 0; cg
< cgs
->nr
; cg
++)
95 a1
= cgs
->index
[cg
+1];
99 for (a
= a0
; a
< a1
; a
++)
101 rvec_inc(cen
, x
[a_mol
+a
]);
103 svmul(1.0/(a1
-a0
), cen
, cen
);
104 for (a
= a0
; a
< a1
; a
++)
106 r2
= distance2(cen
, x
[a_mol
+a
]);
107 if (r2
> r2v2
&& (bLJ
[atom
[a
].type
] ||
121 (atom
[a
].q
!= 0 || atom
[a
].qB
!= 0))
136 a_mol
+= molb
->natoms_mol
;
142 *rvdw1
= std::sqrt(r2v1
);
143 *rvdw2
= std::sqrt(r2v2
);
144 *rcoul1
= std::sqrt(r2c1
);
145 *rcoul2
= std::sqrt(r2c2
);
148 void calc_cgcm(FILE gmx_unused
*fplog
, int cg0
, int cg1
, const t_block
*cgs
,
149 rvec pos
[], rvec cg_cm
[])
151 int icg
, k
, k0
, k1
, d
;
157 fprintf(fplog
, "Calculating centre of geometry for charge groups %d to %d\n",
160 cgindex
= cgs
->index
;
162 /* Compute the center of geometry for all charge groups */
163 for (icg
= cg0
; (icg
< cg1
); icg
++)
170 copy_rvec(pos
[k0
], cg_cm
[icg
]);
177 for (k
= k0
; (k
< k1
); k
++)
179 for (d
= 0; (d
< DIM
); d
++)
184 for (d
= 0; (d
< DIM
); d
++)
186 cg_cm
[icg
][d
] = inv_ncg
*cg
[d
];
192 void put_charge_groups_in_box(FILE gmx_unused
*fplog
, int cg0
, int cg1
,
193 int ePBC
, matrix box
, t_block
*cgs
,
194 rvec pos
[], rvec cg_cm
[])
197 int npbcdim
, icg
, k
, k0
, k1
, d
, e
;
203 if (ePBC
== epbcNONE
)
205 gmx_incons("Calling put_charge_groups_in_box for a system without PBC");
209 fprintf(fplog
, "Putting cgs %d to %d in box\n", cg0
, cg1
);
211 cgindex
= cgs
->index
;
222 bTric
= TRICLINIC(box
);
224 for (icg
= cg0
; (icg
< cg1
); icg
++)
226 /* First compute the center of geometry for this charge group */
233 copy_rvec(pos
[k0
], cg_cm
[icg
]);
240 for (k
= k0
; (k
< k1
); k
++)
242 for (d
= 0; d
< DIM
; d
++)
247 for (d
= 0; d
< DIM
; d
++)
249 cg_cm
[icg
][d
] = inv_ncg
*cg
[d
];
252 /* Now check pbc for this cg */
255 for (d
= npbcdim
-1; d
>= 0; d
--)
257 while (cg_cm
[icg
][d
] < 0)
259 for (e
= d
; e
>= 0; e
--)
261 cg_cm
[icg
][e
] += box
[d
][e
];
262 for (k
= k0
; (k
< k1
); k
++)
264 pos
[k
][e
] += box
[d
][e
];
268 while (cg_cm
[icg
][d
] >= box
[d
][d
])
270 for (e
= d
; e
>= 0; e
--)
272 cg_cm
[icg
][e
] -= box
[d
][e
];
273 for (k
= k0
; (k
< k1
); k
++)
275 pos
[k
][e
] -= box
[d
][e
];
283 for (d
= 0; d
< npbcdim
; d
++)
285 while (cg_cm
[icg
][d
] < 0)
287 cg_cm
[icg
][d
] += box
[d
][d
];
288 for (k
= k0
; (k
< k1
); k
++)
290 pos
[k
][d
] += box
[d
][d
];
293 while (cg_cm
[icg
][d
] >= box
[d
][d
])
295 cg_cm
[icg
][d
] -= box
[d
][d
];
296 for (k
= k0
; (k
< k1
); k
++)
298 pos
[k
][d
] -= box
[d
][d
];
304 for (d
= 0; (d
< npbcdim
); d
++)
306 if ((cg_cm
[icg
][d
] < 0) || (cg_cm
[icg
][d
] >= box
[d
][d
]))
308 gmx_fatal(FARGS
, "cg_cm[%d] = %15f %15f %15f\n"
309 "box = %15f %15f %15f\n",
310 icg
, cg_cm
[icg
][XX
], cg_cm
[icg
][YY
], cg_cm
[icg
][ZZ
],
311 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);