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