src/mdlib/rf_util.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, 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 "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "force.h"
  45 #include "names.h"
  46 #include "vec.h"
  47 #include "physics.h"
  48 #include "copyrite.h"
  49 #include "pbc.h"
  50
  51 real RF_excl_correction(FILE *log,
  52                         const t_forcerec *fr,t_graph *g,
  53                         const t_mdatoms *mdatoms,const t_blocka *excl,
  54                         rvec x[],rvec f[],rvec *fshift,const t_pbc *pbc,
  55                         real lambda,real *dvdlambda)
  56 {
  57   /* Calculate the reaction-field energy correction for this node:
  58    * epsfac q_i q_j (k_rf r_ij^2 - c_rf)
  59    * and force correction for all excluded pairs, including self pairs.
  60    */
  61   int    top,i,j,j1,j2,k,ki;
  62   double q2sumA,q2sumB,ener;
  63   const real *chargeA,*chargeB;
  64   real   ek,ec,L1,qiA,qiB,qqA,qqB,qqL,v;
  65   rvec   dx,df;
  66   atom_id *AA;
  67   ivec   dt;
  68   int    start = mdatoms->start;
  69   int    end   = mdatoms->homenr+start;
  70   int    niat;
  71   gmx_bool   bMolPBC = fr->bMolPBC;
  72
  73   if (fr->n_tpi)
  74     /* For test particle insertion we only correct for the test molecule */
  75     start = mdatoms->nr - fr->n_tpi;
  76
  77   ek = fr->epsfac*fr->k_rf;
  78   ec = fr->epsfac*fr->c_rf;
  79   chargeA = mdatoms->chargeA;
  80   chargeB = mdatoms->chargeB;
  81   AA = excl->a;
  82   ki = CENTRAL;
  83
  84   if (fr->bDomDec)
  85     niat = excl->nr;
  86   else
  87     niat = end;
  88
  89   q2sumA = 0;
  90   q2sumB = 0;
  91   ener = 0;
  92   if (mdatoms->nChargePerturbed == 0) {
  93     for(i=start; i<niat; i++) {
  94       qiA = chargeA[i];
  95       if (i < end)
  96         q2sumA += qiA*qiA;
  97       /* Do the exclusions */
  98       j1  = excl->index[i];
  99       j2  = excl->index[i+1];
 100       for(j=j1; j<j2; j++) {
 101         k = AA[j];
 102         if (k > i) {
 103           qqA = qiA*chargeA[k];
 104           if (qqA != 0) {
 105             if (g) {
 106               rvec_sub(x[i],x[k],dx);
 107               ivec_sub(SHIFT_IVEC(g,i),SHIFT_IVEC(g,k),dt);
 108               ki=IVEC2IS(dt);
 109             } else if (bMolPBC) {
 110               ki = pbc_dx_aiuc(pbc,x[i],x[k],dx);
 111             } else
 112               rvec_sub(x[i],x[k],dx);
 113             ener += qqA*(ek*norm2(dx) - ec);
 114             svmul(-2*qqA*ek,dx,df);
 115             rvec_inc(f[i],df);
 116             rvec_dec(f[k],df);
 117             rvec_inc(fshift[ki],df);
 118             rvec_dec(fshift[CENTRAL],df);
 119           }
 120         }
 121       }
 122     }
 123     ener += -0.5*ec*q2sumA;
 124   } else {
 125     L1 = 1.0 - lambda;
 126     for(i=start; i<niat; i++) {
 127       qiA = chargeA[i];
 128       qiB = chargeB[i];
 129       if (i < end) {
 130         q2sumA += qiA*qiA;
 131         q2sumB += qiB*qiB;
 132       }
 133       /* Do the exclusions */
 134       j1  = excl->index[i];
 135       j2  = excl->index[i+1];
 136       for(j=j1; j<j2; j++) {
 137         k = AA[j];
 138         if (k > i) {
 139           qqA = qiA*chargeA[k];
 140           qqB = qiB*chargeB[k];
 141           if (qqA != 0 || qqB != 0) {
 142             qqL = L1*qqA + lambda*qqB;
 143             if (g) {
 144               rvec_sub(x[i],x[k],dx);
 145               ivec_sub(SHIFT_IVEC(g,i),SHIFT_IVEC(g,k),dt);
 146               ki=IVEC2IS(dt);
 147             } else if (bMolPBC) {
 148               ki = pbc_dx_aiuc(pbc,x[i],x[k],dx);
 149             } else
 150               rvec_sub(x[i],x[k],dx);
 151             v = ek*norm2(dx) - ec;
 152             ener += qqL*v;
 153             svmul(-2*qqL*ek,dx,df);
 154             rvec_inc(f[i],df);
 155             rvec_dec(f[k],df);
 156             rvec_inc(fshift[ki],df);
 157             rvec_dec(fshift[CENTRAL],df);
 158             *dvdlambda += (qqB - qqA)*v;
 159           }
 160         }
 161       }
 162     }
 163     ener += -0.5*ec*(L1*q2sumA + lambda*q2sumB);
 164     *dvdlambda += -0.5*ec*(q2sumB - q2sumA);
 165   }
 166
 167   if (debug)
 168     fprintf(debug,"RF exclusion energy: %g\n",ener);
 169
 170   return ener;
 171 }
 172
 173 void calc_rffac(FILE *fplog,int eel,real eps_r,real eps_rf,real Rc,real Temp,
 174                 real zsq,matrix box,
 175                 real *kappa,real *krf,real *crf)
 176 {
 177   /* Compute constants for Generalized reaction field */
 178   real   k1,k2,I,vol,rmin;
 179
 180   if (EEL_RF(eel)) {
 181     vol     = det(box);
 182     if (eel == eelGRF) {
 183       /* Consistency check */
 184       if (Temp <= 0.0)
 185         gmx_fatal(FARGS,"Temperature is %f while using"
 186                     " Generalized Reaction Field\n",Temp);
 187       /* Ionic strength (only needed for eelGRF */
 188       I       = 0.5*zsq/vol;
 189       *kappa  = sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));
 190     }
 191     else {
 192       I      = 0;
 193       *kappa = 0;
 194     }
 195
 196     /* eps == 0 signals infinite dielectric */
 197     if (eps_rf == 0) {
 198       *krf = 1/(2*Rc*Rc*Rc);
 199     } else {
 200       k1   = 1 + *kappa*Rc;
 201       k2   = eps_rf*sqr((real)(*kappa*Rc));
 202
 203       *krf = ((eps_rf - eps_r)*k1 + 0.5*k2)/((2*eps_rf + eps_r)*k1 + k2)/(Rc*Rc*Rc);
 204     }
 205     *crf   = 1/Rc + *krf*Rc*Rc;
 206     rmin   = pow(*krf*2.0,-1.0/3.0);
 207
 208     if (fplog) {
 209       if (eel == eelGRF) {
 210         please_cite(fplog,"Tironi95a");
 211         fprintf(fplog,"%s:\n"
 212                 "epsRF = %10g, I   = %10g, volume = %10g, kappa  = %10g\n"
 213                 "rc    = %10g, krf = %10g, crf    = %10g, epsfac = %10g\n",
 214                 eel_names[eel],eps_rf,I,vol,*kappa,Rc,*krf,*crf,
 215                 ONE_4PI_EPS0/eps_r);
 216       } else {
 217         fprintf(fplog,"%s:\n"
 218                 "epsRF = %g, rc = %g, krf = %g, crf = %g, epsfac = %g\n",
 219                 eel_names[eel],eps_rf,Rc,*krf,*crf,ONE_4PI_EPS0/eps_r);
 220       }
 221       fprintf(fplog,
 222               "The electrostatics potential has its minimum at r = %g\n",
 223               rmin);
 224     }
 225   }
 226 }
 227
 228 void init_generalized_rf(FILE *fplog,
 229                          const gmx_mtop_t *mtop,const t_inputrec *ir,
 230                          t_forcerec *fr)
 231 {
 232   int  mb,i,j;
 233   real q,zsq,nrdf,T;
 234   const gmx_moltype_t *molt;
 235   const t_block *cgs;
 236
 237   if (ir->efep != efepNO && fplog) {
 238     fprintf(fplog,"\nWARNING: the generalized reaction field constants are determined from topology A only\n\n");
 239   }
 240   zsq = 0.0;
 241   for(mb=0; mb<mtop->nmolblock; mb++) {
 242     molt = &mtop->moltype[mtop->molblock[mb].type];
 243     cgs  = &molt->cgs;
 244     for(i=0; (i<cgs->nr); i++) {
 245       q = 0;
 246       for(j=cgs->index[i]; (j<cgs->index[i+1]); j++) {
 247         q += molt->atoms.atom[j].q;
 248       }
 249       zsq += mtop->molblock[mb].nmol*q*q;
 250     }
 251     fr->zsquare = zsq;
 252   }
 253
 254   T    = 0.0;
 255   nrdf = 0.0;
 256   for(i=0; (i<ir->opts.ngtc); i++) {
 257     nrdf += ir->opts.nrdf[i];
 258     T    += (ir->opts.nrdf[i] * ir->opts.ref_t[i]);
 259   }
 260   if (nrdf == 0) {
 261     gmx_fatal(FARGS,"No degrees of freedom!");
 262   }
 263   fr->temp   = T/nrdf;
 264 }