changed reading hint
[gromacs/adressmacs.git] / src / tools / calcpot.c
blob419001a60f6bda4996ae4eed50e2640ba31a3f20
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 2.0
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
29 static char *SRCID_calcpot_c = "$Id$";
31 #include "vec.h"
32 #include "calcpot.h"
33 #include "nrnb.h"
34 #include "mdebin.h"
35 #include "mshift.h"
36 #include "smalloc.h"
37 #include "force.h"
38 #include "main.h"
39 #include "filenm.h"
40 #include "fatal.h"
41 #include "mdrun.h"
42 #include "ns.h"
43 #include "txtdump.h"
45 static void c_tabpot(real tabscale, real VFtab[],
46 int nri, int iinr[],
47 int shift[],
48 int jindex[], int jjnr[],
49 real pos[],
50 real facel, real charge[],
51 real pot[], real shiftvec[])
53 /* Local variables */
54 const real nul = 0.000000;
56 /* Table stuff */
57 const real one = 1.000000;
58 const real two = 2.000000;
59 real r1,r1t,fijC,eps,eps2,Y,F,Fp,Geps,Heps2,VV,FF;
60 int n0,n1,nnn;
62 /* General and coulomb stuff */
63 int ii,k,n,jnr,ii3,nj0,nj1,is3,j3,ggid;
64 real fxJ,fyJ,fzJ,fxO,fyO,fzO;
65 real ixO,iyO,izO,dxO,dyO,dzO;
66 real txO,tyO,tzO,vcO,fsO,qO,rsqO,rinv1O,rinv2O;
67 real qqO,qj;
68 real jx,jy,jz,shX,shY,shZ,poti;
70 /* Outer loop (over i particles) starts here */
71 for(n=0; (n<nri); n++) {
73 /* Unpack shift vector */
74 is3 = 3*shift[n];
75 shX = shiftvec[is3];
76 shY = shiftvec[is3+1];
77 shZ = shiftvec[is3+2];
79 /* Unpack I particle */
80 ii = iinr[n];
81 ii3 = 3*ii;
83 /* Charge of i particle(s) divided by 4 pi eps0 */
84 qO = facel*charge[ii];
86 /* Bounds for the innerloop */
87 nj0 = jindex[n];
88 nj1 = jindex[n+1];
90 /* Compute shifted i position */
91 ixO = shX + pos[ii3];
92 iyO = shY + pos[ii3+1];
93 izO = shZ + pos[ii3+2];
94 poti = nul;
96 /* Inner loop (over j-particles) starts right here */
97 for(k=nj0; (k<nj1); k++) {
99 /* Unpack neighbourlist */
100 jnr = jjnr[k];
101 j3 = 3*jnr;
102 qj = facel*charge[jnr];
103 jx = pos[j3];
104 jy = pos[j3+1];
105 jz = pos[j3+2];
107 /* First one is for oxygen, with LJ */
108 dxO = ixO - jx;
109 dyO = iyO - jy;
110 dzO = izO - jz;
111 rsqO = dxO*dxO + dyO*dyO + dzO*dzO;
113 /* Doing fast invsqrt */
114 rinv1O = invsqrt(rsqO);
116 /* O block */
117 r1 = one/rinv1O;
118 r1t = r1*tabscale;
119 n0 = r1t;
120 n1 = 12*n0;
121 eps = r1t-n0;
122 eps2 = eps*eps;
123 nnn = n1;
124 Y = VFtab[nnn];
125 F = VFtab[nnn+1];
126 Geps = eps*VFtab[nnn+2];
127 Heps2 = eps2*VFtab[nnn+3];
128 Fp = F+Geps+Heps2;
129 VV = Y+eps*Fp;
131 pot[jnr] += VV*qO;
132 poti += VV*qj;
135 pot[ii] += poti;
139 static void low_calc_pot(FILE *log,int ftype,t_forcerec *fr,
140 rvec x[],t_mdatoms *mdatoms,real pot[])
142 t_nblist *nlist;
144 if (ftype == F_SR)
145 nlist = &fr->nlist_sr[eNL_QQ];
146 else
147 nlist = &fr->nlist_sr[eNL_VDW];
149 c_tabpot(fr->tabscale,fr->VFtab,nlist->nri,nlist->iinr,
150 nlist->shift,nlist->jindex,nlist->jjnr,
151 x[0],fr->epsfac,mdatoms->chargeA,pot,fr->shift_vec[0]);
153 fprintf(log,"There were %d interactions\n",nlist->nrj);
156 void calc_pot(FILE *logf,t_nsborder *nsb,t_commrec *cr,t_groups *grps,
157 t_parm *parm,t_topology *top,rvec x[],t_forcerec *fr,
158 t_mdatoms *mdatoms,real pot[])
160 static bool bFirst=TRUE;
161 static t_nrnb nrnb;
162 static rvec *f;
163 real lam=0,dum=0;
164 rvec box_size;
165 int i,m;
167 /* Calc the force */
168 fprintf(stderr,"Doing single force calculation...\n");
170 if (bFirst) {
171 snew(f, nsb->natoms);
173 bFirst = FALSE;
175 /* Reset long range forces if necessary */
176 if (fr->bTwinRange) {
177 clear_rvecs(nsb->natoms,fr->flr);
178 clear_rvecs(SHIFTS,fr->fshift_lr);
180 if (parm->ir.epc != epcNO)
181 calc_shifts(parm->box,box_size,fr->shift_vec,FALSE);
182 put_charge_groups_in_box(stdlog,0,top->blocks[ebCGS].nr,FALSE,
183 parm->box,box_size,&(top->blocks[ebCGS]),x,
184 fr->shift_vec,fr->cg_cm);
185 /* mk_mshift(stdlog,graph,parm->box,x);*/
186 /* Do the actual neighbour searching and if twin range electrostatics
187 * also do the calculation of long range forces and energies.
190 ns(logf,fr,x,f,parm->box,grps,&(parm->ir.opts),top,mdatoms,cr,
191 &nrnb,nsb,0,lam,&dum);
192 for(m=0; (m<DIM); m++)
193 box_size[m] = parm->box[m][m];
194 for(i=0; (i<mdatoms->nr); i++)
195 pot[i] = 0;
196 if (debug) {
197 pr_rvecs(debug,0,"x",x,mdatoms->nr);
198 pr_rvecs(debug,0,"cgcm",fr->cg_cm,top->blocks[ebCGS].nr);
200 /* electrostatics from any atom to atoms without LJ */
201 low_calc_pot(logf,F_SR,fr,x,mdatoms,pot);
202 /* electrostatics from any atom to atoms with LJ */
203 low_calc_pot(logf,F_LJ,fr,x,mdatoms,pot);
206 void init_calcpot(int nfile,t_filenm fnm[],t_topology *top,
207 rvec **x,t_parm *parm,t_commrec *cr,
208 t_graph **graph,t_mdatoms **mdatoms,
209 t_nsborder *nsb,t_groups *grps,
210 t_forcerec **fr,real **pot,
211 matrix box)
213 real t,t0,lam,lam0,SAfac;
214 bool bTYZ;
215 char *traj,*xtc_traj;
216 rvec *v;
217 t_nrnb nrnb;
218 t_mdebin *mdebin;
219 int fp_ene,m;
220 rvec vcm,box_size;
221 tensor force_vir,shake_vir;
223 /* Initiate */
224 cr->nprocs = 1; cr->pid = 0; cr->left = 0; cr->right = 1;
225 open_log(ftp2fn(efLOG,nfile,fnm),cr);
226 #ifdef CINVSQRT
227 init_lookup_table(stdlog);
228 #endif
230 init_nrnb(&nrnb);
231 init_single(stdlog,parm,ftp2fn(efTPX,nfile,fnm),top,x,&v,mdatoms,nsb);
232 init_md(cr,&(parm->ir),&t,&t0,&lam,&lam0,&SAfac,
233 &nrnb,&bTYZ,top,-1,NULL,&traj,&xtc_traj,&fp_ene,NULL,
234 &mdebin,grps,vcm,force_vir,shake_vir,*mdatoms);
235 init_groups(stdlog,*mdatoms,&(parm->ir.opts),grps);
237 /* Calculate intramolecular shift vectors to make molecules whole again */
238 *graph = mk_graph(&(top->idef),top->atoms.nr,FALSE);
239 mk_mshift(stdlog,*graph,parm->box,*x);
241 /* Turn off watertype optimizations, to ease coding above. */
242 parm->ir.solvent_opt = -1;
244 /* Turn off twin range if appropriate */
245 parm->ir.rvdw = parm->ir.rcoulomb;
246 parm->ir.rlist = parm->ir.rcoulomb;
247 fprintf(stderr,"Using a coulomb cut-off of %g nm\n",parm->ir.rcoulomb);
249 /* Turn off free energy computation */
250 parm->ir.bPert = FALSE;
252 /* Set vanderwaals to shift, to force tables */
253 parm->ir.vdwtype = evdwSHIFT;
254 parm->ir.rvdw_switch = 0.0;
256 /* Initiate forcerecord */
257 *fr = mk_forcerec();
258 init_forcerec(stdlog,*fr,&(parm->ir),&(top->blocks[ebMOLS]),cr,
259 &(top->blocks[ebCGS]),&(top->idef),*mdatoms,
260 nsb,parm->box,FALSE);
262 /* Remove periodicity */
263 for(m=0; (m<DIM); m++)
264 box_size[m] = parm->box[m][m];
265 if (parm->ir.eBox != ebtNONE)
266 do_pbc_first(stdlog,parm,box_size,*fr,*graph,*x);
268 copy_mat(parm->box,box);
270 snew(*pot,nsb->natoms);
272 sfree(v);