src/gmxlib/chargegroup.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "vec.h"
  46 #include "pbc.h"
  47 #include "smalloc.h"
  48 #include "gmx_fatal.h"
  49 #include "chargegroup.h"
  50
  51
  52 void calc_chargegroup_radii(const gmx_mtop_t *mtop,rvec *x,
  53                             real *rvdw1,real *rvdw2,
  54                             real *rcoul1,real *rcoul2)
  55 {
  56     real r2v1,r2v2,r2c1,r2c2,r2;
  57     int  ntype,i,j,mb,m,cg,a_mol,a0,a1,a;
  58     gmx_bool *bLJ;
  59     gmx_molblock_t *molb;
  60     gmx_moltype_t *molt;
  61     t_block *cgs;
  62     t_atom *atom;
  63     rvec cen;
  64
  65     r2v1 = 0;
  66     r2v2 = 0;
  67     r2c1 = 0;
  68     r2c2 = 0;
  69
  70     ntype = mtop->ffparams.atnr;
  71     snew(bLJ,ntype);
  72     for(i=0; i<ntype; i++)
  73     {
  74         bLJ[i] = FALSE;
  75         for(j=0; j<ntype; j++)
  76         {
  77             if (mtop->ffparams.iparams[i*ntype+j].lj.c6  != 0 ||
  78                 mtop->ffparams.iparams[i*ntype+j].lj.c12 != 0)
  79             {
  80                 bLJ[i] = TRUE;
  81             }
  82         }
  83     }
  84
  85     a_mol = 0;
  86     for(mb=0; mb<mtop->nmolblock; mb++)
  87     {
  88         molb = &mtop->molblock[mb];
  89         molt = &mtop->moltype[molb->type];
  90         cgs  = &molt->cgs;
  91         atom = molt->atoms.atom;
  92         for(m=0; m<molb->nmol; m++)
  93         {
  94             for(cg=0; cg<cgs->nr; cg++)
  95             {
  96                 a0 = cgs->index[cg];
  97                 a1 = cgs->index[cg+1];
  98                 if (a1 - a0 > 1)
  99                 {
 100                     clear_rvec(cen);
 101                     for(a=a0; a<a1; a++)
 102                     {
 103                         rvec_inc(cen,x[a_mol+a]);
 104                     }
 105                     svmul(1.0/(a1-a0),cen,cen);
 106                     for(a=a0; a<a1; a++)
 107                     {
 108                         r2 = distance2(cen,x[a_mol+a]);
 109                         if (r2 > r2v2 && (bLJ[atom[a].type ] ||
 110                                           bLJ[atom[a].typeB]))
 111                         {
 112                             if (r2 > r2v1)
 113                             {
 114                                 r2v2 = r2v1;
 115                                 r2v1 = r2;
 116                             }
 117                             else
 118                             {
 119                                 r2v2 = r2;
 120                             }
 121                         }
 122                         if (r2 > r2c2 &&
 123                             (atom[a].q != 0 || atom[a].qB != 0))
 124                         {
 125                             if (r2 > r2c1)
 126                             {
 127                                 r2c2 = r2c1;
 128                                 r2c1 = r2;
 129                             }
 130                             else
 131                             {
 132                                 r2c2 = r2;
 133                             }
 134                         }
 135                     }
 136                 }
 137             }
 138             a_mol += molb->natoms_mol;
 139         }
 140     }
 141
 142     sfree(bLJ);
 143
 144     *rvdw1  = sqrt(r2v1);
 145     *rvdw2  = sqrt(r2v2);
 146     *rcoul1 = sqrt(r2c1);
 147     *rcoul2 = sqrt(r2c2);
 148 }
 149
 150 void calc_cgcm(FILE *fplog,int cg0,int cg1,t_block *cgs,
 151                rvec pos[],rvec cg_cm[])
 152 {
 153     int  icg,k,k0,k1,d;
 154     rvec cg;
 155     real nrcg,inv_ncg;
 156     atom_id *cgindex;
 157
 158 #ifdef DEBUG
 159     fprintf(fplog,"Calculating centre of geometry for charge groups %d to %d\n",
 160             cg0,cg1);
 161 #endif
 162     cgindex = cgs->index;
 163
 164     /* Compute the center of geometry for all charge groups */
 165     for(icg=cg0; (icg<cg1); icg++) {
 166         k0      = cgindex[icg];
 167         k1      = cgindex[icg+1];
 168         nrcg    = k1-k0;
 169         if (nrcg == 1) {
 170             copy_rvec(pos[k0],cg_cm[icg]);
 171         }
 172         else {
 173             inv_ncg = 1.0/nrcg;
 174
 175             clear_rvec(cg);
 176             for(k=k0; (k<k1); k++)  {
 177                 for(d=0; (d<DIM); d++)
 178                     cg[d] += pos[k][d];
 179             }
 180             for(d=0; (d<DIM); d++)
 181                 cg_cm[icg][d] = inv_ncg*cg[d];
 182         }
 183     }
 184 }
 185
 186 void put_charge_groups_in_box(FILE *fplog,int cg0,int cg1,
 187                               int ePBC,matrix box,t_block *cgs,
 188                               rvec pos[],rvec cg_cm[])
 189
 190 {
 191     int  npbcdim,icg,k,k0,k1,d,e;
 192     rvec cg;
 193     real nrcg,inv_ncg;
 194     atom_id *cgindex;
 195     gmx_bool bTric;
 196
 197     if (ePBC == epbcNONE)
 198         gmx_incons("Calling put_charge_groups_in_box for a system without PBC");
 199
 200 #ifdef DEBUG
 201     fprintf(fplog,"Putting cgs %d to %d in box\n",cg0,cg1);
 202 #endif
 203     cgindex = cgs->index;
 204
 205     if (ePBC == epbcXY)
 206         npbcdim = 2;
 207     else
 208         npbcdim = 3;
 209
 210     bTric = TRICLINIC(box);
 211
 212     for(icg=cg0; (icg<cg1); icg++) {
 213         /* First compute the center of geometry for this charge group */
 214         k0      = cgindex[icg];
 215         k1      = cgindex[icg+1];
 216         nrcg    = k1-k0;
 217
 218         if (nrcg == 1) {
 219             copy_rvec(pos[k0],cg_cm[icg]);
 220         } else {
 221             inv_ncg = 1.0/nrcg;
 222
 223             clear_rvec(cg);
 224             for(k=k0; (k<k1); k++)  {
 225                 for(d=0; d<DIM; d++)
 226                     cg[d] += pos[k][d];
 227             }
 228             for(d=0; d<DIM; d++)
 229                 cg_cm[icg][d] = inv_ncg*cg[d];
 230         }
 231         /* Now check pbc for this cg */
 232         if (bTric) {
 233             for(d=npbcdim-1; d>=0; d--) {
 234                 while(cg_cm[icg][d] < 0) {
 235                     for(e=d; e>=0; e--) {
 236                         cg_cm[icg][e] += box[d][e];
 237                         for(k=k0; (k<k1); k++)
 238                             pos[k][e] += box[d][e];
 239                     }
 240                 }
 241                 while(cg_cm[icg][d] >= box[d][d]) {
 242                     for(e=d; e>=0; e--) {
 243                         cg_cm[icg][e] -= box[d][e];
 244                         for(k=k0; (k<k1); k++)
 245                             pos[k][e] -= box[d][e];
 246                     }
 247                 }
 248             }
 249         } else {
 250             for(d=0; d<npbcdim; d++) {
 251                 while(cg_cm[icg][d] < 0) {
 252                     cg_cm[icg][d] += box[d][d];
 253                     for(k=k0; (k<k1); k++)
 254                         pos[k][d] += box[d][d];
 255                 }
 256                 while(cg_cm[icg][d] >= box[d][d]) {
 257                     cg_cm[icg][d] -= box[d][d];
 258                     for(k=k0; (k<k1); k++)
 259                         pos[k][d] -= box[d][d];
 260                 }
 261             }
 262         }
 263 #ifdef DEBUG_PBC
 264         for(d=0; (d<npbcdim); d++) {
 265             if ((cg_cm[icg][d] < 0) || (cg_cm[icg][d] >= box[d][d]))
 266                 gmx_fatal(FARGS,"cg_cm[%d] = %15f  %15f  %15f\n"
 267                           "box = %15f  %15f  %15f\n",
 268                           icg,cg_cm[icg][XX],cg_cm[icg][YY],cg_cm[icg][ZZ],
 269                           box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
 270         }
 271 #endif
 272     }
 273 }