Merge remote-tracking branch 'origin/release-4-6' into HEAD
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_adress_c / nb_kernel131_c_adress.c
1 /*
2  * Copyright (c) Erik Lindahl, David van der Spoel 2003
3  * 
4  * This file is generated automatically at compile time
5  * by the program mknb in the Gromacs distribution.
6  *
7  * Options used when generation this file:
8  * Language:         c
9  * Precision:        single
10  * Threads:          yes
11  * Software invsqrt: no
12  * PowerPC invsqrt:  no
13  * Prefetch forces:  no
14  * Adress kernel:  yes
15  * Comments:         no
16  */
17 #ifdef HAVE_CONFIG_H
18 #include<config.h>
19 #endif
20 #ifdef GMX_THREAD_SHM_FDECOMP
21 #include<thread_mpi.h>
22 #endif
23 #define ALMOST_ZERO 1e-30
24 #define ALMOST_ONE 1-(1e-30)
25 #include<math.h>
26
27 #include "nb_kernel131_adress.h"
28
29
30
31 /*
32  * Gromacs nonbonded kernel nb_kernel131_adress_cg
33  * Coulomb interaction:     Normal Coulomb
34  * VdW interaction:         Tabulated
35  * water optimization:      SPC/TIP3P - other atoms
36  * Calculate forces:        yes
37  */
38 void nb_kernel131_adress_cg(
39                     int *           p_nri,
40                     int *           iinr,
41                     int *           jindex,
42                     int *           jjnr,
43                     int *           shift,
44                     real *         shiftvec,
45                     real *         fshift,
46                     int *           gid,
47                     real *         pos,
48                     real *         faction,
49                     real *         charge,
50                     real *         p_facel,
51                     real *         p_krf,
52                     real *         p_crf,
53                     real *         Vc,
54                     int *           type,
55                     int *           p_ntype,
56                     real *         vdwparam,
57                     real *         Vvdw,
58                     real *         p_tabscale,
59                     real *         VFtab,
60                     real *         invsqrta,
61                     real *         dvda,
62                     real *         p_gbtabscale,
63                     real *         GBtab,
64                     int *           p_nthreads,
65                     int *           count,
66                     void *          mtx,
67                     int *           outeriter,
68                     int *           inneriter,
69                     real           force_cap,
70                     real *         wf)
71 {
72     int           nri,ntype,nthreads;
73     real         facel,krf,crf,tabscale,gbtabscale;
74     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
75     int           nn0,nn1,nouter,ninner;
76     real         shX,shY,shZ;
77     real         fscal,tx,ty,tz;
78     real         rinvsq;
79     real         jq;
80     real         qq,vcoul,vctot;
81     int           nti;
82     int           tj;
83     real         Vvdw6,Vvdwtot;
84     real         Vvdw12;
85     real         r,rt,eps,eps2;
86     int           n0,nnn;
87     real         Y,F,Geps,Heps2,Fp,VV;
88     real         FF;
89     real         fijD,fijR;
90     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
91     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
92     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
93     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
94     real         dx11,dy11,dz11,rsq11,rinv11;
95     real         dx21,dy21,dz21,rsq21,rinv21;
96     real         dx31,dy31,dz31,rsq31,rinv31;
97     real         qO,qH;
98     real         c6,c12;
99     real         weight_cg1, weight_cg2, weight_product;
100     real         hybscal;
101
102     nri              = *p_nri;         
103     ntype            = *p_ntype;       
104     nthreads         = *p_nthreads;    
105     facel            = *p_facel;       
106     krf              = *p_krf;         
107     crf              = *p_crf;         
108     tabscale         = *p_tabscale;    
109     ii               = iinr[0];        
110     qO               = facel*charge[ii];
111     qH               = facel*charge[ii+1];
112     nti              = 2*ntype*type[ii];
113
114     nouter           = 0;              
115     ninner           = 0;              
116     
117     do
118     {
119         #ifdef GMX_THREAD_SHM_FDECOMP
120         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
121         nn0              = *count;         
122         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
123         *count           = nn1;            
124         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
125         if(nn1>nri) nn1=nri;
126         #else
127         nn0 = 0;
128         nn1 = nri;
129         #endif
130         
131         for(n=nn0; (n<nn1); n++)
132         {
133             is3              = 3*shift[n];     
134             shX              = shiftvec[is3];  
135             shY              = shiftvec[is3+1];
136             shZ              = shiftvec[is3+2];
137             nj0              = jindex[n];      
138             nj1              = jindex[n+1];    
139             ii               = iinr[n];        
140             ii3              = 3*ii;           
141             ix1              = shX + pos[ii3+0];
142             iy1              = shY + pos[ii3+1];
143             iz1              = shZ + pos[ii3+2];
144             ix2              = shX + pos[ii3+3];
145             iy2              = shY + pos[ii3+4];
146             iz2              = shZ + pos[ii3+5];
147             ix3              = shX + pos[ii3+6];
148             iy3              = shY + pos[ii3+7];
149             iz3              = shZ + pos[ii3+8];
150             weight_cg1       = wf[ii];         
151             vctot            = 0;              
152             Vvdwtot          = 0;              
153             fix1             = 0;              
154             fiy1             = 0;              
155             fiz1             = 0;              
156             fix2             = 0;              
157             fiy2             = 0;              
158             fiz2             = 0;              
159             fix3             = 0;              
160             fiy3             = 0;              
161             fiz3             = 0;              
162             
163             for(k=nj0; (k<nj1); k++)
164             {
165                 jnr              = jjnr[k];        
166                 weight_cg2       = wf[jnr];        
167                 weight_product   = weight_cg1*weight_cg2;
168                 if (weight_product < ALMOST_ZERO) {
169                        hybscal = 1.0;
170                 }
171                 else if (weight_product >= ALMOST_ONE)
172                 {
173                   /* force is zero, skip this molecule */
174                        continue;
175                 }
176                 else
177                 {
178                    hybscal = 1.0 - weight_product;
179                 }
180                 j3               = 3*jnr;          
181                 jx1              = pos[j3+0];      
182                 jy1              = pos[j3+1];      
183                 jz1              = pos[j3+2];      
184                 dx11             = ix1 - jx1;      
185                 dy11             = iy1 - jy1;      
186                 dz11             = iz1 - jz1;      
187                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
188                 dx21             = ix2 - jx1;      
189                 dy21             = iy2 - jy1;      
190                 dz21             = iz2 - jz1;      
191                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
192                 dx31             = ix3 - jx1;      
193                 dy31             = iy3 - jy1;      
194                 dz31             = iz3 - jz1;      
195                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
196                 rinv11           = 1.0/sqrt(rsq11);
197                 rinv21           = 1.0/sqrt(rsq21);
198                 rinv31           = 1.0/sqrt(rsq31);
199                 jq               = charge[jnr+0];  
200                 qq               = qO*jq;          
201                 tj               = nti+2*type[jnr];
202                 c6               = vdwparam[tj];   
203                 c12              = vdwparam[tj+1]; 
204                 rinvsq           = rinv11*rinv11;  
205                 vcoul            = qq*rinv11;      
206                 vctot            = vctot+vcoul;    
207                 r                = rsq11*rinv11;   
208                 rt               = r*tabscale;     
209                 n0               = rt;             
210                 eps              = rt-n0;          
211                 eps2             = eps*eps;        
212                 nnn              = 8*n0;           
213                 Y                = VFtab[nnn];     
214                 F                = VFtab[nnn+1];   
215                 Geps             = eps*VFtab[nnn+2];
216                 Heps2            = eps2*VFtab[nnn+3];
217                 Fp               = F+Geps+Heps2;   
218                 VV               = Y+eps*Fp;       
219                 FF               = Fp+Geps+2.0*Heps2;
220                 Vvdw6            = c6*VV;          
221                 fijD             = c6*FF;          
222                 nnn              = nnn+4;          
223                 Y                = VFtab[nnn];     
224                 F                = VFtab[nnn+1];   
225                 Geps             = eps*VFtab[nnn+2];
226                 Heps2            = eps2*VFtab[nnn+3];
227                 Fp               = F+Geps+Heps2;   
228                 VV               = Y+eps*Fp;       
229                 FF               = Fp+Geps+2.0*Heps2;
230                 Vvdw12           = c12*VV;         
231                 fijR             = c12*FF;         
232                 Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12;
233                 fscal            = (vcoul)*rinvsq-((fijD+fijR)*tabscale)*rinv11;
234                 fscal *= hybscal;
235                 tx               = fscal*dx11;     
236                 ty               = fscal*dy11;     
237                 tz               = fscal*dz11;     
238                 fix1             = fix1 + tx;      
239                 fiy1             = fiy1 + ty;      
240                 fiz1             = fiz1 + tz;      
241                 fjx1             = faction[j3+0] - tx;
242                 fjy1             = faction[j3+1] - ty;
243                 fjz1             = faction[j3+2] - tz;
244                 qq               = qH*jq;          
245                 rinvsq           = rinv21*rinv21;  
246                 vcoul            = qq*rinv21;      
247                 vctot            = vctot+vcoul;    
248                 fscal            = (vcoul)*rinvsq; 
249                 fscal *= hybscal;
250                 tx               = fscal*dx21;     
251                 ty               = fscal*dy21;     
252                 tz               = fscal*dz21;     
253                 fix2             = fix2 + tx;      
254                 fiy2             = fiy2 + ty;      
255                 fiz2             = fiz2 + tz;      
256                 fjx1             = fjx1 - tx;      
257                 fjy1             = fjy1 - ty;      
258                 fjz1             = fjz1 - tz;      
259                 rinvsq           = rinv31*rinv31;  
260                 vcoul            = qq*rinv31;      
261                 vctot            = vctot+vcoul;    
262                 fscal            = (vcoul)*rinvsq; 
263                 fscal *= hybscal;
264                 tx               = fscal*dx31;     
265                 ty               = fscal*dy31;     
266                 tz               = fscal*dz31;     
267                 fix3             = fix3 + tx;      
268                 fiy3             = fiy3 + ty;      
269                 fiz3             = fiz3 + tz;      
270                 faction[j3+0]    = fjx1 - tx;      
271                 faction[j3+1]    = fjy1 - ty;      
272                 faction[j3+2]    = fjz1 - tz;      
273             }
274             
275             faction[ii3+0]   = faction[ii3+0] + fix1;
276             faction[ii3+1]   = faction[ii3+1] + fiy1;
277             faction[ii3+2]   = faction[ii3+2] + fiz1;
278             faction[ii3+3]   = faction[ii3+3] + fix2;
279             faction[ii3+4]   = faction[ii3+4] + fiy2;
280             faction[ii3+5]   = faction[ii3+5] + fiz2;
281             faction[ii3+6]   = faction[ii3+6] + fix3;
282             faction[ii3+7]   = faction[ii3+7] + fiy3;
283             faction[ii3+8]   = faction[ii3+8] + fiz3;
284             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
285             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
286             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
287             ggid             = gid[n];         
288             Vc[ggid]         = Vc[ggid] + vctot;
289             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
290             ninner           = ninner + nj1 - nj0;
291         }
292         
293         nouter           = nouter + nn1 - nn0;
294     }
295     while (nn1<nri);
296     
297     *outeriter       = nouter;         
298     *inneriter       = ninner;         
299 }
300
301
302
303
304
305 /*
306  * Gromacs nonbonded kernel nb_kernel131_adress_ex
307  * Coulomb interaction:     Normal Coulomb
308  * VdW interaction:         Tabulated
309  * water optimization:      SPC/TIP3P - other atoms
310  * Calculate forces:        yes
311  */
312 void nb_kernel131_adress_ex(
313                     int *           p_nri,
314                     int *           iinr,
315                     int *           jindex,
316                     int *           jjnr,
317                     int *           shift,
318                     real *         shiftvec,
319                     real *         fshift,
320                     int *           gid,
321                     real *         pos,
322                     real *         faction,
323                     real *         charge,
324                     real *         p_facel,
325                     real *         p_krf,
326                     real *         p_crf,
327                     real *         Vc,
328                     int *           type,
329                     int *           p_ntype,
330                     real *         vdwparam,
331                     real *         Vvdw,
332                     real *         p_tabscale,
333                     real *         VFtab,
334                     real *         invsqrta,
335                     real *         dvda,
336                     real *         p_gbtabscale,
337                     real *         GBtab,
338                     int *           p_nthreads,
339                     int *           count,
340                     void *          mtx,
341                     int *           outeriter,
342                     int *           inneriter,
343                     real           force_cap,
344                     real *         wf)
345 {
346     int           nri,ntype,nthreads;
347     real         facel,krf,crf,tabscale,gbtabscale;
348     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
349     int           nn0,nn1,nouter,ninner;
350     real         shX,shY,shZ;
351     real         fscal,tx,ty,tz;
352     real         rinvsq;
353     real         jq;
354     real         qq,vcoul,vctot;
355     int           nti;
356     int           tj;
357     real         Vvdw6,Vvdwtot;
358     real         Vvdw12;
359     real         r,rt,eps,eps2;
360     int           n0,nnn;
361     real         Y,F,Geps,Heps2,Fp,VV;
362     real         FF;
363     real         fijD,fijR;
364     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
365     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
366     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
367     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
368     real         dx11,dy11,dz11,rsq11,rinv11;
369     real         dx21,dy21,dz21,rsq21,rinv21;
370     real         dx31,dy31,dz31,rsq31,rinv31;
371     real         qO,qH;
372     real         c6,c12;
373     real         weight_cg1, weight_cg2, weight_product;
374     real         hybscal;
375
376     nri              = *p_nri;         
377     ntype            = *p_ntype;       
378     nthreads         = *p_nthreads;    
379     facel            = *p_facel;       
380     krf              = *p_krf;         
381     crf              = *p_crf;         
382     tabscale         = *p_tabscale;    
383     ii               = iinr[0];        
384     qO               = facel*charge[ii];
385     qH               = facel*charge[ii+1];
386     nti              = 2*ntype*type[ii];
387
388     nouter           = 0;              
389     ninner           = 0;              
390     
391     do
392     {
393         #ifdef GMX_THREAD_SHM_FDECOMP
394         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
395         nn0              = *count;         
396         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
397         *count           = nn1;            
398         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
399         if(nn1>nri) nn1=nri;
400         #else
401         nn0 = 0;
402         nn1 = nri;
403         #endif
404         
405         for(n=nn0; (n<nn1); n++)
406         {
407             is3              = 3*shift[n];     
408             shX              = shiftvec[is3];  
409             shY              = shiftvec[is3+1];
410             shZ              = shiftvec[is3+2];
411             nj0              = jindex[n];      
412             nj1              = jindex[n+1];    
413             ii               = iinr[n];        
414             ii3              = 3*ii;           
415             ix1              = shX + pos[ii3+0];
416             iy1              = shY + pos[ii3+1];
417             iz1              = shZ + pos[ii3+2];
418             ix2              = shX + pos[ii3+3];
419             iy2              = shY + pos[ii3+4];
420             iz2              = shZ + pos[ii3+5];
421             ix3              = shX + pos[ii3+6];
422             iy3              = shY + pos[ii3+7];
423             iz3              = shZ + pos[ii3+8];
424             weight_cg1       = wf[ii];         
425             vctot            = 0;              
426             Vvdwtot          = 0;              
427             fix1             = 0;              
428             fiy1             = 0;              
429             fiz1             = 0;              
430             fix2             = 0;              
431             fiy2             = 0;              
432             fiz2             = 0;              
433             fix3             = 0;              
434             fiy3             = 0;              
435             fiz3             = 0;              
436             
437             for(k=nj0; (k<nj1); k++)
438             {
439                 jnr              = jjnr[k];        
440                 weight_cg2       = wf[jnr];        
441                 weight_product   = weight_cg1*weight_cg2;
442                 if (weight_product < ALMOST_ZERO) {
443                 /* force is zero, skip this molecule */
444                  continue;
445                 }
446                 else if (weight_product >= ALMOST_ONE)
447                 {
448                        hybscal = 1.0;
449                 }
450                 else
451                 {
452                    hybscal = weight_product;
453                 }
454                 j3               = 3*jnr;          
455                 jx1              = pos[j3+0];      
456                 jy1              = pos[j3+1];      
457                 jz1              = pos[j3+2];      
458                 dx11             = ix1 - jx1;      
459                 dy11             = iy1 - jy1;      
460                 dz11             = iz1 - jz1;      
461                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
462                 dx21             = ix2 - jx1;      
463                 dy21             = iy2 - jy1;      
464                 dz21             = iz2 - jz1;      
465                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
466                 dx31             = ix3 - jx1;      
467                 dy31             = iy3 - jy1;      
468                 dz31             = iz3 - jz1;      
469                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
470                 rinv11           = 1.0/sqrt(rsq11);
471                 rinv21           = 1.0/sqrt(rsq21);
472                 rinv31           = 1.0/sqrt(rsq31);
473                 jq               = charge[jnr+0];  
474                 qq               = qO*jq;          
475                 tj               = nti+2*type[jnr];
476                 c6               = vdwparam[tj];   
477                 c12              = vdwparam[tj+1]; 
478                 rinvsq           = rinv11*rinv11;  
479                 vcoul            = qq*rinv11;      
480                 vctot            = vctot+vcoul;    
481                 r                = rsq11*rinv11;   
482                 rt               = r*tabscale;     
483                 n0               = rt;             
484                 eps              = rt-n0;          
485                 eps2             = eps*eps;        
486                 nnn              = 8*n0;           
487                 Y                = VFtab[nnn];     
488                 F                = VFtab[nnn+1];   
489                 Geps             = eps*VFtab[nnn+2];
490                 Heps2            = eps2*VFtab[nnn+3];
491                 Fp               = F+Geps+Heps2;   
492                 VV               = Y+eps*Fp;       
493                 FF               = Fp+Geps+2.0*Heps2;
494                 Vvdw6            = c6*VV;          
495                 fijD             = c6*FF;          
496                 nnn              = nnn+4;          
497                 Y                = VFtab[nnn];     
498                 F                = VFtab[nnn+1];   
499                 Geps             = eps*VFtab[nnn+2];
500                 Heps2            = eps2*VFtab[nnn+3];
501                 Fp               = F+Geps+Heps2;   
502                 VV               = Y+eps*Fp;       
503                 FF               = Fp+Geps+2.0*Heps2;
504                 Vvdw12           = c12*VV;         
505                 fijR             = c12*FF;         
506                 Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12;
507                 fscal            = (vcoul)*rinvsq-((fijD+fijR)*tabscale)*rinv11;
508                 fscal *= hybscal;
509                 if(force_cap>0 && (fabs(fscal)> force_cap)){
510                 fscal=force_cap*fscal/fabs(fscal);
511                 }
512                 tx               = fscal*dx11;     
513                 ty               = fscal*dy11;     
514                 tz               = fscal*dz11;     
515                 fix1             = fix1 + tx;      
516                 fiy1             = fiy1 + ty;      
517                 fiz1             = fiz1 + tz;      
518                 fjx1             = faction[j3+0] - tx;
519                 fjy1             = faction[j3+1] - ty;
520                 fjz1             = faction[j3+2] - tz;
521                 qq               = qH*jq;          
522                 rinvsq           = rinv21*rinv21;  
523                 vcoul            = qq*rinv21;      
524                 vctot            = vctot+vcoul;    
525                 fscal            = (vcoul)*rinvsq; 
526                 fscal *= hybscal;
527                 if(force_cap>0 && (fabs(fscal)> force_cap)){
528                 fscal=force_cap*fscal/fabs(fscal);
529                 }
530                 tx               = fscal*dx21;     
531                 ty               = fscal*dy21;     
532                 tz               = fscal*dz21;     
533                 fix2             = fix2 + tx;      
534                 fiy2             = fiy2 + ty;      
535                 fiz2             = fiz2 + tz;      
536                 fjx1             = fjx1 - tx;      
537                 fjy1             = fjy1 - ty;      
538                 fjz1             = fjz1 - tz;      
539                 rinvsq           = rinv31*rinv31;  
540                 vcoul            = qq*rinv31;      
541                 vctot            = vctot+vcoul;    
542                 fscal            = (vcoul)*rinvsq; 
543                 fscal *= hybscal;
544                 if(force_cap>0 && (fabs(fscal)> force_cap)){
545                 fscal=force_cap*fscal/fabs(fscal);
546                 }
547                 tx               = fscal*dx31;     
548                 ty               = fscal*dy31;     
549                 tz               = fscal*dz31;     
550                 fix3             = fix3 + tx;      
551                 fiy3             = fiy3 + ty;      
552                 fiz3             = fiz3 + tz;      
553                 faction[j3+0]    = fjx1 - tx;      
554                 faction[j3+1]    = fjy1 - ty;      
555                 faction[j3+2]    = fjz1 - tz;      
556             }
557             
558             faction[ii3+0]   = faction[ii3+0] + fix1;
559             faction[ii3+1]   = faction[ii3+1] + fiy1;
560             faction[ii3+2]   = faction[ii3+2] + fiz1;
561             faction[ii3+3]   = faction[ii3+3] + fix2;
562             faction[ii3+4]   = faction[ii3+4] + fiy2;
563             faction[ii3+5]   = faction[ii3+5] + fiz2;
564             faction[ii3+6]   = faction[ii3+6] + fix3;
565             faction[ii3+7]   = faction[ii3+7] + fiy3;
566             faction[ii3+8]   = faction[ii3+8] + fiz3;
567             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
568             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
569             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
570             ggid             = gid[n];         
571             Vc[ggid]         = Vc[ggid] + vctot;
572             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
573             ninner           = ninner + nj1 - nj0;
574         }
575         
576         nouter           = nouter + nn1 - nn0;
577     }
578     while (nn1<nri);
579     
580     *outeriter       = nouter;         
581     *inneriter       = ninner;         
582 }
583
584