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 * GROningen Mixture of Alchemy and Childrens' Stories
41 #include "types/simple.h"
44 #include "nb_generic_adress.h"
46 #define ALMOST_ZERO 1e-30
47 #define ALMOST_ONE 1-(1e-30)
49 gmx_nb_generic_adress_kernel(t_nblist * nlist,
63 int nri,ntype,table_nelements,ielec,ivdw;
64 real facel,gbtabscale;
65 int n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid,nnn,n0;
70 real qq,vcoul,krsq,vctot;
73 real rt,r,eps,eps2,Y,F,Geps,Heps2,VV,FF,Fp,fijD,fijR;
76 real Vvdw_rep,Vvdw_disp;
77 real ix,iy,iz,fix,fiy,fiz;
79 real dx,dy,dz,rsq,rinv;
80 real c6,c12,cexp1,cexp2,br;
91 real hybscal; /* the multiplicator to the force for hybrid interactions*/
92 gmx_bool bHybrid; /*Are we in the hybrid zone ?*/
97 force_cap = fr->adress_ex_forcecap;
102 /* avoid compiler warnings for cases that cannot happen */
108 /* 3 VdW parameters for buckingham, otherwise 2 */
109 nvdwparam = (nlist->ivdw==2) ? 3 : 2;
110 table_nelements = (ielec==3) ? 4 : 0;
111 table_nelements += (ivdw==3) ? 8 : 0;
113 charge = mdatoms->chargeA;
114 type = mdatoms->typeA;
116 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];
136 iq = facel*charge[ii];
137 nti = nvdwparam*ntype*type[ii];
146 /* TODO: why does this line her not speed up things ?
147 * if ((!bCG) && weight_cg1 < ALMOST_ZERO) continue;
149 for(k=nj0; (k<nj1); k++)
151 jnr = nlist->jjnr[k];
152 weight_cg2 = wf[jnr];
154 weight_product = weight_cg1*weight_cg2;
156 if (weight_product < ALMOST_ZERO)
158 /* if it's a explicit loop, skip this atom */
163 else /* if it's a coarse grained loop, include this atom */
169 else if (weight_product >= ALMOST_ONE)
172 /* if it's a explicit loop, include this atom */
178 else /* if it's a coarse grained loop, skip this atom */
183 /* both have double identity, get hybrid scaling factor */
187 hybscal = weight_product;
191 hybscal = 1.0 - hybscal;
203 rsq = dx*dx+dy*dy+dz*dz;
204 rinv = gmx_invsqrt(rsq);
210 if(ielec==3 || ivdw==3)
217 nnn = table_nelements*n0;
220 /* Coulomb interaction. ielec==0 means no interaction */
228 /* Vanilla cutoff coulomb */
230 fscal = vcoul*rinvsq;
236 vcoul = qq*(rinv+krsq-fr->c_rf);
237 fscal = qq*(rinv-2.0*krsq)*rinvsq;
241 /* Tabulated coulomb */
244 Geps = eps*VFtab[nnn+2];
245 Heps2 = eps2*VFtab[nnn+3];
249 FF = Fp+Geps+2.0*Heps2;
251 fscal = -qq*FF*tabscale*rinv;
256 gmx_fatal(FARGS,"Death & horror! GB generic interaction not implemented.\n");
260 gmx_fatal(FARGS,"Death & horror! No generic coulomb interaction for ielec=%d.\n",ielec);
264 } /* End of coulomb interactions */
266 /* VdW interaction. ivdw==0 means no interaction */
269 tj = nti+nvdwparam*type[jnr];
274 /* Vanilla Lennard-Jones cutoff */
276 c12 = vdwparam[tj+1];
278 rinvsix = rinvsq*rinvsq*rinvsq;
279 Vvdw_disp = c6*rinvsix;
280 Vvdw_rep = c12*rinvsix*rinvsix;
281 fscal += (12.0*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
282 Vvdwtot = Vvdwtot+Vvdw_rep-Vvdw_disp;
288 cexp1 = vdwparam[tj+1];
289 cexp2 = vdwparam[tj+2];
291 rinvsix = rinvsq*rinvsq*rinvsq;
292 Vvdw_disp = c6*rinvsix;
294 Vvdw_rep = cexp1*exp(-br);
295 fscal += (br*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
296 Vvdwtot = Vvdwtot+Vvdw_rep-Vvdw_disp;
302 c12 = vdwparam[tj+1];
306 Geps = eps*VFtab[nnn+2];
307 Heps2 = eps2*VFtab[nnn+3];
310 FF = Fp+Geps+2.0*Heps2;
316 Geps = eps*VFtab[nnn+2];
317 Heps2 = eps2*VFtab[nnn+3];
320 FF = Fp+Geps+2.0*Heps2;
323 fscal += -(fijD+fijR)*tabscale*rinv;
324 Vvdwtot = Vvdwtot + Vvdw_disp + Vvdw_rep;
325 if(!bCG && force_cap>0 && (fabs(fscal)> force_cap))
327 fscal=force_cap*fscal/fabs(fscal);
332 gmx_fatal(FARGS,"Death & horror! No generic VdW interaction for ivdw=%d.\n",ivdw);
335 } /* end VdW interactions */
337 /* force weight is one anyway */
349 f[j3+0] = f[j3+0] - tx;
350 f[j3+1] = f[j3+1] - ty;
351 f[j3+2] = f[j3+2] - tz;
354 f[ii3+0] = f[ii3+0] + fix;
355 f[ii3+1] = f[ii3+1] + fiy;
356 f[ii3+2] = f[ii3+2] + fiz;
357 fshift[is3] = fshift[is3]+fix;
358 fshift[is3+1] = fshift[is3+1]+fiy;
359 fshift[is3+2] = fshift[is3+2]+fiz;
360 ggid = nlist->gid[n];
361 Vc[ggid] = Vc[ggid] + vctot;
362 Vvdw[ggid] = Vvdw[ggid] + Vvdwtot;
365 *outeriter = nlist->nri;
366 *inneriter = nlist->jindex[n];