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 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];
101 rvec_inc(cen,x[a_mol+a]);
103 svmul(1.0/(a1-a0),cen,cen);
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;
144 *rcoul1 = sqrt(r2c1);
145 *rcoul2 = sqrt(r2c2);
148 void calc_cgcm(FILE *fplog,int cg0,int cg1,t_block *cgs,
149 rvec pos[],rvec cg_cm[])
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++) {
168 copy_rvec(pos[k0],cg_cm[icg]);
174 for(k=k0; (k<k1); k++) {
175 for(d=0; (d<DIM); d++)
178 for(d=0; (d<DIM); d++)
179 cg_cm[icg][d] = inv_ncg*cg[d];
184 void put_charge_groups_in_box(FILE *fplog,int cg0,int cg1,
185 int ePBC,matrix box,t_block *cgs,
186 rvec pos[],rvec cg_cm[])
189 int npbcdim,icg,k,k0,k1,d,e;
195 if (ePBC == epbcNONE)
196 gmx_incons("Calling put_charge_groups_in_box for a system without PBC");
199 fprintf(fplog,"Putting cgs %d to %d in box\n",cg0,cg1);
201 cgindex = cgs->index;
208 bTric = TRICLINIC(box);
210 for(icg=cg0; (icg<cg1); icg++) {
211 /* First compute the center of geometry for this charge group */
217 copy_rvec(pos[k0],cg_cm[icg]);
222 for(k=k0; (k<k1); k++) {
227 cg_cm[icg][d] = inv_ncg*cg[d];
229 /* Now check pbc for this cg */
231 for(d=npbcdim-1; d>=0; d--) {
232 while(cg_cm[icg][d] < 0) {
233 for(e=d; e>=0; e--) {
234 cg_cm[icg][e] += box[d][e];
235 for(k=k0; (k<k1); k++)
236 pos[k][e] += box[d][e];
239 while(cg_cm[icg][d] >= box[d][d]) {
240 for(e=d; e>=0; e--) {
241 cg_cm[icg][e] -= box[d][e];
242 for(k=k0; (k<k1); k++)
243 pos[k][e] -= box[d][e];
248 for(d=0; d<npbcdim; d++) {
249 while(cg_cm[icg][d] < 0) {
250 cg_cm[icg][d] += box[d][d];
251 for(k=k0; (k<k1); k++)
252 pos[k][d] += box[d][d];
254 while(cg_cm[icg][d] >= box[d][d]) {
255 cg_cm[icg][d] -= box[d][d];
256 for(k=k0; (k<k1); k++)
257 pos[k][d] -= box[d][d];
262 for(d=0; (d<npbcdim); d++) {
263 if ((cg_cm[icg][d] < 0) || (cg_cm[icg][d] >= box[d][d]))
264 gmx_fatal(FARGS,"cg_cm[%d] = %15f %15f %15f\n"
265 "box = %15f %15f %15f\n",
266 icg,cg_cm[icg][XX],cg_cm[icg][YY],cg_cm[icg][ZZ],
267 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);