src/gromacs/mdlib/rf_util.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROwing Monsters And Cloning Shrimps
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "sysstuff.h"
  40 #include "typedefs.h"
  41 #include "force.h"
  42 #include "names.h"
  43 #include "vec.h"
  44 #include "physics.h"
  45 #include "copyrite.h"
  46 #include "pbc.h"
  47
  48 real RF_excl_correction(const t_forcerec *fr, t_graph *g,
  49                         const t_mdatoms *mdatoms, const t_blocka *excl,
  50                         rvec x[], rvec f[], rvec *fshift, const t_pbc *pbc,
  51                         real lambda, real *dvdlambda)
  52 {
  53     /* Calculate the reaction-field energy correction for this node:
  54      * epsfac q_i q_j (k_rf r_ij^2 - c_rf)
  55      * and force correction for all excluded pairs, including self pairs.
  56      */
  57     int         top, i, j, j1, j2, k, ki;
  58     double      q2sumA, q2sumB, ener;
  59     const real *chargeA, *chargeB;
  60     real        ek, ec, L1, qiA, qiB, qqA, qqB, qqL, v;
  61     rvec        dx, df;
  62     atom_id    *AA;
  63     ivec        dt;
  64     int         start = mdatoms->start;
  65     int         end   = mdatoms->homenr+start;
  66     int         niat;
  67     gmx_bool    bMolPBC = fr->bMolPBC;
  68
  69     if (fr->n_tpi)
  70     {
  71         /* For test particle insertion we only correct for the test molecule */
  72         start = mdatoms->nr - fr->n_tpi;
  73     }
  74
  75     ek      = fr->epsfac*fr->k_rf;
  76     ec      = fr->epsfac*fr->c_rf;
  77     chargeA = mdatoms->chargeA;
  78     chargeB = mdatoms->chargeB;
  79     AA      = excl->a;
  80     ki      = CENTRAL;
  81
  82     if (fr->bDomDec)
  83     {
  84         niat = excl->nr;
  85     }
  86     else
  87     {
  88         niat = end;
  89     }
  90
  91     q2sumA = 0;
  92     q2sumB = 0;
  93     ener   = 0;
  94     if (mdatoms->nChargePerturbed == 0)
  95     {
  96         for (i = start; i < niat; i++)
  97         {
  98             qiA = chargeA[i];
  99             if (i < end)
 100             {
 101                 q2sumA += qiA*qiA;
 102             }
 103             /* Do the exclusions */
 104             j1  = excl->index[i];
 105             j2  = excl->index[i+1];
 106             for (j = j1; j < j2; j++)
 107             {
 108                 k = AA[j];
 109                 if (k > i)
 110                 {
 111                     qqA = qiA*chargeA[k];
 112                     if (qqA != 0)
 113                     {
 114                         if (g)
 115                         {
 116                             rvec_sub(x[i], x[k], dx);
 117                             ivec_sub(SHIFT_IVEC(g, i), SHIFT_IVEC(g, k), dt);
 118                             ki = IVEC2IS(dt);
 119                         }
 120                         else if (bMolPBC)
 121                         {
 122                             ki = pbc_dx_aiuc(pbc, x[i], x[k], dx);
 123                         }
 124                         else
 125                         {
 126                             rvec_sub(x[i], x[k], dx);
 127                         }
 128                         ener += qqA*(ek*norm2(dx) - ec);
 129                         svmul(-2*qqA*ek, dx, df);
 130                         rvec_inc(f[i], df);
 131                         rvec_dec(f[k], df);
 132                         rvec_inc(fshift[ki], df);
 133                         rvec_dec(fshift[CENTRAL], df);
 134                     }
 135                 }
 136             }
 137         }
 138         ener += -0.5*ec*q2sumA;
 139     }
 140     else
 141     {
 142         L1 = 1.0 - lambda;
 143         for (i = start; i < niat; i++)
 144         {
 145             qiA = chargeA[i];
 146             qiB = chargeB[i];
 147             if (i < end)
 148             {
 149                 q2sumA += qiA*qiA;
 150                 q2sumB += qiB*qiB;
 151             }
 152             /* Do the exclusions */
 153             j1  = excl->index[i];
 154             j2  = excl->index[i+1];
 155             for (j = j1; j < j2; j++)
 156             {
 157                 k = AA[j];
 158                 if (k > i)
 159                 {
 160                     qqA = qiA*chargeA[k];
 161                     qqB = qiB*chargeB[k];
 162                     if (qqA != 0 || qqB != 0)
 163                     {
 164                         qqL = L1*qqA + lambda*qqB;
 165                         if (g)
 166                         {
 167                             rvec_sub(x[i], x[k], dx);
 168                             ivec_sub(SHIFT_IVEC(g, i), SHIFT_IVEC(g, k), dt);
 169                             ki = IVEC2IS(dt);
 170                         }
 171                         else if (bMolPBC)
 172                         {
 173                             ki = pbc_dx_aiuc(pbc, x[i], x[k], dx);
 174                         }
 175                         else
 176                         {
 177                             rvec_sub(x[i], x[k], dx);
 178                         }
 179                         v     = ek*norm2(dx) - ec;
 180                         ener += qqL*v;
 181                         svmul(-2*qqL*ek, dx, df);
 182                         rvec_inc(f[i], df);
 183                         rvec_dec(f[k], df);
 184                         rvec_inc(fshift[ki], df);
 185                         rvec_dec(fshift[CENTRAL], df);
 186                         *dvdlambda += (qqB - qqA)*v;
 187                     }
 188                 }
 189             }
 190         }
 191         ener       += -0.5*ec*(L1*q2sumA + lambda*q2sumB);
 192         *dvdlambda += -0.5*ec*(q2sumB - q2sumA);
 193     }
 194
 195     if (debug)
 196     {
 197         fprintf(debug, "RF exclusion energy: %g\n", ener);
 198     }
 199
 200     return ener;
 201 }
 202
 203 void calc_rffac(FILE *fplog, int eel, real eps_r, real eps_rf, real Rc, real Temp,
 204                 real zsq, matrix box,
 205                 real *kappa, real *krf, real *crf)
 206 {
 207     /* Compute constants for Generalized reaction field */
 208     real   k1, k2, I, vol, rmin;
 209
 210     if (EEL_RF(eel))
 211     {
 212         vol     = det(box);
 213         if (eel == eelGRF)
 214         {
 215             /* Consistency check */
 216             if (Temp <= 0.0)
 217             {
 218                 gmx_fatal(FARGS, "Temperature is %f while using"
 219                           " Generalized Reaction Field\n", Temp);
 220             }
 221             /* Ionic strength (only needed for eelGRF */
 222             I       = 0.5*zsq/vol;
 223             *kappa  = sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));
 224         }
 225         else
 226         {
 227             I      = 0;
 228             *kappa = 0;
 229         }
 230
 231         /* eps == 0 signals infinite dielectric */
 232         if (eps_rf == 0)
 233         {
 234             *krf = 1/(2*Rc*Rc*Rc);
 235         }
 236         else
 237         {
 238             k1   = 1 + *kappa*Rc;
 239             k2   = eps_rf*sqr((real)(*kappa*Rc));
 240
 241             *krf = ((eps_rf - eps_r)*k1 + 0.5*k2)/((2*eps_rf + eps_r)*k1 + k2)/(Rc*Rc*Rc);
 242         }
 243         *crf   = 1/Rc + *krf*Rc*Rc;
 244         rmin   = pow(*krf*2.0, -1.0/3.0);
 245
 246         if (fplog)
 247         {
 248             if (eel == eelGRF)
 249             {
 250                 please_cite(fplog, "Tironi95a");
 251                 fprintf(fplog, "%s:\n"
 252                         "epsRF = %10g, I   = %10g, volume = %10g, kappa  = %10g\n"
 253                         "rc    = %10g, krf = %10g, crf    = %10g, epsfac = %10g\n",
 254                         eel_names[eel], eps_rf, I, vol, *kappa, Rc, *krf, *crf,
 255                         ONE_4PI_EPS0/eps_r);
 256             }
 257             else
 258             {
 259                 fprintf(fplog, "%s:\n"
 260                         "epsRF = %g, rc = %g, krf = %g, crf = %g, epsfac = %g\n",
 261                         eel_names[eel], eps_rf, Rc, *krf, *crf, ONE_4PI_EPS0/eps_r);
 262             }
 263             fprintf(fplog,
 264                     "The electrostatics potential has its minimum at r = %g\n",
 265                     rmin);
 266         }
 267     }
 268 }
 269
 270 void init_generalized_rf(FILE *fplog,
 271                          const gmx_mtop_t *mtop, const t_inputrec *ir,
 272                          t_forcerec *fr)
 273 {
 274     int                  mb, i, j;
 275     real                 q, zsq, nrdf, T;
 276     const gmx_moltype_t *molt;
 277     const t_block       *cgs;
 278
 279     if (ir->efep != efepNO && fplog)
 280     {
 281         fprintf(fplog, "\nWARNING: the generalized reaction field constants are determined from topology A only\n\n");
 282     }
 283     zsq = 0.0;
 284     for (mb = 0; mb < mtop->nmolblock; mb++)
 285     {
 286         molt = &mtop->moltype[mtop->molblock[mb].type];
 287         cgs  = &molt->cgs;
 288         for (i = 0; (i < cgs->nr); i++)
 289         {
 290             q = 0;
 291             for (j = cgs->index[i]; (j < cgs->index[i+1]); j++)
 292             {
 293                 q += molt->atoms.atom[j].q;
 294             }
 295             zsq += mtop->molblock[mb].nmol*q*q;
 296         }
 297         fr->zsquare = zsq;
 298     }
 299
 300     T    = 0.0;
 301     nrdf = 0.0;
 302     for (i = 0; (i < ir->opts.ngtc); i++)
 303     {
 304         nrdf += ir->opts.nrdf[i];
 305         T    += (ir->opts.nrdf[i] * ir->opts.ref_t[i]);
 306     }
 307     if (nrdf == 0)
 308     {
 309         gmx_fatal(FARGS, "No degrees of freedom!");
 310     }
 311     fr->temp   = T/nrdf;
 312 }