?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Prepared t_mdatoms for using vector
[gromacs.git] / src / gromacs / mdlib / genborn.cpp
blob7d9cc0f617658a47b0d6d882355f94c056d367a1
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
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,2017, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37
38 #include "gmxpre.h"
39
40 #include "genborn.h"
41
42 #include <string.h>
43
44 #include <cmath>
45
46 #include <algorithm>
47
48 #include "gromacs/domdec/domdec.h"
49 #include "gromacs/domdec/domdec_struct.h"
50 #include "gromacs/fileio/pdbio.h"
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/mdlib/genborn_allvsall.h"
57 #include "gromacs/mdtypes/commrec.h"
58 #include "gromacs/mdtypes/inputrec.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/mdtypes/nblist.h"
61 #include "gromacs/pbcutil/ishift.h"
62 #include "gromacs/pbcutil/mshift.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/utility/fatalerror.h"
66 #include "gromacs/utility/gmxmpi.h"
67 #include "gromacs/utility/smalloc.h"
68
69
70 typedef struct {
71 int shift;
72 int naj;
73 int *aj;
74 int aj_nalloc;
75 } gbtmpnbl_t;
76
77 typedef struct gbtmpnbls {
78 int nlist;
79 gbtmpnbl_t *list;
80 int list_nalloc;
81 } t_gbtmpnbls;
82
83 /* This function is exactly the same as the one in listed-forces/bonded.cpp. The reason
84 * it is copied here is that the bonded gb-interactions are evaluated
85 * not in calc_bonds, but rather in calc_gb_forces
86 */
87 static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx)
88 {
89 if (pbc)
90 {
91 return pbc_dx_aiuc(pbc, xi, xj, dx);
92 }
93 else
94 {
95 rvec_sub(xi, xj, dx);
96 return CENTRAL;
97 }
98 }
99
100 static int init_gb_nblist(int natoms, t_nblist *nl)
101 {
102 nl->maxnri = natoms*4;
103 nl->maxnrj = 0;
104 nl->nri = 0;
105 nl->nrj = 0;
106 nl->iinr = nullptr;
107 nl->gid = nullptr;
108 nl->shift = nullptr;
109 nl->jindex = nullptr;
110 nl->jjnr = nullptr;
111 /*nl->nltype = nltype;*/
112
113 srenew(nl->iinr, nl->maxnri);
114 srenew(nl->gid, nl->maxnri);
115 srenew(nl->shift, nl->maxnri);
116 srenew(nl->jindex, nl->maxnri+1);
117
118 nl->jindex[0] = 0;
119
120 return 0;
121 }
122
123
124 static int init_gb_still(const t_atomtypes *atype, t_idef *idef, t_atoms *atoms,
125 gmx_genborn_t *born, int natoms)
126 {
127
128 int i, j, m, ia, ib;
129 real r, ri, rj, ri2, rj2, r3, r4, ratio, term, h, doffset;
130
131 real *vsol;
132 real *gp;
133
134 snew(vsol, natoms);
135 snew(gp, natoms);
136 snew(born->gpol_still_work, natoms+3);
137
138 doffset = born->gb_doffset;
139
140 for (i = 0; i < natoms; i++)
141 {
142 born->gpol_globalindex[i] = born->vsolv_globalindex[i] =
143 born->gb_radius_globalindex[i] = 0;
144 }
145
146 /* Compute atomic solvation volumes for Still method */
147 for (i = 0; i < natoms; i++)
148 {
149 ri = atype->gb_radius[atoms->atom[i].type];
150 born->gb_radius_globalindex[i] = ri;
151 r3 = ri*ri*ri;
152 born->vsolv_globalindex[i] = (4*M_PI/3)*r3;
153 }
154
155 for (j = 0; j < idef->il[F_GB12].nr; j += 3)
156 {
157 m = idef->il[F_GB12].iatoms[j];
158 ia = idef->il[F_GB12].iatoms[j+1];
159 ib = idef->il[F_GB12].iatoms[j+2];
160
161 r = 1.01*idef->iparams[m].gb.st;
162
163 ri = atype->gb_radius[atoms->atom[ia].type];
164 rj = atype->gb_radius[atoms->atom[ib].type];
165
166 ri2 = ri*ri;
167 rj2 = rj*rj;
168
169 ratio = (rj2-ri2-r*r)/(2*ri*r);
170 h = ri*(1+ratio);
171 term = (M_PI/3.0)*h*h*(3.0*ri-h);
172
173 born->vsolv_globalindex[ia] -= term;
174
175 ratio = (ri2-rj2-r*r)/(2*rj*r);
176 h = rj*(1+ratio);
177 term = (M_PI/3.0)*h*h*(3.0*rj-h);
178
179 born->vsolv_globalindex[ib] -= term;
180 }
181
182 /* Get the self-, 1-2 and 1-3 polarization energies for analytical Still
183 method */
184 /* Self */
185 for (j = 0; j < natoms; j++)
186 {
187 if (born->use_globalindex[j] == 1)
188 {
189 born->gpol_globalindex[j] = -0.5*ONE_4PI_EPS0/
190 (atype->gb_radius[atoms->atom[j].type]-doffset+STILL_P1);
191 }
192 }
193
194 /* 1-2 */
195 for (j = 0; j < idef->il[F_GB12].nr; j += 3)
196 {
197 m = idef->il[F_GB12].iatoms[j];
198 ia = idef->il[F_GB12].iatoms[j+1];
199 ib = idef->il[F_GB12].iatoms[j+2];
200
201 r = idef->iparams[m].gb.st;
202
203 r4 = r*r*r*r;
204
205 born->gpol_globalindex[ia] = born->gpol_globalindex[ia]+
206 STILL_P2*born->vsolv_globalindex[ib]/r4;
207 born->gpol_globalindex[ib] = born->gpol_globalindex[ib]+
208 STILL_P2*born->vsolv_globalindex[ia]/r4;
209 }
210
211 /* 1-3 */
212 for (j = 0; j < idef->il[F_GB13].nr; j += 3)
213 {
214 m = idef->il[F_GB13].iatoms[j];
215 ia = idef->il[F_GB13].iatoms[j+1];
216 ib = idef->il[F_GB13].iatoms[j+2];
217
218 r = idef->iparams[m].gb.st;
219 r4 = r*r*r*r;
220
221 born->gpol_globalindex[ia] = born->gpol_globalindex[ia]+
222 STILL_P3*born->vsolv_globalindex[ib]/r4;
223 born->gpol_globalindex[ib] = born->gpol_globalindex[ib]+
224 STILL_P3*born->vsolv_globalindex[ia]/r4;
225 }
226
227 sfree(vsol);
228 sfree(gp);
229
230 return 0;
231 }
232
233 /* Initialize all GB datastructs and compute polarization energies */
234 int init_gb(gmx_genborn_t **p_born,
235 t_forcerec *fr, const t_inputrec *ir,
236 const gmx_mtop_t *mtop, int gb_algorithm)
237 {
238 int i, jj, natoms;
239 real rai, sk, doffset;
240
241 t_atoms atoms;
242 gmx_genborn_t *born;
243 gmx_localtop_t *localtop;
244
245 natoms = mtop->natoms;
246
247 atoms = gmx_mtop_global_atoms(mtop);
248 localtop = gmx_mtop_generate_local_top(mtop, ir->efep != efepNO);
249
250 snew(born, 1);
251 *p_born = born;
252
253 born->nr = natoms;
254
255 snew(born->drobc, natoms);
256 snew(born->bRad, natoms);
257
258 /* Allocate memory for the global data arrays */
259 snew(born->param_globalindex, natoms+3);
260 snew(born->gpol_globalindex, natoms+3);
261 snew(born->vsolv_globalindex, natoms+3);
262 snew(born->gb_radius_globalindex, natoms+3);
263 snew(born->use_globalindex, natoms+3);
264
265 snew(fr->invsqrta, natoms);
266 snew(fr->dvda, natoms);
267
268 fr->dadx = nullptr;
269 fr->dadx_rawptr = nullptr;
270 fr->nalloc_dadx = 0;
271 born->gpol_still_work = nullptr;
272 born->gpol_hct_work = nullptr;
273
274 /* snew(born->asurf,natoms); */
275 /* snew(born->dasurf,natoms); */
276
277 /* Initialize the gb neighbourlist */
278 snew(fr->gblist, 1);
279 init_gb_nblist(natoms, fr->gblist);
280
281 /* Do the Vsites exclusions (if any) */
282 for (i = 0; i < natoms; i++)
283 {
284 jj = atoms.atom[i].type;
285 if (mtop->atomtypes.gb_radius[atoms.atom[i].type] > 0)
286 {
287 born->use_globalindex[i] = 1;
288 }
289 else
290 {
291 born->use_globalindex[i] = 0;
292 }
293
294 /* If we have a Vsite, put vs_globalindex[i]=0 */
295 if (C6 (fr->nbfp, fr->ntype, jj, jj) == 0 &&
296 C12(fr->nbfp, fr->ntype, jj, jj) == 0 &&
297 atoms.atom[i].q == 0)
298 {
299 born->use_globalindex[i] = 0;
300 }
301 }
302
303 /* Copy algorithm parameters from inputrecord to local structure */
304 born->obc_alpha = ir->gb_obc_alpha;
305 born->obc_beta = ir->gb_obc_beta;
306 born->obc_gamma = ir->gb_obc_gamma;
307 born->gb_doffset = ir->gb_dielectric_offset;
308 born->gb_epsilon_solvent = ir->gb_epsilon_solvent;
309 born->epsilon_r = ir->epsilon_r;
310
311 doffset = born->gb_doffset;
312
313 /* Set the surface tension */
314 born->sa_surface_tension = ir->sa_surface_tension;
315
316 /* If Still model, initialise the polarisation energies */
317 if (gb_algorithm == egbSTILL)
318 {
319 init_gb_still(&(mtop->atomtypes), &(localtop->idef), &atoms,
320 born, natoms);
321 }
322
323
324 /* If HCT/OBC, precalculate the sk*atype->S_hct factors */
325 else if (gb_algorithm == egbHCT || gb_algorithm == egbOBC)
326 {
327
328 snew(born->gpol_hct_work, natoms+3);
329
330 for (i = 0; i < natoms; i++)
331 {
332 if (born->use_globalindex[i] == 1)
333 {
334 rai = mtop->atomtypes.gb_radius[atoms.atom[i].type]-doffset;
335 sk = rai * mtop->atomtypes.S_hct[atoms.atom[i].type];
336 born->param_globalindex[i] = sk;
337 born->gb_radius_globalindex[i] = rai;
338 }
339 else
340 {
341 born->param_globalindex[i] = 0;
342 born->gb_radius_globalindex[i] = 0;
343 }
344 }
345 }
346
347 /* Allocate memory for work arrays for temporary use */
348 snew(born->work, natoms+4);
349 snew(born->count, natoms);
350 snew(born->nblist_work, natoms);
351
352 /* Domain decomposition specific stuff */
353 born->nalloc = 0;
354
355 return 0;
356 }
357
358
359
360 static int
361 calc_gb_rad_still(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top,
362 rvec x[], t_nblist *nl,
363 gmx_genborn_t *born, t_mdatoms *md)
364 {
365 int i, k, n, nj0, nj1, ai, aj;
366 int shift;
367 real shX, shY, shZ;
368 real gpi, dr2, idr4, rvdw, ratio, ccf, theta, term, rai, raj;
369 real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11;
370 real rinv, idr2, idr6, vaj, dccf, cosq, sinq, prod, gpi2;
371 real factor;
372 real vai, prod_ai, icf4, icf6;
373
374 factor = 0.5*ONE_4PI_EPS0;
375 n = 0;
376
377 for (i = 0; i < born->nr; i++)
378 {
379 born->gpol_still_work[i] = 0;
380 }
381
382 for (i = 0; i < nl->nri; i++)
383 {
384 ai = nl->iinr[i];
385
386 nj0 = nl->jindex[i];
387 nj1 = nl->jindex[i+1];
388
389 /* Load shifts for this list */
390 shift = nl->shift[i];
391 shX = fr->shift_vec[shift][0];
392 shY = fr->shift_vec[shift][1];
393 shZ = fr->shift_vec[shift][2];
394
395 gpi = 0;
396
397 rai = top->atomtypes.gb_radius[md->typeA[ai]];
398 vai = born->vsolv[ai];
399 prod_ai = STILL_P4*vai;
400
401 /* Load atom i coordinates, add shift vectors */
402 ix1 = shX + x[ai][0];
403 iy1 = shY + x[ai][1];
404 iz1 = shZ + x[ai][2];
405
406 for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
407 {
408 aj = nl->jjnr[k];
409 jx1 = x[aj][0];
410 jy1 = x[aj][1];
411 jz1 = x[aj][2];
412
413 dx11 = ix1-jx1;
414 dy11 = iy1-jy1;
415 dz11 = iz1-jz1;
416
417 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
418 rinv = gmx::invsqrt(dr2);
419 idr2 = rinv*rinv;
420 idr4 = idr2*idr2;
421 idr6 = idr4*idr2;
422
423 raj = top->atomtypes.gb_radius[md->typeA[aj]];
424
425 rvdw = rai + raj;
426
427 ratio = dr2 / (rvdw * rvdw);
428 vaj = born->vsolv[aj];
429
430 if (ratio > STILL_P5INV)
431 {
432 ccf = 1.0;
433 dccf = 0.0;
434 }
435 else
436 {
437 theta = ratio*STILL_PIP5;
438 cosq = cos(theta);
439 term = 0.5*(1.0-cosq);
440 ccf = term*term;
441 sinq = 1.0 - cosq*cosq;
442 dccf = 2.0*term*sinq*gmx::invsqrt(sinq)*theta;
443 }
444
445 prod = STILL_P4*vaj;
446 icf4 = ccf*idr4;
447 icf6 = (4*ccf-dccf)*idr6;
448 born->gpol_still_work[aj] += prod_ai*icf4;
449 gpi = gpi+prod*icf4;
450
451 /* Save ai->aj and aj->ai chain rule terms */
452 fr->dadx[n++] = prod*icf6;
453 fr->dadx[n++] = prod_ai*icf6;
454 }
455 born->gpol_still_work[ai] += gpi;
456 }
457
458 /* Parallel summations */
459 if (DOMAINDECOMP(cr))
460 {
461 dd_atom_sum_real(cr->dd, born->gpol_still_work);
462 }
463
464 /* Calculate the radii */
465 for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
466 {
467 if (born->use[i] != 0)
468 {
469 gpi = born->gpol[i]+born->gpol_still_work[i];
470 gpi2 = gpi * gpi;
471 born->bRad[i] = factor*gmx::invsqrt(gpi2);
472 fr->invsqrta[i] = gmx::invsqrt(born->bRad[i]);
473 }
474 }
475
476 /* Extra communication required for DD */
477 if (DOMAINDECOMP(cr))
478 {
479 dd_atom_spread_real(cr->dd, born->bRad);
480 dd_atom_spread_real(cr->dd, fr->invsqrta);
481 }
482
483 return 0;
484
485 }
486
487
488 static int
489 calc_gb_rad_hct(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top,
490 rvec x[], t_nblist *nl,
491 gmx_genborn_t *born, t_mdatoms *md)
492 {
493 int i, k, n, ai, aj, nj0, nj1;
494 int shift;
495 real shX, shY, shZ;
496 real rai, raj, dr2, dr, sk, sk_ai, sk2, sk2_ai, lij, uij, diff2, tmp, sum_ai;
497 real rad, min_rad, rinv, rai_inv;
498 real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11;
499 real lij2, uij2, lij3, uij3, t1, t2, t3;
500 real lij_inv, dlij, sk2_rinv, prod, log_term;
501 real doffset, raj_inv, dadx_val;
502 real *gb_radius;
503
504 doffset = born->gb_doffset;
505 gb_radius = born->gb_radius;
506
507 for (i = 0; i < born->nr; i++)
508 {
509 born->gpol_hct_work[i] = 0;
510 }
511
512 /* Keep the compiler happy */
513 n = 0;
514
515 for (i = 0; i < nl->nri; i++)
516 {
517 ai = nl->iinr[i];
518
519 nj0 = nl->jindex[i];
520 nj1 = nl->jindex[i+1];
521
522 /* Load shifts for this list */
523 shift = nl->shift[i];
524 shX = fr->shift_vec[shift][0];
525 shY = fr->shift_vec[shift][1];
526 shZ = fr->shift_vec[shift][2];
527
528 rai = gb_radius[ai];
529 rai_inv = 1.0/rai;
530
531 sk_ai = born->param[ai];
532 sk2_ai = sk_ai*sk_ai;
533
534 /* Load atom i coordinates, add shift vectors */
535 ix1 = shX + x[ai][0];
536 iy1 = shY + x[ai][1];
537 iz1 = shZ + x[ai][2];
538
539 sum_ai = 0;
540
541 for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
542 {
543 aj = nl->jjnr[k];
544
545 jx1 = x[aj][0];
546 jy1 = x[aj][1];
547 jz1 = x[aj][2];
548
549 dx11 = ix1 - jx1;
550 dy11 = iy1 - jy1;
551 dz11 = iz1 - jz1;
552
553 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
554 rinv = gmx::invsqrt(dr2);
555 dr = rinv*dr2;
556
557 sk = born->param[aj];
558 raj = gb_radius[aj];
559
560 /* aj -> ai interaction */
561 if (rai < dr+sk)
562 {
563 lij = 1.0/(dr-sk);
564 dlij = 1.0;
565
566 if (rai > dr-sk)
567 {
568 lij = rai_inv;
569 dlij = 0.0;
570 }
571
572 lij2 = lij*lij;
573 lij3 = lij2*lij;
574
575 uij = 1.0/(dr+sk);
576 uij2 = uij*uij;
577 uij3 = uij2*uij;
578
579 diff2 = uij2-lij2;
580
581 lij_inv = gmx::invsqrt(lij2);
582 sk2 = sk*sk;
583 sk2_rinv = sk2*rinv;
584 prod = 0.25*sk2_rinv;
585
586 log_term = std::log(uij*lij_inv);
587
588 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
589 prod*(-diff2);
590
591 if (rai < sk-dr)
592 {
593 tmp = tmp + 2.0 * (rai_inv-lij);
594 }
595
596 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
597 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
598 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
599
600 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
601 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */
602 /* rb2 is moved to chainrule */
603
604 sum_ai += 0.5*tmp;
605 }
606 else
607 {
608 dadx_val = 0.0;
609 }
610 fr->dadx[n++] = dadx_val;
611
612
613 /* ai -> aj interaction */
614 if (raj < dr + sk_ai)
615 {
616 lij = 1.0/(dr-sk_ai);
617 dlij = 1.0;
618 raj_inv = 1.0/raj;
619
620 if (raj > dr-sk_ai)
621 {
622 lij = raj_inv;
623 dlij = 0.0;
624 }
625
626 lij2 = lij * lij;
627 lij3 = lij2 * lij;
628
629 uij = 1.0/(dr+sk_ai);
630 uij2 = uij * uij;
631 uij3 = uij2 * uij;
632
633 diff2 = uij2-lij2;
634
635 lij_inv = gmx::invsqrt(lij2);
636 sk2 = sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
637 sk2_rinv = sk2*rinv;
638 prod = 0.25 * sk2_rinv;
639
640 /* log_term = table_log(uij*lij_inv,born->log_table,
641 LOG_TABLE_ACCURACY); */
642 log_term = std::log(uij*lij_inv);
643
644 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
645 prod*(-diff2);
646
647 if (raj < sk_ai-dr)
648 {
649 tmp = tmp + 2.0 * (raj_inv-lij);
650 }
651
652 /* duij = 1.0 */
653 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
654 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
655 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
656
657 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
658 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */ /* rb2 is moved to chainrule */
659
660 born->gpol_hct_work[aj] += 0.5*tmp;
661 }
662 else
663 {
664 dadx_val = 0.0;
665 }
666 fr->dadx[n++] = dadx_val;
667 }
668
669 born->gpol_hct_work[ai] += sum_ai;
670 }
671
672 /* Parallel summations */
673 if (DOMAINDECOMP(cr))
674 {
675 dd_atom_sum_real(cr->dd, born->gpol_hct_work);
676 }
677
678 for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
679 {
680 if (born->use[i] != 0)
681 {
682 rai = top->atomtypes.gb_radius[md->typeA[i]]-doffset;
683 sum_ai = 1.0/rai - born->gpol_hct_work[i];
684 min_rad = rai + doffset;
685 rad = 1.0/sum_ai;
686
687 born->bRad[i] = std::max(rad, min_rad);
688 fr->invsqrta[i] = gmx::invsqrt(born->bRad[i]);
689 }
690 }
691
692 /* Extra communication required for DD */
693 if (DOMAINDECOMP(cr))
694 {
695 dd_atom_spread_real(cr->dd, born->bRad);
696 dd_atom_spread_real(cr->dd, fr->invsqrta);
697 }
698
699
700 return 0;
701 }
702
703 static int
704 calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top,
705 rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md)
706 {
707 int i, k, ai, aj, nj0, nj1, n;
708 int shift;
709 real shX, shY, shZ;
710 real rai, raj, dr2, dr, sk, sk2, lij, uij, diff2, tmp, sum_ai;
711 real sum_ai2, sum_ai3, tsum, tchain, rinv, rai_inv, lij_inv, rai_inv2;
712 real log_term, prod, sk2_rinv, sk_ai, sk2_ai;
713 real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11;
714 real lij2, uij2, lij3, uij3, dlij, t1, t2, t3;
715 real doffset, raj_inv, dadx_val;
716 real *gb_radius;
717
718 /* Keep the compiler happy */
719 n = 0;
720
721 doffset = born->gb_doffset;
722 gb_radius = born->gb_radius;
723
724 for (i = 0; i < born->nr; i++)
725 {
726 born->gpol_hct_work[i] = 0;
727 }
728
729 for (i = 0; i < nl->nri; i++)
730 {
731 ai = nl->iinr[i];
732
733 nj0 = nl->jindex[i];
734 nj1 = nl->jindex[i+1];
735
736 /* Load shifts for this list */
737 shift = nl->shift[i];
738 shX = fr->shift_vec[shift][0];
739 shY = fr->shift_vec[shift][1];
740 shZ = fr->shift_vec[shift][2];
741
742 rai = gb_radius[ai];
743 rai_inv = 1.0/rai;
744
745 sk_ai = born->param[ai];
746 sk2_ai = sk_ai*sk_ai;
747
748 /* Load atom i coordinates, add shift vectors */
749 ix1 = shX + x[ai][0];
750 iy1 = shY + x[ai][1];
751 iz1 = shZ + x[ai][2];
752
753 sum_ai = 0;
754
755 for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
756 {
757 aj = nl->jjnr[k];
758
759 jx1 = x[aj][0];
760 jy1 = x[aj][1];
761 jz1 = x[aj][2];
762
763 dx11 = ix1 - jx1;
764 dy11 = iy1 - jy1;
765 dz11 = iz1 - jz1;
766
767 dr2 = dx11*dx11+dy11*dy11+dz11*dz11;
768 rinv = gmx::invsqrt(dr2);
769 dr = dr2*rinv;
770
771 /* sk is precalculated in init_gb() */
772 sk = born->param[aj];
773 raj = gb_radius[aj];
774
775 /* aj -> ai interaction */
776 if (rai < dr+sk)
777 {
778 lij = 1.0/(dr-sk);
779 dlij = 1.0;
780
781 if (rai > dr-sk)
782 {
783 lij = rai_inv;
784 dlij = 0.0;
785 }
786
787 uij = 1.0/(dr+sk);
788 lij2 = lij * lij;
789 lij3 = lij2 * lij;
790 uij2 = uij * uij;
791 uij3 = uij2 * uij;
792
793 diff2 = uij2-lij2;
794
795 lij_inv = gmx::invsqrt(lij2);
796 sk2 = sk*sk;
797 sk2_rinv = sk2*rinv;
798 prod = 0.25*sk2_rinv;
799
800 log_term = std::log(uij*lij_inv);
801
802 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
803
804 if (rai < sk-dr)
805 {
806 tmp = tmp + 2.0 * (rai_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
816 sum_ai += 0.5*tmp;
817 }
818 else
819 {
820 dadx_val = 0.0;
821 }
822 fr->dadx[n++] = dadx_val;
823
824 /* ai -> aj interaction */
825 if (raj < dr + sk_ai)
826 {
827 lij = 1.0/(dr-sk_ai);
828 dlij = 1.0;
829 raj_inv = 1.0/raj;
830
831 if (raj > dr-sk_ai)
832 {
833 lij = raj_inv;
834 dlij = 0.0;
835 }
836
837 lij2 = lij * lij;
838 lij3 = lij2 * lij;
839
840 uij = 1.0/(dr+sk_ai);
841 uij2 = uij * uij;
842 uij3 = uij2 * uij;
843
844 diff2 = uij2-lij2;
845
846 lij_inv = gmx::invsqrt(lij2);
847 sk2 = sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
848 sk2_rinv = sk2*rinv;
849 prod = 0.25 * sk2_rinv;
850
851 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
852 log_term = std::log(uij*lij_inv);
853
854 tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
855
856 if (raj < sk_ai-dr)
857 {
858 tmp = tmp + 2.0 * (raj_inv-lij);
859 }
860
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
867 born->gpol_hct_work[aj] += 0.5*tmp;
868
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
881 /* Parallel summations */
882 if (DOMAINDECOMP(cr))
883 {
884 dd_atom_sum_real(cr->dd, born->gpol_hct_work);
885 }
886
887 for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
888 {
889 if (born->use[i] != 0)
890 {
891 rai = top->atomtypes.gb_radius[md->typeA[i]];
892 rai_inv2 = 1.0/rai;
893 rai = rai-doffset;
894 rai_inv = 1.0/rai;
895 sum_ai = rai * born->gpol_hct_work[i];
896 sum_ai2 = sum_ai * sum_ai;
897 sum_ai3 = sum_ai2 * sum_ai;
898
899 tsum = tanh(born->obc_alpha*sum_ai-born->obc_beta*sum_ai2+born->obc_gamma*sum_ai3);
900 born->bRad[i] = rai_inv - tsum*rai_inv2;
901 born->bRad[i] = 1.0 / born->bRad[i];
902
903 fr->invsqrta[i] = gmx::invsqrt(born->bRad[i]);
904
905 tchain = rai * (born->obc_alpha-2*born->obc_beta*sum_ai+3*born->obc_gamma*sum_ai2);
906 born->drobc[i] = (1.0-tsum*tsum)*tchain*rai_inv2;
907 }
908 }
909
910 /* Extra (local) communication required for DD */
911 if (DOMAINDECOMP(cr))
912 {
913 dd_atom_spread_real(cr->dd, born->bRad);
914 dd_atom_spread_real(cr->dd, fr->invsqrta);
915 dd_atom_spread_real(cr->dd, born->drobc);
916 }
917
918 return 0;
919
920 }
921
922
923
924 int calc_gb_rad(t_commrec *cr, t_forcerec *fr, t_inputrec *ir, gmx_localtop_t *top,
925 rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md, t_nrnb *nrnb)
926 {
927 int cnt;
928 int ndadx;
929
930 if (fr->bAllvsAll && fr->dadx == nullptr)
931 {
932 /* We might need up to 8 atoms of padding before and after,
933 * and another 4 units to guarantee SSE alignment.
934 */
935 fr->nalloc_dadx = 2*(md->homenr+12)*(md->nr/2+1+12);
936 snew(fr->dadx_rawptr, fr->nalloc_dadx);
937 fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
938 }
939 else
940 {
941 /* In the SSE-enabled gb-loops, when writing to dadx, we
942 * always write 2*4 elements at a time, even in the case with only
943 * 1-3 j particles, where we only really need to write 2*(1-3)
944 * elements. This is because we want dadx to be aligned to a 16-
945 * byte boundary, and being able to use _mm_store/load_ps
946 */
947 ndadx = 2 * (nl->nrj + 3*nl->nri);
948
949 /* First, reallocate the dadx array, we need 3 extra for SSE */
950 if (ndadx + 3 > fr->nalloc_dadx)
951 {
952 fr->nalloc_dadx = over_alloc_large(ndadx) + 3;
953 srenew(fr->dadx_rawptr, fr->nalloc_dadx);
954 fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
955 }
956 }
957
958 if (fr->bAllvsAll)
959 {
960 cnt = md->homenr*(md->nr/2+1);
961
962 if (ir->gb_algorithm == egbSTILL)
963 {
964 genborn_allvsall_calc_still_radii(fr, md, born, top, x[0], &fr->AllvsAll_workgb);
965 /* 13 flops in outer loop, 47 flops in inner loop */
966 inc_nrnb(nrnb, eNR_BORN_AVA_RADII_STILL, md->homenr*13+cnt*47);
967 }
968 else if (ir->gb_algorithm == egbHCT || ir->gb_algorithm == egbOBC)
969 {
970 genborn_allvsall_calc_hct_obc_radii(fr, md, born, ir->gb_algorithm, top, x[0], &fr->AllvsAll_workgb);
971 /* 24 flops in outer loop, 183 in inner */
972 inc_nrnb(nrnb, eNR_BORN_AVA_RADII_HCT_OBC, md->homenr*24+cnt*183);
973 }
974 else
975 {
976 gmx_fatal(FARGS, "Bad gb algorithm for all-vs-all interactions");
977 }
978 return 0;
979 }
980
981 /* Switch for determining which algorithm to use for Born radii calculation */
982 #if GMX_DOUBLE
983
984 switch (ir->gb_algorithm)
985 {
986 case egbSTILL:
987 calc_gb_rad_still(cr, fr, top, x, nl, born, md);
988 break;
989 case egbHCT:
990 calc_gb_rad_hct(cr, fr, top, x, nl, born, md);
991 break;
992 case egbOBC:
993 calc_gb_rad_obc(cr, fr, top, x, nl, born, md);
994 break;
995
996 default:
997 gmx_fatal(FARGS, "Unknown double precision algorithm for Born radii calculation: %d", ir->gb_algorithm);
998 }
999
1000 #else
1001
1002 switch (ir->gb_algorithm)
1003 {
1004 case egbSTILL:
1005 calc_gb_rad_still(cr, fr, top, x, nl, born, md);
1006 break;
1007 case egbHCT:
1008 calc_gb_rad_hct(cr, fr, top, x, nl, born, md);
1009 break;
1010 case egbOBC:
1011 calc_gb_rad_obc(cr, fr, top, x, nl, born, md);
1012 break;
1013
1014 default:
1015 gmx_fatal(FARGS, "Unknown algorithm for Born radii calculation: %d", ir->gb_algorithm);
1016 }
1017
1018 #endif /* Double or single precision */
1019
1020 if (fr->bAllvsAll == FALSE)
1021 {
1022 switch (ir->gb_algorithm)
1023 {
1024 case egbSTILL:
1025 /* 17 flops per outer loop iteration, 47 flops per inner loop */
1026 inc_nrnb(nrnb, eNR_BORN_RADII_STILL, nl->nri*17+nl->nrj*47);
1027 break;
1028 case egbHCT:
1029 case egbOBC:
1030 /* 61 (assuming 10 for tanh) flops for outer loop iteration, 183 flops per inner loop */
1031 inc_nrnb(nrnb, eNR_BORN_RADII_HCT_OBC, nl->nri*61+nl->nrj*183);
1032 break;
1033
1034 default:
1035 break;
1036 }
1037 }
1038
1039 return 0;
1040 }
1041
1042
1043
1044 real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1045 real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
1046 real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph)
1047 {
1048 int i, j, n0, m, nnn, ai, aj;
1049 int ki;
1050
1051 real isai, isaj;
1052 real r, rsq11;
1053 real rinv11, iq;
1054 real isaprod, qq, gbscale, gbtabscale, Y, F, Geps, Heps2, Fp, VV, FF, rt, eps, eps2;
1055 real vgb, fgb, fijC, dvdatmp, fscal;
1056 real vctot;
1057
1058 rvec dx;
1059 ivec dt;
1060
1061 t_iatom *forceatoms;
1062
1063 /* Scale the electrostatics by gb_epsilon_solvent */
1064 facel = facel * ((1.0/epsilon_r) - 1.0/gb_epsilon_solvent);
1065
1066 gbtabscale = *p_gbtabscale;
1067 vctot = 0.0;
1068
1069 for (j = F_GB12; j <= F_GB14; j++)
1070 {
1071 forceatoms = idef->il[j].iatoms;
1072
1073 for (i = 0; i < idef->il[j].nr; )
1074 {
1075 /* To avoid reading in the interaction type, we just increment i to pass over
1076 * the types in the forceatoms array, this saves some memory accesses
1077 */
1078 i++;
1079 ai = forceatoms[i++];
1080 aj = forceatoms[i++];
1081
1082 ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx);
1083 rsq11 = iprod(dx, dx);
1084
1085 isai = invsqrta[ai];
1086 iq = (-1)*facel*charge[ai];
1087
1088 rinv11 = gmx::invsqrt(rsq11);
1089 isaj = invsqrta[aj];
1090 isaprod = isai*isaj;
1091 qq = isaprod*iq*charge[aj];
1092 gbscale = isaprod*gbtabscale;
1093 r = rsq11*rinv11;
1094 rt = r*gbscale;
1095 n0 = static_cast<int>(rt);
1096 eps = rt-n0;
1097 eps2 = eps*eps;
1098 nnn = 4*n0;
1099 Y = GBtab[nnn];
1100 F = GBtab[nnn+1];
1101 Geps = eps*GBtab[nnn+2];
1102 Heps2 = eps2*GBtab[nnn+3];
1103 Fp = F+Geps+Heps2;
1104 VV = Y+eps*Fp;
1105 FF = Fp+Geps+2.0*Heps2;
1106 vgb = qq*VV;
1107 fijC = qq*FF*gbscale;
1108 dvdatmp = -(vgb+fijC*r)*0.5;
1109 dvda[aj] = dvda[aj] + dvdatmp*isaj*isaj;
1110 dvda[ai] = dvda[ai] + dvdatmp*isai*isai;
1111 vctot = vctot + vgb;
1112 fgb = -(fijC)*rinv11;
1113
1114 if (graph)
1115 {
1116 ivec_sub(SHIFT_IVEC(graph, ai), SHIFT_IVEC(graph, aj), dt);
1117 ki = IVEC2IS(dt);
1118 }
1119
1120 for (m = 0; (m < DIM); m++) /* 15 */
1121 {
1122 fscal = fgb*dx[m];
1123 f[ai][m] += fscal;
1124 f[aj][m] -= fscal;
1125 fshift[ki][m] += fscal;
1126 fshift[CENTRAL][m] -= fscal;
1127 }
1128 }
1129 }
1130
1131 return vctot;
1132 }
1133
1134 static real calc_gb_selfcorrections(t_commrec *cr, int natoms,
1135 real *charge, gmx_genborn_t *born, real *dvda, double facel)
1136 {
1137 int i, ai, at0, at1;
1138 real rai, e, derb, q, q2, fi, rai_inv, vtot;
1139
1140 if (DOMAINDECOMP(cr))
1141 {
1142 at0 = 0;
1143 at1 = cr->dd->nat_home;
1144 }
1145 else
1146 {
1147 at0 = 0;
1148 at1 = natoms;
1149
1150 }
1151
1152 /* Scale the electrostatics by gb_epsilon_solvent */
1153 facel = facel * ((1.0/born->epsilon_r) - 1.0/born->gb_epsilon_solvent);
1154
1155 vtot = 0.0;
1156
1157 /* Apply self corrections */
1158 for (i = at0; i < at1; i++)
1159 {
1160 ai = i;
1161
1162 if (born->use[ai] == 1)
1163 {
1164 rai = born->bRad[ai];
1165 rai_inv = 1.0/rai;
1166 q = charge[ai];
1167 q2 = q*q;
1168 fi = facel*q2;
1169 e = fi*rai_inv;
1170 derb = 0.5*e*rai_inv*rai_inv;
1171 dvda[ai] += derb*rai;
1172 vtot -= 0.5*e;
1173 }
1174 }
1175
1176 return vtot;
1177
1178 }
1179
1180 static real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr, int natoms, gmx_genborn_t *born, gmx_localtop_t *top,
1181 real *dvda, t_mdatoms *md)
1182 {
1183 int ai, i, at0, at1;
1184 real e, es, rai, term, probe, tmp, factor;
1185 real rbi_inv, rbi_inv2;
1186
1187 if (DOMAINDECOMP(cr))
1188 {
1189 at0 = 0;
1190 at1 = cr->dd->nat_home;
1191 }
1192 else
1193 {
1194 at0 = 0;
1195 at1 = natoms;
1196 }
1197
1198 /* factor is the surface tension */
1199 factor = born->sa_surface_tension;
1200
1201 es = 0;
1202 probe = 0.14;
1203 term = M_PI*4;
1204
1205 for (i = at0; i < at1; i++)
1206 {
1207 ai = i;
1208
1209 if (born->use[ai] == 1)
1210 {
1211 rai = top->atomtypes.gb_radius[md->typeA[ai]];
1212 rbi_inv = fr->invsqrta[ai];
1213 rbi_inv2 = rbi_inv * rbi_inv;
1214 tmp = (rai*rbi_inv2)*(rai*rbi_inv2);
1215 tmp = tmp*tmp*tmp;
1216 e = factor*term*(rai+probe)*(rai+probe)*tmp;
1217 dvda[ai] = dvda[ai] - 6*e*rbi_inv2;
1218 es = es + e;
1219 }
1220 }
1221
1222 return es;
1223 }
1224
1225
1226
1227 static real calc_gb_chainrule(int natoms, t_nblist *nl, real *dadx, real *dvda, rvec x[], rvec t[], rvec fshift[],
1228 rvec shift_vec[], int gb_algorithm, gmx_genborn_t *born)
1229 {
1230 int i, k, n, ai, aj, nj0, nj1, n0, n1;
1231 int shift;
1232 real shX, shY, shZ;
1233 real fgb, rbi, fix1, fiy1, fiz1;
1234 real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11;
1235 real tx, ty, tz, rbai, rbaj, fgb_ai;
1236 real *rb;
1237
1238 n = 0;
1239 rb = born->work;
1240
1241 n0 = 0;
1242 n1 = natoms;
1243
1244 if (gb_algorithm == egbSTILL)
1245 {
1246 for (i = n0; i < n1; i++)
1247 {
1248 rbi = born->bRad[i];
1249 rb[i] = (2 * rbi * rbi * dvda[i])/ONE_4PI_EPS0;
1250 }
1251 }
1252 else if (gb_algorithm == egbHCT)
1253 {
1254 for (i = n0; i < n1; i++)
1255 {
1256 rbi = born->bRad[i];
1257 rb[i] = rbi * rbi * dvda[i];
1258 }
1259 }
1260 else if (gb_algorithm == egbOBC)
1261 {
1262 for (i = n0; i < n1; i++)
1263 {
1264 rbi = born->bRad[i];
1265 rb[i] = rbi * rbi * born->drobc[i] * dvda[i];
1266 }
1267 }
1268
1269 for (i = 0; i < nl->nri; i++)
1270 {
1271 ai = nl->iinr[i];
1272
1273 nj0 = nl->jindex[i];
1274 nj1 = nl->jindex[i+1];
1275
1276 /* Load shifts for this list */
1277 shift = nl->shift[i];
1278 shX = shift_vec[shift][0];
1279 shY = shift_vec[shift][1];
1280 shZ = shift_vec[shift][2];
1281
1282 /* Load atom i coordinates, add shift vectors */
1283 ix1 = shX + x[ai][0];
1284 iy1 = shY + x[ai][1];
1285 iz1 = shZ + x[ai][2];
1286
1287 fix1 = 0;
1288 fiy1 = 0;
1289 fiz1 = 0;
1290
1291 rbai = rb[ai];
1292
1293 for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
1294 {
1295 aj = nl->jjnr[k];
1296
1297 jx1 = x[aj][0];
1298 jy1 = x[aj][1];
1299 jz1 = x[aj][2];
1300
1301 dx11 = ix1 - jx1;
1302 dy11 = iy1 - jy1;
1303 dz11 = iz1 - jz1;
1304
1305 rbaj = rb[aj];
1306
1307 fgb = rbai*dadx[n++];
1308 fgb_ai = rbaj*dadx[n++];
1309
1310 /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1311 fgb = fgb + fgb_ai;
1312
1313 tx = fgb * dx11;
1314 ty = fgb * dy11;
1315 tz = fgb * dz11;
1316
1317 fix1 = fix1 + tx;
1318 fiy1 = fiy1 + ty;
1319 fiz1 = fiz1 + tz;
1320
1321 /* Update force on atom aj */
1322 t[aj][0] = t[aj][0] - tx;
1323 t[aj][1] = t[aj][1] - ty;
1324 t[aj][2] = t[aj][2] - tz;
1325 }
1326
1327 /* Update force and shift forces on atom ai */
1328 t[ai][0] = t[ai][0] + fix1;
1329 t[ai][1] = t[ai][1] + fiy1;
1330 t[ai][2] = t[ai][2] + fiz1;
1331
1332 fshift[shift][0] = fshift[shift][0] + fix1;
1333 fshift[shift][1] = fshift[shift][1] + fiy1;
1334 fshift[shift][2] = fshift[shift][2] + fiz1;
1335
1336 }
1337
1338 return 0;
1339 }
1340
1341
1342 void
1343 calc_gb_forces(t_commrec *cr, t_mdatoms *md, gmx_genborn_t *born, gmx_localtop_t *top,
1344 rvec x[], rvec f[], t_forcerec *fr, t_idef *idef, int gb_algorithm, int sa_algorithm, t_nrnb *nrnb,
1345 const t_pbc *pbc, const t_graph *graph, gmx_enerdata_t *enerd)
1346 {
1347 int cnt;
1348
1349 /* PBC or not? */
1350 const t_pbc *pbc_null;
1351
1352 if (fr->bMolPBC)
1353 {
1354 pbc_null = pbc;
1355 }
1356 else
1357 {
1358 pbc_null = nullptr;
1359 }
1360
1361 if (sa_algorithm == esaAPPROX)
1362 {
1363 /* Do a simple ACE type approximation for the non-polar solvation */
1364 enerd->term[F_NPSOLVATION] += calc_gb_nonpolar(cr, fr, born->nr, born, top, fr->dvda, md);
1365 }
1366
1367 /* Calculate the bonded GB-interactions using either table or analytical formula */
1368 enerd->term[F_GBPOL] += gb_bonds_tab(x, f, fr->fshift, md->chargeA, &(fr->gbtabscale),
1369 fr->invsqrta, fr->dvda, fr->gbtab->data, idef, born->epsilon_r, born->gb_epsilon_solvent, fr->ic->epsfac, pbc_null, graph);
1370
1371 /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1372 enerd->term[F_GBPOL] += calc_gb_selfcorrections(cr, born->nr, md->chargeA, born, fr->dvda, fr->ic->epsfac);
1373
1374 /* If parallel, sum the derivative of the potential w.r.t the born radii */
1375 if (DOMAINDECOMP(cr))
1376 {
1377 dd_atom_sum_real(cr->dd, fr->dvda);
1378 dd_atom_spread_real(cr->dd, fr->dvda);
1379 }
1380
1381 if (fr->bAllvsAll)
1382 {
1383 genborn_allvsall_calc_chainrule(fr, md, born, x[0], f[0], gb_algorithm, fr->AllvsAll_workgb);
1384 cnt = md->homenr*(md->nr/2+1);
1385 /* 9 flops for outer loop, 15 for inner */
1386 inc_nrnb(nrnb, eNR_BORN_AVA_CHAINRULE, md->homenr*9+cnt*15);
1387 return;
1388 }
1389
1390 calc_gb_chainrule(fr->natoms_force, fr->gblist, fr->dadx, fr->dvda,
1391 x, f, fr->fshift, fr->shift_vec, gb_algorithm, born);
1392
1393 if (!fr->bAllvsAll)
1394 {
1395 /* 9 flops for outer loop, 15 for inner */
1396 inc_nrnb(nrnb, eNR_BORN_CHAINRULE, fr->gblist->nri*9+fr->gblist->nrj*15);
1397 }
1398 }
1399
1400 static void add_j_to_gblist(gbtmpnbl_t *list, int aj)
1401 {
1402 if (list->naj >= list->aj_nalloc)
1403 {
1404 list->aj_nalloc = over_alloc_large(list->naj+1);
1405 srenew(list->aj, list->aj_nalloc);
1406 }
1407
1408 list->aj[list->naj++] = aj;
1409 }
1410
1411 static gbtmpnbl_t *find_gbtmplist(struct gbtmpnbls *lists, int shift)
1412 {
1413 int ind, i;
1414
1415 /* Search the list with the same shift, if there is one */
1416 ind = 0;
1417 while (ind < lists->nlist && shift != lists->list[ind].shift)
1418 {
1419 ind++;
1420 }
1421 if (ind == lists->nlist)
1422 {
1423 if (lists->nlist == lists->list_nalloc)
1424 {
1425 lists->list_nalloc++;
1426 srenew(lists->list, lists->list_nalloc);
1427 for (i = lists->nlist; i < lists->list_nalloc; i++)
1428 {
1429 lists->list[i].aj = nullptr;
1430 lists->list[i].aj_nalloc = 0;
1431 }
1432
1433 }
1434
1435 lists->list[lists->nlist].shift = shift;
1436 lists->list[lists->nlist].naj = 0;
1437 lists->nlist++;
1438 }
1439
1440 return &lists->list[ind];
1441 }
1442
1443 static void add_bondeds_to_gblist(t_ilist *il,
1444 gmx_bool bMolPBC, t_pbc *pbc, t_graph *g, rvec *x,
1445 struct gbtmpnbls *nls)
1446 {
1447 int ind, j, ai, aj, found;
1448 rvec dx;
1449 ivec dt;
1450 gbtmpnbl_t *list;
1451
1452 for (ind = 0; ind < il->nr; ind += 3)
1453 {
1454 ai = il->iatoms[ind+1];
1455 aj = il->iatoms[ind+2];
1456
1457 int shift = CENTRAL;
1458 if (g != nullptr)
1459 {
1460 rvec_sub(x[ai], x[aj], dx);
1461 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
1462 shift = IVEC2IS(dt);
1463 }
1464 else if (bMolPBC)
1465 {
1466 shift = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
1467 }
1468
1469 /* Find the list for this shift or create one */
1470 list = find_gbtmplist(&nls[ai], shift);
1471
1472 found = 0;
1473
1474 /* So that we do not add the same bond twice.
1475 * This happens with some constraints between 1-3 atoms
1476 * that are in the bond-list but should not be in the GB nb-list */
1477 for (j = 0; j < list->naj; j++)
1478 {
1479 if (list->aj[j] == aj)
1480 {
1481 found = 1;
1482 }
1483 }
1484
1485 if (found == 0)
1486 {
1487 if (ai == aj)
1488 {
1489 gmx_incons("ai == aj");
1490 }
1491
1492 add_j_to_gblist(list, aj);
1493 }
1494 }
1495 }
1496
1497
1498 int make_gb_nblist(t_commrec *cr, int gb_algorithm,
1499 rvec x[], matrix box,
1500 t_forcerec *fr, t_idef *idef, t_graph *graph, gmx_genborn_t *born)
1501 {
1502 int i, j, k, n, nj0, nj1, ai, shift, s;
1503 t_nblist *nblist;
1504 t_pbc pbc;
1505
1506 struct gbtmpnbls *nls;
1507 gbtmpnbl_t *list = nullptr;
1508
1509 set_pbc(&pbc, fr->ePBC, box);
1510 nls = born->nblist_work;
1511
1512 for (i = 0; i < born->nr; i++)
1513 {
1514 nls[i].nlist = 0;
1515 }
1516
1517 if (fr->bMolPBC)
1518 {
1519 set_pbc_dd(&pbc, fr->ePBC, cr->dd->nc, TRUE, box);
1520 }
1521
1522 switch (gb_algorithm)
1523 {
1524 case egbHCT:
1525 case egbOBC:
1526 /* Loop over 1-2, 1-3 and 1-4 interactions */
1527 for (j = F_GB12; j <= F_GB14; j++)
1528 {
1529 add_bondeds_to_gblist(&idef->il[j], fr->bMolPBC, &pbc, graph, x, nls);
1530 }
1531 break;
1532 case egbSTILL:
1533 /* Loop over 1-4 interactions */
1534 add_bondeds_to_gblist(&idef->il[F_GB14], fr->bMolPBC, &pbc, graph, x, nls);
1535 break;
1536 default:
1537 gmx_incons("Unknown GB algorithm");
1538 }
1539
1540 /* Loop over the VDWQQ and VDW nblists to set up the nonbonded part of the GB list */
1541 for (n = 0; (n < fr->nnblists); n++)
1542 {
1543 for (i = 0; (i < eNL_NR); i++)
1544 {
1545 nblist = &(fr->nblists[n].nlist_sr[i]);
1546
1547 if (nblist->nri > 0 && (i == eNL_VDWQQ || i == eNL_QQ))
1548 {
1549 for (j = 0; j < nblist->nri; j++)
1550 {
1551 ai = nblist->iinr[j];
1552 shift = nblist->shift[j];
1553
1554 /* Find the list for this shift or create one */
1555 list = find_gbtmplist(&nls[ai], shift);
1556
1557 nj0 = nblist->jindex[j];
1558 nj1 = nblist->jindex[j+1];
1559
1560 /* Add all the j-atoms in the non-bonded list to the GB list */
1561 for (k = nj0; k < nj1; k++)
1562 {
1563 add_j_to_gblist(list, nblist->jjnr[k]);
1564 }
1565 }
1566 }
1567 }
1568 }
1569
1570 /* Zero out some counters */
1571 fr->gblist->nri = 0;
1572 fr->gblist->nrj = 0;
1573
1574 fr->gblist->jindex[0] = fr->gblist->nri;
1575
1576 for (i = 0; i < fr->natoms_force; i++)
1577 {
1578 for (s = 0; s < nls[i].nlist; s++)
1579 {
1580 list = &nls[i].list[s];
1581
1582 /* Only add those atoms that actually have neighbours */
1583 if (born->use[i] != 0)
1584 {
1585 fr->gblist->iinr[fr->gblist->nri] = i;
1586 fr->gblist->shift[fr->gblist->nri] = list->shift;
1587 fr->gblist->nri++;
1588
1589 for (k = 0; k < list->naj; k++)
1590 {
1591 /* Memory allocation for jjnr */
1592 if (fr->gblist->nrj >= fr->gblist->maxnrj)
1593 {
1594 fr->gblist->maxnrj += over_alloc_large(fr->gblist->maxnrj);
1595
1596 if (debug)
1597 {
1598 fprintf(debug, "Increasing GB neighbourlist j size to %d\n", fr->gblist->maxnrj);
1599 }
1600
1601 srenew(fr->gblist->jjnr, fr->gblist->maxnrj);
1602 }
1603
1604 /* Put in list */
1605 if (i == list->aj[k])
1606 {
1607 gmx_incons("i == list->aj[k]");
1608 }
1609 fr->gblist->jjnr[fr->gblist->nrj++] = list->aj[k];
1610 }
1611
1612 fr->gblist->jindex[fr->gblist->nri] = fr->gblist->nrj;
1613 }
1614 }
1615 }
1616
1617 return 0;
1618 }
1619
1620 void make_local_gb(const t_commrec *cr, gmx_genborn_t *born, int gb_algorithm)
1621 {
1622 int i, at0, at1;
1623 gmx_domdec_t *dd = nullptr;
1624
1625 if (DOMAINDECOMP(cr))
1626 {
1627 dd = cr->dd;
1628 at0 = 0;
1629 at1 = dd->nat_tot;
1630 }
1631 else
1632 {
1633 /* Single node, just copy pointers and return */
1634 if (gb_algorithm == egbSTILL)
1635 {
1636 born->gpol = born->gpol_globalindex;
1637 born->vsolv = born->vsolv_globalindex;
1638 born->gb_radius = born->gb_radius_globalindex;
1639 }
1640 else
1641 {
1642 born->param = born->param_globalindex;
1643 born->gb_radius = born->gb_radius_globalindex;
1644 }
1645
1646 born->use = born->use_globalindex;
1647
1648 return;
1649 }
1650
1651 /* Reallocation of local arrays if necessary */
1652 /* fr->natoms_force is equal to dd->nat_tot */
1653 if (DOMAINDECOMP(cr) && dd->nat_tot > born->nalloc)
1654 {
1655 int nalloc;
1656
1657 nalloc = dd->nat_tot;
1658
1659 /* Arrays specific to different gb algorithms */
1660 if (gb_algorithm == egbSTILL)
1661 {
1662 srenew(born->gpol, nalloc+3);
1663 srenew(born->vsolv, nalloc+3);
1664 srenew(born->gb_radius, nalloc+3);
1665 for (i = born->nalloc; (i < nalloc+3); i++)
1666 {
1667 born->gpol[i] = 0;
1668 born->vsolv[i] = 0;
1669 born->gb_radius[i] = 0;
1670 }
1671 }
1672 else
1673 {
1674 srenew(born->param, nalloc+3);
1675 srenew(born->gb_radius, nalloc+3);
1676 for (i = born->nalloc; (i < nalloc+3); i++)
1677 {
1678 born->param[i] = 0;
1679 born->gb_radius[i] = 0;
1680 }
1681 }
1682
1683 /* All gb-algorithms use the array for vsites exclusions */
1684 srenew(born->use, nalloc+3);
1685 for (i = born->nalloc; (i < nalloc+3); i++)
1686 {
1687 born->use[i] = 0;
1688 }
1689
1690 born->nalloc = nalloc;
1691 }
1692
1693 /* With dd, copy algorithm specific arrays */
1694 if (gb_algorithm == egbSTILL)
1695 {
1696 for (i = at0; i < at1; i++)
1697 {
1698 born->gpol[i] = born->gpol_globalindex[dd->gatindex[i]];
1699 born->vsolv[i] = born->vsolv_globalindex[dd->gatindex[i]];
1700 born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
1701 born->use[i] = born->use_globalindex[dd->gatindex[i]];
1702 }
1703 }
1704 else
1705 {
1706 for (i = at0; i < at1; i++)
1707 {
1708 born->param[i] = born->param_globalindex[dd->gatindex[i]];
1709 born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
1710 born->use[i] = born->use_globalindex[dd->gatindex[i]];
1711 }
1712 }
1713 }
The ultimate molecular dynamics simulation package