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-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,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.
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.
46 #include "gromacs/domdec/domdec.h"
47 #include "gromacs/domdec/domdec_struct.h"
48 #include "gromacs/gmxlib/network.h"
49 #include "gromacs/gmxlib/nrnb.h"
50 #include "gromacs/math/functions.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdlib/gmx_omp_nthreads.h"
53 #include "gromacs/mdtypes/commrec.h"
54 #include "gromacs/pbcutil/ishift.h"
55 #include "gromacs/pbcutil/mshift.h"
56 #include "gromacs/pbcutil/pbc.h"
57 #include "gromacs/topology/ifunc.h"
58 #include "gromacs/topology/mtop_util.h"
59 #include "gromacs/utility/exceptions.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/gmxassert.h"
62 #include "gromacs/utility/gmxomp.h"
63 #include "gromacs/utility/smalloc.h"
66 /* The strategy used here for assigning virtual sites to (thread-)tasks
69 * We divide the atom range that vsites operate on (natoms_local with DD,
70 * 0 - last atom involved in vsites without DD) equally over all threads.
72 * Vsites in the local range constructed from atoms in the local range
73 * and/or other vsites that are fully local are assigned to a simple,
76 * Vsites that are not assigned after using the above criterion get assigned
77 * to a so called "interdependent" thread task when none of the constructing
78 * atoms is a vsite. These tasks are called interdependent, because one task
79 * accesses atoms assigned to a different task/thread.
80 * Note that this option is turned off with large (local) atom counts
81 * to avoid high memory usage.
83 * Any remaining vsites are assigned to a separate master thread task.
88 static void init_ilist(t_ilist
*ilist
)
90 for (int i
= 0; i
< F_NRE
; i
++)
94 ilist
[i
].iatoms
= nullptr;
98 /*! \brief List of atom indices belonging to a task */
100 //! List of atom indices
101 std::vector
<int> atom
;
104 /*! \brief Data structure for thread tasks that use constructing atoms outside their own atom range */
105 struct InterdependentTask
107 //! The interaction lists, only vsite entries are used
108 t_ilist ilist
[F_NRE
];
109 //! Thread/task-local force buffer
110 std::vector
<RVec
> force
;
111 //! The atom indices of the vsites of our task
112 std::vector
<int> vsite
;
113 //! Flags if elements in force are spread to or not
114 std::vector
<bool> use
;
115 //! The number of entries set to true in use
117 //! Array of atoms indices, size nthreads, covering all nuse set elements in use
118 std::vector
<AtomIndex
> atomIndex
;
119 //! List of tasks (force blocks) this task spread forces to
120 std::vector
<int> spreadTask
;
121 //! List of tasks that write to this tasks force block range
122 std::vector
<int> reduceTask
;
131 /*! \brief Vsite thread task data structure */
133 //! Start of atom range of this task
135 //! End of atom range of this task
137 //! The interaction lists, only vsite entries are used
138 t_ilist ilist
[F_NRE
];
139 //! Local fshift accumulation buffer
141 //! Local virial dx*df accumulation buffer
143 //! Tells if interdependent task idTask should be used (in addition to the rest of this task), this bool has the same value on all threads
144 bool useInterdependentTask
;
145 //! Data for vsites that involve constructing atoms in the atom range of other threads/tasks
146 InterdependentTask idTask
;
153 clear_rvecs(SHIFTS
, fshift
);
155 useInterdependentTask
= false;
160 /* The start and end values of for the vsite indices in the ftype enum.
161 * The validity of these values is checked in init_vsite.
162 * This is used to avoid loops over all ftypes just to get the vsite entries.
163 * (We should replace the fixed ilist array by only the used entries.)
165 static const int c_ftypeVsiteStart
= F_VSITE2
;
166 static const int c_ftypeVsiteEnd
= F_VSITEN
+ 1;
169 /* Returns the sum of the vsite ilist sizes over all vsite types */
170 static int vsiteIlistNrCount(const t_ilist
*ilist
)
173 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
175 nr
+= ilist
[ftype
].nr
;
181 static int pbc_rvec_sub(const t_pbc
*pbc
, const rvec xi
, const rvec xj
, rvec dx
)
185 return pbc_dx_aiuc(pbc
, xi
, xj
, dx
);
189 rvec_sub(xi
, xj
, dx
);
194 /* Vsite construction routines */
196 static void constr_vsite2(const rvec xi
, const rvec xj
, rvec x
, real a
, const t_pbc
*pbc
)
204 pbc_dx_aiuc(pbc
, xj
, xi
, dx
);
205 x
[XX
] = xi
[XX
] + a
*dx
[XX
];
206 x
[YY
] = xi
[YY
] + a
*dx
[YY
];
207 x
[ZZ
] = xi
[ZZ
] + a
*dx
[ZZ
];
211 x
[XX
] = b
*xi
[XX
] + a
*xj
[XX
];
212 x
[YY
] = b
*xi
[YY
] + a
*xj
[YY
];
213 x
[ZZ
] = b
*xi
[ZZ
] + a
*xj
[ZZ
];
217 /* TOTAL: 10 flops */
220 static void constr_vsite3(const rvec xi
, const rvec xj
, const rvec xk
, rvec x
, real a
, real b
,
230 pbc_dx_aiuc(pbc
, xj
, xi
, dxj
);
231 pbc_dx_aiuc(pbc
, xk
, xi
, dxk
);
232 x
[XX
] = xi
[XX
] + a
*dxj
[XX
] + b
*dxk
[XX
];
233 x
[YY
] = xi
[YY
] + a
*dxj
[YY
] + b
*dxk
[YY
];
234 x
[ZZ
] = xi
[ZZ
] + a
*dxj
[ZZ
] + b
*dxk
[ZZ
];
238 x
[XX
] = c
*xi
[XX
] + a
*xj
[XX
] + b
*xk
[XX
];
239 x
[YY
] = c
*xi
[YY
] + a
*xj
[YY
] + b
*xk
[YY
];
240 x
[ZZ
] = c
*xi
[ZZ
] + a
*xj
[ZZ
] + b
*xk
[ZZ
];
244 /* TOTAL: 17 flops */
247 static void constr_vsite3FD(const rvec xi
, const rvec xj
, const rvec xk
, rvec x
, real a
, real b
,
253 pbc_rvec_sub(pbc
, xj
, xi
, xij
);
254 pbc_rvec_sub(pbc
, xk
, xj
, xjk
);
257 /* temp goes from i to a point on the line jk */
258 temp
[XX
] = xij
[XX
] + a
*xjk
[XX
];
259 temp
[YY
] = xij
[YY
] + a
*xjk
[YY
];
260 temp
[ZZ
] = xij
[ZZ
] + a
*xjk
[ZZ
];
263 c
= b
*gmx::invsqrt(iprod(temp
, temp
));
266 x
[XX
] = xi
[XX
] + c
*temp
[XX
];
267 x
[YY
] = xi
[YY
] + c
*temp
[YY
];
268 x
[ZZ
] = xi
[ZZ
] + c
*temp
[ZZ
];
271 /* TOTAL: 34 flops */
274 static void constr_vsite3FAD(const rvec xi
, const rvec xj
, const rvec xk
, rvec x
, real a
, real b
, const t_pbc
*pbc
)
277 real a1
, b1
, c1
, invdij
;
279 pbc_rvec_sub(pbc
, xj
, xi
, xij
);
280 pbc_rvec_sub(pbc
, xk
, xj
, xjk
);
283 invdij
= gmx::invsqrt(iprod(xij
, xij
));
284 c1
= invdij
* invdij
* iprod(xij
, xjk
);
285 xp
[XX
] = xjk
[XX
] - c1
*xij
[XX
];
286 xp
[YY
] = xjk
[YY
] - c1
*xij
[YY
];
287 xp
[ZZ
] = xjk
[ZZ
] - c1
*xij
[ZZ
];
289 b1
= b
*gmx::invsqrt(iprod(xp
, xp
));
292 x
[XX
] = xi
[XX
] + a1
*xij
[XX
] + b1
*xp
[XX
];
293 x
[YY
] = xi
[YY
] + a1
*xij
[YY
] + b1
*xp
[YY
];
294 x
[ZZ
] = xi
[ZZ
] + a1
*xij
[ZZ
] + b1
*xp
[ZZ
];
297 /* TOTAL: 63 flops */
300 static void constr_vsite3OUT(const rvec xi
, const rvec xj
, const rvec xk
, rvec x
,
301 real a
, real b
, real c
, const t_pbc
*pbc
)
305 pbc_rvec_sub(pbc
, xj
, xi
, xij
);
306 pbc_rvec_sub(pbc
, xk
, xi
, xik
);
307 cprod(xij
, xik
, temp
);
310 x
[XX
] = xi
[XX
] + a
*xij
[XX
] + b
*xik
[XX
] + c
*temp
[XX
];
311 x
[YY
] = xi
[YY
] + a
*xij
[YY
] + b
*xik
[YY
] + c
*temp
[YY
];
312 x
[ZZ
] = xi
[ZZ
] + a
*xij
[ZZ
] + b
*xik
[ZZ
] + c
*temp
[ZZ
];
315 /* TOTAL: 33 flops */
318 static void constr_vsite4FD(const rvec xi
, const rvec xj
, const rvec xk
, const rvec xl
, rvec x
,
319 real a
, real b
, real c
, const t_pbc
*pbc
)
321 rvec xij
, xjk
, xjl
, temp
;
324 pbc_rvec_sub(pbc
, xj
, xi
, xij
);
325 pbc_rvec_sub(pbc
, xk
, xj
, xjk
);
326 pbc_rvec_sub(pbc
, xl
, xj
, xjl
);
329 /* temp goes from i to a point on the plane jkl */
330 temp
[XX
] = xij
[XX
] + a
*xjk
[XX
] + b
*xjl
[XX
];
331 temp
[YY
] = xij
[YY
] + a
*xjk
[YY
] + b
*xjl
[YY
];
332 temp
[ZZ
] = xij
[ZZ
] + a
*xjk
[ZZ
] + b
*xjl
[ZZ
];
335 d
= c
*gmx::invsqrt(iprod(temp
, temp
));
338 x
[XX
] = xi
[XX
] + d
*temp
[XX
];
339 x
[YY
] = xi
[YY
] + d
*temp
[YY
];
340 x
[ZZ
] = xi
[ZZ
] + d
*temp
[ZZ
];
343 /* TOTAL: 43 flops */
346 static void constr_vsite4FDN(const rvec xi
, const rvec xj
, const rvec xk
, const rvec xl
, rvec x
,
347 real a
, real b
, real c
, const t_pbc
*pbc
)
349 rvec xij
, xik
, xil
, ra
, rb
, rja
, rjb
, rm
;
352 pbc_rvec_sub(pbc
, xj
, xi
, xij
);
353 pbc_rvec_sub(pbc
, xk
, xi
, xik
);
354 pbc_rvec_sub(pbc
, xl
, xi
, xil
);
367 rvec_sub(ra
, xij
, rja
);
368 rvec_sub(rb
, xij
, rjb
);
374 d
= c
*gmx::invsqrt(norm2(rm
));
377 x
[XX
] = xi
[XX
] + d
*rm
[XX
];
378 x
[YY
] = xi
[YY
] + d
*rm
[YY
];
379 x
[ZZ
] = xi
[ZZ
] + d
*rm
[ZZ
];
382 /* TOTAL: 47 flops */
386 static int constr_vsiten(const t_iatom
*ia
, const t_iparams ip
[],
387 rvec
*x
, const t_pbc
*pbc
)
394 n3
= 3*ip
[ia
[0]].vsiten
.n
;
397 copy_rvec(x
[ai
], x1
);
399 for (int i
= 3; i
< n3
; i
+= 3)
402 a
= ip
[ia
[i
]].vsiten
.a
;
405 pbc_dx_aiuc(pbc
, x
[ai
], x1
, dx
);
409 rvec_sub(x
[ai
], x1
, dx
);
411 dsum
[XX
] += a
*dx
[XX
];
412 dsum
[YY
] += a
*dx
[YY
];
413 dsum
[ZZ
] += a
*dx
[ZZ
];
417 x
[av
][XX
] = x1
[XX
] + dsum
[XX
];
418 x
[av
][YY
] = x1
[YY
] + dsum
[YY
];
419 x
[av
][ZZ
] = x1
[ZZ
] + dsum
[ZZ
];
425 static void construct_vsites_thread(const gmx_vsite_t
*vsite
,
428 const t_iparams ip
[], const t_ilist ilist
[],
429 const t_pbc
*pbc_null
)
433 const t_pbc
*pbc_null2
;
445 bPBCAll
= (pbc_null
!= nullptr && !vsite
->bHaveChargeGroups
);
449 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
451 if (ilist
[ftype
].nr
== 0)
457 int nra
= interaction_function
[ftype
].nratoms
;
459 int nr
= ilist
[ftype
].nr
;
461 const t_iatom
*ia
= ilist
[ftype
].iatoms
;
465 pbc_null2
= pbc_null
;
467 else if (pbc_null
!= nullptr)
469 vsite_pbc
= vsite
->vsite_pbc_loc
[ftype
- c_ftypeVsiteStart
];
472 for (int i
= 0; i
< nr
; )
475 /* The vsite and constructing atoms */
478 /* Constants for constructing vsites */
479 real a1
= ip
[tp
].vsite
.a
;
480 /* Check what kind of pbc we need to use */
485 /* No charge groups, vsite follows its own pbc */
487 copy_rvec(x
[avsite
], xpbc
);
489 else if (vsite_pbc
!= nullptr)
491 pbc_atom
= vsite_pbc
[i
/(1 + nra
)];
496 /* We need to copy the coordinates here,
497 * single for single atom cg's pbc_atom
498 * is the vsite itself.
500 copy_rvec(x
[pbc_atom
], xpbc
);
502 pbc_null2
= pbc_null
;
513 /* Copy the old position */
515 copy_rvec(x
[avsite
], xv
);
517 /* Construct the vsite depending on type */
524 constr_vsite2(x
[ai
], x
[aj
], x
[avsite
], a1
, pbc_null2
);
530 constr_vsite3(x
[ai
], x
[aj
], x
[ak
], x
[avsite
], a1
, b1
, pbc_null2
);
536 constr_vsite3FD(x
[ai
], x
[aj
], x
[ak
], x
[avsite
], a1
, b1
, pbc_null2
);
542 constr_vsite3FAD(x
[ai
], x
[aj
], x
[ak
], x
[avsite
], a1
, b1
, pbc_null2
);
549 constr_vsite3OUT(x
[ai
], x
[aj
], x
[ak
], x
[avsite
], a1
, b1
, c1
, pbc_null2
);
557 constr_vsite4FD(x
[ai
], x
[aj
], x
[ak
], x
[al
], x
[avsite
], a1
, b1
, c1
,
566 constr_vsite4FDN(x
[ai
], x
[aj
], x
[ak
], x
[al
], x
[avsite
], a1
, b1
, c1
,
570 inc
= constr_vsiten(ia
, ip
, x
, pbc_null2
);
573 gmx_fatal(FARGS
, "No such vsite type %d in %s, line %d",
574 ftype
, __FILE__
, __LINE__
);
579 /* Match the pbc of this vsite to the rest of its charge group */
581 int ishift
= pbc_dx_aiuc(pbc_null
, x
[avsite
], xpbc
, dx
);
582 if (ishift
!= CENTRAL
)
584 rvec_add(xpbc
, dx
, x
[avsite
]);
589 /* Calculate velocity of vsite... */
591 rvec_sub(x
[avsite
], xv
, vv
);
592 svmul(inv_dt
, vv
, v
[avsite
]);
595 /* Increment loop variables */
603 void construct_vsites(const gmx_vsite_t
*vsite
,
606 const t_iparams ip
[], const t_ilist ilist
[],
607 int ePBC
, gmx_bool bMolPBC
,
608 t_commrec
*cr
, matrix box
)
610 t_pbc pbc
, *pbc_null
;
612 gmx_bool bDomDec
= cr
&& DOMAINDECOMP(cr
);
614 /* We only need to do pbc when we have inter-cg vsites */
615 if (ePBC
!= epbcNONE
&& (bDomDec
|| bMolPBC
) && vsite
->n_intercg_vsite
)
617 /* This is wasting some CPU time as we now do this multiple times
621 clear_ivec(null_ivec
);
622 pbc_null
= set_pbc_dd(&pbc
, ePBC
,
623 DOMAINDECOMP(cr
) ? cr
->dd
->nc
: null_ivec
,
633 dd_move_x_vsites(cr
->dd
, box
, x
);
636 if (vsite
->nthreads
== 1)
638 construct_vsites_thread(vsite
,
645 #pragma omp parallel num_threads(vsite->nthreads)
649 int th
= gmx_omp_get_thread_num();
650 construct_vsites_thread(vsite
,
652 ip
, vsite
->tData
[th
]->ilist
,
654 if (vsite
->tData
[th
]->useInterdependentTask
)
656 /* Here we don't need a barrier (unlike the spreading),
657 * since both tasks only construct vsites from particles,
658 * or local vsites, not from non-local vsites.
660 construct_vsites_thread(vsite
,
662 ip
, vsite
->tData
[th
]->idTask
.ilist
,
666 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
668 /* Now we can construct the vsites that might depend on other vsites */
669 construct_vsites_thread(vsite
,
671 ip
, vsite
->tData
[vsite
->nthreads
]->ilist
,
676 static void spread_vsite2(const t_iatom ia
[], real a
,
678 rvec f
[], rvec fshift
[],
679 const t_pbc
*pbc
, const t_graph
*g
)
690 svmul(1 - a
, f
[av
], fi
);
691 svmul( a
, f
[av
], fj
);
700 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, av
), di
);
702 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), di
);
707 siv
= pbc_dx_aiuc(pbc
, x
[ai
], x
[av
], dx
);
708 sij
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], dx
);
716 if (fshift
&& (siv
!= CENTRAL
|| sij
!= CENTRAL
))
718 rvec_inc(fshift
[siv
], f
[av
]);
719 rvec_dec(fshift
[CENTRAL
], fi
);
720 rvec_dec(fshift
[sij
], fj
);
723 /* TOTAL: 13 flops */
726 void construct_vsites_mtop(gmx_vsite_t
*vsite
,
727 gmx_mtop_t
*mtop
, rvec x
[])
730 for (int mb
= 0; mb
< mtop
->nmolblock
; mb
++)
732 const gmx_molblock_t
*molb
= &mtop
->molblock
[mb
];
733 const gmx_moltype_t
*molt
= &mtop
->moltype
[molb
->type
];
734 for (int mol
= 0; mol
< molb
->nmol
; mol
++)
736 construct_vsites(vsite
, x
+as
, 0.0, nullptr,
737 mtop
->ffparams
.iparams
, molt
->ilist
,
738 epbcNONE
, TRUE
, nullptr, nullptr);
739 as
+= molt
->atoms
.nr
;
744 static void spread_vsite3(const t_iatom ia
[], real a
, real b
,
746 rvec f
[], rvec fshift
[],
747 const t_pbc
*pbc
, const t_graph
*g
)
759 svmul(1 - a
- b
, f
[av
], fi
);
760 svmul( a
, f
[av
], fj
);
761 svmul( b
, f
[av
], fk
);
771 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, ia
[1]), di
);
773 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), di
);
775 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, ak
), di
);
780 siv
= pbc_dx_aiuc(pbc
, x
[ai
], x
[av
], dx
);
781 sij
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], dx
);
782 sik
= pbc_dx_aiuc(pbc
, x
[ai
], x
[ak
], dx
);
791 if (fshift
&& (siv
!= CENTRAL
|| sij
!= CENTRAL
|| sik
!= CENTRAL
))
793 rvec_inc(fshift
[siv
], f
[av
]);
794 rvec_dec(fshift
[CENTRAL
], fi
);
795 rvec_dec(fshift
[sij
], fj
);
796 rvec_dec(fshift
[sik
], fk
);
799 /* TOTAL: 20 flops */
802 static void spread_vsite3FD(const t_iatom ia
[], real a
, real b
,
804 rvec f
[], rvec fshift
[],
805 gmx_bool VirCorr
, matrix dxdf
,
806 const t_pbc
*pbc
, const t_graph
*g
)
808 real c
, invl
, fproj
, a1
;
809 rvec xvi
, xij
, xjk
, xix
, fv
, temp
;
810 t_iatom av
, ai
, aj
, ak
;
818 copy_rvec(f
[av
], fv
);
820 sji
= pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], xij
);
821 skj
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], xjk
);
824 /* xix goes from i to point x on the line jk */
825 xix
[XX
] = xij
[XX
]+a
*xjk
[XX
];
826 xix
[YY
] = xij
[YY
]+a
*xjk
[YY
];
827 xix
[ZZ
] = xij
[ZZ
]+a
*xjk
[ZZ
];
830 invl
= gmx::invsqrt(iprod(xix
, xix
));
834 fproj
= iprod(xix
, fv
)*invl
*invl
; /* = (xix . f)/(xix . xix) */
836 temp
[XX
] = c
*(fv
[XX
]-fproj
*xix
[XX
]);
837 temp
[YY
] = c
*(fv
[YY
]-fproj
*xix
[YY
]);
838 temp
[ZZ
] = c
*(fv
[ZZ
]-fproj
*xix
[ZZ
]);
841 /* c is already calculated in constr_vsite3FD
842 storing c somewhere will save 26 flops! */
845 f
[ai
][XX
] += fv
[XX
] - temp
[XX
];
846 f
[ai
][YY
] += fv
[YY
] - temp
[YY
];
847 f
[ai
][ZZ
] += fv
[ZZ
] - temp
[ZZ
];
848 f
[aj
][XX
] += a1
*temp
[XX
];
849 f
[aj
][YY
] += a1
*temp
[YY
];
850 f
[aj
][ZZ
] += a1
*temp
[ZZ
];
851 f
[ak
][XX
] += a
*temp
[XX
];
852 f
[ak
][YY
] += a
*temp
[YY
];
853 f
[ak
][ZZ
] += a
*temp
[ZZ
];
858 ivec_sub(SHIFT_IVEC(g
, ia
[1]), SHIFT_IVEC(g
, ai
), di
);
860 ivec_sub(SHIFT_IVEC(g
, aj
), SHIFT_IVEC(g
, ai
), di
);
862 ivec_sub(SHIFT_IVEC(g
, ak
), SHIFT_IVEC(g
, aj
), di
);
867 svi
= pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xvi
);
874 if (fshift
&& (svi
!= CENTRAL
|| sji
!= CENTRAL
|| skj
!= CENTRAL
))
876 rvec_dec(fshift
[svi
], fv
);
877 fshift
[CENTRAL
][XX
] += fv
[XX
] - (1 + a
)*temp
[XX
];
878 fshift
[CENTRAL
][YY
] += fv
[YY
] - (1 + a
)*temp
[YY
];
879 fshift
[CENTRAL
][ZZ
] += fv
[ZZ
] - (1 + a
)*temp
[ZZ
];
880 fshift
[ sji
][XX
] += temp
[XX
];
881 fshift
[ sji
][YY
] += temp
[YY
];
882 fshift
[ sji
][ZZ
] += temp
[ZZ
];
883 fshift
[ skj
][XX
] += a
*temp
[XX
];
884 fshift
[ skj
][YY
] += a
*temp
[YY
];
885 fshift
[ skj
][ZZ
] += a
*temp
[ZZ
];
890 /* When VirCorr=TRUE, the virial for the current forces is not
891 * calculated from the redistributed forces. This means that
892 * the effect of non-linear virtual site constructions on the virial
893 * needs to be added separately. This contribution can be calculated
894 * in many ways, but the simplest and cheapest way is to use
895 * the first constructing atom ai as a reference position in space:
896 * subtract (xv-xi)*fv and add (xj-xi)*fj + (xk-xi)*fk.
900 pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xiv
);
902 for (int i
= 0; i
< DIM
; i
++)
904 for (int j
= 0; j
< DIM
; j
++)
906 /* As xix is a linear combination of j and k, use that here */
907 dxdf
[i
][j
] += -xiv
[i
]*fv
[j
] + xix
[i
]*temp
[j
];
912 /* TOTAL: 61 flops */
915 static void spread_vsite3FAD(const t_iatom ia
[], real a
, real b
,
917 rvec f
[], rvec fshift
[],
918 gmx_bool VirCorr
, matrix dxdf
,
919 const t_pbc
*pbc
, const t_graph
*g
)
921 rvec xvi
, xij
, xjk
, xperp
, Fpij
, Fppp
, fv
, f1
, f2
, f3
;
922 real a1
, b1
, c1
, c2
, invdij
, invdij2
, invdp
, fproj
;
923 t_iatom av
, ai
, aj
, ak
;
924 int svi
, sji
, skj
, d
;
931 copy_rvec(f
[ia
[1]], fv
);
933 sji
= pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], xij
);
934 skj
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], xjk
);
937 invdij
= gmx::invsqrt(iprod(xij
, xij
));
938 invdij2
= invdij
* invdij
;
939 c1
= iprod(xij
, xjk
) * invdij2
;
940 xperp
[XX
] = xjk
[XX
] - c1
*xij
[XX
];
941 xperp
[YY
] = xjk
[YY
] - c1
*xij
[YY
];
942 xperp
[ZZ
] = xjk
[ZZ
] - c1
*xij
[ZZ
];
943 /* xperp in plane ijk, perp. to ij */
944 invdp
= gmx::invsqrt(iprod(xperp
, xperp
));
949 /* a1, b1 and c1 are already calculated in constr_vsite3FAD
950 storing them somewhere will save 45 flops! */
952 fproj
= iprod(xij
, fv
)*invdij2
;
953 svmul(fproj
, xij
, Fpij
); /* proj. f on xij */
954 svmul(iprod(xperp
, fv
)*invdp
*invdp
, xperp
, Fppp
); /* proj. f on xperp */
955 svmul(b1
*fproj
, xperp
, f3
);
958 rvec_sub(fv
, Fpij
, f1
); /* f1 = f - Fpij */
959 rvec_sub(f1
, Fppp
, f2
); /* f2 = f - Fpij - Fppp */
960 for (d
= 0; (d
< DIM
); d
++)
968 f
[ai
][XX
] += fv
[XX
] - f1
[XX
] + c1
*f2
[XX
] + f3
[XX
];
969 f
[ai
][YY
] += fv
[YY
] - f1
[YY
] + c1
*f2
[YY
] + f3
[YY
];
970 f
[ai
][ZZ
] += fv
[ZZ
] - f1
[ZZ
] + c1
*f2
[ZZ
] + f3
[ZZ
];
971 f
[aj
][XX
] += f1
[XX
] - c2
*f2
[XX
] - f3
[XX
];
972 f
[aj
][YY
] += f1
[YY
] - c2
*f2
[YY
] - f3
[YY
];
973 f
[aj
][ZZ
] += f1
[ZZ
] - c2
*f2
[ZZ
] - f3
[ZZ
];
981 ivec_sub(SHIFT_IVEC(g
, ia
[1]), SHIFT_IVEC(g
, ai
), di
);
983 ivec_sub(SHIFT_IVEC(g
, aj
), SHIFT_IVEC(g
, ai
), di
);
985 ivec_sub(SHIFT_IVEC(g
, ak
), SHIFT_IVEC(g
, aj
), di
);
990 svi
= pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xvi
);
997 if (fshift
&& (svi
!= CENTRAL
|| sji
!= CENTRAL
|| skj
!= CENTRAL
))
999 rvec_dec(fshift
[svi
], fv
);
1000 fshift
[CENTRAL
][XX
] += fv
[XX
] - f1
[XX
] - (1-c1
)*f2
[XX
] + f3
[XX
];
1001 fshift
[CENTRAL
][YY
] += fv
[YY
] - f1
[YY
] - (1-c1
)*f2
[YY
] + f3
[YY
];
1002 fshift
[CENTRAL
][ZZ
] += fv
[ZZ
] - f1
[ZZ
] - (1-c1
)*f2
[ZZ
] + f3
[ZZ
];
1003 fshift
[ sji
][XX
] += f1
[XX
] - c1
*f2
[XX
] - f3
[XX
];
1004 fshift
[ sji
][YY
] += f1
[YY
] - c1
*f2
[YY
] - f3
[YY
];
1005 fshift
[ sji
][ZZ
] += f1
[ZZ
] - c1
*f2
[ZZ
] - f3
[ZZ
];
1006 fshift
[ skj
][XX
] += f2
[XX
];
1007 fshift
[ skj
][YY
] += f2
[YY
];
1008 fshift
[ skj
][ZZ
] += f2
[ZZ
];
1016 pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xiv
);
1018 for (i
= 0; i
< DIM
; i
++)
1020 for (j
= 0; j
< DIM
; j
++)
1022 /* Note that xik=xij+xjk, so we have to add xij*f2 */
1025 + xij
[i
]*(f1
[j
] + (1 - c2
)*f2
[j
] - f3
[j
])
1031 /* TOTAL: 113 flops */
1034 static void spread_vsite3OUT(const t_iatom ia
[], real a
, real b
, real c
,
1036 rvec f
[], rvec fshift
[],
1037 gmx_bool VirCorr
, matrix dxdf
,
1038 const t_pbc
*pbc
, const t_graph
*g
)
1040 rvec xvi
, xij
, xik
, fv
, fj
, fk
;
1051 sji
= pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], xij
);
1052 ski
= pbc_rvec_sub(pbc
, x
[ak
], x
[ai
], xik
);
1055 copy_rvec(f
[av
], fv
);
1062 fj
[XX
] = a
*fv
[XX
] - xik
[ZZ
]*cfy
+ xik
[YY
]*cfz
;
1063 fj
[YY
] = xik
[ZZ
]*cfx
+ a
*fv
[YY
] - xik
[XX
]*cfz
;
1064 fj
[ZZ
] = -xik
[YY
]*cfx
+ xik
[XX
]*cfy
+ a
*fv
[ZZ
];
1066 fk
[XX
] = b
*fv
[XX
] + xij
[ZZ
]*cfy
- xij
[YY
]*cfz
;
1067 fk
[YY
] = -xij
[ZZ
]*cfx
+ b
*fv
[YY
] + xij
[XX
]*cfz
;
1068 fk
[ZZ
] = xij
[YY
]*cfx
- xij
[XX
]*cfy
+ b
*fv
[ZZ
];
1071 f
[ai
][XX
] += fv
[XX
] - fj
[XX
] - fk
[XX
];
1072 f
[ai
][YY
] += fv
[YY
] - fj
[YY
] - fk
[YY
];
1073 f
[ai
][ZZ
] += fv
[ZZ
] - fj
[ZZ
] - fk
[ZZ
];
1074 rvec_inc(f
[aj
], fj
);
1075 rvec_inc(f
[ak
], fk
);
1080 ivec_sub(SHIFT_IVEC(g
, ia
[1]), SHIFT_IVEC(g
, ai
), di
);
1082 ivec_sub(SHIFT_IVEC(g
, aj
), SHIFT_IVEC(g
, ai
), di
);
1084 ivec_sub(SHIFT_IVEC(g
, ak
), SHIFT_IVEC(g
, ai
), di
);
1089 svi
= pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xvi
);
1096 if (fshift
&& (svi
!= CENTRAL
|| sji
!= CENTRAL
|| ski
!= CENTRAL
))
1098 rvec_dec(fshift
[svi
], fv
);
1099 fshift
[CENTRAL
][XX
] += fv
[XX
] - fj
[XX
] - fk
[XX
];
1100 fshift
[CENTRAL
][YY
] += fv
[YY
] - fj
[YY
] - fk
[YY
];
1101 fshift
[CENTRAL
][ZZ
] += fv
[ZZ
] - fj
[ZZ
] - fk
[ZZ
];
1102 rvec_inc(fshift
[sji
], fj
);
1103 rvec_inc(fshift
[ski
], fk
);
1110 pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xiv
);
1112 for (int i
= 0; i
< DIM
; i
++)
1114 for (int j
= 0; j
< DIM
; j
++)
1116 dxdf
[i
][j
] += -xiv
[i
]*fv
[j
] + xij
[i
]*fj
[j
] + xik
[i
]*fk
[j
];
1121 /* TOTAL: 54 flops */
1124 static void spread_vsite4FD(const t_iatom ia
[], real a
, real b
, real c
,
1126 rvec f
[], rvec fshift
[],
1127 gmx_bool VirCorr
, matrix dxdf
,
1128 const t_pbc
*pbc
, const t_graph
*g
)
1130 real d
, invl
, fproj
, a1
;
1131 rvec xvi
, xij
, xjk
, xjl
, xix
, fv
, temp
;
1132 int av
, ai
, aj
, ak
, al
;
1134 int svi
, sji
, skj
, slj
, m
;
1142 sji
= pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], xij
);
1143 skj
= pbc_rvec_sub(pbc
, x
[ak
], x
[aj
], xjk
);
1144 slj
= pbc_rvec_sub(pbc
, x
[al
], x
[aj
], xjl
);
1147 /* xix goes from i to point x on the plane jkl */
1148 for (m
= 0; m
< DIM
; m
++)
1150 xix
[m
] = xij
[m
] + a
*xjk
[m
] + b
*xjl
[m
];
1154 invl
= gmx::invsqrt(iprod(xix
, xix
));
1156 /* 4 + ?10? flops */
1158 copy_rvec(f
[av
], fv
);
1160 fproj
= iprod(xix
, fv
)*invl
*invl
; /* = (xix . f)/(xix . xix) */
1162 for (m
= 0; m
< DIM
; m
++)
1164 temp
[m
] = d
*(fv
[m
] - fproj
*xix
[m
]);
1168 /* c is already calculated in constr_vsite3FD
1169 storing c somewhere will save 35 flops! */
1172 for (m
= 0; m
< DIM
; m
++)
1174 f
[ai
][m
] += fv
[m
] - temp
[m
];
1175 f
[aj
][m
] += a1
*temp
[m
];
1176 f
[ak
][m
] += a
*temp
[m
];
1177 f
[al
][m
] += b
*temp
[m
];
1183 ivec_sub(SHIFT_IVEC(g
, ia
[1]), SHIFT_IVEC(g
, ai
), di
);
1185 ivec_sub(SHIFT_IVEC(g
, aj
), SHIFT_IVEC(g
, ai
), di
);
1187 ivec_sub(SHIFT_IVEC(g
, ak
), SHIFT_IVEC(g
, aj
), di
);
1189 ivec_sub(SHIFT_IVEC(g
, al
), SHIFT_IVEC(g
, aj
), di
);
1194 svi
= pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xvi
);
1202 (svi
!= CENTRAL
|| sji
!= CENTRAL
|| skj
!= CENTRAL
|| slj
!= CENTRAL
))
1204 rvec_dec(fshift
[svi
], fv
);
1205 for (m
= 0; m
< DIM
; m
++)
1207 fshift
[CENTRAL
][m
] += fv
[m
] - (1 + a
+ b
)*temp
[m
];
1208 fshift
[ sji
][m
] += temp
[m
];
1209 fshift
[ skj
][m
] += a
*temp
[m
];
1210 fshift
[ slj
][m
] += b
*temp
[m
];
1219 pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xiv
);
1221 for (i
= 0; i
< DIM
; i
++)
1223 for (j
= 0; j
< DIM
; j
++)
1225 dxdf
[i
][j
] += -xiv
[i
]*fv
[j
] + xix
[i
]*temp
[j
];
1230 /* TOTAL: 77 flops */
1234 static void spread_vsite4FDN(const t_iatom ia
[], real a
, real b
, real c
,
1236 rvec f
[], rvec fshift
[],
1237 gmx_bool VirCorr
, matrix dxdf
,
1238 const t_pbc
*pbc
, const t_graph
*g
)
1240 rvec xvi
, xij
, xik
, xil
, ra
, rb
, rja
, rjb
, rab
, rm
, rt
;
1241 rvec fv
, fj
, fk
, fl
;
1245 int av
, ai
, aj
, ak
, al
;
1246 int svi
, sij
, sik
, sil
;
1248 /* DEBUG: check atom indices */
1255 copy_rvec(f
[av
], fv
);
1257 sij
= pbc_rvec_sub(pbc
, x
[aj
], x
[ai
], xij
);
1258 sik
= pbc_rvec_sub(pbc
, x
[ak
], x
[ai
], xik
);
1259 sil
= pbc_rvec_sub(pbc
, x
[al
], x
[ai
], xil
);
1272 rvec_sub(ra
, xij
, rja
);
1273 rvec_sub(rb
, xij
, rjb
);
1274 rvec_sub(rb
, ra
, rab
);
1277 cprod(rja
, rjb
, rm
);
1280 invrm
= gmx::invsqrt(norm2(rm
));
1281 denom
= invrm
*invrm
;
1284 cfx
= c
*invrm
*fv
[XX
];
1285 cfy
= c
*invrm
*fv
[YY
];
1286 cfz
= c
*invrm
*fv
[ZZ
];
1297 fj
[XX
] = ( -rm
[XX
]*rt
[XX
]) * cfx
+ ( rab
[ZZ
]-rm
[YY
]*rt
[XX
]) * cfy
+ (-rab
[YY
]-rm
[ZZ
]*rt
[XX
]) * cfz
;
1298 fj
[YY
] = (-rab
[ZZ
]-rm
[XX
]*rt
[YY
]) * cfx
+ ( -rm
[YY
]*rt
[YY
]) * cfy
+ ( rab
[XX
]-rm
[ZZ
]*rt
[YY
]) * cfz
;
1299 fj
[ZZ
] = ( rab
[YY
]-rm
[XX
]*rt
[ZZ
]) * cfx
+ (-rab
[XX
]-rm
[YY
]*rt
[ZZ
]) * cfy
+ ( -rm
[ZZ
]*rt
[ZZ
]) * cfz
;
1310 fk
[XX
] = ( -rm
[XX
]*rt
[XX
]) * cfx
+ (-a
*rjb
[ZZ
]-rm
[YY
]*rt
[XX
]) * cfy
+ ( a
*rjb
[YY
]-rm
[ZZ
]*rt
[XX
]) * cfz
;
1311 fk
[YY
] = ( a
*rjb
[ZZ
]-rm
[XX
]*rt
[YY
]) * cfx
+ ( -rm
[YY
]*rt
[YY
]) * cfy
+ (-a
*rjb
[XX
]-rm
[ZZ
]*rt
[YY
]) * cfz
;
1312 fk
[ZZ
] = (-a
*rjb
[YY
]-rm
[XX
]*rt
[ZZ
]) * cfx
+ ( a
*rjb
[XX
]-rm
[YY
]*rt
[ZZ
]) * cfy
+ ( -rm
[ZZ
]*rt
[ZZ
]) * cfz
;
1323 fl
[XX
] = ( -rm
[XX
]*rt
[XX
]) * cfx
+ ( b
*rja
[ZZ
]-rm
[YY
]*rt
[XX
]) * cfy
+ (-b
*rja
[YY
]-rm
[ZZ
]*rt
[XX
]) * cfz
;
1324 fl
[YY
] = (-b
*rja
[ZZ
]-rm
[XX
]*rt
[YY
]) * cfx
+ ( -rm
[YY
]*rt
[YY
]) * cfy
+ ( b
*rja
[XX
]-rm
[ZZ
]*rt
[YY
]) * cfz
;
1325 fl
[ZZ
] = ( b
*rja
[YY
]-rm
[XX
]*rt
[ZZ
]) * cfx
+ (-b
*rja
[XX
]-rm
[YY
]*rt
[ZZ
]) * cfy
+ ( -rm
[ZZ
]*rt
[ZZ
]) * cfz
;
1328 f
[ai
][XX
] += fv
[XX
] - fj
[XX
] - fk
[XX
] - fl
[XX
];
1329 f
[ai
][YY
] += fv
[YY
] - fj
[YY
] - fk
[YY
] - fl
[YY
];
1330 f
[ai
][ZZ
] += fv
[ZZ
] - fj
[ZZ
] - fk
[ZZ
] - fl
[ZZ
];
1331 rvec_inc(f
[aj
], fj
);
1332 rvec_inc(f
[ak
], fk
);
1333 rvec_inc(f
[al
], fl
);
1338 ivec_sub(SHIFT_IVEC(g
, av
), SHIFT_IVEC(g
, ai
), di
);
1340 ivec_sub(SHIFT_IVEC(g
, aj
), SHIFT_IVEC(g
, ai
), di
);
1342 ivec_sub(SHIFT_IVEC(g
, ak
), SHIFT_IVEC(g
, ai
), di
);
1344 ivec_sub(SHIFT_IVEC(g
, al
), SHIFT_IVEC(g
, ai
), di
);
1349 svi
= pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xvi
);
1356 if (fshift
&& (svi
!= CENTRAL
|| sij
!= CENTRAL
|| sik
!= CENTRAL
|| sil
!= CENTRAL
))
1358 rvec_dec(fshift
[svi
], fv
);
1359 fshift
[CENTRAL
][XX
] += fv
[XX
] - fj
[XX
] - fk
[XX
] - fl
[XX
];
1360 fshift
[CENTRAL
][YY
] += fv
[YY
] - fj
[YY
] - fk
[YY
] - fl
[YY
];
1361 fshift
[CENTRAL
][ZZ
] += fv
[ZZ
] - fj
[ZZ
] - fk
[ZZ
] - fl
[ZZ
];
1362 rvec_inc(fshift
[sij
], fj
);
1363 rvec_inc(fshift
[sik
], fk
);
1364 rvec_inc(fshift
[sil
], fl
);
1372 pbc_rvec_sub(pbc
, x
[av
], x
[ai
], xiv
);
1374 for (i
= 0; i
< DIM
; i
++)
1376 for (j
= 0; j
< DIM
; j
++)
1378 dxdf
[i
][j
] += -xiv
[i
]*fv
[j
] + xij
[i
]*fj
[j
] + xik
[i
]*fk
[j
] + xil
[i
]*fl
[j
];
1383 /* Total: 207 flops (Yuck!) */
1387 static int spread_vsiten(const t_iatom ia
[], const t_iparams ip
[],
1389 rvec f
[], rvec fshift
[],
1390 const t_pbc
*pbc
, const t_graph
*g
)
1398 n3
= 3*ip
[ia
[0]].vsiten
.n
;
1400 copy_rvec(x
[av
], xv
);
1402 for (i
= 0; i
< n3
; i
+= 3)
1407 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, av
), di
);
1412 siv
= pbc_dx_aiuc(pbc
, x
[ai
], xv
, dx
);
1418 a
= ip
[ia
[i
]].vsiten
.a
;
1419 svmul(a
, f
[av
], fi
);
1420 rvec_inc(f
[ai
], fi
);
1421 if (fshift
&& siv
!= CENTRAL
)
1423 rvec_inc(fshift
[siv
], fi
);
1424 rvec_dec(fshift
[CENTRAL
], fi
);
1433 static int vsite_count(const t_ilist
*ilist
, int ftype
)
1435 if (ftype
== F_VSITEN
)
1437 return ilist
[ftype
].nr
/3;
1441 return ilist
[ftype
].nr
/(1 + interaction_function
[ftype
].nratoms
);
1445 static void spread_vsite_f_thread(const gmx_vsite_t
*vsite
,
1447 rvec f
[], rvec
*fshift
,
1448 gmx_bool VirCorr
, matrix dxdf
,
1449 t_iparams ip
[], const t_ilist ilist
[],
1450 const t_graph
*g
, const t_pbc
*pbc_null
)
1453 const t_pbc
*pbc_null2
;
1454 const int *vsite_pbc
;
1456 bPBCAll
= (pbc_null
!= nullptr && !vsite
->bHaveChargeGroups
);
1458 /* this loop goes backwards to be able to build *
1459 * higher type vsites from lower types */
1460 pbc_null2
= nullptr;
1461 vsite_pbc
= nullptr;
1462 for (int ftype
= c_ftypeVsiteEnd
- 1; ftype
>= c_ftypeVsiteStart
; ftype
--)
1464 if (ilist
[ftype
].nr
== 0)
1470 int nra
= interaction_function
[ftype
].nratoms
;
1472 int nr
= ilist
[ftype
].nr
;
1474 const t_iatom
*ia
= ilist
[ftype
].iatoms
;
1478 pbc_null2
= pbc_null
;
1480 else if (pbc_null
!= nullptr)
1482 vsite_pbc
= vsite
->vsite_pbc_loc
[ftype
- c_ftypeVsiteStart
];
1485 for (int i
= 0; i
< nr
; )
1487 if (vsite_pbc
!= nullptr)
1489 if (vsite_pbc
[i
/(1 + nra
)] > -2)
1491 pbc_null2
= pbc_null
;
1495 pbc_null2
= nullptr;
1501 /* Constants for constructing */
1503 a1
= ip
[tp
].vsite
.a
;
1504 /* Construct the vsite depending on type */
1508 spread_vsite2(ia
, a1
, x
, f
, fshift
, pbc_null2
, g
);
1511 b1
= ip
[tp
].vsite
.b
;
1512 spread_vsite3(ia
, a1
, b1
, x
, f
, fshift
, pbc_null2
, g
);
1515 b1
= ip
[tp
].vsite
.b
;
1516 spread_vsite3FD(ia
, a1
, b1
, x
, f
, fshift
, VirCorr
, dxdf
, pbc_null2
, g
);
1519 b1
= ip
[tp
].vsite
.b
;
1520 spread_vsite3FAD(ia
, a1
, b1
, x
, f
, fshift
, VirCorr
, dxdf
, pbc_null2
, g
);
1523 b1
= ip
[tp
].vsite
.b
;
1524 c1
= ip
[tp
].vsite
.c
;
1525 spread_vsite3OUT(ia
, a1
, b1
, c1
, x
, f
, fshift
, VirCorr
, dxdf
, pbc_null2
, g
);
1528 b1
= ip
[tp
].vsite
.b
;
1529 c1
= ip
[tp
].vsite
.c
;
1530 spread_vsite4FD(ia
, a1
, b1
, c1
, x
, f
, fshift
, VirCorr
, dxdf
, pbc_null2
, g
);
1533 b1
= ip
[tp
].vsite
.b
;
1534 c1
= ip
[tp
].vsite
.c
;
1535 spread_vsite4FDN(ia
, a1
, b1
, c1
, x
, f
, fshift
, VirCorr
, dxdf
, pbc_null2
, g
);
1538 inc
= spread_vsiten(ia
, ip
, x
, f
, fshift
, pbc_null2
, g
);
1541 gmx_fatal(FARGS
, "No such vsite type %d in %s, line %d",
1542 ftype
, __FILE__
, __LINE__
);
1544 clear_rvec(f
[ia
[1]]);
1546 /* Increment loop variables */
1554 /*! \brief Clears the task force buffer elements that are written by task idTask */
1555 static void clearTaskForceBufferUsedElements(InterdependentTask
*idTask
)
1557 int ntask
= idTask
->spreadTask
.size();
1558 for (int ti
= 0; ti
< ntask
; ti
++)
1560 const AtomIndex
*atomList
= &idTask
->atomIndex
[idTask
->spreadTask
[ti
]];
1561 int natom
= atomList
->atom
.size();
1562 RVec
*force
= idTask
->force
.data();
1563 for (int i
= 0; i
< natom
; i
++)
1565 clear_rvec(force
[atomList
->atom
[i
]]);
1570 void spread_vsite_f(const gmx_vsite_t
*vsite
,
1571 const rvec
* gmx_restrict x
,
1572 rvec
* gmx_restrict f
, rvec
* gmx_restrict fshift
,
1573 gmx_bool VirCorr
, matrix vir
,
1574 t_nrnb
*nrnb
, const t_idef
*idef
,
1575 int ePBC
, gmx_bool bMolPBC
, const t_graph
*g
, const matrix box
,
1578 t_pbc pbc
, *pbc_null
;
1580 /* We only need to do pbc when we have inter-cg vsites */
1581 if ((DOMAINDECOMP(cr
) || bMolPBC
) && vsite
->n_intercg_vsite
)
1583 /* This is wasting some CPU time as we now do this multiple times
1586 pbc_null
= set_pbc_dd(&pbc
, ePBC
, cr
->dd
? cr
->dd
->nc
: nullptr, FALSE
, box
);
1593 if (DOMAINDECOMP(cr
))
1595 dd_clear_f_vsites(cr
->dd
, f
);
1598 if (vsite
->nthreads
== 1)
1602 clear_mat(vsite
->tData
[0]->dxdf
);
1604 spread_vsite_f_thread(vsite
,
1606 VirCorr
, vsite
->tData
[0]->dxdf
,
1607 idef
->iparams
, idef
->il
,
1612 /* First spread the vsites that might depend on non-local vsites */
1615 clear_mat(vsite
->tData
[vsite
->nthreads
]->dxdf
);
1617 spread_vsite_f_thread(vsite
,
1619 VirCorr
, vsite
->tData
[vsite
->nthreads
]->dxdf
,
1621 vsite
->tData
[vsite
->nthreads
]->ilist
,
1624 #pragma omp parallel num_threads(vsite->nthreads)
1628 int thread
= gmx_omp_get_thread_num();
1629 VsiteThread
*tData
= vsite
->tData
[thread
];
1632 if (thread
== 0 || fshift
== nullptr)
1638 fshift_t
= tData
->fshift
;
1640 for (int i
= 0; i
< SHIFTS
; i
++)
1642 clear_rvec(fshift_t
[i
]);
1647 clear_mat(tData
->dxdf
);
1650 if (tData
->useInterdependentTask
)
1652 /* Spread the vsites that spread outside our local range.
1653 * This is done using a thread-local force buffer force.
1654 * First we need to copy the input vsite forces to force.
1656 InterdependentTask
*idTask
= &tData
->idTask
;
1658 /* Clear the buffer elements set by our task during
1659 * the last call to spread_vsite_f.
1661 clearTaskForceBufferUsedElements(idTask
);
1663 int nvsite
= idTask
->vsite
.size();
1664 for (int i
= 0; i
< nvsite
; i
++)
1666 copy_rvec(f
[idTask
->vsite
[i
]],
1667 idTask
->force
[idTask
->vsite
[i
]]);
1669 spread_vsite_f_thread(vsite
,
1670 x
, as_rvec_array(idTask
->force
.data()), fshift_t
,
1671 VirCorr
, tData
->dxdf
,
1673 tData
->idTask
.ilist
,
1676 /* We need a barrier before reducing forces below
1677 * that have been produced by a different thread above.
1681 /* Loop over all thread task and reduce forces they
1682 * produced on atoms that fall in our range.
1683 * Note that atomic reduction would be a simpler solution,
1684 * but that might not have good support on all platforms.
1686 int ntask
= idTask
->reduceTask
.size();
1687 for (int ti
= 0; ti
< ntask
; ti
++)
1689 const InterdependentTask
*idt_foreign
= &vsite
->tData
[idTask
->reduceTask
[ti
]]->idTask
;
1690 const AtomIndex
*atomList
= &idt_foreign
->atomIndex
[thread
];
1691 const RVec
*f_foreign
= idt_foreign
->force
.data();
1693 int natom
= atomList
->atom
.size();
1694 for (int i
= 0; i
< natom
; i
++)
1696 int ind
= atomList
->atom
[i
];
1697 rvec_inc(f
[ind
], f_foreign
[ind
]);
1698 /* Clearing of f_foreign is done at the next step */
1701 /* Clear the vsite forces, both in f and force */
1702 for (int i
= 0; i
< nvsite
; i
++)
1704 int ind
= tData
->idTask
.vsite
[i
];
1706 clear_rvec(tData
->idTask
.force
[ind
]);
1710 /* Spread the vsites that spread locally only */
1711 spread_vsite_f_thread(vsite
,
1713 VirCorr
, tData
->dxdf
,
1718 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
1721 if (fshift
!= nullptr)
1723 for (int th
= 1; th
< vsite
->nthreads
; th
++)
1725 for (int i
= 0; i
< SHIFTS
; i
++)
1727 rvec_inc(fshift
[i
], vsite
->tData
[th
]->fshift
[i
]);
1735 for (int th
= 0; th
< (vsite
->nthreads
== 1 ? 1 : vsite
->nthreads
+ 1); th
++)
1737 for (int i
= 0; i
< DIM
; i
++)
1739 for (int j
= 0; j
< DIM
; j
++)
1741 vir
[i
][j
] += -0.5*vsite
->tData
[th
]->dxdf
[i
][j
];
1747 if (DOMAINDECOMP(cr
))
1749 dd_move_f_vsites(cr
->dd
, f
, fshift
);
1752 inc_nrnb(nrnb
, eNR_VSITE2
, vsite_count(idef
->il
, F_VSITE2
));
1753 inc_nrnb(nrnb
, eNR_VSITE3
, vsite_count(idef
->il
, F_VSITE3
));
1754 inc_nrnb(nrnb
, eNR_VSITE3FD
, vsite_count(idef
->il
, F_VSITE3FD
));
1755 inc_nrnb(nrnb
, eNR_VSITE3FAD
, vsite_count(idef
->il
, F_VSITE3FAD
));
1756 inc_nrnb(nrnb
, eNR_VSITE3OUT
, vsite_count(idef
->il
, F_VSITE3OUT
));
1757 inc_nrnb(nrnb
, eNR_VSITE4FD
, vsite_count(idef
->il
, F_VSITE4FD
));
1758 inc_nrnb(nrnb
, eNR_VSITE4FDN
, vsite_count(idef
->il
, F_VSITE4FDN
));
1759 inc_nrnb(nrnb
, eNR_VSITEN
, vsite_count(idef
->il
, F_VSITEN
));
1762 static int *atom2cg(t_block
*cgs
)
1766 snew(a2cg
, cgs
->index
[cgs
->nr
]);
1767 for (int cg
= 0; cg
< cgs
->nr
; cg
++)
1769 for (int i
= cgs
->index
[cg
]; i
< cgs
->index
[cg
+1]; i
++)
1778 int count_intercg_vsites(const gmx_mtop_t
*mtop
)
1780 gmx_molblock_t
*molb
;
1781 gmx_moltype_t
*molt
;
1783 int n_intercg_vsite
;
1785 n_intercg_vsite
= 0;
1786 for (int mb
= 0; mb
< mtop
->nmolblock
; mb
++)
1788 molb
= &mtop
->molblock
[mb
];
1789 molt
= &mtop
->moltype
[molb
->type
];
1791 a2cg
= atom2cg(&molt
->cgs
);
1792 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
1794 int nral
= NRAL(ftype
);
1795 t_ilist
*il
= &molt
->ilist
[ftype
];
1796 const t_iatom
*ia
= il
->iatoms
;
1797 for (int i
= 0; i
< il
->nr
; i
+= 1 + nral
)
1799 int cg
= a2cg
[ia
[1+i
]];
1800 for (int a
= 1; a
< nral
; a
++)
1802 if (a2cg
[ia
[1+a
]] != cg
)
1804 n_intercg_vsite
+= molb
->nmol
;
1813 return n_intercg_vsite
;
1816 static int **get_vsite_pbc(const t_iparams
*iparams
, const t_ilist
*ilist
,
1817 const t_atom
*atom
, const t_mdatoms
*md
,
1818 const t_block
*cgs
, const int *a2cg
)
1822 /* Make an array that tells if the pbc of an atom is set */
1823 snew(pbc_set
, cgs
->index
[cgs
->nr
]);
1824 /* PBC is set for all non vsites */
1825 for (int a
= 0; a
< cgs
->index
[cgs
->nr
]; a
++)
1827 if ((atom
&& atom
[a
].ptype
!= eptVSite
) ||
1828 (md
&& md
->ptype
[a
] != eptVSite
))
1835 snew(vsite_pbc
, F_VSITEN
-F_VSITE2
+1);
1837 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
1840 int nral
= NRAL(ftype
);
1841 const t_ilist
*il
= &ilist
[ftype
];
1842 const t_iatom
*ia
= il
->iatoms
;
1845 snew(vsite_pbc
[ftype
-F_VSITE2
], il
->nr
/(1 + nral
));
1846 vsite_pbc_f
= vsite_pbc
[ftype
-F_VSITE2
];
1851 int vsi
= i
/(1 + nral
);
1852 t_iatom vsite
= ia
[i
+1];
1853 int cg_v
= a2cg
[vsite
];
1854 /* A value of -2 signals that this vsite and its contructing
1855 * atoms are all within the same cg, so no pbc is required.
1857 vsite_pbc_f
[vsi
] = -2;
1858 /* Check if constructing atoms are outside the vsite's cg */
1860 if (ftype
== F_VSITEN
)
1862 nc3
= 3*iparams
[ia
[i
]].vsiten
.n
;
1863 for (int j
= 0; j
< nc3
; j
+= 3)
1865 if (a2cg
[ia
[i
+j
+2]] != cg_v
)
1867 vsite_pbc_f
[vsi
] = -1;
1873 for (int a
= 1; a
< nral
; a
++)
1875 if (a2cg
[ia
[i
+1+a
]] != cg_v
)
1877 vsite_pbc_f
[vsi
] = -1;
1881 if (vsite_pbc_f
[vsi
] == -1)
1883 /* Check if this is the first processed atom of a vsite only cg */
1884 gmx_bool bViteOnlyCG_and_FirstAtom
= TRUE
;
1885 for (int a
= cgs
->index
[cg_v
]; a
< cgs
->index
[cg_v
+1]; a
++)
1887 /* Non-vsites already have pbc set, so simply check for pbc_set */
1890 bViteOnlyCG_and_FirstAtom
= FALSE
;
1894 if (bViteOnlyCG_and_FirstAtom
)
1896 /* First processed atom of a vsite only charge group.
1897 * The pbc of the input coordinates to construct_vsites
1898 * should be preserved.
1900 vsite_pbc_f
[vsi
] = vsite
;
1902 else if (cg_v
!= a2cg
[ia
[1+i
+1]])
1904 /* This vsite has a different charge group index
1905 * than it's first constructing atom
1906 * and the charge group has more than one atom,
1907 * search for the first normal particle
1908 * or vsite that already had its pbc defined.
1909 * If nothing is found, use full pbc for this vsite.
1911 for (int a
= cgs
->index
[cg_v
]; a
< cgs
->index
[cg_v
+1]; a
++)
1913 if (a
!= vsite
&& pbc_set
[a
])
1915 vsite_pbc_f
[vsi
] = a
;
1918 fprintf(debug
, "vsite %d match pbc with atom %d\n",
1926 fprintf(debug
, "vsite atom %d cg %d - %d pbc atom %d\n",
1927 vsite
+1, cgs
->index
[cg_v
]+1, cgs
->index
[cg_v
+1],
1928 vsite_pbc_f
[vsi
]+1);
1932 if (ftype
== F_VSITEN
)
1934 /* The other entries in vsite_pbc_f are not used for center vsites */
1942 /* This vsite now has its pbc defined */
1954 gmx_vsite_t
*init_vsite(const gmx_mtop_t
*mtop
, t_commrec
*cr
,
1955 gmx_bool bSerial_NoPBC
)
1958 gmx_moltype_t
*molt
;
1960 /* check if there are vsites */
1962 for (int ftype
= 0; ftype
< F_NRE
; ftype
++)
1964 if (interaction_function
[ftype
].flags
& IF_VSITE
)
1966 GMX_ASSERT(ftype
>= c_ftypeVsiteStart
&& ftype
< c_ftypeVsiteEnd
, "c_ftypeVsiteStart and/or c_ftypeVsiteEnd do not have correct values");
1968 nvsite
+= gmx_mtop_ftype_count(mtop
, ftype
);
1972 GMX_ASSERT(ftype
< c_ftypeVsiteStart
|| ftype
>= c_ftypeVsiteEnd
, "c_ftypeVsiteStart and/or c_ftypeVsiteEnd do not have correct values");
1983 vsite
->n_intercg_vsite
= count_intercg_vsites(mtop
);
1985 vsite
->bHaveChargeGroups
= (ncg_mtop(mtop
) < mtop
->natoms
);
1987 /* If we don't have charge groups, the vsite follows its own pbc */
1988 if (!bSerial_NoPBC
&&
1989 vsite
->bHaveChargeGroups
&&
1990 vsite
->n_intercg_vsite
> 0 && DOMAINDECOMP(cr
))
1992 vsite
->nvsite_pbc_molt
= mtop
->nmoltype
;
1993 snew(vsite
->vsite_pbc_molt
, vsite
->nvsite_pbc_molt
);
1994 for (int mt
= 0; mt
< mtop
->nmoltype
; mt
++)
1996 molt
= &mtop
->moltype
[mt
];
1997 /* Make an atom to charge group index */
1998 int *a2cg
= atom2cg(&molt
->cgs
);
1999 vsite
->vsite_pbc_molt
[mt
] = get_vsite_pbc(mtop
->ffparams
.iparams
,
2001 molt
->atoms
.atom
, nullptr,
2006 snew(vsite
->vsite_pbc_loc_nalloc
, c_ftypeVsiteEnd
- c_ftypeVsiteStart
);
2007 snew(vsite
->vsite_pbc_loc
, c_ftypeVsiteEnd
- c_ftypeVsiteStart
);
2012 vsite
->nthreads
= 1;
2016 vsite
->nthreads
= gmx_omp_nthreads_get(emntVSITE
);
2020 /* We need one extra thread data structure for the overlap vsites */
2021 snew(vsite
->tData
, vsite
->nthreads
+ 1);
2022 #pragma omp parallel for num_threads(vsite->nthreads) schedule(static)
2023 for (int thread
= 0; thread
< vsite
->nthreads
; thread
++)
2027 vsite
->tData
[thread
] = new VsiteThread
;
2029 InterdependentTask
*idTask
= &vsite
->tData
[thread
]->idTask
;
2031 idTask
->atomIndex
.resize(vsite
->nthreads
);
2033 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2035 if (vsite
->nthreads
> 1)
2037 vsite
->tData
[vsite
->nthreads
] = new VsiteThread
;
2041 vsite
->taskIndex
= nullptr;
2042 vsite
->taskIndexNalloc
= 0;
2047 static gmx_inline
void flagAtom(InterdependentTask
*idTask
, int atom
,
2048 int thread
, int nthread
, int natperthread
)
2050 if (!idTask
->use
[atom
])
2052 idTask
->use
[atom
] = true;
2053 thread
= atom
/natperthread
;
2054 /* Assign all non-local atom force writes to thread 0 */
2055 if (thread
>= nthread
)
2059 idTask
->atomIndex
[thread
].atom
.push_back(atom
);
2063 /*\brief Here we try to assign all vsites that are in our local range.
2065 * Our task local atom range is tData->rangeStart - tData->rangeEnd.
2066 * Vsites that depend only on local atoms, as indicated by taskIndex[]==thread,
2067 * are assigned to task tData->ilist. Vsites that depend on non-local atoms
2068 * but not on other vsites are assigned to task tData->id_task.ilist.
2069 * taskIndex[] is set for all vsites in our range, either to our local tasks
2070 * or to the single last task as taskIndex[]=2*nthreads.
2072 static void assignVsitesToThread(VsiteThread
*tData
,
2077 const t_ilist
*ilist
,
2078 const t_iparams
*ip
,
2079 const unsigned short *ptype
)
2081 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
2083 tData
->ilist
[ftype
].nr
= 0;
2084 tData
->idTask
.ilist
[ftype
].nr
= 0;
2086 int nral1
= 1 + NRAL(ftype
);
2088 t_iatom
*iat
= ilist
[ftype
].iatoms
;
2089 for (int i
= 0; i
< ilist
[ftype
].nr
; )
2091 if (ftype
== F_VSITEN
)
2093 /* The 3 below is from 1+NRAL(ftype)=3 */
2094 inc
= ip
[iat
[i
]].vsiten
.n
*3;
2097 if (iat
[1 + i
] < tData
->rangeStart
||
2098 iat
[1 + i
] >= tData
->rangeEnd
)
2100 /* This vsite belongs to a different thread */
2105 /* We would like to assign this vsite to task thread,
2106 * but it might depend on atoms outside the atom range of thread
2107 * or on another vsite not assigned to task thread.
2110 if (ftype
!= F_VSITEN
)
2112 for (int j
= i
+ 2; j
< i
+ nral1
; j
++)
2114 /* Do a range check to avoid a harmless race on taskIndex */
2115 if (iat
[j
] < tData
->rangeStart
||
2116 iat
[j
] >= tData
->rangeEnd
||
2117 taskIndex
[iat
[j
]] != thread
)
2119 if (!tData
->useInterdependentTask
||
2120 ptype
[iat
[j
]] == eptVSite
)
2122 /* At least one constructing atom is a vsite
2123 * that is not assigned to the same thread.
2124 * Put this vsite into a separate task.
2130 /* There are constructing atoms outside our range,
2131 * put this vsite into a second task to be executed
2132 * on the same thread. During construction no barrier
2133 * is needed between the two tasks on the same thread.
2134 * During spreading we need to run this task with
2135 * an additional thread-local intermediate force buffer
2136 * (or atomic reduction) and a barrier between the two
2139 task
= nthread
+ thread
;
2145 for (int j
= i
+ 2; j
< i
+ inc
; j
+= 3)
2147 /* Do a range check to avoid a harmless race on taskIndex */
2148 if (iat
[j
] < tData
->rangeStart
||
2149 iat
[j
] >= tData
->rangeEnd
||
2150 taskIndex
[iat
[j
]] != thread
)
2152 GMX_ASSERT(ptype
[iat
[j
]] != eptVSite
, "A vsite to be assigned in assignVsitesToThread has a vsite as a constructing atom that does not belong to our task, such vsites should be assigned to the single 'master' task");
2154 task
= nthread
+ thread
;
2159 /* Update this vsite's thread index entry */
2160 taskIndex
[iat
[1+i
]] = task
;
2162 if (task
== thread
|| task
== nthread
+ thread
)
2164 /* Copy this vsite to the thread data struct of thread */
2168 il_task
= &tData
->ilist
[ftype
];
2172 il_task
= &tData
->idTask
.ilist
[ftype
];
2174 /* Ensure we have sufficient memory allocated */
2175 if (il_task
->nr
+ inc
> il_task
->nalloc
)
2177 il_task
->nalloc
= over_alloc_large(il_task
->nr
+ inc
);
2178 srenew(il_task
->iatoms
, il_task
->nalloc
);
2180 /* Copy the vsite data to the thread-task local array */
2181 for (int j
= i
; j
< i
+ inc
; j
++)
2183 il_task
->iatoms
[il_task
->nr
++] = iat
[j
];
2185 if (task
== nthread
+ thread
)
2187 /* This vsite write outside our own task force block.
2188 * Put it into the interdependent task list and flag
2189 * the atoms involved for reduction.
2191 tData
->idTask
.vsite
.push_back(iat
[i
+ 1]);
2192 if (ftype
!= F_VSITEN
)
2194 for (int j
= i
+ 2; j
< i
+ nral1
; j
++)
2196 flagAtom(&tData
->idTask
, iat
[j
],
2197 thread
, nthread
, natperthread
);
2202 for (int j
= i
+ 2; j
< i
+ inc
; j
+= 3)
2204 flagAtom(&tData
->idTask
, iat
[j
],
2205 thread
, nthread
, natperthread
);
2216 /*! \brief Assign all vsites with taskIndex[]==task to task tData */
2217 static void assignVsitesToSingleTask(VsiteThread
*tData
,
2219 const int *taskIndex
,
2220 const t_ilist
*ilist
,
2221 const t_iparams
*ip
)
2223 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
2225 tData
->ilist
[ftype
].nr
= 0;
2226 tData
->idTask
.ilist
[ftype
].nr
= 0;
2228 int nral1
= 1 + NRAL(ftype
);
2230 t_iatom
*iat
= ilist
[ftype
].iatoms
;
2231 t_ilist
*il_task
= &tData
->ilist
[ftype
];
2233 for (int i
= 0; i
< ilist
[ftype
].nr
; )
2235 if (ftype
== F_VSITEN
)
2237 /* The 3 below is from 1+NRAL(ftype)=3 */
2238 inc
= ip
[iat
[i
]].vsiten
.n
*3;
2240 /* Check if the vsite is assigned to our task */
2241 if (taskIndex
[iat
[1 + i
]] == task
)
2243 /* Ensure we have sufficient memory allocated */
2244 if (il_task
->nr
+ inc
> il_task
->nalloc
)
2246 il_task
->nalloc
= over_alloc_large(il_task
->nr
+ inc
);
2247 srenew(il_task
->iatoms
, il_task
->nalloc
);
2249 /* Copy the vsite data to the thread-task local array */
2250 for (int j
= i
; j
< i
+ inc
; j
++)
2252 il_task
->iatoms
[il_task
->nr
++] = iat
[j
];
2261 void split_vsites_over_threads(const t_ilist
*ilist
,
2262 const t_iparams
*ip
,
2263 const t_mdatoms
*mdatoms
,
2264 gmx_bool bLimitRange
,
2267 int vsite_atom_range
, natperthread
;
2269 if (vsite
->nthreads
== 1)
2275 /* The current way of distributing the vsites over threads in primitive.
2276 * We divide the atom range 0 - natoms_in_vsite uniformly over threads,
2277 * without taking into account how the vsites are distributed.
2278 * Without domain decomposition we bLimitRange=TRUE and we at least
2279 * tighten the upper bound of the range (useful for common systems
2280 * such as a vsite-protein in 3-site water).
2281 * With domain decomposition, as long as the vsites are distributed
2282 * uniformly in each domain along the major dimension, usually x,
2283 * it will also perform well.
2287 vsite_atom_range
= -1;
2288 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
2291 if (ftype
!= F_VSITEN
)
2293 int nral1
= 1 + NRAL(ftype
);
2294 const t_iatom
*iat
= ilist
[ftype
].iatoms
;
2295 for (int i
= 0; i
< ilist
[ftype
].nr
; i
+= nral1
)
2297 for (int j
= i
+ 1; j
< i
+ nral1
; j
++)
2299 vsite_atom_range
= std::max(vsite_atom_range
, iat
[j
]);
2307 const t_iatom
*iat
= ilist
[ftype
].iatoms
;
2310 while (i
< ilist
[ftype
].nr
)
2312 /* The 3 below is from 1+NRAL(ftype)=3 */
2313 vs_ind_end
= i
+ ip
[iat
[i
]].vsiten
.n
*3;
2315 vsite_atom_range
= std::max(vsite_atom_range
, iat
[i
+1]);
2316 while (i
< vs_ind_end
)
2318 vsite_atom_range
= std::max(vsite_atom_range
, iat
[i
+2]);
2326 natperthread
= (vsite_atom_range
+ vsite
->nthreads
- 1)/vsite
->nthreads
;
2330 /* Any local or not local atom could be involved in virtual sites.
2331 * But since we usually have very few non-local virtual sites
2332 * (only non-local vsites that depend on local vsites),
2333 * we distribute the local atom range equally over the threads.
2334 * When assigning vsites to threads, we should take care that the last
2335 * threads also covers the non-local range.
2337 vsite_atom_range
= mdatoms
->nr
;
2338 natperthread
= (mdatoms
->homenr
+ vsite
->nthreads
- 1)/vsite
->nthreads
;
2343 fprintf(debug
, "virtual site thread dist: natoms %d, range %d, natperthread %d\n", mdatoms
->nr
, vsite_atom_range
, natperthread
);
2346 /* To simplify the vsite assignment, we make an index which tells us
2347 * to which task particles, both non-vsites and vsites, are assigned.
2349 if (mdatoms
->nr
> vsite
->taskIndexNalloc
)
2351 vsite
->taskIndexNalloc
= over_alloc_large(mdatoms
->nr
);
2352 srenew(vsite
->taskIndex
, vsite
->taskIndexNalloc
);
2355 /* Initialize the task index array. Here we assign the non-vsite
2356 * particles to task=thread, so we easily figure out if vsites
2357 * depend on local and/or non-local particles in assignVsitesToThread.
2359 int *taskIndex
= vsite
->taskIndex
;
2362 for (int i
= 0; i
< mdatoms
->nr
; i
++)
2364 if (mdatoms
->ptype
[i
] == eptVSite
)
2366 /* vsites are not assigned to a task yet */
2371 /* assign non-vsite particles to task thread */
2372 taskIndex
[i
] = thread
;
2374 if (i
== (thread
+ 1)*natperthread
&& thread
< vsite
->nthreads
)
2381 #pragma omp parallel num_threads(vsite->nthreads)
2385 int thread
= gmx_omp_get_thread_num();
2386 VsiteThread
*tData
= vsite
->tData
[thread
];
2388 /* Clear the buffer use flags that were set before */
2389 if (tData
->useInterdependentTask
)
2391 InterdependentTask
*idTask
= &tData
->idTask
;
2393 /* To avoid an extra OpenMP barrier in spread_vsite_f,
2394 * we clear the force buffer at the next step,
2395 * so we need to do it here as well.
2397 clearTaskForceBufferUsedElements(idTask
);
2399 idTask
->vsite
.resize(0);
2400 for (int t
= 0; t
< vsite
->nthreads
; t
++)
2402 AtomIndex
*atomIndex
= &idTask
->atomIndex
[t
];
2403 int natom
= atomIndex
->atom
.size();
2404 for (int i
= 0; i
< natom
; i
++)
2406 idTask
->use
[atomIndex
->atom
[i
]] = false;
2408 atomIndex
->atom
.resize(0);
2413 /* To avoid large f_buf allocations of #threads*vsite_atom_range
2414 * we don't use task2 with more than 200000 atoms. This doesn't
2415 * affect performance, since with such a large range relatively few
2416 * vsites will end up in the separate task.
2417 * Note that useTask2 should be the same for all threads.
2419 tData
->useInterdependentTask
= (vsite_atom_range
<= 200000);
2420 if (tData
->useInterdependentTask
)
2422 size_t natoms_use_in_vsites
= vsite_atom_range
;
2423 InterdependentTask
*idTask
= &tData
->idTask
;
2424 /* To avoid resizing and re-clearing every nstlist steps,
2425 * we never down size the force buffer.
2427 if (natoms_use_in_vsites
> idTask
->force
.size() ||
2428 natoms_use_in_vsites
> idTask
->use
.size())
2430 idTask
->force
.resize(natoms_use_in_vsites
, { 0, 0, 0 });
2431 idTask
->use
.resize(natoms_use_in_vsites
, false);
2435 /* Assign all vsites that can execute independently on threads */
2436 tData
->rangeStart
= thread
*natperthread
;
2437 if (thread
< vsite
->nthreads
- 1)
2439 tData
->rangeEnd
= (thread
+ 1)*natperthread
;
2443 /* The last thread should cover up to the end of the range */
2444 tData
->rangeEnd
= mdatoms
->nr
;
2446 assignVsitesToThread(tData
,
2447 thread
, vsite
->nthreads
,
2450 ilist
, ip
, mdatoms
->ptype
);
2452 if (tData
->useInterdependentTask
)
2454 /* In the worst case, all tasks write to force ranges of
2455 * all other tasks, leading to #tasks^2 scaling (this is only
2456 * the overhead, the actual flops remain constant).
2457 * But in most cases there is far less coupling. To improve
2458 * scaling at high thread counts we therefore construct
2459 * an index to only loop over the actually affected tasks.
2461 InterdependentTask
*idTask
= &tData
->idTask
;
2463 /* Ensure assignVsitesToThread finished on other threads */
2466 idTask
->spreadTask
.resize(0);
2467 idTask
->reduceTask
.resize(0);
2468 for (int t
= 0; t
< vsite
->nthreads
; t
++)
2470 /* Do we write to the force buffer of task t? */
2471 if (idTask
->atomIndex
[t
].atom
.size() > 0)
2473 idTask
->spreadTask
.push_back(t
);
2475 /* Does task t write to our force buffer? */
2476 if (vsite
->tData
[t
]->idTask
.atomIndex
[thread
].atom
.size() > 0)
2478 idTask
->reduceTask
.push_back(t
);
2483 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2485 /* Assign all remaining vsites, that will have taskIndex[]=2*vsite->nthreads,
2486 * to a single task that will not run in parallel with other tasks.
2488 assignVsitesToSingleTask(vsite
->tData
[vsite
->nthreads
],
2493 if (debug
&& vsite
->nthreads
> 1)
2495 fprintf(debug
, "virtual site useInterdependentTask %d, nuse:\n",
2496 vsite
->tData
[0]->useInterdependentTask
);
2497 for (int th
= 0; th
< vsite
->nthreads
+ 1; th
++)
2499 fprintf(debug
, " %4d", vsite
->tData
[th
]->idTask
.nuse
);
2501 fprintf(debug
, "\n");
2503 for (int ftype
= c_ftypeVsiteStart
; ftype
< c_ftypeVsiteEnd
; ftype
++)
2505 if (ilist
[ftype
].nr
> 0)
2507 fprintf(debug
, "%-20s thread dist:",
2508 interaction_function
[ftype
].longname
);
2509 for (int th
= 0; th
< vsite
->nthreads
+ 1; th
++)
2511 fprintf(debug
, " %4d %4d ",
2512 vsite
->tData
[th
]->ilist
[ftype
].nr
,
2513 vsite
->tData
[th
]->idTask
.ilist
[ftype
].nr
);
2515 fprintf(debug
, "\n");
2521 int nrOrig
= vsiteIlistNrCount(ilist
);
2523 for (int th
= 0; th
< vsite
->nthreads
+ 1; th
++)
2526 vsiteIlistNrCount(vsite
->tData
[th
]->ilist
) +
2527 vsiteIlistNrCount(vsite
->tData
[th
]->idTask
.ilist
);
2529 GMX_ASSERT(nrThreaded
== nrOrig
, "The number of virtual sites assigned to all thread task has to match the total number of virtual sites");
2533 void set_vsite_top(gmx_vsite_t
*vsite
, gmx_localtop_t
*top
, t_mdatoms
*md
,
2536 if (vsite
->n_intercg_vsite
> 0)
2538 if (vsite
->bHaveChargeGroups
)
2540 /* Make an atom to charge group index */
2541 int *a2cg
= atom2cg(&top
->cgs
);
2542 vsite
->vsite_pbc_loc
= get_vsite_pbc(top
->idef
.iparams
,
2543 top
->idef
.il
, nullptr, md
,
2550 if (vsite
->nthreads
> 1)
2552 if (vsite
->bHaveChargeGroups
)
2554 gmx_fatal(FARGS
, "The combination of threading, virtual sites and charge groups is not implemented");
2557 split_vsites_over_threads(top
->idef
.il
, top
->idef
.iparams
,
2558 md
, !DOMAINDECOMP(cr
), vsite
);