2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gromacs/domdec/domdec.h"
49 #include "gromacs/fileio/pdbio.h"
50 #include "gromacs/gmxlib/network.h"
51 #include "gromacs/legacyheaders/names.h"
52 #include "gromacs/legacyheaders/nrnb.h"
53 #include "gromacs/legacyheaders/types/commrec.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/mdlib/genborn_allvsall.h"
57 #include "gromacs/pbcutil/ishift.h"
58 #include "gromacs/pbcutil/mshift.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/topology/mtop_util.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/gmxmpi.h"
63 #include "gromacs/utility/smalloc.h"
73 typedef struct gbtmpnbls
{
79 /* This function is exactly the same as the one in listed-forces/bonded.cpp. The reason
80 * it is copied here is that the bonded gb-interactions are evaluated
81 * not in calc_bonds, but rather in calc_gb_forces
83 static int pbc_rvec_sub(const t_pbc
*pbc
, const rvec xi
, const rvec xj
, rvec dx
)
87 return pbc_dx_aiuc(pbc
, xi
, xj
, dx
);
96 static int init_gb_nblist(int natoms
, t_nblist
*nl
)
98 nl
->maxnri
= natoms
*4;
107 /*nl->nltype = nltype;*/
109 srenew(nl
->iinr
, nl
->maxnri
);
110 srenew(nl
->gid
, nl
->maxnri
);
111 srenew(nl
->shift
, nl
->maxnri
);
112 srenew(nl
->jindex
, nl
->maxnri
+1);
120 static int init_gb_still(const t_atomtypes
*atype
, t_idef
*idef
, t_atoms
*atoms
,
121 gmx_genborn_t
*born
, int natoms
)
125 real r
, ri
, rj
, ri2
, rj2
, r3
, r4
, ratio
, term
, h
, doffset
;
132 snew(born
->gpol_still_work
, natoms
+3);
134 doffset
= born
->gb_doffset
;
136 for (i
= 0; i
< natoms
; i
++)
138 born
->gpol_globalindex
[i
] = born
->vsolv_globalindex
[i
] =
139 born
->gb_radius_globalindex
[i
] = 0;
142 /* Compute atomic solvation volumes for Still method */
143 for (i
= 0; i
< natoms
; i
++)
145 ri
= atype
->gb_radius
[atoms
->atom
[i
].type
];
146 born
->gb_radius_globalindex
[i
] = ri
;
148 born
->vsolv_globalindex
[i
] = (4*M_PI
/3)*r3
;
151 for (j
= 0; j
< idef
->il
[F_GB12
].nr
; j
+= 3)
153 m
= idef
->il
[F_GB12
].iatoms
[j
];
154 ia
= idef
->il
[F_GB12
].iatoms
[j
+1];
155 ib
= idef
->il
[F_GB12
].iatoms
[j
+2];
157 r
= 1.01*idef
->iparams
[m
].gb
.st
;
159 ri
= atype
->gb_radius
[atoms
->atom
[ia
].type
];
160 rj
= atype
->gb_radius
[atoms
->atom
[ib
].type
];
165 ratio
= (rj2
-ri2
-r
*r
)/(2*ri
*r
);
167 term
= (M_PI
/3.0)*h
*h
*(3.0*ri
-h
);
169 born
->vsolv_globalindex
[ia
] -= term
;
171 ratio
= (ri2
-rj2
-r
*r
)/(2*rj
*r
);
173 term
= (M_PI
/3.0)*h
*h
*(3.0*rj
-h
);
175 born
->vsolv_globalindex
[ib
] -= term
;
178 /* Get the self-, 1-2 and 1-3 polarization energies for analytical Still
181 for (j
= 0; j
< natoms
; j
++)
183 if (born
->use_globalindex
[j
] == 1)
185 born
->gpol_globalindex
[j
] = -0.5*ONE_4PI_EPS0
/
186 (atype
->gb_radius
[atoms
->atom
[j
].type
]-doffset
+STILL_P1
);
191 for (j
= 0; j
< idef
->il
[F_GB12
].nr
; j
+= 3)
193 m
= idef
->il
[F_GB12
].iatoms
[j
];
194 ia
= idef
->il
[F_GB12
].iatoms
[j
+1];
195 ib
= idef
->il
[F_GB12
].iatoms
[j
+2];
197 r
= idef
->iparams
[m
].gb
.st
;
201 born
->gpol_globalindex
[ia
] = born
->gpol_globalindex
[ia
]+
202 STILL_P2
*born
->vsolv_globalindex
[ib
]/r4
;
203 born
->gpol_globalindex
[ib
] = born
->gpol_globalindex
[ib
]+
204 STILL_P2
*born
->vsolv_globalindex
[ia
]/r4
;
208 for (j
= 0; j
< idef
->il
[F_GB13
].nr
; j
+= 3)
210 m
= idef
->il
[F_GB13
].iatoms
[j
];
211 ia
= idef
->il
[F_GB13
].iatoms
[j
+1];
212 ib
= idef
->il
[F_GB13
].iatoms
[j
+2];
214 r
= idef
->iparams
[m
].gb
.st
;
217 born
->gpol_globalindex
[ia
] = born
->gpol_globalindex
[ia
]+
218 STILL_P3
*born
->vsolv_globalindex
[ib
]/r4
;
219 born
->gpol_globalindex
[ib
] = born
->gpol_globalindex
[ib
]+
220 STILL_P3
*born
->vsolv_globalindex
[ia
]/r4
;
229 /* Initialize all GB datastructs and compute polarization energies */
230 int init_gb(gmx_genborn_t
**p_born
,
231 t_forcerec
*fr
, const t_inputrec
*ir
,
232 const gmx_mtop_t
*mtop
, int gb_algorithm
)
235 real rai
, sk
, doffset
;
239 gmx_localtop_t
*localtop
;
241 natoms
= mtop
->natoms
;
243 atoms
= gmx_mtop_global_atoms(mtop
);
244 localtop
= gmx_mtop_generate_local_top(mtop
, ir
);
251 snew(born
->drobc
, natoms
);
252 snew(born
->bRad
, natoms
);
254 /* Allocate memory for the global data arrays */
255 snew(born
->param_globalindex
, natoms
+3);
256 snew(born
->gpol_globalindex
, natoms
+3);
257 snew(born
->vsolv_globalindex
, natoms
+3);
258 snew(born
->gb_radius_globalindex
, natoms
+3);
259 snew(born
->use_globalindex
, natoms
+3);
261 snew(fr
->invsqrta
, natoms
);
262 snew(fr
->dvda
, natoms
);
265 fr
->dadx_rawptr
= NULL
;
267 born
->gpol_still_work
= NULL
;
268 born
->gpol_hct_work
= NULL
;
270 /* snew(born->asurf,natoms); */
271 /* snew(born->dasurf,natoms); */
273 /* Initialize the gb neighbourlist */
274 init_gb_nblist(natoms
, &(fr
->gblist
));
276 /* Do the Vsites exclusions (if any) */
277 for (i
= 0; i
< natoms
; i
++)
279 jj
= atoms
.atom
[i
].type
;
280 if (mtop
->atomtypes
.gb_radius
[atoms
.atom
[i
].type
] > 0)
282 born
->use_globalindex
[i
] = 1;
286 born
->use_globalindex
[i
] = 0;
289 /* If we have a Vsite, put vs_globalindex[i]=0 */
290 if (C6 (fr
->nbfp
, fr
->ntype
, jj
, jj
) == 0 &&
291 C12(fr
->nbfp
, fr
->ntype
, jj
, jj
) == 0 &&
292 atoms
.atom
[i
].q
== 0)
294 born
->use_globalindex
[i
] = 0;
298 /* Copy algorithm parameters from inputrecord to local structure */
299 born
->obc_alpha
= ir
->gb_obc_alpha
;
300 born
->obc_beta
= ir
->gb_obc_beta
;
301 born
->obc_gamma
= ir
->gb_obc_gamma
;
302 born
->gb_doffset
= ir
->gb_dielectric_offset
;
303 born
->gb_epsilon_solvent
= ir
->gb_epsilon_solvent
;
304 born
->epsilon_r
= ir
->epsilon_r
;
306 doffset
= born
->gb_doffset
;
308 /* Set the surface tension */
309 born
->sa_surface_tension
= ir
->sa_surface_tension
;
311 /* If Still model, initialise the polarisation energies */
312 if (gb_algorithm
== egbSTILL
)
314 init_gb_still(&(mtop
->atomtypes
), &(localtop
->idef
), &atoms
,
319 /* If HCT/OBC, precalculate the sk*atype->S_hct factors */
320 else if (gb_algorithm
== egbHCT
|| gb_algorithm
== egbOBC
)
323 snew(born
->gpol_hct_work
, natoms
+3);
325 for (i
= 0; i
< natoms
; i
++)
327 if (born
->use_globalindex
[i
] == 1)
329 rai
= mtop
->atomtypes
.gb_radius
[atoms
.atom
[i
].type
]-doffset
;
330 sk
= rai
* mtop
->atomtypes
.S_hct
[atoms
.atom
[i
].type
];
331 born
->param_globalindex
[i
] = sk
;
332 born
->gb_radius_globalindex
[i
] = rai
;
336 born
->param_globalindex
[i
] = 0;
337 born
->gb_radius_globalindex
[i
] = 0;
342 /* Allocate memory for work arrays for temporary use */
343 snew(born
->work
, natoms
+4);
344 snew(born
->count
, natoms
);
345 snew(born
->nblist_work
, natoms
);
347 /* Domain decomposition specific stuff */
356 calc_gb_rad_still(t_commrec
*cr
, t_forcerec
*fr
, gmx_localtop_t
*top
,
357 rvec x
[], t_nblist
*nl
,
358 gmx_genborn_t
*born
, t_mdatoms
*md
)
360 int i
, k
, n
, nj0
, nj1
, ai
, aj
;
363 real gpi
, dr2
, idr4
, rvdw
, ratio
, ccf
, theta
, term
, rai
, raj
;
364 real ix1
, iy1
, iz1
, jx1
, jy1
, jz1
, dx11
, dy11
, dz11
;
365 real rinv
, idr2
, idr6
, vaj
, dccf
, cosq
, sinq
, prod
, gpi2
;
367 real vai
, prod_ai
, icf4
, icf6
;
369 factor
= 0.5*ONE_4PI_EPS0
;
372 for (i
= 0; i
< born
->nr
; i
++)
374 born
->gpol_still_work
[i
] = 0;
377 for (i
= 0; i
< nl
->nri
; i
++)
382 nj1
= nl
->jindex
[i
+1];
384 /* Load shifts for this list */
385 shift
= nl
->shift
[i
];
386 shX
= fr
->shift_vec
[shift
][0];
387 shY
= fr
->shift_vec
[shift
][1];
388 shZ
= fr
->shift_vec
[shift
][2];
392 rai
= top
->atomtypes
.gb_radius
[md
->typeA
[ai
]];
393 vai
= born
->vsolv
[ai
];
394 prod_ai
= STILL_P4
*vai
;
396 /* Load atom i coordinates, add shift vectors */
397 ix1
= shX
+ x
[ai
][0];
398 iy1
= shY
+ x
[ai
][1];
399 iz1
= shZ
+ x
[ai
][2];
401 for (k
= nj0
; k
< nj1
&& nl
->jjnr
[k
] >= 0; k
++)
412 dr2
= dx11
*dx11
+dy11
*dy11
+dz11
*dz11
;
413 rinv
= gmx_invsqrt(dr2
);
418 raj
= top
->atomtypes
.gb_radius
[md
->typeA
[aj
]];
422 ratio
= dr2
/ (rvdw
* rvdw
);
423 vaj
= born
->vsolv
[aj
];
425 if (ratio
> STILL_P5INV
)
432 theta
= ratio
*STILL_PIP5
;
434 term
= 0.5*(1.0-cosq
);
436 sinq
= 1.0 - cosq
*cosq
;
437 dccf
= 2.0*term
*sinq
*gmx_invsqrt(sinq
)*theta
;
442 icf6
= (4*ccf
-dccf
)*idr6
;
443 born
->gpol_still_work
[aj
] += prod_ai
*icf4
;
446 /* Save ai->aj and aj->ai chain rule terms */
447 fr
->dadx
[n
++] = prod
*icf6
;
448 fr
->dadx
[n
++] = prod_ai
*icf6
;
450 born
->gpol_still_work
[ai
] += gpi
;
453 /* Parallel summations */
454 if (DOMAINDECOMP(cr
))
456 dd_atom_sum_real(cr
->dd
, born
->gpol_still_work
);
459 /* Calculate the radii */
460 for (i
= 0; i
< fr
->natoms_force
; i
++) /* PELA born->nr */
462 if (born
->use
[i
] != 0)
464 gpi
= born
->gpol
[i
]+born
->gpol_still_work
[i
];
466 born
->bRad
[i
] = factor
*gmx_invsqrt(gpi2
);
467 fr
->invsqrta
[i
] = gmx_invsqrt(born
->bRad
[i
]);
471 /* Extra communication required for DD */
472 if (DOMAINDECOMP(cr
))
474 dd_atom_spread_real(cr
->dd
, born
->bRad
);
475 dd_atom_spread_real(cr
->dd
, fr
->invsqrta
);
484 calc_gb_rad_hct(t_commrec
*cr
, t_forcerec
*fr
, gmx_localtop_t
*top
,
485 rvec x
[], t_nblist
*nl
,
486 gmx_genborn_t
*born
, t_mdatoms
*md
)
488 int i
, k
, n
, ai
, aj
, nj0
, nj1
;
491 real rai
, raj
, dr2
, dr
, sk
, sk_ai
, sk2
, sk2_ai
, lij
, uij
, diff2
, tmp
, sum_ai
;
492 real rad
, min_rad
, rinv
, rai_inv
;
493 real ix1
, iy1
, iz1
, jx1
, jy1
, jz1
, dx11
, dy11
, dz11
;
494 real lij2
, uij2
, lij3
, uij3
, t1
, t2
, t3
;
495 real lij_inv
, dlij
, sk2_rinv
, prod
, log_term
;
496 real doffset
, raj_inv
, dadx_val
;
499 doffset
= born
->gb_doffset
;
500 gb_radius
= born
->gb_radius
;
502 for (i
= 0; i
< born
->nr
; i
++)
504 born
->gpol_hct_work
[i
] = 0;
507 /* Keep the compiler happy */
510 for (i
= 0; i
< nl
->nri
; i
++)
515 nj1
= nl
->jindex
[i
+1];
517 /* Load shifts for this list */
518 shift
= nl
->shift
[i
];
519 shX
= fr
->shift_vec
[shift
][0];
520 shY
= fr
->shift_vec
[shift
][1];
521 shZ
= fr
->shift_vec
[shift
][2];
526 sk_ai
= born
->param
[ai
];
527 sk2_ai
= sk_ai
*sk_ai
;
529 /* Load atom i coordinates, add shift vectors */
530 ix1
= shX
+ x
[ai
][0];
531 iy1
= shY
+ x
[ai
][1];
532 iz1
= shZ
+ x
[ai
][2];
536 for (k
= nj0
; k
< nj1
&& nl
->jjnr
[k
] >= 0; k
++)
548 dr2
= dx11
*dx11
+dy11
*dy11
+dz11
*dz11
;
549 rinv
= gmx_invsqrt(dr2
);
552 sk
= born
->param
[aj
];
555 /* aj -> ai interaction */
576 lij_inv
= gmx_invsqrt(lij2
);
579 prod
= 0.25*sk2_rinv
;
581 log_term
= std::log(uij
*lij_inv
);
583 tmp
= lij
-uij
+ 0.25*dr
*diff2
+ (0.5*rinv
)*log_term
+
588 tmp
= tmp
+ 2.0 * (rai_inv
-lij
);
591 t1
= 0.5*lij2
+ prod
*lij3
- 0.25*(lij
*rinv
+lij3
*dr
);
592 t2
= -0.5*uij2
- 0.25*sk2_rinv
*uij3
+ 0.25*(uij
*rinv
+uij3
*dr
);
593 t3
= 0.125*(1.0+sk2_rinv
*rinv
)*(-diff2
)+0.25*log_term
*rinv
*rinv
;
595 dadx_val
= (dlij
*t1
+t2
+t3
)*rinv
; /* rb2 is moved to chainrule */
596 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */
597 /* rb2 is moved to chainrule */
605 fr
->dadx
[n
++] = dadx_val
;
608 /* ai -> aj interaction */
609 if (raj
< dr
+ sk_ai
)
611 lij
= 1.0/(dr
-sk_ai
);
624 uij
= 1.0/(dr
+sk_ai
);
630 lij_inv
= gmx_invsqrt(lij2
);
631 sk2
= sk2_ai
; /* sk2_ai = sk_ai * sk_ai in i loop above */
633 prod
= 0.25 * sk2_rinv
;
635 /* log_term = table_log(uij*lij_inv,born->log_table,
636 LOG_TABLE_ACCURACY); */
637 log_term
= std::log(uij
*lij_inv
);
639 tmp
= lij
-uij
+ 0.25*dr
*diff2
+ (0.5*rinv
)*log_term
+
644 tmp
= tmp
+ 2.0 * (raj_inv
-lij
);
648 t1
= 0.5*lij2
+ prod
*lij3
- 0.25*(lij
*rinv
+lij3
*dr
);
649 t2
= -0.5*uij2
- 0.25*sk2_rinv
*uij3
+ 0.25*(uij
*rinv
+uij3
*dr
);
650 t3
= 0.125*(1.0+sk2_rinv
*rinv
)*(-diff2
)+0.25*log_term
*rinv
*rinv
;
652 dadx_val
= (dlij
*t1
+t2
+t3
)*rinv
; /* rb2 is moved to chainrule */
653 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */ /* rb2 is moved to chainrule */
655 born
->gpol_hct_work
[aj
] += 0.5*tmp
;
661 fr
->dadx
[n
++] = dadx_val
;
664 born
->gpol_hct_work
[ai
] += sum_ai
;
667 /* Parallel summations */
668 if (DOMAINDECOMP(cr
))
670 dd_atom_sum_real(cr
->dd
, born
->gpol_hct_work
);
673 for (i
= 0; i
< fr
->natoms_force
; i
++) /* PELA born->nr */
675 if (born
->use
[i
] != 0)
677 rai
= top
->atomtypes
.gb_radius
[md
->typeA
[i
]]-doffset
;
678 sum_ai
= 1.0/rai
- born
->gpol_hct_work
[i
];
679 min_rad
= rai
+ doffset
;
682 born
->bRad
[i
] = std::max(rad
, min_rad
);
683 fr
->invsqrta
[i
] = gmx_invsqrt(born
->bRad
[i
]);
687 /* Extra communication required for DD */
688 if (DOMAINDECOMP(cr
))
690 dd_atom_spread_real(cr
->dd
, born
->bRad
);
691 dd_atom_spread_real(cr
->dd
, fr
->invsqrta
);
699 calc_gb_rad_obc(t_commrec
*cr
, t_forcerec
*fr
, gmx_localtop_t
*top
,
700 rvec x
[], t_nblist
*nl
, gmx_genborn_t
*born
, t_mdatoms
*md
)
702 int i
, k
, ai
, aj
, nj0
, nj1
, n
;
705 real rai
, raj
, dr2
, dr
, sk
, sk2
, lij
, uij
, diff2
, tmp
, sum_ai
;
706 real sum_ai2
, sum_ai3
, tsum
, tchain
, rinv
, rai_inv
, lij_inv
, rai_inv2
;
707 real log_term
, prod
, sk2_rinv
, sk_ai
, sk2_ai
;
708 real ix1
, iy1
, iz1
, jx1
, jy1
, jz1
, dx11
, dy11
, dz11
;
709 real lij2
, uij2
, lij3
, uij3
, dlij
, t1
, t2
, t3
;
710 real doffset
, raj_inv
, dadx_val
;
713 /* Keep the compiler happy */
716 doffset
= born
->gb_doffset
;
717 gb_radius
= born
->gb_radius
;
719 for (i
= 0; i
< born
->nr
; i
++)
721 born
->gpol_hct_work
[i
] = 0;
724 for (i
= 0; i
< nl
->nri
; i
++)
729 nj1
= nl
->jindex
[i
+1];
731 /* Load shifts for this list */
732 shift
= nl
->shift
[i
];
733 shX
= fr
->shift_vec
[shift
][0];
734 shY
= fr
->shift_vec
[shift
][1];
735 shZ
= fr
->shift_vec
[shift
][2];
740 sk_ai
= born
->param
[ai
];
741 sk2_ai
= sk_ai
*sk_ai
;
743 /* Load atom i coordinates, add shift vectors */
744 ix1
= shX
+ x
[ai
][0];
745 iy1
= shY
+ x
[ai
][1];
746 iz1
= shZ
+ x
[ai
][2];
750 for (k
= nj0
; k
< nj1
&& nl
->jjnr
[k
] >= 0; k
++)
762 dr2
= dx11
*dx11
+dy11
*dy11
+dz11
*dz11
;
763 rinv
= gmx_invsqrt(dr2
);
766 /* sk is precalculated in init_gb() */
767 sk
= born
->param
[aj
];
770 /* aj -> ai interaction */
790 lij_inv
= gmx_invsqrt(lij2
);
793 prod
= 0.25*sk2_rinv
;
795 log_term
= std::log(uij
*lij_inv
);
797 tmp
= lij
-uij
+ 0.25*dr
*diff2
+ (0.5*rinv
)*log_term
+ prod
*(-diff2
);
801 tmp
= tmp
+ 2.0 * (rai_inv
-lij
);
805 t1
= 0.5*lij2
+ prod
*lij3
- 0.25*(lij
*rinv
+lij3
*dr
);
806 t2
= -0.5*uij2
- 0.25*sk2_rinv
*uij3
+ 0.25*(uij
*rinv
+uij3
*dr
);
807 t3
= 0.125*(1.0+sk2_rinv
*rinv
)*(-diff2
)+0.25*log_term
*rinv
*rinv
;
809 dadx_val
= (dlij
*t1
+t2
+t3
)*rinv
; /* rb2 is moved to chainrule */
817 fr
->dadx
[n
++] = dadx_val
;
819 /* ai -> aj interaction */
820 if (raj
< dr
+ sk_ai
)
822 lij
= 1.0/(dr
-sk_ai
);
835 uij
= 1.0/(dr
+sk_ai
);
841 lij_inv
= gmx_invsqrt(lij2
);
842 sk2
= sk2_ai
; /* sk2_ai = sk_ai * sk_ai in i loop above */
844 prod
= 0.25 * sk2_rinv
;
846 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
847 log_term
= std::log(uij
*lij_inv
);
849 tmp
= lij
-uij
+ 0.25*dr
*diff2
+ (0.5*rinv
)*log_term
+ prod
*(-diff2
);
853 tmp
= tmp
+ 2.0 * (raj_inv
-lij
);
856 t1
= 0.5*lij2
+ prod
*lij3
- 0.25*(lij
*rinv
+lij3
*dr
);
857 t2
= -0.5*uij2
- 0.25*sk2_rinv
*uij3
+ 0.25*(uij
*rinv
+uij3
*dr
);
858 t3
= 0.125*(1.0+sk2_rinv
*rinv
)*(-diff2
)+0.25*log_term
*rinv
*rinv
;
860 dadx_val
= (dlij
*t1
+t2
+t3
)*rinv
; /* rb2 is moved to chainrule */
862 born
->gpol_hct_work
[aj
] += 0.5*tmp
;
869 fr
->dadx
[n
++] = dadx_val
;
872 born
->gpol_hct_work
[ai
] += sum_ai
;
876 /* Parallel summations */
877 if (DOMAINDECOMP(cr
))
879 dd_atom_sum_real(cr
->dd
, born
->gpol_hct_work
);
882 for (i
= 0; i
< fr
->natoms_force
; i
++) /* PELA born->nr */
884 if (born
->use
[i
] != 0)
886 rai
= top
->atomtypes
.gb_radius
[md
->typeA
[i
]];
890 sum_ai
= rai
* born
->gpol_hct_work
[i
];
891 sum_ai2
= sum_ai
* sum_ai
;
892 sum_ai3
= sum_ai2
* sum_ai
;
894 tsum
= tanh(born
->obc_alpha
*sum_ai
-born
->obc_beta
*sum_ai2
+born
->obc_gamma
*sum_ai3
);
895 born
->bRad
[i
] = rai_inv
- tsum
*rai_inv2
;
896 born
->bRad
[i
] = 1.0 / born
->bRad
[i
];
898 fr
->invsqrta
[i
] = gmx_invsqrt(born
->bRad
[i
]);
900 tchain
= rai
* (born
->obc_alpha
-2*born
->obc_beta
*sum_ai
+3*born
->obc_gamma
*sum_ai2
);
901 born
->drobc
[i
] = (1.0-tsum
*tsum
)*tchain
*rai_inv2
;
905 /* Extra (local) communication required for DD */
906 if (DOMAINDECOMP(cr
))
908 dd_atom_spread_real(cr
->dd
, born
->bRad
);
909 dd_atom_spread_real(cr
->dd
, fr
->invsqrta
);
910 dd_atom_spread_real(cr
->dd
, born
->drobc
);
919 int calc_gb_rad(t_commrec
*cr
, t_forcerec
*fr
, t_inputrec
*ir
, gmx_localtop_t
*top
,
920 rvec x
[], t_nblist
*nl
, gmx_genborn_t
*born
, t_mdatoms
*md
, t_nrnb
*nrnb
)
925 if (fr
->bAllvsAll
&& fr
->dadx
== NULL
)
927 /* We might need up to 8 atoms of padding before and after,
928 * and another 4 units to guarantee SSE alignment.
930 fr
->nalloc_dadx
= 2*(md
->homenr
+12)*(md
->nr
/2+1+12);
931 snew(fr
->dadx_rawptr
, fr
->nalloc_dadx
);
932 fr
->dadx
= (real
*) (((size_t) fr
->dadx_rawptr
+ 16) & (~((size_t) 15)));
936 /* In the SSE-enabled gb-loops, when writing to dadx, we
937 * always write 2*4 elements at a time, even in the case with only
938 * 1-3 j particles, where we only really need to write 2*(1-3)
939 * elements. This is because we want dadx to be aligned to a 16-
940 * byte boundary, and being able to use _mm_store/load_ps
942 ndadx
= 2 * (nl
->nrj
+ 3*nl
->nri
);
944 /* First, reallocate the dadx array, we need 3 extra for SSE */
945 if (ndadx
+ 3 > fr
->nalloc_dadx
)
947 fr
->nalloc_dadx
= over_alloc_large(ndadx
) + 3;
948 srenew(fr
->dadx_rawptr
, fr
->nalloc_dadx
);
949 fr
->dadx
= (real
*) (((size_t) fr
->dadx_rawptr
+ 16) & (~((size_t) 15)));
955 cnt
= md
->homenr
*(md
->nr
/2+1);
957 if (ir
->gb_algorithm
== egbSTILL
)
959 genborn_allvsall_calc_still_radii(fr
, md
, born
, top
, x
[0], &fr
->AllvsAll_workgb
);
960 /* 13 flops in outer loop, 47 flops in inner loop */
961 inc_nrnb(nrnb
, eNR_BORN_AVA_RADII_STILL
, md
->homenr
*13+cnt
*47);
963 else if (ir
->gb_algorithm
== egbHCT
|| ir
->gb_algorithm
== egbOBC
)
965 genborn_allvsall_calc_hct_obc_radii(fr
, md
, born
, ir
->gb_algorithm
, top
, x
[0], &fr
->AllvsAll_workgb
);
966 /* 24 flops in outer loop, 183 in inner */
967 inc_nrnb(nrnb
, eNR_BORN_AVA_RADII_HCT_OBC
, md
->homenr
*24+cnt
*183);
971 gmx_fatal(FARGS
, "Bad gb algorithm for all-vs-all interactions");
976 /* Switch for determining which algorithm to use for Born radii calculation */
979 switch (ir
->gb_algorithm
)
982 calc_gb_rad_still(cr
, fr
, top
, x
, nl
, born
, md
);
985 calc_gb_rad_hct(cr
, fr
, top
, x
, nl
, born
, md
);
988 calc_gb_rad_obc(cr
, fr
, top
, x
, nl
, born
, md
);
992 gmx_fatal(FARGS
, "Unknown double precision algorithm for Born radii calculation: %d", ir
->gb_algorithm
);
997 switch (ir
->gb_algorithm
)
1000 calc_gb_rad_still(cr
, fr
, top
, x
, nl
, born
, md
);
1003 calc_gb_rad_hct(cr
, fr
, top
, x
, nl
, born
, md
);
1006 calc_gb_rad_obc(cr
, fr
, top
, x
, nl
, born
, md
);
1010 gmx_fatal(FARGS
, "Unknown algorithm for Born radii calculation: %d", ir
->gb_algorithm
);
1013 #endif /* Double or single precision */
1015 if (fr
->bAllvsAll
== FALSE
)
1017 switch (ir
->gb_algorithm
)
1020 /* 17 flops per outer loop iteration, 47 flops per inner loop */
1021 inc_nrnb(nrnb
, eNR_BORN_RADII_STILL
, nl
->nri
*17+nl
->nrj
*47);
1025 /* 61 (assuming 10 for tanh) flops for outer loop iteration, 183 flops per inner loop */
1026 inc_nrnb(nrnb
, eNR_BORN_RADII_HCT_OBC
, nl
->nri
*61+nl
->nrj
*183);
1039 real
gb_bonds_tab(rvec x
[], rvec f
[], rvec fshift
[], real
*charge
, real
*p_gbtabscale
,
1040 real
*invsqrta
, real
*dvda
, real
*GBtab
, t_idef
*idef
, real epsilon_r
,
1041 real gb_epsilon_solvent
, real facel
, const t_pbc
*pbc
, const t_graph
*graph
)
1043 int i
, j
, n0
, m
, nnn
, ai
, aj
;
1049 real isaprod
, qq
, gbscale
, gbtabscale
, Y
, F
, Geps
, Heps2
, Fp
, VV
, FF
, rt
, eps
, eps2
;
1050 real vgb
, fgb
, fijC
, dvdatmp
, fscal
;
1056 t_iatom
*forceatoms
;
1058 /* Scale the electrostatics by gb_epsilon_solvent */
1059 facel
= facel
* ((1.0/epsilon_r
) - 1.0/gb_epsilon_solvent
);
1061 gbtabscale
= *p_gbtabscale
;
1064 for (j
= F_GB12
; j
<= F_GB14
; j
++)
1066 forceatoms
= idef
->il
[j
].iatoms
;
1068 for (i
= 0; i
< idef
->il
[j
].nr
; )
1070 /* To avoid reading in the interaction type, we just increment i to pass over
1071 * the types in the forceatoms array, this saves some memory accesses
1074 ai
= forceatoms
[i
++];
1075 aj
= forceatoms
[i
++];
1077 ki
= pbc_rvec_sub(pbc
, x
[ai
], x
[aj
], dx
);
1078 rsq11
= iprod(dx
, dx
);
1080 isai
= invsqrta
[ai
];
1081 iq
= (-1)*facel
*charge
[ai
];
1083 rinv11
= gmx_invsqrt(rsq11
);
1084 isaj
= invsqrta
[aj
];
1085 isaprod
= isai
*isaj
;
1086 qq
= isaprod
*iq
*charge
[aj
];
1087 gbscale
= isaprod
*gbtabscale
;
1090 n0
= static_cast<int>(rt
);
1096 Geps
= eps
*GBtab
[nnn
+2];
1097 Heps2
= eps2
*GBtab
[nnn
+3];
1100 FF
= Fp
+Geps
+2.0*Heps2
;
1102 fijC
= qq
*FF
*gbscale
;
1103 dvdatmp
= -(vgb
+fijC
*r
)*0.5;
1104 dvda
[aj
] = dvda
[aj
] + dvdatmp
*isaj
*isaj
;
1105 dvda
[ai
] = dvda
[ai
] + dvdatmp
*isai
*isai
;
1106 vctot
= vctot
+ vgb
;
1107 fgb
= -(fijC
)*rinv11
;
1111 ivec_sub(SHIFT_IVEC(graph
, ai
), SHIFT_IVEC(graph
, aj
), dt
);
1115 for (m
= 0; (m
< DIM
); m
++) /* 15 */
1120 fshift
[ki
][m
] += fscal
;
1121 fshift
[CENTRAL
][m
] -= fscal
;
1129 real
calc_gb_selfcorrections(t_commrec
*cr
, int natoms
,
1130 real
*charge
, gmx_genborn_t
*born
, real
*dvda
, double facel
)
1132 int i
, ai
, at0
, at1
;
1133 real rai
, e
, derb
, q
, q2
, fi
, rai_inv
, vtot
;
1135 if (DOMAINDECOMP(cr
))
1138 at1
= cr
->dd
->nat_home
;
1147 /* Scale the electrostatics by gb_epsilon_solvent */
1148 facel
= facel
* ((1.0/born
->epsilon_r
) - 1.0/born
->gb_epsilon_solvent
);
1152 /* Apply self corrections */
1153 for (i
= at0
; i
< at1
; i
++)
1157 if (born
->use
[ai
] == 1)
1159 rai
= born
->bRad
[ai
];
1165 derb
= 0.5*e
*rai_inv
*rai_inv
;
1166 dvda
[ai
] += derb
*rai
;
1175 real
calc_gb_nonpolar(t_commrec
*cr
, t_forcerec
*fr
, int natoms
, gmx_genborn_t
*born
, gmx_localtop_t
*top
,
1176 real
*dvda
, t_mdatoms
*md
)
1178 int ai
, i
, at0
, at1
;
1179 real e
, es
, rai
, term
, probe
, tmp
, factor
;
1180 real rbi_inv
, rbi_inv2
;
1182 if (DOMAINDECOMP(cr
))
1185 at1
= cr
->dd
->nat_home
;
1193 /* factor is the surface tension */
1194 factor
= born
->sa_surface_tension
;
1200 for (i
= at0
; i
< at1
; i
++)
1204 if (born
->use
[ai
] == 1)
1206 rai
= top
->atomtypes
.gb_radius
[md
->typeA
[ai
]];
1207 rbi_inv
= fr
->invsqrta
[ai
];
1208 rbi_inv2
= rbi_inv
* rbi_inv
;
1209 tmp
= (rai
*rbi_inv2
)*(rai
*rbi_inv2
);
1211 e
= factor
*term
*(rai
+probe
)*(rai
+probe
)*tmp
;
1212 dvda
[ai
] = dvda
[ai
] - 6*e
*rbi_inv2
;
1222 real
calc_gb_chainrule(int natoms
, t_nblist
*nl
, real
*dadx
, real
*dvda
, rvec x
[], rvec t
[], rvec fshift
[],
1223 rvec shift_vec
[], int gb_algorithm
, gmx_genborn_t
*born
)
1225 int i
, k
, n
, ai
, aj
, nj0
, nj1
, n0
, n1
;
1228 real fgb
, rbi
, fix1
, fiy1
, fiz1
;
1229 real ix1
, iy1
, iz1
, jx1
, jy1
, jz1
, dx11
, dy11
, dz11
;
1230 real tx
, ty
, tz
, rbai
, rbaj
, fgb_ai
;
1239 if (gb_algorithm
== egbSTILL
)
1241 for (i
= n0
; i
< n1
; i
++)
1243 rbi
= born
->bRad
[i
];
1244 rb
[i
] = (2 * rbi
* rbi
* dvda
[i
])/ONE_4PI_EPS0
;
1247 else if (gb_algorithm
== egbHCT
)
1249 for (i
= n0
; i
< n1
; i
++)
1251 rbi
= born
->bRad
[i
];
1252 rb
[i
] = rbi
* rbi
* dvda
[i
];
1255 else if (gb_algorithm
== egbOBC
)
1257 for (i
= n0
; i
< n1
; i
++)
1259 rbi
= born
->bRad
[i
];
1260 rb
[i
] = rbi
* rbi
* born
->drobc
[i
] * dvda
[i
];
1264 for (i
= 0; i
< nl
->nri
; i
++)
1268 nj0
= nl
->jindex
[i
];
1269 nj1
= nl
->jindex
[i
+1];
1271 /* Load shifts for this list */
1272 shift
= nl
->shift
[i
];
1273 shX
= shift_vec
[shift
][0];
1274 shY
= shift_vec
[shift
][1];
1275 shZ
= shift_vec
[shift
][2];
1277 /* Load atom i coordinates, add shift vectors */
1278 ix1
= shX
+ x
[ai
][0];
1279 iy1
= shY
+ x
[ai
][1];
1280 iz1
= shZ
+ x
[ai
][2];
1288 for (k
= nj0
; k
< nj1
&& nl
->jjnr
[k
] >= 0; k
++)
1302 fgb
= rbai
*dadx
[n
++];
1303 fgb_ai
= rbaj
*dadx
[n
++];
1305 /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1316 /* Update force on atom aj */
1317 t
[aj
][0] = t
[aj
][0] - tx
;
1318 t
[aj
][1] = t
[aj
][1] - ty
;
1319 t
[aj
][2] = t
[aj
][2] - tz
;
1322 /* Update force and shift forces on atom ai */
1323 t
[ai
][0] = t
[ai
][0] + fix1
;
1324 t
[ai
][1] = t
[ai
][1] + fiy1
;
1325 t
[ai
][2] = t
[ai
][2] + fiz1
;
1327 fshift
[shift
][0] = fshift
[shift
][0] + fix1
;
1328 fshift
[shift
][1] = fshift
[shift
][1] + fiy1
;
1329 fshift
[shift
][2] = fshift
[shift
][2] + fiz1
;
1338 calc_gb_forces(t_commrec
*cr
, t_mdatoms
*md
, gmx_genborn_t
*born
, gmx_localtop_t
*top
,
1339 rvec x
[], rvec f
[], t_forcerec
*fr
, t_idef
*idef
, int gb_algorithm
, int sa_algorithm
, t_nrnb
*nrnb
,
1340 const t_pbc
*pbc
, const t_graph
*graph
, gmx_enerdata_t
*enerd
)
1345 const t_pbc
*pbc_null
;
1356 if (sa_algorithm
== esaAPPROX
)
1358 /* Do a simple ACE type approximation for the non-polar solvation */
1359 enerd
->term
[F_NPSOLVATION
] += calc_gb_nonpolar(cr
, fr
, born
->nr
, born
, top
, fr
->dvda
, md
);
1362 /* Calculate the bonded GB-interactions using either table or analytical formula */
1363 enerd
->term
[F_GBPOL
] += gb_bonds_tab(x
, f
, fr
->fshift
, md
->chargeA
, &(fr
->gbtabscale
),
1364 fr
->invsqrta
, fr
->dvda
, fr
->gbtab
.data
, idef
, born
->epsilon_r
, born
->gb_epsilon_solvent
, fr
->epsfac
, pbc_null
, graph
);
1366 /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1367 enerd
->term
[F_GBPOL
] += calc_gb_selfcorrections(cr
, born
->nr
, md
->chargeA
, born
, fr
->dvda
, fr
->epsfac
);
1369 /* If parallel, sum the derivative of the potential w.r.t the born radii */
1370 if (DOMAINDECOMP(cr
))
1372 dd_atom_sum_real(cr
->dd
, fr
->dvda
);
1373 dd_atom_spread_real(cr
->dd
, fr
->dvda
);
1378 genborn_allvsall_calc_chainrule(fr
, md
, born
, x
[0], f
[0], gb_algorithm
, fr
->AllvsAll_workgb
);
1379 cnt
= md
->homenr
*(md
->nr
/2+1);
1380 /* 9 flops for outer loop, 15 for inner */
1381 inc_nrnb(nrnb
, eNR_BORN_AVA_CHAINRULE
, md
->homenr
*9+cnt
*15);
1385 calc_gb_chainrule(fr
->natoms_force
, &(fr
->gblist
), fr
->dadx
, fr
->dvda
,
1386 x
, f
, fr
->fshift
, fr
->shift_vec
, gb_algorithm
, born
);
1390 /* 9 flops for outer loop, 15 for inner */
1391 inc_nrnb(nrnb
, eNR_BORN_CHAINRULE
, fr
->gblist
.nri
*9+fr
->gblist
.nrj
*15);
1395 static void add_j_to_gblist(gbtmpnbl_t
*list
, int aj
)
1397 if (list
->naj
>= list
->aj_nalloc
)
1399 list
->aj_nalloc
= over_alloc_large(list
->naj
+1);
1400 srenew(list
->aj
, list
->aj_nalloc
);
1403 list
->aj
[list
->naj
++] = aj
;
1406 static gbtmpnbl_t
*find_gbtmplist(struct gbtmpnbls
*lists
, int shift
)
1410 /* Search the list with the same shift, if there is one */
1412 while (ind
< lists
->nlist
&& shift
!= lists
->list
[ind
].shift
)
1416 if (ind
== lists
->nlist
)
1418 if (lists
->nlist
== lists
->list_nalloc
)
1420 lists
->list_nalloc
++;
1421 srenew(lists
->list
, lists
->list_nalloc
);
1422 for (i
= lists
->nlist
; i
< lists
->list_nalloc
; i
++)
1424 lists
->list
[i
].aj
= NULL
;
1425 lists
->list
[i
].aj_nalloc
= 0;
1430 lists
->list
[lists
->nlist
].shift
= shift
;
1431 lists
->list
[lists
->nlist
].naj
= 0;
1435 return &lists
->list
[ind
];
1438 static void add_bondeds_to_gblist(t_ilist
*il
,
1439 gmx_bool bMolPBC
, t_pbc
*pbc
, t_graph
*g
, rvec
*x
,
1440 struct gbtmpnbls
*nls
)
1442 int ind
, j
, ai
, aj
, found
;
1447 for (ind
= 0; ind
< il
->nr
; ind
+= 3)
1449 ai
= il
->iatoms
[ind
+1];
1450 aj
= il
->iatoms
[ind
+2];
1452 int shift
= CENTRAL
;
1455 rvec_sub(x
[ai
], x
[aj
], dx
);
1456 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
1457 shift
= IVEC2IS(dt
);
1461 shift
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], dx
);
1464 /* Find the list for this shift or create one */
1465 list
= find_gbtmplist(&nls
[ai
], shift
);
1469 /* So that we do not add the same bond twice.
1470 * This happens with some constraints between 1-3 atoms
1471 * that are in the bond-list but should not be in the GB nb-list */
1472 for (j
= 0; j
< list
->naj
; j
++)
1474 if (list
->aj
[j
] == aj
)
1484 gmx_incons("ai == aj");
1487 add_j_to_gblist(list
, aj
);
1493 int make_gb_nblist(t_commrec
*cr
, int gb_algorithm
,
1494 rvec x
[], matrix box
,
1495 t_forcerec
*fr
, t_idef
*idef
, t_graph
*graph
, gmx_genborn_t
*born
)
1497 int i
, j
, k
, n
, nj0
, nj1
, ai
, shift
, s
;
1501 struct gbtmpnbls
*nls
;
1502 gbtmpnbl_t
*list
= NULL
;
1504 set_pbc(&pbc
, fr
->ePBC
, box
);
1505 nls
= born
->nblist_work
;
1507 for (i
= 0; i
< born
->nr
; i
++)
1514 set_pbc_dd(&pbc
, fr
->ePBC
, cr
->dd
, TRUE
, box
);
1517 switch (gb_algorithm
)
1521 /* Loop over 1-2, 1-3 and 1-4 interactions */
1522 for (j
= F_GB12
; j
<= F_GB14
; j
++)
1524 add_bondeds_to_gblist(&idef
->il
[j
], fr
->bMolPBC
, &pbc
, graph
, x
, nls
);
1528 /* Loop over 1-4 interactions */
1529 add_bondeds_to_gblist(&idef
->il
[F_GB14
], fr
->bMolPBC
, &pbc
, graph
, x
, nls
);
1532 gmx_incons("Unknown GB algorithm");
1535 /* Loop over the VDWQQ and VDW nblists to set up the nonbonded part of the GB list */
1536 for (n
= 0; (n
< fr
->nnblists
); n
++)
1538 for (i
= 0; (i
< eNL_NR
); i
++)
1540 nblist
= &(fr
->nblists
[n
].nlist_sr
[i
]);
1542 if (nblist
->nri
> 0 && (i
== eNL_VDWQQ
|| i
== eNL_QQ
))
1544 for (j
= 0; j
< nblist
->nri
; j
++)
1546 ai
= nblist
->iinr
[j
];
1547 shift
= nblist
->shift
[j
];
1549 /* Find the list for this shift or create one */
1550 list
= find_gbtmplist(&nls
[ai
], shift
);
1552 nj0
= nblist
->jindex
[j
];
1553 nj1
= nblist
->jindex
[j
+1];
1555 /* Add all the j-atoms in the non-bonded list to the GB list */
1556 for (k
= nj0
; k
< nj1
; k
++)
1558 add_j_to_gblist(list
, nblist
->jjnr
[k
]);
1565 /* Zero out some counters */
1569 fr
->gblist
.jindex
[0] = fr
->gblist
.nri
;
1571 for (i
= 0; i
< fr
->natoms_force
; i
++)
1573 for (s
= 0; s
< nls
[i
].nlist
; s
++)
1575 list
= &nls
[i
].list
[s
];
1577 /* Only add those atoms that actually have neighbours */
1578 if (born
->use
[i
] != 0)
1580 fr
->gblist
.iinr
[fr
->gblist
.nri
] = i
;
1581 fr
->gblist
.shift
[fr
->gblist
.nri
] = list
->shift
;
1584 for (k
= 0; k
< list
->naj
; k
++)
1586 /* Memory allocation for jjnr */
1587 if (fr
->gblist
.nrj
>= fr
->gblist
.maxnrj
)
1589 fr
->gblist
.maxnrj
+= over_alloc_large(fr
->gblist
.maxnrj
);
1593 fprintf(debug
, "Increasing GB neighbourlist j size to %d\n", fr
->gblist
.maxnrj
);
1596 srenew(fr
->gblist
.jjnr
, fr
->gblist
.maxnrj
);
1600 if (i
== list
->aj
[k
])
1602 gmx_incons("i == list->aj[k]");
1604 fr
->gblist
.jjnr
[fr
->gblist
.nrj
++] = list
->aj
[k
];
1607 fr
->gblist
.jindex
[fr
->gblist
.nri
] = fr
->gblist
.nrj
;
1615 void make_local_gb(const t_commrec
*cr
, gmx_genborn_t
*born
, int gb_algorithm
)
1618 gmx_domdec_t
*dd
= NULL
;
1620 if (DOMAINDECOMP(cr
))
1628 /* Single node, just copy pointers and return */
1629 if (gb_algorithm
== egbSTILL
)
1631 born
->gpol
= born
->gpol_globalindex
;
1632 born
->vsolv
= born
->vsolv_globalindex
;
1633 born
->gb_radius
= born
->gb_radius_globalindex
;
1637 born
->param
= born
->param_globalindex
;
1638 born
->gb_radius
= born
->gb_radius_globalindex
;
1641 born
->use
= born
->use_globalindex
;
1646 /* Reallocation of local arrays if necessary */
1647 /* fr->natoms_force is equal to dd->nat_tot */
1648 if (DOMAINDECOMP(cr
) && dd
->nat_tot
> born
->nalloc
)
1652 nalloc
= dd
->nat_tot
;
1654 /* Arrays specific to different gb algorithms */
1655 if (gb_algorithm
== egbSTILL
)
1657 srenew(born
->gpol
, nalloc
+3);
1658 srenew(born
->vsolv
, nalloc
+3);
1659 srenew(born
->gb_radius
, nalloc
+3);
1660 for (i
= born
->nalloc
; (i
< nalloc
+3); i
++)
1664 born
->gb_radius
[i
] = 0;
1669 srenew(born
->param
, nalloc
+3);
1670 srenew(born
->gb_radius
, nalloc
+3);
1671 for (i
= born
->nalloc
; (i
< nalloc
+3); i
++)
1674 born
->gb_radius
[i
] = 0;
1678 /* All gb-algorithms use the array for vsites exclusions */
1679 srenew(born
->use
, nalloc
+3);
1680 for (i
= born
->nalloc
; (i
< nalloc
+3); i
++)
1685 born
->nalloc
= nalloc
;
1688 /* With dd, copy algorithm specific arrays */
1689 if (gb_algorithm
== egbSTILL
)
1691 for (i
= at0
; i
< at1
; i
++)
1693 born
->gpol
[i
] = born
->gpol_globalindex
[dd
->gatindex
[i
]];
1694 born
->vsolv
[i
] = born
->vsolv_globalindex
[dd
->gatindex
[i
]];
1695 born
->gb_radius
[i
] = born
->gb_radius_globalindex
[dd
->gatindex
[i
]];
1696 born
->use
[i
] = born
->use_globalindex
[dd
->gatindex
[i
]];
1701 for (i
= at0
; i
< at1
; i
++)
1703 born
->param
[i
] = born
->param_globalindex
[dd
->gatindex
[i
]];
1704 born
->gb_radius
[i
] = born
->gb_radius_globalindex
[dd
->gatindex
[i
]];
1705 born
->use
[i
] = born
->use_globalindex
[dd
->gatindex
[i
]];