src/tools/calcpot.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 <stddef.h>
  43
  44 #include "types/commrec.h"
  45 #include "vec.h"
  46 #include "calcpot.h"
  47 #include "nrnb.h"
  48 #include "mdebin.h"
  49 #include "mshift.h"
  50 #include "smalloc.h"
  51 #include "force.h"
  52 #include "main.h"
  53 #include "filenm.h"
  54 #include "gmx_fatal.h"
  55 #include "mdrun.h"
  56 #include "ns.h"
  57 #include "txtdump.h"
  58 #include "mdatoms.h"
  59 #include "main.h"
  60 #include "mtop_util.h"
  61 #include "chargegroup.h"
  62
  63 static void c_tabpot(real tabscale,   real VFtab[],
  64                      int  nri,        int  iinr[],
  65                      int  shift[],
  66                      int  jindex[],   int  jjnr[],
  67                      real pos[],
  68                      real facel,      real charge[],
  69                      real pot[],      real shiftvec[])
  70 {
  71   /* Local variables */
  72   const real nul = 0.000000;
  73
  74   /* Table stuff */
  75   const real one = 1.000000;
  76   const real two = 2.000000;
  77   real r1,r1t,fijC,eps,eps2,Y,F,Fp,Geps,Heps2,VV,FF;
  78   int  n0,n1,nnn;
  79
  80   /* General and coulomb stuff */
  81   int  ii,k,n,jnr,ii3,nj0,nj1,is3,j3,ggid;
  82   real fxJ,fyJ,fzJ,fxO,fyO,fzO;
  83   real ixO,iyO,izO,dxO,dyO,dzO;
  84   real txO,tyO,tzO,vcO,fsO,qO,rsqO,rinv1O,rinv2O;
  85   real qqO,qj;
  86   real jx,jy,jz,shX,shY,shZ,poti;
  87
  88   /* Outer loop (over i particles) starts here */
  89   for(n=0; (n<nri); n++) {
  90
  91     /* Unpack shift vector */
  92     is3               = 3*shift[n];
  93     shX               = shiftvec[is3];
  94     shY               = shiftvec[is3+1];
  95     shZ               = shiftvec[is3+2];
  96
  97     /* Unpack I particle */
  98     ii                = iinr[n];
  99     ii3               = 3*ii;
 100
 101     /* Charge of i particle(s) divided by 4 pi eps0 */
 102     qO                = facel*charge[ii];
 103
 104     /* Bounds for the innerloop */
 105     nj0               = jindex[n];
 106     nj1               = jindex[n+1];
 107
 108     /* Compute shifted i position */
 109     ixO               = shX + pos[ii3];
 110     iyO               = shY + pos[ii3+1];
 111     izO               = shZ + pos[ii3+2];
 112     poti              = nul;
 113
 114     /* Inner loop (over j-particles) starts right here */
 115     for(k=nj0; (k<nj1); k++) {
 116
 117       /* Unpack neighbourlist */
 118       jnr               = jjnr[k];
 119       j3                = 3*jnr;
 120       qj                = facel*charge[jnr];
 121       jx                = pos[j3];
 122       jy                = pos[j3+1];
 123       jz                = pos[j3+2];
 124
 125       /* First one is for oxygen, with LJ */
 126       dxO               = ixO - jx;
 127       dyO               = iyO - jy;
 128       dzO               = izO - jz;
 129       rsqO              = dxO*dxO + dyO*dyO + dzO*dzO;
 130
 131       /* Doing fast gmx_invsqrt */
 132       rinv1O            = gmx_invsqrt(rsqO);
 133
 134       /* O block */
 135       r1                = one/rinv1O;
 136       r1t               = r1*tabscale;
 137       n0                = r1t;
 138       n1                = 12*n0;
 139       eps               = r1t-n0;
 140       eps2              = eps*eps;
 141       nnn               = n1;
 142       Y                 = VFtab[nnn];
 143       F                 = VFtab[nnn+1];
 144       Geps              = eps*VFtab[nnn+2];
 145       Heps2             = eps2*VFtab[nnn+3];
 146       Fp                = F+Geps+Heps2;
 147       VV                = Y+eps*Fp;
 148
 149       pot[jnr]         += VV*qO;
 150       poti             += VV*qj;
 151
 152     }
 153     pot[ii] += poti;
 154   }
 155 }
 156
 157 static void low_calc_pot(FILE *log,int nl_type,t_forcerec *fr,
 158                          rvec x[],t_mdatoms *mdatoms,real pot[])
 159 {
 160   t_nblist *nlist;
 161
 162   nlist = &fr->nblists[0].nlist_sr[nl_type];
 163
 164   c_tabpot(fr->nblists[0].table_elec_vdw.scale,fr->nblists[0].table_elec_vdw.data,
 165            nlist->nri,nlist->iinr,
 166            nlist->shift,nlist->jindex,nlist->jjnr,
 167            x[0],fr->epsfac,mdatoms->chargeA,pot,fr->shift_vec[0]);
 168
 169   fprintf(log,"There were %d interactions\n",nlist->nrj);
 170 }
 171
 172 void calc_pot(FILE *logf,t_commrec *cr,
 173               gmx_mtop_t *mtop,
 174               t_inputrec *inputrec,gmx_localtop_t *top,rvec x[],
 175               t_forcerec *fr,gmx_enerdata_t *enerd,
 176               t_mdatoms *mdatoms,real pot[],matrix box,t_graph *graph)
 177 {
 178   static t_nrnb      nrnb;
 179   real        lam[efptNR],dum[efptNR];
 180   rvec        box_size;
 181   int         i,m;
 182
 183   /* Calc the force */
 184   fprintf(stderr,"Doing single force calculation...\n");
 185
 186   if (inputrec->ePBC != epbcNONE)
 187     calc_shifts(box,fr->shift_vec);
 188
 189   put_charge_groups_in_box(logf,0,top->cgs.nr,fr->ePBC,box,&(top->cgs),
 190                            x,fr->cg_cm);
 191   if (graph)
 192     mk_mshift(logf,graph,fr->ePBC,box,x);
 193   /* Do the actual neighbour searching and if twin range electrostatics
 194    * also do the calculation of long range forces and energies.
 195    */
 196   ns(logf,fr,x,box,&mtop->groups,&(inputrec->opts),top,mdatoms,cr,
 197      &nrnb,&lam[0],&dum[0],&enerd->grpp,TRUE,FALSE);
 198   for(m=0; (m<DIM); m++)
 199     box_size[m] = box[m][m];
 200   for(i=0; (i<mdatoms->nr); i++)
 201     pot[i] = 0;
 202   if (debug) {
 203     pr_rvecs(debug,0,"x",x,mdatoms->nr);
 204     pr_rvecs(debug,0,"cgcm",fr->cg_cm,top->cgs.nr);
 205   }
 206   /* electrostatics from any atom to atoms without LJ */
 207   low_calc_pot(logf,eNL_QQ,fr,x,mdatoms,pot);
 208     /* electrostatics from any atom to atoms without charge */
 209   low_calc_pot(logf,eNL_VDW,fr,x,mdatoms,pot);
 210   /* electrostatics from any atom to atoms with LJ */
 211   low_calc_pot(logf,eNL_VDWQQ,fr,x,mdatoms,pot);
 212 }
 213
 214 FILE *init_calcpot(const char *log,const char *tpx,const char *table,
 215                    gmx_mtop_t *mtop,gmx_localtop_t *top,
 216                    t_inputrec *inputrec,t_commrec **cr,
 217                    t_graph **graph,t_mdatoms **mdatoms,
 218                    t_forcerec **fr,
 219                    gmx_enerdata_t *enerd,
 220                    real **pot,
 221                    matrix box,rvec **x, const output_env_t oenv)
 222 {
 223   gmx_localtop_t *ltop;
 224   double   t,t0;
 225   real     lam[efptNR];
 226   int      fep_state;
 227   gmx_bool     bNEMD,bSA;
 228   int      traj=0,xtc_traj=0;
 229   t_state  *state;
 230   rvec     mutot;
 231   t_nrnb   nrnb;
 232   int      m;
 233   rvec     box_size;
 234   tensor   force_vir,shake_vir;
 235   FILE     *fplog;
 236
 237   /* Initiate */
 238   *cr = init_par(NULL,NULL);
 239   gmx_log_open(log,*cr,FALSE,0,&fplog);
 240
 241   init_nrnb(&nrnb);
 242   snew(state,1);
 243   read_tpx_state(tpx,inputrec,state, NULL, mtop);
 244   set_state_entries(state,inputrec,1);
 245   if (fplog)
 246       pr_inputrec(fplog,0,"Input Parameters",inputrec,FALSE);
 247
 248
 249
 250   if (inputrec->efep) {
 251     fprintf(stderr,"WARNING: turning of free energy, will use lambda=0\n");
 252     inputrec->efep = 0;
 253   }
 254
 255   clear_rvec(mutot);
 256   init_md(fplog,*cr,inputrec,oenv,&t,&t0,lam,&fep_state,NULL,
 257           &nrnb,mtop,NULL,-1,NULL,NULL,NULL,
 258           force_vir,shake_vir,mutot,&bSA,NULL,NULL,0);
 259
 260   init_enerdata(mtop->groups.grps[egcENER].nr,0,enerd);
 261
 262   ltop = gmx_mtop_generate_local_top(mtop,inputrec);
 263   *top = *ltop;
 264   sfree(ltop);
 265
 266   *mdatoms = init_mdatoms(fplog,mtop,FALSE);
 267   atoms2md(mtop,inputrec,0,NULL,0,mtop->natoms,*mdatoms);
 268
 269   if (inputrec->ePBC == epbcXYZ) {
 270     /* Calculate intramolecular shift vectors to make molecules whole again */
 271     *graph = mk_graph(fplog,&(top->idef),0,mtop->natoms,FALSE,FALSE);
 272     mk_mshift(fplog,*graph,inputrec->ePBC,state->box,state->x);
 273   } else {
 274     *graph = NULL;
 275   }
 276
 277   /* Turn off twin range if appropriate */
 278   inputrec->rvdw  = inputrec->rcoulomb;
 279   inputrec->rlist = inputrec->rcoulomb;
 280   fprintf(stderr,"Using a coulomb cut-off of %g nm\n",inputrec->rcoulomb);
 281
 282   /* Turn off free energy computation */
 283   inputrec->efep = 0;
 284
 285   /* Set vanderwaals to shift, to force tables */
 286   inputrec->vdwtype     = evdwSHIFT;
 287   inputrec->rvdw_switch = 0.0;
 288   inputrec->eDispCorr   = edispcNO;
 289
 290   /* Initiate forcerecord */
 291   *fr = mk_forcerec();
 292   init_forcerec(fplog,oenv,*fr,NULL,inputrec,mtop,*cr,
 293                 state->box,FALSE,table,NULL,table,NULL,NULL,TRUE,-1);
 294
 295   /* Remove periodicity */
 296   for(m=0; (m<DIM); m++)
 297     box_size[m] = state->box[m][m];
 298   if (inputrec->ePBC != epbcNONE)
 299     do_pbc_first(fplog,state->box,*fr,*graph,state->x);
 300
 301   copy_mat(state->box,box);
 302   *x = state->x;
 303   state->x = NULL;
 304   done_state(state);
 305   sfree(state);
 306
 307   snew(*pot,mtop->natoms);
 308
 309   return fplog;
 310 }