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