src/tools/calcpot.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include "vec.h"
  34 #include "calcpot.h"
  35 #include "nrnb.h"
  36 #include "mdebin.h"
  37 #include "mshift.h"
  38 #include "smalloc.h"
  39 #include "force.h"
  40 #include "main.h"
  41 #include "filenm.h"
  42 #include "fatal.h"
  43 #include "mdrun.h"
  44 #include "ns.h"
  45 #include "txtdump.h"
  46
  47 static void c_tabpot(real tabscale,   real VFtab[],
  48                      int  nri,        int  iinr[],
  49                      int  shift[],
  50                      int  jindex[],   int  jjnr[],
  51                      real pos[],
  52                      real facel,      real charge[],
  53                      real pot[],      real shiftvec[])
  54 {
  55   /* Local variables */
  56   const real nul = 0.000000;
  57
  58   /* Table stuff */
  59   const real one = 1.000000;
  60   const real two = 2.000000;
  61   real r1,r1t,fijC,eps,eps2,Y,F,Fp,Geps,Heps2,VV,FF;
  62   int  n0,n1,nnn;
  63
  64   /* General and coulomb stuff */
  65   int  ii,k,n,jnr,ii3,nj0,nj1,is3,j3,ggid;
  66   real fxJ,fyJ,fzJ,fxO,fyO,fzO;
  67   real ixO,iyO,izO,dxO,dyO,dzO;
  68   real txO,tyO,tzO,vcO,fsO,qO,rsqO,rinv1O,rinv2O;
  69   real qqO,qj;
  70   real jx,jy,jz,shX,shY,shZ,poti;
  71
  72   /* Outer loop (over i particles) starts here */
  73   for(n=0; (n<nri); n++) {
  74
  75     /* Unpack shift vector */
  76     is3               = 3*shift[n];
  77     shX               = shiftvec[is3];
  78     shY               = shiftvec[is3+1];
  79     shZ               = shiftvec[is3+2];
  80
  81     /* Unpack I particle */
  82     ii                = iinr[n];
  83     ii3               = 3*ii;
  84
  85     /* Charge of i particle(s) divided by 4 pi eps0 */
  86     qO                = facel*charge[ii];
  87
  88     /* Bounds for the innerloop */
  89     nj0               = jindex[n];
  90     nj1               = jindex[n+1];
  91
  92     /* Compute shifted i position */
  93     ixO               = shX + pos[ii3];
  94     iyO               = shY + pos[ii3+1];
  95     izO               = shZ + pos[ii3+2];
  96     poti              = nul;
  97
  98     /* Inner loop (over j-particles) starts right here */
  99     for(k=nj0; (k<nj1); k++) {
 100
 101       /* Unpack neighbourlist */
 102       jnr               = jjnr[k];
 103       j3                = 3*jnr;
 104       qj                = facel*charge[jnr];
 105       jx                = pos[j3];
 106       jy                = pos[j3+1];
 107       jz                = pos[j3+2];
 108
 109       /* First one is for oxygen, with LJ */
 110       dxO               = ixO - jx;
 111       dyO               = iyO - jy;
 112       dzO               = izO - jz;
 113       rsqO              = dxO*dxO + dyO*dyO + dzO*dzO;
 114
 115       /* Doing fast invsqrt */
 116       rinv1O            = invsqrt(rsqO);
 117
 118       /* O block */
 119       r1                = one/rinv1O;
 120       r1t               = r1*tabscale;
 121       n0                = r1t;
 122       n1                = 12*n0;
 123       eps               = r1t-n0;
 124       eps2              = eps*eps;
 125       nnn               = n1;
 126       Y                 = VFtab[nnn];
 127       F                 = VFtab[nnn+1];
 128       Geps              = eps*VFtab[nnn+2];
 129       Heps2             = eps2*VFtab[nnn+3];
 130       Fp                = F+Geps+Heps2;
 131       VV                = Y+eps*Fp;
 132
 133       pot[jnr]         += VV*qO;
 134       poti             += VV*qj;
 135
 136     }
 137     pot[ii] += poti;
 138   }
 139 }
 140
 141 static void low_calc_pot(FILE *log,int nl_type,t_forcerec *fr,
 142                          rvec x[],t_mdatoms *mdatoms,real pot[])
 143 {
 144   t_nblist *nlist;
 145
 146   nlist = &fr->nlist_sr[nl_type];
 147
 148   c_tabpot(fr->tabscale,fr->coulvdwtab,nlist->nri,nlist->iinr,
 149            nlist->shift,nlist->jindex,nlist->jjnr,
 150            x[0],fr->epsfac,mdatoms->chargeA,pot,fr->shift_vec[0]);
 151
 152   fprintf(log,"There were %d interactions\n",nlist->nrj);
 153 }
 154
 155 void calc_pot(FILE *logf,t_nsborder *nsb,t_commrec *cr,t_groups *grps,
 156               t_parm *parm,t_topology *top,rvec x[],t_forcerec *fr,
 157               t_mdatoms *mdatoms,real pot[])
 158 {
 159   static bool        bFirst=TRUE;
 160   static t_nrnb      nrnb;
 161   static rvec        *f;
 162   real        lam=0,dum=0;
 163   rvec        box_size;
 164   int         i,m;
 165
 166   /* Calc the force */
 167   fprintf(stderr,"Doing single force calculation...\n");
 168
 169   if (bFirst) {
 170     snew(f,   nsb->natoms);
 171
 172     bFirst = FALSE;
 173   }
 174   /* Reset long range forces if necessary */
 175   if (fr->bTwinRange) {
 176     clear_rvecs(nsb->natoms,fr->f_twin);
 177     clear_rvecs(SHIFTS,fr->fshift_twin);
 178   }
 179   if (parm->ir.epc != epcNO)
 180     calc_shifts(parm->box,box_size,fr->shift_vec);
 181   put_charge_groups_in_box(stdlog,0,top->blocks[ebCGS].nr,
 182                            parm->box,box_size,&(top->blocks[ebCGS]),x,
 183                            fr->cg_cm);
 184   /* mk_mshift(stdlog,graph,parm->box,x);*/
 185   /* Do the actual neighbour searching and if twin range electrostatics
 186    * also do the calculation of long range forces and energies.
 187    */
 188
 189   ns(logf,fr,x,f,parm->box,grps,&(parm->ir.opts),top,mdatoms,cr,
 190      &nrnb,nsb,0,lam,&dum);
 191   for(m=0; (m<DIM); m++)
 192     box_size[m] = parm->box[m][m];
 193   for(i=0; (i<mdatoms->nr); i++)
 194     pot[i] = 0;
 195   if (debug) {
 196     pr_rvecs(debug,0,"x",x,mdatoms->nr);
 197     pr_rvecs(debug,0,"cgcm",fr->cg_cm,top->blocks[ebCGS].nr);
 198   }
 199   /* electrostatics from any atom to atoms without LJ */
 200   low_calc_pot(logf,eNL_QQ,fr,x,mdatoms,pot);
 201     /* electrostatics from any atom to atoms without charge */
 202   low_calc_pot(logf,eNL_VDW,fr,x,mdatoms,pot);
 203   /* electrostatics from any atom to atoms with LJ */
 204   low_calc_pot(logf,eNL_VDWQQ,fr,x,mdatoms,pot);
 205 }
 206
 207 void init_calcpot(int nfile,t_filenm fnm[],t_topology *top,
 208                   rvec **x,t_parm *parm,t_commrec *cr,
 209                   t_graph **graph,t_mdatoms **mdatoms,
 210                   t_nsborder *nsb,t_groups *grps,
 211                   t_forcerec **fr,real **pot,
 212                   matrix box)
 213 {
 214   real     t,t0,lam,lam0;
 215   bool     bTYZ,bNEMD,bSA;
 216   char     *traj,*xtc_traj;
 217   rvec     *v,mutot;
 218   t_nrnb   nrnb;
 219   t_mdebin *mdebin;
 220   t_vcm    *vcm=NULL;
 221   int      fp_ene,m;
 222   rvec     box_size;
 223   tensor   force_vir,pme_vir,shake_vir;
 224
 225   /* Initiate */
 226   cr->nnodes = 1; cr->nodeid    = 0; cr->left   = 0; cr->right  = 1;
 227   cr->nthreads = 1 ; cr->threadid = 0;
 228   open_log(ftp2fn(efLOG,nfile,fnm),cr);
 229
 230   if (parm->ir.efep) {
 231     fprintf(stderr,"WARNING: turning of free energy, will use lambda=0\n");
 232     parm->ir.efep = 0;
 233   }
 234
 235   init_nrnb(&nrnb);
 236   init_single(stdlog,parm,ftp2fn(efTPX,nfile,fnm),top,x,&v,mdatoms,nsb);
 237   init_md(cr,&(parm->ir),parm->box,&t,&t0,&lam,&lam0,
 238           &nrnb,&bTYZ,top,-1,NULL,&traj,&xtc_traj,&fp_ene,NULL,
 239           &mdebin,grps,force_vir,pme_vir,
 240           shake_vir,*mdatoms,mutot,&bNEMD,&bSA,&vcm,nsb);
 241   init_groups(stdlog,*mdatoms,&(parm->ir.opts),grps);
 242
 243   /* Calculate intramolecular shift vectors to make molecules whole again */
 244   *graph = mk_graph(&(top->idef),top->atoms.nr,FALSE,FALSE);
 245   mk_mshift(stdlog,*graph,parm->box,*x);
 246
 247   /* Turn off twin range if appropriate */
 248   parm->ir.rvdw  = parm->ir.rcoulomb;
 249   parm->ir.rlist = parm->ir.rcoulomb;
 250   fprintf(stderr,"Using a coulomb cut-off of %g nm\n",parm->ir.rcoulomb);
 251
 252   /* Turn off free energy computation */
 253   parm->ir.efep = 0;
 254
 255   /* Set vanderwaals to shift, to force tables */
 256   parm->ir.vdwtype     = evdwSHIFT;
 257   parm->ir.rvdw_switch = 0.0;
 258
 259   /* Initiate forcerecord */
 260   *fr = mk_forcerec();
 261   init_forcerec(stdlog,*fr,&(parm->ir),top,cr,*mdatoms,
 262                 nsb,parm->box,FALSE,NULL,TRUE);
 263
 264   /* Remove periodicity */
 265   for(m=0; (m<DIM); m++)
 266     box_size[m] = parm->box[m][m];
 267   if (parm->ir.ePBC != epbcNONE)
 268     do_pbc_first(stdlog,parm,box_size,*fr,*graph,*x);
 269
 270   copy_mat(parm->box,box);
 271
 272   snew(*pot,nsb->natoms);
 273
 274   sfree(v);
 275 }