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