1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROningen Mixture of Alchemy and Childrens' Stories
46 #include "gmx_fatal.h"
47 #include "chargegroup.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 mb,m,cg,a_mol,a0,a1,a;
68 ip = mtop->ffparams.iparams;
71 for(mb=0; mb<mtop->nmolblock; mb++)
73 molb = &mtop->molblock[mb];
74 molt = &mtop->moltype[molb->type];
76 atom = molt->atoms.atom;
77 for(m=0; m<molb->nmol; m++)
79 for(cg=0; cg<cgs->nr; cg++)
82 a1 = cgs->index[cg+1];
88 rvec_inc(cen,x[a_mol+a]);
90 svmul(1.0/(a1-a0),cen,cen);
93 r2 = distance2(cen,x[a_mol+a]);
95 (ip[atom[a].type ].lj.c6 != 0 ||
96 ip[atom[a].type ].lj.c12 != 0 ||
97 ip[atom[a].typeB].lj.c6 != 0 ||
98 ip[atom[a].typeB].lj.c12 != 0))
111 (atom[a].q != 0 || atom[a].qB != 0))
126 a_mol += molb->natoms_mol;
132 *rcoul1 = sqrt(r2c1);
133 *rcoul2 = sqrt(r2c2);
136 void calc_cgcm(FILE *fplog,int cg0,int cg1,t_block *cgs,
137 rvec pos[],rvec cg_cm[])
145 fprintf(fplog,"Calculating centre of geometry for charge groups %d to %d\n",
148 cgindex = cgs->index;
150 /* Compute the center of geometry for all charge groups */
151 for(icg=cg0; (icg<cg1); icg++) {
156 copy_rvec(pos[k0],cg_cm[icg]);
162 for(k=k0; (k<k1); k++) {
163 for(d=0; (d<DIM); d++)
166 for(d=0; (d<DIM); d++)
167 cg_cm[icg][d] = inv_ncg*cg[d];
172 void put_charge_groups_in_box(FILE *fplog,int cg0,int cg1,
173 int ePBC,matrix box,t_block *cgs,
174 rvec pos[],rvec cg_cm[])
177 int npbcdim,icg,k,k0,k1,d,e;
183 if (ePBC == epbcNONE)
184 gmx_incons("Calling put_charge_groups_in_box for a system without PBC");
187 fprintf(fplog,"Putting cgs %d to %d in box\n",cg0,cg1);
189 cgindex = cgs->index;
196 bTric = TRICLINIC(box);
198 for(icg=cg0; (icg<cg1); icg++) {
199 /* First compute the center of geometry for this charge group */
205 copy_rvec(pos[k0],cg_cm[icg]);
210 for(k=k0; (k<k1); k++) {
215 cg_cm[icg][d] = inv_ncg*cg[d];
217 /* Now check pbc for this cg */
219 for(d=npbcdim-1; d>=0; d--) {
220 while(cg_cm[icg][d] < 0) {
221 for(e=d; e>=0; e--) {
222 cg_cm[icg][e] += box[d][e];
223 for(k=k0; (k<k1); k++)
224 pos[k][e] += box[d][e];
227 while(cg_cm[icg][d] >= box[d][d]) {
228 for(e=d; e>=0; e--) {
229 cg_cm[icg][e] -= box[d][e];
230 for(k=k0; (k<k1); k++)
231 pos[k][e] -= box[d][e];
236 for(d=0; d<npbcdim; d++) {
237 while(cg_cm[icg][d] < 0) {
238 cg_cm[icg][d] += box[d][d];
239 for(k=k0; (k<k1); k++)
240 pos[k][d] += box[d][d];
242 while(cg_cm[icg][d] >= box[d][d]) {
243 cg_cm[icg][d] -= box[d][d];
244 for(k=k0; (k<k1); k++)
245 pos[k][d] -= box[d][d];
250 for(d=0; (d<npbcdim); d++) {
251 if ((cg_cm[icg][d] < 0) || (cg_cm[icg][d] >= box[d][d]))
252 gmx_fatal(FARGS,"cg_cm[%d] = %15f %15f %15f\n"
253 "box = %15f %15f %15f\n",
254 icg,cg_cm[icg][XX],cg_cm[icg][YY],cg_cm[icg][ZZ],
255 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);