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/simple.h"
47 #include "nb_generic_cg.h"
48 #include "nonbonded.h"
49 #include "nb_kernel.h"
53 gmx_nb_generic_cg_kernel(t_nblist * nlist,
58 nb_kernel_data_t * kernel_data,
61 int nri,ntype,table_nelements,ielec,ivdw;
62 real facel,gbtabscale;
63 int n,is3,i3,k,nj0,nj1,j3,ggid,nnn,n0;
64 int ai0,ai1,ai,aj0,aj1,aj;
69 real qq,vcoul,krsq,vctot;
72 real rt,r,eps,eps2,Y,F,Geps,Heps2,VV,FF,Fp,fijD,fijR;
75 real Vvdw_rep,Vvdw_disp;
76 real ix,iy,iz,fix,fiy,fiz;
78 real dx,dy,dz,rsq,rinv;
79 real c6,c12,cexp1,cexp2,br;
99 fshift = fr->fshift[0];
100 Vc = kernel_data->energygrp_elec;
101 Vvdw = kernel_data->energygrp_vdw;
102 tabscale = kernel_data->table_elec_vdw->scale;
103 VFtab = kernel_data->table_elec_vdw->data;
105 /* avoid compiler warnings for cases that cannot happen */
111 /* 3 VdW parameters for buckingham, otherwise 2 */
112 nvdwparam = (nlist->ivdw==2) ? 3 : 2;
113 table_nelements = (ielec==3) ? 4 : 0;
114 table_nelements += (ivdw==3) ? 8 : 0;
116 charge = mdatoms->chargeA;
117 type = mdatoms->typeA;
119 shiftvec = fr->shift_vec[0];
123 for(n=0; (n<nlist->nri); n++)
125 is3 = 3*nlist->shift[n];
127 shY = shiftvec[is3+1];
128 shZ = shiftvec[is3+2];
129 nj0 = nlist->jindex[n];
130 nj1 = nlist->jindex[n+1];
131 ai0 = nlist->iinr[n];
132 ai1 = nlist->iinr_end[n];
139 for(k=nj0; (k<nj1); k++)
141 aj0 = nlist->jjnr[k];
142 aj1 = nlist->jjnr_end[k];
143 excl = &nlist->excl[k*MAX_CGCGSIZE];
145 for(ai=ai0; (ai<ai1); ai++)
151 iq = facel*charge[ai];
152 nti = nvdwparam*ntype*type[ai];
154 /* Note that this code currently calculates
155 * all LJ and Coulomb interactions,
156 * even if the LJ parameters or charges are zero.
157 * If required, this can be optimized.
160 for(aj=aj0; (aj<aj1); aj++)
162 /* Check if this interaction is excluded */
163 if (excl[aj-aj0] & (1<<(ai-ai0)))
175 rsq = dx*dx+dy*dy+dz*dz;
176 rinv = gmx_invsqrt(rsq);
180 if (ielec==3 || ivdw==3)
187 nnn = table_nelements*n0;
190 /* Coulomb interaction. ielec==0 means no interaction */
198 /* Vanilla cutoff coulomb */
200 fscal = vcoul*rinvsq;
206 vcoul = qq*(rinv+krsq-fr->c_rf);
207 fscal = qq*(rinv-2.0*krsq)*rinvsq;
211 /* Tabulated coulomb */
214 Geps = eps*VFtab[nnn+2];
215 Heps2 = eps2*VFtab[nnn+3];
219 FF = Fp+Geps+2.0*Heps2;
221 fscal = -qq*FF*tabscale*rinv;
226 gmx_fatal(FARGS,"Death & horror! GB generic interaction not implemented.\n");
230 gmx_fatal(FARGS,"Death & horror! No generic coulomb interaction for ielec=%d.\n",ielec);
234 } /* End of coulomb interactions */
237 /* VdW interaction. ivdw==0 means no interaction */
240 tj = nti+nvdwparam*type[aj];
245 /* Vanilla Lennard-Jones cutoff */
247 c12 = vdwparam[tj+1];
249 rinvsix = rinvsq*rinvsq*rinvsq;
250 Vvdw_disp = c6*rinvsix;
251 Vvdw_rep = c12*rinvsix*rinvsix;
252 fscal += (12.0*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
253 Vvdwtot = Vvdwtot+Vvdw_rep-Vvdw_disp;
259 cexp1 = vdwparam[tj+1];
260 cexp2 = vdwparam[tj+2];
262 rinvsix = rinvsq*rinvsq*rinvsq;
263 Vvdw_disp = c6*rinvsix;
265 Vvdw_rep = cexp1*exp(-br);
266 fscal += (br*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
267 Vvdwtot = Vvdwtot+Vvdw_rep-Vvdw_disp;
273 c12 = vdwparam[tj+1];
277 Geps = eps*VFtab[nnn+2];
278 Heps2 = eps2*VFtab[nnn+3];
281 FF = Fp+Geps+2.0*Heps2;
287 Geps = eps*VFtab[nnn+2];
288 Heps2 = eps2*VFtab[nnn+3];
291 FF = Fp+Geps+2.0*Heps2;
294 fscal += -(fijD+fijR)*tabscale*rinv;
295 Vvdwtot = Vvdwtot + Vvdw_disp + Vvdw_rep;
299 gmx_fatal(FARGS,"Death & horror! No generic VdW interaction for ivdw=%d.\n",ivdw);
302 } /* end VdW interactions */
322 fshift[is3+1] += fiy;
323 fshift[is3+2] += fiz;
324 ggid = nlist->gid[n];
326 Vvdw[ggid] += Vvdwtot;
328 /* Estimate flops, average for generic cg kernel:
329 * 12 flops per outer iteration
330 * 100 flops per inner iteration
332 inc_nrnb(nrnb,eNR_NBKERNEL_FREE_ENERGY,nlist->nri*12 + nlist->jindex[n]*100);