3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
48 #include "gmx_fatal.h"
54 #include "mtop_util.h"
55 #include "chargegroup.h"
57 static void c_tabpot(real tabscale, real VFtab[],
60 int jindex[], int jjnr[],
62 real facel, real charge[],
63 real pot[], real shiftvec[])
66 const real nul = 0.000000;
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;
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;
80 real jx,jy,jz,shX,shY,shZ,poti;
82 /* Outer loop (over i particles) starts here */
83 for(n=0; (n<nri); n++) {
85 /* Unpack shift vector */
88 shY = shiftvec[is3+1];
89 shZ = shiftvec[is3+2];
91 /* Unpack I particle */
95 /* Charge of i particle(s) divided by 4 pi eps0 */
96 qO = facel*charge[ii];
98 /* Bounds for the innerloop */
102 /* Compute shifted i position */
103 ixO = shX + pos[ii3];
104 iyO = shY + pos[ii3+1];
105 izO = shZ + pos[ii3+2];
108 /* Inner loop (over j-particles) starts right here */
109 for(k=nj0; (k<nj1); k++) {
111 /* Unpack neighbourlist */
114 qj = facel*charge[jnr];
119 /* First one is for oxygen, with LJ */
123 rsqO = dxO*dxO + dyO*dyO + dzO*dzO;
125 /* Doing fast gmx_invsqrt */
126 rinv1O = gmx_invsqrt(rsqO);
138 Geps = eps*VFtab[nnn+2];
139 Heps2 = eps2*VFtab[nnn+3];
151 static void low_calc_pot(FILE *log,int nl_type,t_forcerec *fr,
152 rvec x[],t_mdatoms *mdatoms,real pot[])
156 nlist = &fr->nblists[0].nlist_sr[nl_type];
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]);
163 fprintf(log,"There were %d interactions\n",nlist->nrj);
166 void calc_pot(FILE *logf,t_commrec *cr,
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)
178 fprintf(stderr,"Doing single force calculation...\n");
180 if (inputrec->ePBC != epbcNONE)
181 calc_shifts(box,fr->shift_vec);
183 put_charge_groups_in_box(logf,0,top->cgs.nr,fr->ePBC,box,&(top->cgs),
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.
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++)
198 pr_rvecs(debug,0,"x",x,mdatoms->nr);
199 pr_rvecs(debug,0,"cgcm",fr->cg_cm,top->cgs.nr);
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);
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,
214 gmx_enerdata_t *enerd,
216 matrix box,rvec **x, const output_env_t oenv)
218 gmx_localtop_t *ltop;
222 int traj=0,xtc_traj=0;
228 tensor force_vir,shake_vir;
232 *cr = init_cr_nopar();
233 gmx_log_open(log,*cr,FALSE,0,&fplog);
237 read_tpx_state(tpx,inputrec,state, NULL, mtop);
238 set_state_entries(state,inputrec,1);
240 pr_inputrec(fplog,0,"Input Parameters",inputrec,FALSE);
244 if (inputrec->efep) {
245 fprintf(stderr,"WARNING: turning of free energy, will use lambda=0\n");
250 init_md(fplog,*cr,inputrec,oenv,&t,&t0,&lam,&lam0,
251 &nrnb,mtop,NULL,-1,NULL,NULL,NULL,
252 force_vir,shake_vir,mutot,&bSA,NULL,NULL,0);
254 init_enerdata(mtop->groups.grps[egcENER].nr,0,enerd);
256 ltop = gmx_mtop_generate_local_top(mtop,inputrec);
260 *mdatoms = init_mdatoms(fplog,mtop,FALSE);
261 atoms2md(mtop,inputrec,0,NULL,0,mtop->natoms,*mdatoms);
263 if (inputrec->ePBC == epbcXYZ) {
264 /* Calculate intramolecular shift vectors to make molecules whole again */
265 *graph = mk_graph(fplog,&(top->idef),0,mtop->natoms,FALSE,FALSE);
266 mk_mshift(fplog,*graph,inputrec->ePBC,state->box,state->x);
271 /* Turn off twin range if appropriate */
272 inputrec->rvdw = inputrec->rcoulomb;
273 inputrec->rlist = inputrec->rcoulomb;
274 fprintf(stderr,"Using a coulomb cut-off of %g nm\n",inputrec->rcoulomb);
276 /* Turn off free energy computation */
279 /* Set vanderwaals to shift, to force tables */
280 inputrec->vdwtype = evdwSHIFT;
281 inputrec->rvdw_switch = 0.0;
282 inputrec->eDispCorr = edispcNO;
284 /* Initiate forcerecord */
286 init_forcerec(fplog,oenv,*fr,NULL,inputrec,mtop,*cr,
287 state->box,FALSE,table,table,NULL,TRUE,-1);
289 /* Remove periodicity */
290 for(m=0; (m<DIM); m++)
291 box_size[m] = state->box[m][m];
292 if (inputrec->ePBC != epbcNONE)
293 do_pbc_first(fplog,state->box,*fr,*graph,state->x);
295 copy_mat(state->box,box);
301 snew(*pot,mtop->natoms);