2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "types/commrec.h"
54 #include "gmx_fatal.h"
60 #include "mtop_util.h"
61 #include "chargegroup.h"
63 static void c_tabpot(real tabscale, real VFtab[],
66 int jindex[], int jjnr[],
68 real facel, real charge[],
69 real pot[], real shiftvec[])
72 const real nul = 0.000000;
75 const real one = 1.000000;
76 const real two = 2.000000;
77 real r1,r1t,fijC,eps,eps2,Y,F,Fp,Geps,Heps2,VV,FF;
80 /* General and coulomb stuff */
81 int ii,k,n,jnr,ii3,nj0,nj1,is3,j3,ggid;
82 real fxJ,fyJ,fzJ,fxO,fyO,fzO;
83 real ixO,iyO,izO,dxO,dyO,dzO;
84 real txO,tyO,tzO,vcO,fsO,qO,rsqO,rinv1O,rinv2O;
86 real jx,jy,jz,shX,shY,shZ,poti;
88 /* Outer loop (over i particles) starts here */
89 for(n=0; (n<nri); n++) {
91 /* Unpack shift vector */
94 shY = shiftvec[is3+1];
95 shZ = shiftvec[is3+2];
97 /* Unpack I particle */
101 /* Charge of i particle(s) divided by 4 pi eps0 */
102 qO = facel*charge[ii];
104 /* Bounds for the innerloop */
108 /* Compute shifted i position */
109 ixO = shX + pos[ii3];
110 iyO = shY + pos[ii3+1];
111 izO = shZ + pos[ii3+2];
114 /* Inner loop (over j-particles) starts right here */
115 for(k=nj0; (k<nj1); k++) {
117 /* Unpack neighbourlist */
120 qj = facel*charge[jnr];
125 /* First one is for oxygen, with LJ */
129 rsqO = dxO*dxO + dyO*dyO + dzO*dzO;
131 /* Doing fast gmx_invsqrt */
132 rinv1O = gmx_invsqrt(rsqO);
144 Geps = eps*VFtab[nnn+2];
145 Heps2 = eps2*VFtab[nnn+3];
157 static void low_calc_pot(FILE *log,int nl_type,t_forcerec *fr,
158 rvec x[],t_mdatoms *mdatoms,real pot[])
162 nlist = &fr->nblists[0].nlist_sr[nl_type];
164 c_tabpot(fr->nblists[0].table_elec_vdw.scale,fr->nblists[0].table_elec_vdw.data,
165 nlist->nri,nlist->iinr,
166 nlist->shift,nlist->jindex,nlist->jjnr,
167 x[0],fr->epsfac,mdatoms->chargeA,pot,fr->shift_vec[0]);
169 fprintf(log,"There were %d interactions\n",nlist->nrj);
172 void calc_pot(FILE *logf,t_commrec *cr,
174 t_inputrec *inputrec,gmx_localtop_t *top,rvec x[],
175 t_forcerec *fr,gmx_enerdata_t *enerd,
176 t_mdatoms *mdatoms,real pot[],matrix box,t_graph *graph)
179 real lam[efptNR],dum[efptNR];
184 fprintf(stderr,"Doing single force calculation...\n");
186 if (inputrec->ePBC != epbcNONE)
187 calc_shifts(box,fr->shift_vec);
189 put_charge_groups_in_box(logf,0,top->cgs.nr,fr->ePBC,box,&(top->cgs),
192 mk_mshift(logf,graph,fr->ePBC,box,x);
193 /* Do the actual neighbour searching and if twin range electrostatics
194 * also do the calculation of long range forces and energies.
196 ns(logf,fr,x,box,&mtop->groups,&(inputrec->opts),top,mdatoms,cr,
197 &nrnb,&lam[0],&dum[0],&enerd->grpp,TRUE,FALSE);
198 for(m=0; (m<DIM); m++)
199 box_size[m] = box[m][m];
200 for(i=0; (i<mdatoms->nr); i++)
203 pr_rvecs(debug,0,"x",x,mdatoms->nr);
204 pr_rvecs(debug,0,"cgcm",fr->cg_cm,top->cgs.nr);
206 /* electrostatics from any atom to atoms without LJ */
207 low_calc_pot(logf,eNL_QQ,fr,x,mdatoms,pot);
208 /* electrostatics from any atom to atoms without charge */
209 low_calc_pot(logf,eNL_VDW,fr,x,mdatoms,pot);
210 /* electrostatics from any atom to atoms with LJ */
211 low_calc_pot(logf,eNL_VDWQQ,fr,x,mdatoms,pot);
214 FILE *init_calcpot(const char *log,const char *tpx,const char *table,
215 gmx_mtop_t *mtop,gmx_localtop_t *top,
216 t_inputrec *inputrec,t_commrec **cr,
217 t_graph **graph,t_mdatoms **mdatoms,
219 gmx_enerdata_t *enerd,
221 matrix box,rvec **x, const output_env_t oenv)
223 gmx_localtop_t *ltop;
228 int traj=0,xtc_traj=0;
234 tensor force_vir,shake_vir;
238 *cr = init_par(NULL,NULL);
239 gmx_log_open(log,*cr,FALSE,0,&fplog);
243 read_tpx_state(tpx,inputrec,state, NULL, mtop);
244 set_state_entries(state,inputrec,1);
246 pr_inputrec(fplog,0,"Input Parameters",inputrec,FALSE);
250 if (inputrec->efep) {
251 fprintf(stderr,"WARNING: turning of free energy, will use lambda=0\n");
256 init_md(fplog,*cr,inputrec,oenv,&t,&t0,lam,&fep_state,NULL,
257 &nrnb,mtop,NULL,-1,NULL,NULL,NULL,
258 force_vir,shake_vir,mutot,&bSA,NULL,NULL,0);
260 init_enerdata(mtop->groups.grps[egcENER].nr,0,enerd);
262 ltop = gmx_mtop_generate_local_top(mtop,inputrec);
266 *mdatoms = init_mdatoms(fplog,mtop,FALSE);
267 atoms2md(mtop,inputrec,0,NULL,0,mtop->natoms,*mdatoms);
269 if (inputrec->ePBC == epbcXYZ) {
270 /* Calculate intramolecular shift vectors to make molecules whole again */
271 *graph = mk_graph(fplog,&(top->idef),0,mtop->natoms,FALSE,FALSE);
272 mk_mshift(fplog,*graph,inputrec->ePBC,state->box,state->x);
277 /* Turn off twin range if appropriate */
278 inputrec->rvdw = inputrec->rcoulomb;
279 inputrec->rlist = inputrec->rcoulomb;
280 fprintf(stderr,"Using a coulomb cut-off of %g nm\n",inputrec->rcoulomb);
282 /* Turn off free energy computation */
285 /* Set vanderwaals to shift, to force tables */
286 inputrec->vdwtype = evdwSHIFT;
287 inputrec->rvdw_switch = 0.0;
288 inputrec->eDispCorr = edispcNO;
290 /* Initiate forcerecord */
292 init_forcerec(fplog,oenv,*fr,NULL,inputrec,mtop,*cr,
293 state->box,FALSE,table,NULL,table,NULL,NULL,TRUE,-1);
295 /* Remove periodicity */
296 for(m=0; (m<DIM); m++)
297 box_size[m] = state->box[m][m];
298 if (inputrec->ePBC != epbcNONE)
299 do_pbc_first(fplog,state->box,*fr,*graph,state->x);
301 copy_mat(state->box,box);
307 snew(*pot,mtop->natoms);