?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Updated grompp.mdp for tutor/water
[gromacs.git] / src / mdlib / genborn.c
blob6fa2d7b6b048d625d086c9fc3a216e630e9047d8
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
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-2008, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
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.
19 *
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
34 */
35
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <math.h>
41 #include <string.h>
42
43 #include "typedefs.h"
44 #include "smalloc.h"
45 #include "genborn.h"
46 #include "vec.h"
47 #include "grompp.h"
48 #include "pdbio.h"
49 #include "names.h"
50 #include "physics.h"
51 #include "partdec.h"
52 #include "domdec.h"
53 #include "network.h"
54 #include "gmx_fatal.h"
55 #include "mtop_util.h"
56 #include "pbc.h"
57 #include "nrnb.h"
58 #include "bondf.h"
59
60 #ifdef GMX_LIB_MPI
61 #include <mpi.h>
62 #endif
63 #ifdef GMX_THREADS
64 #include "tmpi.h"
65 #endif
66
67 #ifdef GMX_DOUBLE
68 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_SSE2) )
69 #include "genborn_sse2_double.h"
70 #include "genborn_allvsall_sse2_double.h"
71 #endif
72 #else
73 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) )
74 #include "genborn_sse2_single.h"
75 #include "genborn_allvsall_sse2_single.h"
76 #endif /* GMX_SSE */
77 #endif /* GMX_DOUBLE */
78
79 #include "genborn_allvsall.h"
80
81 /*#define DISABLE_SSE*/
82
83 typedef struct {
84 int shift;
85 int naj;
86 int *aj;
87 int aj_nalloc;
88 } gbtmpnbl_t;
89
90 typedef struct gbtmpnbls {
91 int nlist;
92 gbtmpnbl_t *list;
93 int list_nalloc;
94 } t_gbtmpnbls;
95
96 /* This function is exactly the same as the one in bondfree.c. The reason
97 * it is copied here is that the bonded gb-interactions are evaluated
98 * not in calc_bonds, but rather in calc_gb_forces
99 */
100 static int pbc_rvec_sub(const t_pbc *pbc,const rvec xi,const rvec xj,rvec dx)
101 {
102 if (pbc) {
103 return pbc_dx_aiuc(pbc,xi,xj,dx);
104 }
105 else {
106 rvec_sub(xi,xj,dx);
107 return CENTRAL;
108 }
109 }
110
111 int init_gb_nblist(int natoms, t_nblist *nl)
112 {
113 nl->maxnri = natoms*4;
114 nl->maxnrj = 0;
115 nl->maxlen = 0;
116 nl->nri = 0;
117 nl->nrj = 0;
118 nl->iinr = NULL;
119 nl->gid = NULL;
120 nl->shift = NULL;
121 nl->jindex = NULL;
122 nl->jjnr = NULL;
123 /*nl->nltype = nltype;*/
124
125 srenew(nl->iinr, nl->maxnri);
126 srenew(nl->gid, nl->maxnri);
127 srenew(nl->shift, nl->maxnri);
128 srenew(nl->jindex, nl->maxnri+1);
129
130 nl->jindex[0] = 0;
131
132 return 0;
133 }
134
135 int print_nblist(int natoms, t_nblist *nl)
136 {
137 int i,k,ai,aj,nj0,nj1;
138
139 printf("genborn.c: print_nblist, natoms=%d\n",nl->nri);
140
141 for(i=0;i<nl->nri;i++)
142 {
143 ai=nl->iinr[i];
144 nj0=nl->jindex[i];
145 nj1=nl->jindex[i+1];
146
147 for(k=nj0;k<nj1;k++)
148 {
149 aj=nl->jjnr[k];
150 printf("ai=%d, aj=%d\n",ai,aj);
151 }
152 }
153
154 return 0;
155 }
156
157
158 void gb_pd_send(t_commrec *cr, real *send_data, int nr)
159 {
160 #ifdef GMX_MPI
161 int i,cur;
162 int *index,*sendc,*disp;
163
164 snew(sendc,cr->nnodes);
165 snew(disp,cr->nnodes);
166
167 index = pd_index(cr);
168 cur = cr->nodeid;
169
170 /* Setup count/index arrays */
171 for(i=0;i<cr->nnodes;i++)
172 {
173 sendc[i] = index[i+1]-index[i];
174 disp[i] = index[i];
175 }
176
177 /* Do communication */
178 MPI_Gatherv(send_data+index[cur],sendc[cur],GMX_MPI_REAL,send_data,sendc,
179 disp,GMX_MPI_REAL,0,cr->mpi_comm_mygroup);
180 MPI_Bcast(send_data,nr,GMX_MPI_REAL,0,cr->mpi_comm_mygroup);
181
182 #endif
183 }
184
185
186 int init_gb_plist(t_params *p_list)
187 {
188 p_list->nr = 0;
189 p_list->param = NULL;
190
191 return 0;
192 }
193
194
195
196 int init_gb_still(const t_commrec *cr, t_forcerec *fr,
197 const t_atomtypes *atype, t_idef *idef, t_atoms *atoms,
198 gmx_genborn_t *born,int natoms)
199 {
200
201 int i,j,i1,i2,k,m,nbond,nang,ia,ib,ic,id,nb,idx,idx2,at;
202 int iam,ibm;
203 int at0,at1;
204 real length,angle;
205 real r,ri,rj,ri2,ri3,rj2,r2,r3,r4,rk,ratio,term,h,doffset;
206 real p1,p2,p3,factor,cosine,rab,rbc;
207
208 real *vsol;
209 real *gp;
210
211 snew(vsol,natoms);
212 snew(gp,natoms);
213 snew(born->gpol_still_work,natoms+3);
214
215 if(PAR(cr))
216 {
217 if(PARTDECOMP(cr))
218 {
219 pd_at_range(cr,&at0,&at1);
220
221 for(i=0;i<natoms;i++)
222 {
223 vsol[i] = gp[i] = 0;
224 }
225 }
226 else
227 {
228 at0 = 0;
229 at1 = natoms;
230 }
231 }
232 else
233 {
234 at0 = 0;
235 at1 = natoms;
236 }
237
238 doffset = born->gb_doffset;
239
240 for(i=0;i<natoms;i++)
241 {
242 born->gpol_globalindex[i]=born->vsolv_globalindex[i]=
243 born->gb_radius_globalindex[i]=0;
244 }
245
246 /* Compute atomic solvation volumes for Still method */
247 for(i=0;i<natoms;i++)
248 {
249 ri=atype->gb_radius[atoms->atom[i].type];
250 born->gb_radius_globalindex[i] = ri;
251 r3=ri*ri*ri;
252 born->vsolv_globalindex[i]=(4*M_PI/3)*r3;
253 }
254
255 for(j=0;j<idef->il[F_GB12].nr;j+=3)
256 {
257 m=idef->il[F_GB12].iatoms[j];
258 ia=idef->il[F_GB12].iatoms[j+1];
259 ib=idef->il[F_GB12].iatoms[j+2];
260
261 r=1.01*idef->iparams[m].gb.st;
262
263 ri = atype->gb_radius[atoms->atom[ia].type];
264 rj = atype->gb_radius[atoms->atom[ib].type];
265
266 ri2 = ri*ri;
267 ri3 = ri2*ri;
268 rj2 = rj*rj;
269
270 ratio = (rj2-ri2-r*r)/(2*ri*r);
271 h = ri*(1+ratio);
272 term = (M_PI/3.0)*h*h*(3.0*ri-h);
273
274 if(PARTDECOMP(cr))
275 {
276 vsol[ia]+=term;
277 }
278 else
279 {
280 born->vsolv_globalindex[ia] -= term;
281 }
282
283 ratio = (ri2-rj2-r*r)/(2*rj*r);
284 h = rj*(1+ratio);
285 term = (M_PI/3.0)*h*h*(3.0*rj-h);
286
287 if(PARTDECOMP(cr))
288 {
289 vsol[ib]+=term;
290 }
291 else
292 {
293 born->vsolv_globalindex[ib] -= term;
294 }
295 }
296
297 if(PARTDECOMP(cr))
298 {
299 gmx_sum(natoms,vsol,cr);
300
301 for(i=0;i<natoms;i++)
302 {
303 born->vsolv_globalindex[i]=born->vsolv_globalindex[i]-vsol[i];
304 }
305 }
306
307 /* Get the self-, 1-2 and 1-3 polarization energies for analytical Still
308 method */
309 /* Self */
310 for(j=0;j<natoms;j++)
311 {
312 if(born->use_globalindex[j]==1)
313 {
314 born->gpol_globalindex[j]=-0.5*ONE_4PI_EPS0/
315 (atype->gb_radius[atoms->atom[j].type]-doffset+STILL_P1);
316 }
317 }
318
319 /* 1-2 */
320 for(j=0;j<idef->il[F_GB12].nr;j+=3)
321 {
322 m=idef->il[F_GB12].iatoms[j];
323 ia=idef->il[F_GB12].iatoms[j+1];
324 ib=idef->il[F_GB12].iatoms[j+2];
325
326 r=idef->iparams[m].gb.st;
327
328 r4=r*r*r*r;
329
330 if(PARTDECOMP(cr))
331 {
332 gp[ia]+=STILL_P2*born->vsolv_globalindex[ib]/r4;
333 gp[ib]+=STILL_P2*born->vsolv_globalindex[ia]/r4;
334 }
335 else
336 {
337 born->gpol_globalindex[ia]=born->gpol_globalindex[ia]+
338 STILL_P2*born->vsolv_globalindex[ib]/r4;
339 born->gpol_globalindex[ib]=born->gpol_globalindex[ib]+
340 STILL_P2*born->vsolv_globalindex[ia]/r4;
341 }
342 }
343
344 /* 1-3 */
345 for(j=0;j<idef->il[F_GB13].nr;j+=3)
346 {
347 m=idef->il[F_GB13].iatoms[j];
348 ia=idef->il[F_GB13].iatoms[j+1];
349 ib=idef->il[F_GB13].iatoms[j+2];
350
351 r=idef->iparams[m].gb.st;
352 r4=r*r*r*r;
353
354 if(PARTDECOMP(cr))
355 {
356 gp[ia]+=STILL_P3*born->vsolv[ib]/r4;
357 gp[ib]+=STILL_P3*born->vsolv[ia]/r4;
358 }
359 else
360 {
361 born->gpol_globalindex[ia]=born->gpol_globalindex[ia]+
362 STILL_P3*born->vsolv_globalindex[ib]/r4;
363 born->gpol_globalindex[ib]=born->gpol_globalindex[ib]+
364 STILL_P3*born->vsolv_globalindex[ia]/r4;
365 }
366 }
367
368 if(PARTDECOMP(cr))
369 {
370 gmx_sum(natoms,gp,cr);
371
372 for(i=0;i<natoms;i++)
373 {
374 born->gpol_globalindex[i]=born->gpol_globalindex[i]+gp[i];
375 }
376 }
377
378 sfree(vsol);
379 sfree(gp);
380
381 return 0;
382 }
383
384 /* Initialize all GB datastructs and compute polarization energies */
385 int init_gb(gmx_genborn_t **p_born,
386 const t_commrec *cr, t_forcerec *fr, const t_inputrec *ir,
387 const gmx_mtop_t *mtop, real rgbradii, int gb_algorithm)
388 {
389 int i,j,m,ai,aj,jj,natoms,nalloc;
390 real rai,sk,p,doffset;
391
392 t_atoms atoms;
393 gmx_genborn_t *born;
394 gmx_localtop_t *localtop;
395
396 natoms = mtop->natoms;
397
398 atoms = gmx_mtop_global_atoms(mtop);
399 localtop = gmx_mtop_generate_local_top(mtop,ir);
400
401 snew(born,1);
402 *p_born = born;
403
404 born->nr = natoms;
405
406 snew(born->drobc, natoms);
407 snew(born->bRad, natoms);
408
409 /* Allocate memory for the global data arrays */
410 snew(born->param_globalindex, natoms+3);
411 snew(born->gpol_globalindex, natoms+3);
412 snew(born->vsolv_globalindex, natoms+3);
413 snew(born->gb_radius_globalindex, natoms+3);
414 snew(born->use_globalindex, natoms+3);
415
416 snew(fr->invsqrta, natoms);
417 snew(fr->dvda, natoms);
418
419 fr->dadx = NULL;
420 fr->dadx_rawptr = NULL;
421 fr->nalloc_dadx = 0;
422 born->gpol_still_work = NULL;
423 born->gpol_hct_work = NULL;
424
425 /* snew(born->asurf,natoms); */
426 /* snew(born->dasurf,natoms); */
427
428 /* Initialize the gb neighbourlist */
429 init_gb_nblist(natoms,&(fr->gblist));
430
431 /* Do the Vsites exclusions (if any) */
432 for(i=0;i<natoms;i++)
433 {
434 jj = atoms.atom[i].type;
435 if (mtop->atomtypes.gb_radius[atoms.atom[i].type] > 0)
436 {
437 born->use_globalindex[i] = 1;
438 }
439 else
440 {
441 born->use_globalindex[i] = 0;
442 }
443
444 /* If we have a Vsite, put vs_globalindex[i]=0 */
445 if (C6 (fr->nbfp,fr->ntype,jj,jj) == 0 &&
446 C12(fr->nbfp,fr->ntype,jj,jj) == 0 &&
447 atoms.atom[i].q == 0)
448 {
449 born->use_globalindex[i]=0;
450 }
451 }
452
453 /* Copy algorithm parameters from inputrecord to local structure */
454 born->obc_alpha = ir->gb_obc_alpha;
455 born->obc_beta = ir->gb_obc_beta;
456 born->obc_gamma = ir->gb_obc_gamma;
457 born->gb_doffset = ir->gb_dielectric_offset;
458 born->gb_epsilon_solvent = ir->gb_epsilon_solvent;
459 born->epsilon_r = ir->epsilon_r;
460
461 doffset = born->gb_doffset;
462
463 /* Set the surface tension */
464 born->sa_surface_tension = ir->sa_surface_tension;
465
466 /* If Still model, initialise the polarisation energies */
467 if(gb_algorithm==egbSTILL)
468 {
469 init_gb_still(cr, fr,&(mtop->atomtypes), &(localtop->idef), &atoms,
470 born, natoms);
471 }
472
473
474 /* If HCT/OBC, precalculate the sk*atype->S_hct factors */
475 else if(gb_algorithm==egbHCT || gb_algorithm==egbOBC)
476 {
477
478 snew(born->gpol_hct_work, natoms+3);
479
480 for(i=0;i<natoms;i++)
481 {
482 if(born->use_globalindex[i]==1)
483 {
484 rai = mtop->atomtypes.gb_radius[atoms.atom[i].type]-doffset;
485 sk = rai * mtop->atomtypes.S_hct[atoms.atom[i].type];
486 born->param_globalindex[i] = sk;
487 born->gb_radius_globalindex[i] = rai;
488 }
489 else
490 {
491 born->param_globalindex[i] = 0;
492 born->gb_radius_globalindex[i] = 0;
493 }
494 }
495 }
496
497 /* Allocate memory for work arrays for temporary use */
498 snew(born->work,natoms+4);
499 snew(born->count,natoms);
500 snew(born->nblist_work,natoms);
501
502 /* Domain decomposition specific stuff */
503 born->nalloc = 0;
504
505 return 0;
506 }
507
508
509
510 static int
511 calc_gb_rad_still(t_commrec *cr, t_forcerec *fr,int natoms, gmx_localtop_t *top,
512 const t_atomtypes *atype, rvec x[], t_nblist *nl,
513 gmx_genborn_t *born,t_mdatoms *md)
514 {
515 int i,k,n,nj0,nj1,ai,aj,type;
516 int shift;
517 real shX,shY,shZ;
518 real gpi,dr,dr2,dr4,idr4,rvdw,ratio,ccf,theta,term,rai,raj;
519 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
520 real rinv,idr2,idr6,vaj,dccf,cosq,sinq,prod,gpi2;
521 real factor;
522 real vai, prod_ai, icf4,icf6;
523
524 factor = 0.5*ONE_4PI_EPS0;
525 n = 0;
526
527 for(i=0;i<born->nr;i++)
528 {
529 born->gpol_still_work[i]=0;
530 }
531
532 for(i=0;i<nl->nri;i++ )
533 {
534 ai = nl->iinr[i];
535
536 nj0 = nl->jindex[i];
537 nj1 = nl->jindex[i+1];
538
539 /* Load shifts for this list */
540 shift = nl->shift[i];
541 shX = fr->shift_vec[shift][0];
542 shY = fr->shift_vec[shift][1];
543 shZ = fr->shift_vec[shift][2];
544
545 gpi = 0;
546
547 rai = top->atomtypes.gb_radius[md->typeA[ai]];
548 vai = born->vsolv[ai];
549 prod_ai = STILL_P4*vai;
550
551 /* Load atom i coordinates, add shift vectors */
552 ix1 = shX + x[ai][0];
553 iy1 = shY + x[ai][1];
554 iz1 = shZ + x[ai][2];
555
556 for(k=nj0;k<nj1;k++)
557 {
558 aj = nl->jjnr[k];
559 jx1 = x[aj][0];
560 jy1 = x[aj][1];
561 jz1 = x[aj][2];
562
563 dx11 = ix1-jx1;
564 dy11 = iy1-jy1;
565 dz11 = iz1-jz1;
566
567 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
568 rinv = gmx_invsqrt(dr2);
569 idr2 = rinv*rinv;
570 idr4 = idr2*idr2;
571 idr6 = idr4*idr2;
572
573 raj = top->atomtypes.gb_radius[md->typeA[aj]];
574
575 rvdw = rai + raj;
576
577 ratio = dr2 / (rvdw * rvdw);
578 vaj = born->vsolv[aj];
579
580 if(ratio>STILL_P5INV)
581 {
582 ccf=1.0;
583 dccf=0.0;
584 }
585 else
586 {
587 theta = ratio*STILL_PIP5;
588 cosq = cos(theta);
589 term = 0.5*(1.0-cosq);
590 ccf = term*term;
591 sinq = 1.0 - cosq*cosq;
592 dccf = 2.0*term*sinq*gmx_invsqrt(sinq)*theta;
593 }
594
595 prod = STILL_P4*vaj;
596 icf4 = ccf*idr4;
597 icf6 = (4*ccf-dccf)*idr6;
598
599 born->gpol_still_work[aj] += prod_ai*icf4;
600 gpi = gpi+prod*icf4;
601
602 /* Save ai->aj and aj->ai chain rule terms */
603 fr->dadx[n++] = prod*icf6;
604 fr->dadx[n++] = prod_ai*icf6;
605 }
606 born->gpol_still_work[ai]+=gpi;
607 }
608
609 /* Parallel summations */
610 if(PARTDECOMP(cr))
611 {
612 gmx_sum(natoms, born->gpol_still_work, cr);
613 }
614 else if(DOMAINDECOMP(cr))
615 {
616 dd_atom_sum_real(cr->dd, born->gpol_still_work);
617 }
618
619 /* Calculate the radii */
620 for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
621 {
622 if(born->use[i] != 0)
623 {
624
625 gpi = born->gpol[i]+born->gpol_still_work[i];
626 gpi2 = gpi * gpi;
627 born->bRad[i] = factor*gmx_invsqrt(gpi2);
628 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
629 }
630 }
631
632 /* Extra communication required for DD */
633 if(DOMAINDECOMP(cr))
634 {
635 dd_atom_spread_real(cr->dd, born->bRad);
636 dd_atom_spread_real(cr->dd, fr->invsqrta);
637 }
638
639 return 0;
640
641 }
642
643
644 static int
645 calc_gb_rad_hct(t_commrec *cr,t_forcerec *fr,int natoms, gmx_localtop_t *top,
646 const t_atomtypes *atype, rvec x[], t_nblist *nl,
647 gmx_genborn_t *born,t_mdatoms *md)
648 {
649 int i,k,n,ai,aj,nj0,nj1,at0,at1;
650 int shift;
651 real shX,shY,shZ;
652 real rai,raj,gpi,dr2,dr,sk,sk_ai,sk2,sk2_ai,lij,uij,diff2,tmp,sum_ai;
653 real rad,min_rad,rinv,rai_inv;
654 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
655 real lij2, uij2, lij3, uij3, t1,t2,t3;
656 real lij_inv,dlij,duij,sk2_rinv,prod,log_term;
657 real doffset,raj_inv,dadx_val;
658 real *gb_radius;
659
660 doffset = born->gb_doffset;
661 gb_radius = born->gb_radius;
662
663 for(i=0;i<born->nr;i++)
664 {
665 born->gpol_hct_work[i] = 0;
666 }
667
668 /* Keep the compiler happy */
669 n = 0;
670 prod = 0;
671
672 for(i=0;i<nl->nri;i++)
673 {
674 ai = nl->iinr[i];
675
676 nj0 = nl->jindex[i];
677 nj1 = nl->jindex[i+1];
678
679 /* Load shifts for this list */
680 shift = nl->shift[i];
681 shX = fr->shift_vec[shift][0];
682 shY = fr->shift_vec[shift][1];
683 shZ = fr->shift_vec[shift][2];
684
685 rai = gb_radius[ai];
686 rai_inv = 1.0/rai;
687
688 sk_ai = born->param[ai];
689 sk2_ai = sk_ai*sk_ai;
690
691 /* Load atom i coordinates, add shift vectors */
692 ix1 = shX + x[ai][0];
693 iy1 = shY + x[ai][1];
694 iz1 = shZ + x[ai][2];
695
696 sum_ai = 0;
697
698 for(k=nj0;k<nj1;k++)
699 {
700 aj = nl->jjnr[k];
701
702 jx1 = x[aj][0];
703 jy1 = x[aj][1];
704 jz1 = x[aj][2];
705
706 dx11 = ix1 - jx1;
707 dy11 = iy1 - jy1;
708 dz11 = iz1 - jz1;
709
710 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
711 rinv = gmx_invsqrt(dr2);
712 dr = rinv*dr2;
713
714 sk = born->param[aj];
715 raj = gb_radius[aj];
716
717 /* aj -> ai interaction */
718 if(rai < dr+sk)
719 {
720 lij = 1.0/(dr-sk);
721 dlij = 1.0;
722
723 if(rai>dr-sk)
724 {
725 lij = rai_inv;
726 dlij = 0.0;
727 }
728
729 lij2 = lij*lij;
730 lij3 = lij2*lij;
731
732 uij = 1.0/(dr+sk);
733 uij2 = uij*uij;
734 uij3 = uij2*uij;
735
736 diff2 = uij2-lij2;
737
738 lij_inv = gmx_invsqrt(lij2);
739 sk2 = sk*sk;
740 sk2_rinv = sk2*rinv;
741 prod = 0.25*sk2_rinv;
742
743 log_term = log(uij*lij_inv);
744
745 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
746 prod*(-diff2);
747
748 if(rai<sk-dr)
749 {
750 tmp = tmp + 2.0 * (rai_inv-lij);
751 }
752
753 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
754 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
755 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
756
757 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
758 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */
759 /* rb2 is moved to chainrule */
760
761 sum_ai += 0.5*tmp;
762 }
763 else
764 {
765 dadx_val = 0.0;
766 }
767 fr->dadx[n++] = dadx_val;
768
769
770 /* ai -> aj interaction */
771 if(raj < dr + sk_ai)
772 {
773 lij = 1.0/(dr-sk_ai);
774 dlij = 1.0;
775 raj_inv = 1.0/raj;
776
777 if(raj>dr-sk_ai)
778 {
779 lij = raj_inv;
780 dlij = 0.0;
781 }
782
783 lij2 = lij * lij;
784 lij3 = lij2 * lij;
785
786 uij = 1.0/(dr+sk_ai);
787 uij2 = uij * uij;
788 uij3 = uij2 * uij;
789
790 diff2 = uij2-lij2;
791
792 lij_inv = gmx_invsqrt(lij2);
793 sk2 = sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
794 sk2_rinv = sk2*rinv;
795 prod = 0.25 * sk2_rinv;
796
797 /* log_term = table_log(uij*lij_inv,born->log_table,
798 LOG_TABLE_ACCURACY); */
799 log_term = log(uij*lij_inv);
800
801 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
802 prod*(-diff2);
803
804 if(raj<sk_ai-dr)
805 {
806 tmp = tmp + 2.0 * (raj_inv-lij);
807 }
808
809 /* duij = 1.0 */
810 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
811 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
812 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
813
814 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
815 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */ /* rb2 is moved to chainrule */
816
817 born->gpol_hct_work[aj] += 0.5*tmp;
818 }
819 else
820 {
821 dadx_val = 0.0;
822 }
823 fr->dadx[n++] = dadx_val;
824 }
825
826 born->gpol_hct_work[ai] += sum_ai;
827 }
828
829 /* Parallel summations */
830 if(PARTDECOMP(cr))
831 {
832 gmx_sum(natoms, born->gpol_hct_work, cr);
833 }
834 else if(DOMAINDECOMP(cr))
835 {
836 dd_atom_sum_real(cr->dd, born->gpol_hct_work);
837 }
838
839 for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
840 {
841 if(born->use[i] != 0)
842 {
843 rai = top->atomtypes.gb_radius[md->typeA[i]]-doffset;
844 sum_ai = 1.0/rai - born->gpol_hct_work[i];
845 min_rad = rai + doffset;
846 rad = 1.0/sum_ai;
847
848 born->bRad[i] = rad > min_rad ? rad : min_rad;
849 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
850 }
851 }
852
853 /* Extra communication required for DD */
854 if(DOMAINDECOMP(cr))
855 {
856 dd_atom_spread_real(cr->dd, born->bRad);
857 dd_atom_spread_real(cr->dd, fr->invsqrta);
858 }
859
860
861 return 0;
862 }
863
864 static int
865 calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, int natoms, gmx_localtop_t *top,
866 const t_atomtypes *atype, rvec x[], t_nblist *nl, gmx_genborn_t *born,t_mdatoms *md)
867 {
868 int i,k,ai,aj,nj0,nj1,n,at0,at1;
869 int shift;
870 real shX,shY,shZ;
871 real rai,raj,gpi,dr2,dr,sk,sk2,lij,uij,diff2,tmp,sum_ai;
872 real rad, min_rad,sum_ai2,sum_ai3,tsum,tchain,rinv,rai_inv,lij_inv,rai_inv2;
873 real log_term,prod,sk2_rinv,sk_ai,sk2_ai;
874 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
875 real lij2,uij2,lij3,uij3,dlij,duij,t1,t2,t3;
876 real doffset,raj_inv,dadx_val;
877 real *gb_radius;
878
879 /* Keep the compiler happy */
880 n = 0;
881 prod = 0;
882 raj = 0;
883
884 doffset = born->gb_doffset;
885 gb_radius = born->gb_radius;
886
887 for(i=0;i<born->nr;i++)
888 {
889 born->gpol_hct_work[i] = 0;
890 }
891
892 for(i=0;i<nl->nri;i++)
893 {
894 ai = nl->iinr[i];
895
896 nj0 = nl->jindex[i];
897 nj1 = nl->jindex[i+1];
898
899 /* Load shifts for this list */
900 shift = nl->shift[i];
901 shX = fr->shift_vec[shift][0];
902 shY = fr->shift_vec[shift][1];
903 shZ = fr->shift_vec[shift][2];
904
905 rai = gb_radius[ai];
906 rai_inv = 1.0/rai;
907
908 sk_ai = born->param[ai];
909 sk2_ai = sk_ai*sk_ai;
910
911 /* Load atom i coordinates, add shift vectors */
912 ix1 = shX + x[ai][0];
913 iy1 = shY + x[ai][1];
914 iz1 = shZ + x[ai][2];
915
916 sum_ai = 0;
917
918 for(k=nj0;k<nj1;k++)
919 {
920 aj = nl->jjnr[k];
921
922 jx1 = x[aj][0];
923 jy1 = x[aj][1];
924 jz1 = x[aj][2];
925
926 dx11 = ix1 - jx1;
927 dy11 = iy1 - jy1;
928 dz11 = iz1 - jz1;
929
930 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
931 rinv = gmx_invsqrt(dr2);
932 dr = dr2*rinv;
933
934 /* sk is precalculated in init_gb() */
935 sk = born->param[aj];
936 raj = gb_radius[aj];
937
938 /* aj -> ai interaction */
939 if(rai < dr+sk)
940 {
941 lij = 1.0/(dr-sk);
942 dlij = 1.0;
943
944 if(rai>dr-sk)
945 {
946 lij = rai_inv;
947 dlij = 0.0;
948 }
949
950 uij = 1.0/(dr+sk);
951 lij2 = lij * lij;
952 lij3 = lij2 * lij;
953 uij2 = uij * uij;
954 uij3 = uij2 * uij;
955
956 diff2 = uij2-lij2;
957
958 lij_inv = gmx_invsqrt(lij2);
959 sk2 = sk*sk;
960 sk2_rinv = sk2*rinv;
961 prod = 0.25*sk2_rinv;
962
963 log_term = log(uij*lij_inv);
964
965 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
966
967 if(rai < sk-dr)
968 {
969 tmp = tmp + 2.0 * (rai_inv-lij);
970 }
971
972 /* duij = 1.0; */
973 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
974 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
975 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
976
977 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
978
979 sum_ai += 0.5*tmp;
980 }
981 else
982 {
983 dadx_val = 0.0;
984 }
985 fr->dadx[n++] = dadx_val;
986
987 /* ai -> aj interaction */
988 if(raj < dr + sk_ai)
989 {
990 lij = 1.0/(dr-sk_ai);
991 dlij = 1.0;
992 raj_inv = 1.0/raj;
993
994 if(raj>dr-sk_ai)
995 {
996 lij = raj_inv;
997 dlij = 0.0;
998 }
999
1000 lij2 = lij * lij;
1001 lij3 = lij2 * lij;
1002
1003 uij = 1.0/(dr+sk_ai);
1004 uij2 = uij * uij;
1005 uij3 = uij2 * uij;
1006
1007 diff2 = uij2-lij2;
1008
1009 lij_inv = gmx_invsqrt(lij2);
1010 sk2 = sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
1011 sk2_rinv = sk2*rinv;
1012 prod = 0.25 * sk2_rinv;
1013
1014 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
1015 log_term = log(uij*lij_inv);
1016
1017 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
1018
1019 if(raj<sk_ai-dr)
1020 {
1021 tmp = tmp + 2.0 * (raj_inv-lij);
1022 }
1023
1024 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
1025 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
1026 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
1027
1028 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
1029
1030 born->gpol_hct_work[aj] += 0.5*tmp;
1031
1032 }
1033 else
1034 {
1035 dadx_val = 0.0;
1036 }
1037 fr->dadx[n++] = dadx_val;
1038
1039 }
1040 born->gpol_hct_work[ai] += sum_ai;
1041
1042 }
1043
1044 /* Parallel summations */
1045 if(PARTDECOMP(cr))
1046 {
1047 gmx_sum(natoms, born->gpol_hct_work, cr);
1048 }
1049 else if(DOMAINDECOMP(cr))
1050 {
1051 dd_atom_sum_real(cr->dd, born->gpol_hct_work);
1052 }
1053
1054 for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
1055 {
1056 if(born->use[i] != 0)
1057 {
1058 rai = top->atomtypes.gb_radius[md->typeA[i]];
1059 rai_inv2 = 1.0/rai;
1060 rai = rai-doffset;
1061 rai_inv = 1.0/rai;
1062 sum_ai = rai * born->gpol_hct_work[i];
1063 sum_ai2 = sum_ai * sum_ai;
1064 sum_ai3 = sum_ai2 * sum_ai;
1065
1066 tsum = tanh(born->obc_alpha*sum_ai-born->obc_beta*sum_ai2+born->obc_gamma*sum_ai3);
1067 born->bRad[i] = rai_inv - tsum*rai_inv2;
1068 born->bRad[i] = 1.0 / born->bRad[i];
1069
1070 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
1071
1072 tchain = rai * (born->obc_alpha-2*born->obc_beta*sum_ai+3*born->obc_gamma*sum_ai2);
1073 born->drobc[i] = (1.0-tsum*tsum)*tchain*rai_inv2;
1074 }
1075 }
1076
1077 /* Extra (local) communication required for DD */
1078 if(DOMAINDECOMP(cr))
1079 {
1080 dd_atom_spread_real(cr->dd, born->bRad);
1081 dd_atom_spread_real(cr->dd, fr->invsqrta);
1082 dd_atom_spread_real(cr->dd, born->drobc);
1083 }
1084
1085 return 0;
1086
1087 }
1088
1089
1090
1091 int calc_gb_rad(t_commrec *cr, t_forcerec *fr, t_inputrec *ir,gmx_localtop_t *top,
1092 const t_atomtypes *atype, rvec x[], t_nblist *nl, gmx_genborn_t *born,t_mdatoms *md,t_nrnb *nrnb)
1093 {
1094 real *p;
1095 int cnt;
1096 int ndadx;
1097
1098 if(fr->bAllvsAll && fr->dadx==NULL)
1099 {
1100 /* We might need up to 8 atoms of padding before and after,
1101 * and another 4 units to guarantee SSE alignment.
1102 */
1103 fr->nalloc_dadx = 2*(md->homenr+12)*(md->nr/2+1+12);
1104 snew(fr->dadx_rawptr,fr->nalloc_dadx);
1105 fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
1106 }
1107 else
1108 {
1109 /* In the SSE-enabled gb-loops, when writing to dadx, we
1110 * always write 2*4 elements at a time, even in the case with only
1111 * 1-3 j particles, where we only really need to write 2*(1-3)
1112 * elements. This is because we want dadx to be aligned to a 16-
1113 * byte boundary, and being able to use _mm_store/load_ps
1114 */
1115 ndadx = 2 * (nl->nrj + 3*nl->nri);
1116
1117 /* First, reallocate the dadx array, we need 3 extra for SSE */
1118 if (ndadx + 3 > fr->nalloc_dadx)
1119 {
1120 fr->nalloc_dadx = over_alloc_large(ndadx) + 3;
1121 srenew(fr->dadx_rawptr,fr->nalloc_dadx);
1122 fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
1123 }
1124 }
1125
1126 if(fr->bAllvsAll)
1127 {
1128 cnt = md->homenr*(md->nr/2+1);
1129
1130 if(ir->gb_algorithm==egbSTILL)
1131 {
1132 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || (!defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1133 if(fr->UseOptimizedKernels)
1134 {
1135 genborn_allvsall_calc_still_radii_sse2_single(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1136 }
1137 else
1138 {
1139 genborn_allvsall_calc_still_radii(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1140 }
1141 #elif ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || (defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1142 if(fr->UseOptimizedKernels)
1143 {
1144 genborn_allvsall_calc_still_radii_sse2_double(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1145 }
1146 else
1147 {
1148 genborn_allvsall_calc_still_radii(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1149 }
1150 #else
1151 genborn_allvsall_calc_still_radii(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1152 #endif
1153 inc_nrnb(nrnb,eNR_BORN_AVA_RADII_STILL,cnt);
1154 }
1155 else if(ir->gb_algorithm==egbHCT || ir->gb_algorithm==egbOBC)
1156 {
1157 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || (!defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1158 if(fr->UseOptimizedKernels)
1159 {
1160 genborn_allvsall_calc_hct_obc_radii_sse2_single(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1161 }
1162 else
1163 {
1164 genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1165 }
1166 #elif ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || (defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1167 if(fr->UseOptimizedKernels)
1168 {
1169 genborn_allvsall_calc_hct_obc_radii_sse2_double(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1170 }
1171 else
1172 {
1173 genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1174 }
1175 #else
1176 genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1177 #endif
1178 inc_nrnb(nrnb,eNR_BORN_AVA_RADII_HCT_OBC,cnt);
1179 }
1180 else
1181 {
1182 gmx_fatal(FARGS,"Bad gb algorithm for all-vs-all interactions");
1183 }
1184 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,md->homenr);
1185
1186 return 0;
1187 }
1188
1189 /* Switch for determining which algorithm to use for Born radii calculation */
1190 #ifdef GMX_DOUBLE
1191
1192 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_SSE2) )
1193 /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1194 switch(ir->gb_algorithm)
1195 {
1196 case egbSTILL:
1197 if(fr->UseOptimizedKernels)
1198 {
1199 calc_gb_rad_still_sse2_double(cr,fr,md->nr,top, atype, x[0], nl, born);
1200 }
1201 else
1202 {
1203 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1204 }
1205 break;
1206 case egbHCT:
1207 if(fr->UseOptimizedKernels)
1208 {
1209 calc_gb_rad_hct_obc_sse2_double(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1210 }
1211 else
1212 {
1213 calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1214 }
1215 break;
1216 case egbOBC:
1217 if(fr->UseOptimizedKernels)
1218 {
1219 calc_gb_rad_hct_obc_sse2_double(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1220 }
1221 else
1222 {
1223 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1224 }
1225 break;
1226
1227 default:
1228 gmx_fatal(FARGS, "Unknown double precision sse-enabled algorithm for Born radii calculation: %d",ir->gb_algorithm);
1229 }
1230 #else
1231 switch(ir->gb_algorithm)
1232 {
1233 case egbSTILL:
1234 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1235 break;
1236 case egbHCT:
1237 calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1238 break;
1239 case egbOBC:
1240 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1241 break;
1242
1243 default:
1244 gmx_fatal(FARGS, "Unknown double precision algorithm for Born radii calculation: %d",ir->gb_algorithm);
1245 }
1246
1247 #endif
1248
1249 #else
1250
1251 #if (!defined DISABLE_SSE && ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_SSE2) ) )
1252 /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1253 switch(ir->gb_algorithm)
1254 {
1255 case egbSTILL:
1256 if(fr->UseOptimizedKernels)
1257 {
1258 calc_gb_rad_still_sse2_single(cr,fr,born->nr,top, atype, x[0], nl, born);
1259 }
1260 else
1261 {
1262 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1263 }
1264 break;
1265 case egbHCT:
1266 if(fr->UseOptimizedKernels)
1267 {
1268 calc_gb_rad_hct_obc_sse2_single(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1269 }
1270 else
1271 {
1272 calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1273 }
1274 break;
1275
1276 case egbOBC:
1277 if(fr->UseOptimizedKernels)
1278 {
1279 calc_gb_rad_hct_obc_sse2_single(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1280 }
1281 else
1282 {
1283 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1284 }
1285 break;
1286
1287 default:
1288 gmx_fatal(FARGS, "Unknown sse-enabled algorithm for Born radii calculation: %d",ir->gb_algorithm);
1289 }
1290
1291 #else
1292 switch(ir->gb_algorithm)
1293 {
1294 case egbSTILL:
1295 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1296 break;
1297 case egbHCT:
1298 calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1299 break;
1300 case egbOBC:
1301 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1302 break;
1303
1304 default:
1305 gmx_fatal(FARGS, "Unknown algorithm for Born radii calculation: %d",ir->gb_algorithm);
1306 }
1307
1308 #endif /* Single precision sse */
1309
1310 #endif /* Double or single precision */
1311
1312 if(fr->bAllvsAll==FALSE)
1313 {
1314 switch(ir->gb_algorithm)
1315 {
1316 case egbSTILL:
1317 inc_nrnb(nrnb,eNR_BORN_RADII_STILL,nl->nrj);
1318 break;
1319 case egbHCT:
1320 case egbOBC:
1321 inc_nrnb(nrnb,eNR_BORN_RADII_HCT_OBC,nl->nrj);
1322 break;
1323
1324 default:
1325 break;
1326 }
1327 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,nl->nri);
1328 }
1329
1330 return 0;
1331 }
1332
1333
1334
1335 real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1336 real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
1337 real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph)
1338 {
1339 int i,j,n0,m,nnn,type,ai,aj;
1340 int ki;
1341
1342 real isai,isaj;
1343 real r,rsq11;
1344 real rinv11,iq;
1345 real isaprod,qq,gbscale,gbtabscale,Y,F,Geps,Heps2,Fp,VV,FF,rt,eps,eps2;
1346 real vgb,fgb,vcoul,fijC,dvdatmp,fscal,dvdaj;
1347 real vctot;
1348
1349 rvec dx;
1350 ivec dt;
1351
1352 t_iatom *forceatoms;
1353
1354 /* Scale the electrostatics by gb_epsilon_solvent */
1355 facel = facel * ((1.0/epsilon_r) - 1.0/gb_epsilon_solvent);
1356
1357 gbtabscale=*p_gbtabscale;
1358 vctot = 0.0;
1359
1360 for(j=F_GB12;j<=F_GB14;j++)
1361 {
1362 forceatoms = idef->il[j].iatoms;
1363
1364 for(i=0;i<idef->il[j].nr; )
1365 {
1366 /* To avoid reading in the interaction type, we just increment i to pass over
1367 * the types in the forceatoms array, this saves some memory accesses
1368 */
1369 i++;
1370 ai = forceatoms[i++];
1371 aj = forceatoms[i++];
1372
1373 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx);
1374 rsq11 = iprod(dx,dx);
1375
1376 isai = invsqrta[ai];
1377 iq = (-1)*facel*charge[ai];
1378
1379 rinv11 = gmx_invsqrt(rsq11);
1380 isaj = invsqrta[aj];
1381 isaprod = isai*isaj;
1382 qq = isaprod*iq*charge[aj];
1383 gbscale = isaprod*gbtabscale;
1384 r = rsq11*rinv11;
1385 rt = r*gbscale;
1386 n0 = rt;
1387 eps = rt-n0;
1388 eps2 = eps*eps;
1389 nnn = 4*n0;
1390 Y = GBtab[nnn];
1391 F = GBtab[nnn+1];
1392 Geps = eps*GBtab[nnn+2];
1393 Heps2 = eps2*GBtab[nnn+3];
1394 Fp = F+Geps+Heps2;
1395 VV = Y+eps*Fp;
1396 FF = Fp+Geps+2.0*Heps2;
1397 vgb = qq*VV;
1398 fijC = qq*FF*gbscale;
1399 dvdatmp = -(vgb+fijC*r)*0.5;
1400 dvda[aj] = dvda[aj] + dvdatmp*isaj*isaj;
1401 dvda[ai] = dvda[ai] + dvdatmp*isai*isai;
1402 vctot = vctot + vgb;
1403 fgb = -(fijC)*rinv11;
1404
1405 if (graph) {
1406 ivec_sub(SHIFT_IVEC(graph,ai),SHIFT_IVEC(graph,aj),dt);
1407 ki=IVEC2IS(dt);
1408 }
1409
1410 for (m=0; (m<DIM); m++) { /* 15 */
1411 fscal=fgb*dx[m];
1412 f[ai][m]+=fscal;
1413 f[aj][m]-=fscal;
1414 fshift[ki][m]+=fscal;
1415 fshift[CENTRAL][m]-=fscal;
1416 }
1417 }
1418 }
1419
1420 return vctot;
1421 }
1422
1423 real calc_gb_selfcorrections(t_commrec *cr, int natoms,
1424 real *charge, gmx_genborn_t *born, real *dvda, t_mdatoms *md, double facel)
1425 {
1426 int i,ai,at0,at1;
1427 real rai,e,derb,q,q2,fi,rai_inv,vtot;
1428
1429 if(PARTDECOMP(cr))
1430 {
1431 pd_at_range(cr,&at0,&at1);
1432 }
1433 else if(DOMAINDECOMP(cr))
1434 {
1435 at0=0;
1436 at1=cr->dd->nat_home;
1437 }
1438 else
1439 {
1440 at0=0;
1441 at1=natoms;
1442
1443 }
1444
1445 /* Scale the electrostatics by gb_epsilon_solvent */
1446 facel = facel * ((1.0/born->epsilon_r) - 1.0/born->gb_epsilon_solvent);
1447
1448 vtot=0.0;
1449
1450 /* Apply self corrections */
1451 for(i=at0;i<at1;i++)
1452 {
1453 ai = i;
1454
1455 if(born->use[ai]==1)
1456 {
1457 rai = born->bRad[ai];
1458 rai_inv = 1.0/rai;
1459 q = charge[ai];
1460 q2 = q*q;
1461 fi = facel*q2;
1462 e = fi*rai_inv;
1463 derb = 0.5*e*rai_inv*rai_inv;
1464 dvda[ai] += derb*rai;
1465 vtot -= 0.5*e;
1466 }
1467 }
1468
1469 return vtot;
1470
1471 }
1472
1473 real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr,int natoms,gmx_genborn_t *born, gmx_localtop_t *top,
1474 const t_atomtypes *atype, real *dvda,int gb_algorithm, t_mdatoms *md)
1475 {
1476 int ai,i,at0,at1;
1477 real e,es,rai,rbi,term,probe,tmp,factor;
1478 real rbi_inv,rbi_inv2;
1479
1480 /* To keep the compiler happy */
1481 factor=0;
1482
1483 if(PARTDECOMP(cr))
1484 {
1485 pd_at_range(cr,&at0,&at1);
1486 }
1487 else if(DOMAINDECOMP(cr))
1488 {
1489 at0 = 0;
1490 at1 = cr->dd->nat_home;
1491 }
1492 else
1493 {
1494 at0=0;
1495 at1=natoms;
1496 }
1497
1498 /* factor is the surface tension */
1499 factor = born->sa_surface_tension;
1500 /*
1501
1502 // The surface tension factor is 0.0049 for Still model, 0.0054 for HCT/OBC
1503 if(gb_algorithm==egbSTILL)
1504 {
1505 factor=0.0049*100*CAL2JOULE;
1506 }
1507 else
1508 {
1509 factor=0.0054*100*CAL2JOULE;
1510 }
1511 */
1512 /* if(gb_algorithm==egbHCT || gb_algorithm==egbOBC) */
1513
1514 es = 0;
1515 probe = 0.14;
1516 term = M_PI*4;
1517
1518 for(i=at0;i<at1;i++)
1519 {
1520 ai = i;
1521
1522 if(born->use[ai]==1)
1523 {
1524 rai = top->atomtypes.gb_radius[md->typeA[ai]];
1525 rbi_inv = fr->invsqrta[ai];
1526 rbi_inv2 = rbi_inv * rbi_inv;
1527 tmp = (rai*rbi_inv2)*(rai*rbi_inv2);
1528 tmp = tmp*tmp*tmp;
1529 e = factor*term*(rai+probe)*(rai+probe)*tmp;
1530 dvda[ai] = dvda[ai] - 6*e*rbi_inv2;
1531 es = es + e;
1532 }
1533 }
1534
1535 return es;
1536 }
1537
1538
1539
1540 real calc_gb_chainrule(int natoms, t_nblist *nl, real *dadx, real *dvda, rvec x[], rvec t[], rvec fshift[],
1541 rvec shift_vec[], int gb_algorithm, gmx_genborn_t *born, t_mdatoms *md)
1542 {
1543 int i,k,n,ai,aj,nj0,nj1,n0,n1;
1544 int shift;
1545 real shX,shY,shZ;
1546 real fgb,fij,rb2,rbi,fix1,fiy1,fiz1;
1547 real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11,rsq11;
1548 real rinv11,tx,ty,tz,rbai,rbaj,fgb_ai;
1549 real *rb;
1550 volatile int idx;
1551
1552 n = 0;
1553 rb = born->work;
1554
1555 n0 = 0;
1556 n1 = natoms;
1557
1558 if(gb_algorithm==egbSTILL)
1559 {
1560 for(i=n0;i<n1;i++)
1561 {
1562 rbi = born->bRad[i];
1563 rb[i] = (2 * rbi * rbi * dvda[i])/ONE_4PI_EPS0;
1564 }
1565 }
1566 else if(gb_algorithm==egbHCT)
1567 {
1568 for(i=n0;i<n1;i++)
1569 {
1570 rbi = born->bRad[i];
1571 rb[i] = rbi * rbi * dvda[i];
1572 }
1573 }
1574 else if(gb_algorithm==egbOBC)
1575 {
1576 for(i=n0;i<n1;i++)
1577 {
1578 rbi = born->bRad[i];
1579 rb[i] = rbi * rbi * born->drobc[i] * dvda[i];
1580 }
1581 }
1582
1583 for(i=0;i<nl->nri;i++)
1584 {
1585 ai = nl->iinr[i];
1586
1587 nj0 = nl->jindex[i];
1588 nj1 = nl->jindex[i+1];
1589
1590 /* Load shifts for this list */
1591 shift = nl->shift[i];
1592 shX = shift_vec[shift][0];
1593 shY = shift_vec[shift][1];
1594 shZ = shift_vec[shift][2];
1595
1596 /* Load atom i coordinates, add shift vectors */
1597 ix1 = shX + x[ai][0];
1598 iy1 = shY + x[ai][1];
1599 iz1 = shZ + x[ai][2];
1600
1601 fix1 = 0;
1602 fiy1 = 0;
1603 fiz1 = 0;
1604
1605 rbai = rb[ai];
1606
1607 for(k=nj0;k<nj1;k++)
1608 {
1609 aj = nl->jjnr[k];
1610
1611 jx1 = x[aj][0];
1612 jy1 = x[aj][1];
1613 jz1 = x[aj][2];
1614
1615 dx11 = ix1 - jx1;
1616 dy11 = iy1 - jy1;
1617 dz11 = iz1 - jz1;
1618
1619 rbaj = rb[aj];
1620
1621 fgb = rbai*dadx[n++];
1622 fgb_ai = rbaj*dadx[n++];
1623
1624 /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1625 fgb = fgb + fgb_ai;
1626
1627 tx = fgb * dx11;
1628 ty = fgb * dy11;
1629 tz = fgb * dz11;
1630
1631 fix1 = fix1 + tx;
1632 fiy1 = fiy1 + ty;
1633 fiz1 = fiz1 + tz;
1634
1635 /* Update force on atom aj */
1636 t[aj][0] = t[aj][0] - tx;
1637 t[aj][1] = t[aj][1] - ty;
1638 t[aj][2] = t[aj][2] - tz;
1639 }
1640
1641 /* Update force and shift forces on atom ai */
1642 t[ai][0] = t[ai][0] + fix1;
1643 t[ai][1] = t[ai][1] + fiy1;
1644 t[ai][2] = t[ai][2] + fiz1;
1645
1646 fshift[shift][0] = fshift[shift][0] + fix1;
1647 fshift[shift][1] = fshift[shift][1] + fiy1;
1648 fshift[shift][2] = fshift[shift][2] + fiz1;
1649
1650 }
1651
1652 return 0;
1653 }
1654
1655
1656 void
1657 calc_gb_forces(t_commrec *cr, t_mdatoms *md, gmx_genborn_t *born, gmx_localtop_t *top, const t_atomtypes *atype,
1658 rvec x[], rvec f[], t_forcerec *fr, t_idef *idef, int gb_algorithm, int sa_algorithm, t_nrnb *nrnb, gmx_bool bRad,
1659 const t_pbc *pbc, const t_graph *graph, gmx_enerdata_t *enerd)
1660 {
1661 real v=0;
1662 int cnt;
1663 int i;
1664
1665 /* PBC or not? */
1666 const t_pbc *pbc_null;
1667
1668 if (fr->bMolPBC)
1669 pbc_null = pbc;
1670 else
1671 pbc_null = NULL;
1672
1673 if(sa_algorithm == esaAPPROX)
1674 {
1675 /* Do a simple ACE type approximation for the non-polar solvation */
1676 enerd->term[F_NPSOLVATION] += calc_gb_nonpolar(cr, fr,born->nr, born, top, atype, fr->dvda, gb_algorithm,md);
1677 }
1678
1679 /* Calculate the bonded GB-interactions using either table or analytical formula */
1680 enerd->term[F_GBPOL] += gb_bonds_tab(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1681 fr->invsqrta,fr->dvda,fr->gbtab.tab,idef,born->epsilon_r,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1682
1683 /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1684 enerd->term[F_GBPOL] += calc_gb_selfcorrections(cr,born->nr,md->chargeA, born, fr->dvda, md, fr->epsfac);
1685
1686 /* If parallel, sum the derivative of the potential w.r.t the born radii */
1687 if(PARTDECOMP(cr))
1688 {
1689 gmx_sum(md->nr,fr->dvda, cr);
1690 }
1691 else if(DOMAINDECOMP(cr))
1692 {
1693 dd_atom_sum_real(cr->dd,fr->dvda);
1694 dd_atom_spread_real(cr->dd,fr->dvda);
1695 }
1696
1697 if(fr->bAllvsAll)
1698 {
1699 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || (!defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1700 if(fr->UseOptimizedKernels)
1701 {
1702 genborn_allvsall_calc_chainrule_sse2_single(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1703 }
1704 else
1705 {
1706 genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1707 }
1708 #elif ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || (defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1709 if(fr->UseOptimizedKernels)
1710 {
1711 genborn_allvsall_calc_chainrule_sse2_double(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1712 }
1713 else
1714 {
1715 genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1716 }
1717 #else
1718 genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1719 #endif
1720 cnt = md->homenr*(md->nr/2+1);
1721 inc_nrnb(nrnb,eNR_BORN_AVA_CHAINRULE,cnt);
1722 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,md->homenr);
1723 return;
1724 }
1725
1726 #ifdef GMX_DOUBLE
1727
1728 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || (defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1729 if(fr->UseOptimizedKernels)
1730 {
1731 calc_gb_chainrule_sse2_double(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1732 x[0], f[0], fr->fshift[0], fr->shift_vec[0],
1733 gb_algorithm, born, md);
1734 }
1735 else
1736 {
1737 calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1738 x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1739 }
1740 #else
1741 calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1742 x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1743 #endif
1744
1745 #else
1746
1747 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || (!defined(GMX_DOUBLE) && defined(GMX_SSE2)) )
1748 /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1749 if(fr->UseOptimizedKernels)
1750 {
1751 calc_gb_chainrule_sse2_single(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1752 x[0], f[0], fr->fshift[0], fr->shift_vec[0],
1753 gb_algorithm, born, md);
1754 }
1755 else
1756 {
1757 calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1758 x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1759 }
1760
1761 #else
1762 /* Calculate the forces due to chain rule terms with non sse code */
1763 calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1764 x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1765 #endif
1766 #endif
1767
1768 if(!fr->bAllvsAll)
1769 {
1770 inc_nrnb(nrnb,eNR_BORN_CHAINRULE,fr->gblist.nrj);
1771 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,fr->gblist.nri);
1772
1773 }
1774 }
1775
1776 static void add_j_to_gblist(gbtmpnbl_t *list,int aj)
1777 {
1778 if (list->naj >= list->aj_nalloc)
1779 {
1780 list->aj_nalloc = over_alloc_large(list->naj+1);
1781 srenew(list->aj,list->aj_nalloc);
1782 }
1783
1784 list->aj[list->naj++] = aj;
1785 }
1786
1787 static gbtmpnbl_t *find_gbtmplist(struct gbtmpnbls *lists,int shift)
1788 {
1789 int ind,i;
1790
1791 /* Search the list with the same shift, if there is one */
1792 ind = 0;
1793 while (ind < lists->nlist && shift != lists->list[ind].shift)
1794 {
1795 ind++;
1796 }
1797 if (ind == lists->nlist)
1798 {
1799 if (lists->nlist == lists->list_nalloc)
1800 {
1801 lists->list_nalloc++;
1802 srenew(lists->list,lists->list_nalloc);
1803 for(i=lists->nlist; i<lists->list_nalloc; i++)
1804 {
1805 lists->list[i].aj = NULL;
1806 lists->list[i].aj_nalloc = 0;
1807 }
1808
1809 }
1810
1811 lists->list[lists->nlist].shift = shift;
1812 lists->list[lists->nlist].naj = 0;
1813 lists->nlist++;
1814 }
1815
1816 return &lists->list[ind];
1817 }
1818
1819 static void add_bondeds_to_gblist(t_ilist *il,
1820 gmx_bool bMolPBC,t_pbc *pbc,t_graph *g,rvec *x,
1821 struct gbtmpnbls *nls)
1822 {
1823 int ind,j,ai,aj,shift,found;
1824 rvec dx;
1825 ivec dt;
1826 gbtmpnbl_t *list;
1827
1828 shift = CENTRAL;
1829 for(ind=0; ind<il->nr; ind+=3)
1830 {
1831 ai = il->iatoms[ind+1];
1832 aj = il->iatoms[ind+2];
1833
1834 shift = CENTRAL;
1835 if (g != NULL)
1836 {
1837 rvec_sub(x[ai],x[aj],dx);
1838 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1839 shift = IVEC2IS(dt);
1840 }
1841 else if (bMolPBC)
1842 {
1843 shift = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
1844 }
1845
1846 /* Find the list for this shift or create one */
1847 list = find_gbtmplist(&nls[ai],shift);
1848
1849 found=0;
1850
1851 /* So that we do not add the same bond twice.
1852 * This happens with some constraints between 1-3 atoms
1853 * that are in the bond-list but should not be in the GB nb-list */
1854 for(j=0;j<list->naj;j++)
1855 {
1856 if (list->aj[j] == aj)
1857 {
1858 found = 1;
1859 }
1860 }
1861
1862 if (found == 0)
1863 {
1864 if(ai == aj)
1865 {
1866 gmx_incons("ai == aj");
1867 }
1868
1869 add_j_to_gblist(list,aj);
1870 }
1871 }
1872 }
1873
1874 static int
1875 compare_int (const void * a, const void * b)
1876 {
1877 return ( *(int*)a - *(int*)b );
1878 }
1879
1880
1881
1882 int make_gb_nblist(t_commrec *cr, int gb_algorithm, real gbcut,
1883 rvec x[], matrix box,
1884 t_forcerec *fr, t_idef *idef, t_graph *graph, gmx_genborn_t *born)
1885 {
1886 int i,l,ii,j,k,n,nj0,nj1,ai,aj,at0,at1,found,shift,s;
1887 int apa;
1888 t_nblist *nblist;
1889 t_pbc pbc;
1890
1891 struct gbtmpnbls *nls;
1892 gbtmpnbl_t *list =NULL;
1893
1894 set_pbc(&pbc,fr->ePBC,box);
1895 nls = born->nblist_work;
1896
1897 for(i=0;i<born->nr;i++)
1898 {
1899 nls[i].nlist = 0;
1900 }
1901
1902 if (fr->bMolPBC)
1903 {
1904 set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);
1905 }
1906
1907 switch (gb_algorithm)
1908 {
1909 case egbHCT:
1910 case egbOBC:
1911 /* Loop over 1-2, 1-3 and 1-4 interactions */
1912 for(j=F_GB12;j<=F_GB14;j++)
1913 {
1914 add_bondeds_to_gblist(&idef->il[j],fr->bMolPBC,&pbc,graph,x,nls);
1915 }
1916 break;
1917 case egbSTILL:
1918 /* Loop over 1-4 interactions */
1919 add_bondeds_to_gblist(&idef->il[F_GB14],fr->bMolPBC,&pbc,graph,x,nls);
1920 break;
1921 default:
1922 gmx_incons("Unknown GB algorithm");
1923 }
1924
1925 /* Loop over the VDWQQ and VDW nblists to set up the nonbonded part of the GB list */
1926 for(n=0; (n<fr->nnblists); n++)
1927 {
1928 for(i=0; (i<eNL_NR); i++)
1929 {
1930 nblist=&(fr->nblists[n].nlist_sr[i]);
1931
1932 if (nblist->nri > 0 && (i==eNL_VDWQQ || i==eNL_QQ))
1933 {
1934 for(j=0;j<nblist->nri;j++)
1935 {
1936 ai = nblist->iinr[j];
1937 shift = nblist->shift[j];
1938
1939 /* Find the list for this shift or create one */
1940 list = find_gbtmplist(&nls[ai],shift);
1941
1942 nj0 = nblist->jindex[j];
1943 nj1 = nblist->jindex[j+1];
1944
1945 /* Add all the j-atoms in the non-bonded list to the GB list */
1946 for(k=nj0;k<nj1;k++)
1947 {
1948 add_j_to_gblist(list,nblist->jjnr[k]);
1949 }
1950 }
1951 }
1952 }
1953 }
1954
1955 /* Zero out some counters */
1956 fr->gblist.nri=0;
1957 fr->gblist.nrj=0;
1958
1959 fr->gblist.jindex[0] = fr->gblist.nri;
1960
1961 for(i=0;i<fr->natoms_force;i++)
1962 {
1963 for(s=0; s<nls[i].nlist; s++)
1964 {
1965 list = &nls[i].list[s];
1966
1967 /* Only add those atoms that actually have neighbours */
1968 if (born->use[i] != 0)
1969 {
1970 fr->gblist.iinr[fr->gblist.nri] = i;
1971 fr->gblist.shift[fr->gblist.nri] = list->shift;
1972 fr->gblist.nri++;
1973
1974 for(k=0; k<list->naj; k++)
1975 {
1976 /* Memory allocation for jjnr */
1977 if(fr->gblist.nrj >= fr->gblist.maxnrj)
1978 {
1979 fr->gblist.maxnrj += over_alloc_large(fr->gblist.maxnrj);
1980
1981 if (debug)
1982 {
1983 fprintf(debug,"Increasing GB neighbourlist j size to %d\n",fr->gblist.maxnrj);
1984 }
1985
1986 srenew(fr->gblist.jjnr,fr->gblist.maxnrj);
1987 }
1988
1989 /* Put in list */
1990 if(i == list->aj[k])
1991 {
1992 gmx_incons("i == list->aj[k]");
1993 }
1994 fr->gblist.jjnr[fr->gblist.nrj++] = list->aj[k];
1995 }
1996
1997 fr->gblist.jindex[fr->gblist.nri] = fr->gblist.nrj;
1998 }
1999 }
2000 }
2001
2002
2003 #ifdef SORT_GB_LIST
2004 for(i=0;i<fr->gblist.nri;i++)
2005 {
2006 nj0 = fr->gblist.jindex[i];
2007 nj1 = fr->gblist.jindex[i+1];
2008 ai = fr->gblist.iinr[i];
2009
2010 /* Temporary fix */
2011 for(j=nj0;j<nj1;j++)
2012 {
2013 if(fr->gblist.jjnr[j]<ai)
2014 fr->gblist.jjnr[j]+=fr->natoms_force;
2015 }
2016 qsort(fr->gblist.jjnr+nj0,nj1-nj0,sizeof(int),compare_int);
2017 /* Fix back */
2018 for(j=nj0;j<nj1;j++)
2019 {
2020 if(fr->gblist.jjnr[j]>=fr->natoms_force)
2021 fr->gblist.jjnr[j]-=fr->natoms_force;
2022 }
2023
2024 }
2025 #endif
2026
2027 return 0;
2028 }
2029
2030 void make_local_gb(const t_commrec *cr, gmx_genborn_t *born, int gb_algorithm)
2031 {
2032 int i,at0,at1;
2033 gmx_domdec_t *dd=NULL;
2034
2035 if(DOMAINDECOMP(cr))
2036 {
2037 dd = cr->dd;
2038 at0 = 0;
2039 at1 = dd->nat_tot;
2040 }
2041 else
2042 {
2043 /* Single node or particle decomp (global==local), just copy pointers and return */
2044 if(gb_algorithm==egbSTILL)
2045 {
2046 born->gpol = born->gpol_globalindex;
2047 born->vsolv = born->vsolv_globalindex;
2048 born->gb_radius = born->gb_radius_globalindex;
2049 }
2050 else
2051 {
2052 born->param = born->param_globalindex;
2053 born->gb_radius = born->gb_radius_globalindex;
2054 }
2055
2056 born->use = born->use_globalindex;
2057
2058 return;
2059 }
2060
2061 /* Reallocation of local arrays if necessary */
2062 /* fr->natoms_force is equal to dd->nat_tot */
2063 if (DOMAINDECOMP(cr) && dd->nat_tot > born->nalloc)
2064 {
2065 int nalloc;
2066
2067 nalloc = dd->nat_tot;
2068
2069 /* Arrays specific to different gb algorithms */
2070 if (gb_algorithm == egbSTILL)
2071 {
2072 srenew(born->gpol, nalloc+3);
2073 srenew(born->vsolv, nalloc+3);
2074 srenew(born->gb_radius, nalloc+3);
2075 for(i=born->nalloc; (i<nalloc+3); i++)
2076 {
2077 born->gpol[i] = 0;
2078 born->vsolv[i] = 0;
2079 born->gb_radius[i] = 0;
2080 }
2081 }
2082 else
2083 {
2084 srenew(born->param, nalloc+3);
2085 srenew(born->gb_radius, nalloc+3);
2086 for(i=born->nalloc; (i<nalloc+3); i++)
2087 {
2088 born->param[i] = 0;
2089 born->gb_radius[i] = 0;
2090 }
2091 }
2092
2093 /* All gb-algorithms use the array for vsites exclusions */
2094 srenew(born->use, nalloc+3);
2095 for(i=born->nalloc; (i<nalloc+3); i++)
2096 {
2097 born->use[i] = 0;
2098 }
2099
2100 born->nalloc = nalloc;
2101 }
2102
2103 /* With dd, copy algorithm specific arrays */
2104 if(gb_algorithm==egbSTILL)
2105 {
2106 for(i=at0;i<at1;i++)
2107 {
2108 born->gpol[i] = born->gpol_globalindex[dd->gatindex[i]];
2109 born->vsolv[i] = born->vsolv_globalindex[dd->gatindex[i]];
2110 born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2111 born->use[i] = born->use_globalindex[dd->gatindex[i]];
2112 }
2113 }
2114 else
2115 {
2116 for(i=at0;i<at1;i++)
2117 {
2118 born->param[i] = born->param_globalindex[dd->gatindex[i]];
2119 born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2120 born->use[i] = born->use_globalindex[dd->gatindex[i]];
2121 }
2122 }
2123 }
2124
The ultimate molecular dynamics simulation package