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  mb,m,cg,a_mol,a0,a1,a;
  56     t_iparams *ip;
  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     ip = mtop->ffparams.iparams;
  69
  70     a_mol = 0;
  71     for(mb=0; mb<mtop->nmolblock; mb++)
  72     {
  73         molb = &mtop->molblock[mb];
  74         molt = &mtop->moltype[molb->type];
  75         cgs  = &molt->cgs;
  76         atom = molt->atoms.atom;
  77         for(m=0; m<molb->nmol; m++)
  78         {
  79             for(cg=0; cg<cgs->nr; cg++)
  80             {
  81                 a0 = cgs->index[cg];
  82                 a1 = cgs->index[cg+1];
  83                 if (a1 - a0 > 1)
  84                 {
  85                     clear_rvec(cen);
  86                     for(a=a0; a<a1; a++)
  87                     {
  88                         rvec_inc(cen,x[a_mol+a]);
  89                     }
  90                     svmul(1.0/(a1-a0),cen,cen);
  91                     for(a=a0; a<a1; a++)
  92                     {
  93                         r2 = distance2(cen,x[a_mol+a]);
  94                         if (r2 > r2v2 &&
  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))
  99                         {
 100                             if (r2 > r2v1)
 101                             {
 102                                 r2v2 = r2v1;
 103                                 r2v1 = r2;
 104                             }
 105                             else
 106                             {
 107                                 r2v2 = r2;
 108                             }
 109                         }
 110                         if (r2 > r2c2 &&
 111                             (atom[a].q != 0 || atom[a].qB != 0))
 112                         {
 113                             if (r2 > r2c1)
 114                             {
 115                                 r2c2 = r2c1;
 116                                 r2c1 = r2;
 117                             }
 118                             else
 119                             {
 120                                 r2c2 = r2;
 121                             }
 122                         }
 123                     }
 124                 }
 125             }
 126             a_mol += molb->natoms_mol;
 127         }
 128     }
 129
 130     *rvdw1  = sqrt(r2v1);
 131     *rvdw2  = sqrt(r2v2);
 132     *rcoul1 = sqrt(r2c1);
 133     *rcoul2 = sqrt(r2c2);
 134 }
 135
 136 void calc_cgcm(FILE *fplog,int cg0,int cg1,t_block *cgs,
 137                rvec pos[],rvec cg_cm[])
 138 {
 139     int  icg,k,k0,k1,d;
 140     rvec cg;
 141     real nrcg,inv_ncg;
 142     atom_id *cgindex;
 143
 144 #ifdef DEBUG
 145     fprintf(fplog,"Calculating centre of geometry for charge groups %d to %d\n",
 146             cg0,cg1);
 147 #endif
 148     cgindex = cgs->index;
 149
 150     /* Compute the center of geometry for all charge groups */
 151     for(icg=cg0; (icg<cg1); icg++) {
 152         k0      = cgindex[icg];
 153         k1      = cgindex[icg+1];
 154         nrcg    = k1-k0;
 155         if (nrcg == 1) {
 156             copy_rvec(pos[k0],cg_cm[icg]);
 157         }
 158         else {
 159             inv_ncg = 1.0/nrcg;
 160
 161             clear_rvec(cg);
 162             for(k=k0; (k<k1); k++)  {
 163                 for(d=0; (d<DIM); d++)
 164                     cg[d] += pos[k][d];
 165             }
 166             for(d=0; (d<DIM); d++)
 167                 cg_cm[icg][d] = inv_ncg*cg[d];
 168         }
 169     }
 170 }
 171
 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[])
 175
 176 {
 177     int  npbcdim,icg,k,k0,k1,d,e;
 178     rvec cg;
 179     real nrcg,inv_ncg;
 180     atom_id *cgindex;
 181     bool bTric;
 182
 183     if (ePBC == epbcNONE)
 184         gmx_incons("Calling put_charge_groups_in_box for a system without PBC");
 185
 186 #ifdef DEBUG
 187     fprintf(fplog,"Putting cgs %d to %d in box\n",cg0,cg1);
 188 #endif
 189     cgindex = cgs->index;
 190
 191     if (ePBC == epbcXY)
 192         npbcdim = 2;
 193     else
 194         npbcdim = 3;
 195
 196     bTric = TRICLINIC(box);
 197
 198     for(icg=cg0; (icg<cg1); icg++) {
 199         /* First compute the center of geometry for this charge group */
 200         k0      = cgindex[icg];
 201         k1      = cgindex[icg+1];
 202         nrcg    = k1-k0;
 203
 204         if (nrcg == 1) {
 205             copy_rvec(pos[k0],cg_cm[icg]);
 206         } else {
 207             inv_ncg = 1.0/nrcg;
 208
 209             clear_rvec(cg);
 210             for(k=k0; (k<k1); k++)  {
 211                 for(d=0; d<DIM; d++)
 212                     cg[d] += pos[k][d];
 213             }
 214             for(d=0; d<DIM; d++)
 215                 cg_cm[icg][d] = inv_ncg*cg[d];
 216         }
 217         /* Now check pbc for this cg */
 218         if (bTric) {
 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];
 225                     }
 226                 }
 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];
 232                     }
 233                 }
 234             }
 235         } else {
 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];
 241                 }
 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];
 246                 }
 247             }
 248         }
 249 #ifdef DEBUG_PBC
 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]);
 256         }
 257 #endif
 258     }
 259 }