src/gmxlib/chargegroup.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  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.
  15
  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.
  20  *
  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.
  27  *
  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.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include "sysstuff.h"
  42 #include "typedefs.h"
  43 #include "vec.h"
  44 #include "pbc.h"
  45 #include "smalloc.h"
  46 #include "gmx_fatal.h"
  47 #include "chargegroup.h"
  48
  49
  50 void calc_chargegroup_radii(const gmx_mtop_t *mtop,rvec *x,
  51                             real *rvdw1,real *rvdw2,
  52                             real *rcoul1,real *rcoul2)
  53 {
  54     real r2v1,r2v2,r2c1,r2c2,r2;
  55     int  ntype,i,j,mb,m,cg,a_mol,a0,a1,a;
  56     gmx_bool *bLJ;
  57     gmx_molblock_t *molb;
  58     gmx_moltype_t *molt;
  59     t_block *cgs;
  60     t_atom *atom;
  61     rvec cen;
  62
  63     r2v1 = 0;
  64     r2v2 = 0;
  65     r2c1 = 0;
  66     r2c2 = 0;
  67
  68     ntype = mtop->ffparams.atnr;
  69     snew(bLJ,ntype);
  70     for(i=0; i<ntype; i++)
  71     {
  72         bLJ[i] = FALSE;
  73         for(j=0; j<ntype; j++)
  74         {
  75             if (mtop->ffparams.iparams[i*ntype+j].lj.c6  != 0 ||
  76                 mtop->ffparams.iparams[i*ntype+j].lj.c12 != 0)
  77             {
  78                 bLJ[i] = TRUE;
  79             }
  80         }
  81     }
  82
  83     a_mol = 0;
  84     for(mb=0; mb<mtop->nmolblock; mb++)
  85     {
  86         molb = &mtop->molblock[mb];
  87         molt = &mtop->moltype[molb->type];
  88         cgs  = &molt->cgs;
  89         atom = molt->atoms.atom;
  90         for(m=0; m<molb->nmol; m++)
  91         {
  92             for(cg=0; cg<cgs->nr; cg++)
  93             {
  94                 a0 = cgs->index[cg];
  95                 a1 = cgs->index[cg+1];
  96                 if (a1 - a0 > 1)
  97                 {
  98                     clear_rvec(cen);
  99                     for(a=a0; a<a1; a++)
 100                     {
 101                         rvec_inc(cen,x[a_mol+a]);
 102                     }
 103                     svmul(1.0/(a1-a0),cen,cen);
 104                     for(a=a0; a<a1; a++)
 105                     {
 106                         r2 = distance2(cen,x[a_mol+a]);
 107                         if (r2 > r2v2 && (bLJ[atom[a].type ] ||
 108                                           bLJ[atom[a].typeB]))
 109                         {
 110                             if (r2 > r2v1)
 111                             {
 112                                 r2v2 = r2v1;
 113                                 r2v1 = r2;
 114                             }
 115                             else
 116                             {
 117                                 r2v2 = r2;
 118                             }
 119                         }
 120                         if (r2 > r2c2 &&
 121                             (atom[a].q != 0 || atom[a].qB != 0))
 122                         {
 123                             if (r2 > r2c1)
 124                             {
 125                                 r2c2 = r2c1;
 126                                 r2c1 = r2;
 127                             }
 128                             else
 129                             {
 130                                 r2c2 = r2;
 131                             }
 132                         }
 133                     }
 134                 }
 135             }
 136             a_mol += molb->natoms_mol;
 137         }
 138     }
 139
 140     sfree(bLJ);
 141
 142     *rvdw1  = sqrt(r2v1);
 143     *rvdw2  = sqrt(r2v2);
 144     *rcoul1 = sqrt(r2c1);
 145     *rcoul2 = sqrt(r2c2);
 146 }
 147
 148 void calc_cgcm(FILE *fplog,int cg0,int cg1,t_block *cgs,
 149                rvec pos[],rvec cg_cm[])
 150 {
 151     int  icg,k,k0,k1,d;
 152     rvec cg;
 153     real nrcg,inv_ncg;
 154     atom_id *cgindex;
 155
 156 #ifdef DEBUG
 157     fprintf(fplog,"Calculating centre of geometry for charge groups %d to %d\n",
 158             cg0,cg1);
 159 #endif
 160     cgindex = cgs->index;
 161
 162     /* Compute the center of geometry for all charge groups */
 163     for(icg=cg0; (icg<cg1); icg++) {
 164         k0      = cgindex[icg];
 165         k1      = cgindex[icg+1];
 166         nrcg    = k1-k0;
 167         if (nrcg == 1) {
 168             copy_rvec(pos[k0],cg_cm[icg]);
 169         }
 170         else {
 171             inv_ncg = 1.0/nrcg;
 172
 173             clear_rvec(cg);
 174             for(k=k0; (k<k1); k++)  {
 175                 for(d=0; (d<DIM); d++)
 176                     cg[d] += pos[k][d];
 177             }
 178             for(d=0; (d<DIM); d++)
 179                 cg_cm[icg][d] = inv_ncg*cg[d];
 180         }
 181     }
 182 }
 183
 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[])
 187
 188 {
 189     int  npbcdim,icg,k,k0,k1,d,e;
 190     rvec cg;
 191     real nrcg,inv_ncg;
 192     atom_id *cgindex;
 193     gmx_bool bTric;
 194
 195     if (ePBC == epbcNONE)
 196         gmx_incons("Calling put_charge_groups_in_box for a system without PBC");
 197
 198 #ifdef DEBUG
 199     fprintf(fplog,"Putting cgs %d to %d in box\n",cg0,cg1);
 200 #endif
 201     cgindex = cgs->index;
 202
 203     if (ePBC == epbcXY)
 204         npbcdim = 2;
 205     else
 206         npbcdim = 3;
 207
 208     bTric = TRICLINIC(box);
 209
 210     for(icg=cg0; (icg<cg1); icg++) {
 211         /* First compute the center of geometry for this charge group */
 212         k0      = cgindex[icg];
 213         k1      = cgindex[icg+1];
 214         nrcg    = k1-k0;
 215
 216         if (nrcg == 1) {
 217             copy_rvec(pos[k0],cg_cm[icg]);
 218         } else {
 219             inv_ncg = 1.0/nrcg;
 220
 221             clear_rvec(cg);
 222             for(k=k0; (k<k1); k++)  {
 223                 for(d=0; d<DIM; d++)
 224                     cg[d] += pos[k][d];
 225             }
 226             for(d=0; d<DIM; d++)
 227                 cg_cm[icg][d] = inv_ncg*cg[d];
 228         }
 229         /* Now check pbc for this cg */
 230         if (bTric) {
 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];
 237                     }
 238                 }
 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];
 244                     }
 245                 }
 246             }
 247         } else {
 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];
 253                 }
 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];
 258                 }
 259             }
 260         }
 261 #ifdef DEBUG_PBC
 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]);
 268         }
 269 #endif
 270     }
 271 }