From f86bed5e44271e1807d95829087c454db0fe164e Mon Sep 17 00:00:00 2001
From: Berk Hess <hess@kth.se>
Date: Tue, 19 Aug 2014 10:10:48 +0200
Subject: [PATCH] Fixed bugs in vsiteN with OpenMP

Fixes #1579.

Change-Id: I42d234f4ad6a94e8f7b6b8236ea119860dd9f7ab
---
 include/vsite.h    |  1 +
 src/mdlib/domdec.c |  3 ++-
 src/mdlib/vsite.c  | 52 +++++++++++++++++++++++++++++++++-------------
 3 files changed, 41 insertions(+), 15 deletions(-)

diff --git a/include/vsite.h b/include/vsite.h
index edeb426c1c..ecff615006 100644
--- a/include/vsite.h
+++ b/include/vsite.h
@@ -131,6 +131,7 @@ gmx_vsite_t *init_vsite(gmx_mtop_t *mtop, t_commrec *cr,
  */
 
 void split_vsites_over_threads(const t_ilist   *ilist,
+                               const t_iparams *ip,
                                const t_mdatoms *mdatoms,
                                gmx_bool         bLimitRange,
                                gmx_vsite_t     *vsite);
diff --git a/src/mdlib/domdec.c b/src/mdlib/domdec.c
index d4be717774..ac8627bb57 100644
--- a/src/mdlib/domdec.c
+++ b/src/mdlib/domdec.c
@@ -9808,7 +9808,8 @@ void dd_partition_system(FILE                *fplog,
     if (vsite != NULL)
     {
         /* Now we have updated mdatoms, we can do the last vsite bookkeeping */
-        split_vsites_over_threads(top_local->idef.il, mdatoms, FALSE, vsite);
+        split_vsites_over_threads(top_local->idef.il, top_local->idef.iparams,
+                                  mdatoms, FALSE, vsite);
     }
 
     if (shellfc)
diff --git a/src/mdlib/vsite.c b/src/mdlib/vsite.c
index a2d541ea25..cf10fe12b1 100644
--- a/src/mdlib/vsite.c
+++ b/src/mdlib/vsite.c
@@ -1988,6 +1988,7 @@ static void prepare_vsite_thread(const t_ilist      *ilist,
 }
 
 void split_vsites_over_threads(const t_ilist   *ilist,
+                               const t_iparams *ip,
                                const t_mdatoms *mdatoms,
                                gmx_bool         bLimitRange,
                                gmx_vsite_t     *vsite)
@@ -2026,16 +2027,38 @@ void split_vsites_over_threads(const t_ilist   *ilist,
         vsite_atom_range = -1;
         for (ftype = 0; ftype < F_NRE; ftype++)
         {
-            if ((interaction_function[ftype].flags & IF_VSITE) &&
-                ftype != F_VSITEN)
+            if (interaction_function[ftype].flags & IF_VSITE)
             {
-                nral1 = 1 + NRAL(ftype);
-                iat   = ilist[ftype].iatoms;
-                for (i = 0; i < ilist[ftype].nr; i += nral1)
+                if (ftype != F_VSITEN)
                 {
-                    for (j = i+1; j < i+nral1; j++)
+                    nral1 = 1 + NRAL(ftype);
+                    iat   = ilist[ftype].iatoms;
+                    for (i = 0; i < ilist[ftype].nr; i += nral1)
                     {
-                        vsite_atom_range = max(vsite_atom_range, iat[j]);
+                        for (j = i + 1; j < i + nral1; j++)
+                        {
+                            vsite_atom_range = max(vsite_atom_range, iat[j]);
+                        }
+                    }
+                }
+                else
+                {
+                    int vs_ind_end;
+
+                    iat = ilist[ftype].iatoms;
+
+                    i = 0;
+                    while (i < ilist[ftype].nr)
+                    {
+                        /* The 3 below is from 1+NRAL(ftype)=3 */
+                        vs_ind_end = i + ip[iat[i]].vsiten.n*3;
+
+                        vsite_atom_range = max(vsite_atom_range, iat[i+1]);
+                        while (i < vs_ind_end)
+                        {
+                            vsite_atom_range = max(vsite_atom_range, iat[i+2]);
+                            i += 3;
+                        }
                     }
                 }
             }
@@ -2083,8 +2106,7 @@ void split_vsites_over_threads(const t_ilist   *ilist,
 
     for (ftype = 0; ftype < F_NRE; ftype++)
     {
-        if ((interaction_function[ftype].flags & IF_VSITE) &&
-            ftype != F_VSITEN)
+        if (interaction_function[ftype].flags & IF_VSITE)
         {
             nral1 = 1 + NRAL(ftype);
             inc   = nral1;
@@ -2098,7 +2120,7 @@ void split_vsites_over_threads(const t_ilist   *ilist,
                  */
                 if (ftype != F_VSITEN)
                 {
-                    for (j = i+2; j < i+nral1; j++)
+                    for (j = i + 2; j < i + nral1; j++)
                     {
                         if (th_ind[iat[j]] != th)
                         {
@@ -2111,8 +2133,9 @@ void split_vsites_over_threads(const t_ilist   *ilist,
                 }
                 else
                 {
-                    inc = iat[i];
-                    for (j = i+2; j < i+inc; j += 3)
+                    /* The 3 below is from 1+NRAL(ftype)=3 */
+                    inc = ip[iat[i]].vsiten.n*3;
+                    for (j = i + 2; j < i + inc; j += 3)
                     {
                         if (th_ind[iat[j]] != th)
                         {
@@ -2122,7 +2145,7 @@ void split_vsites_over_threads(const t_ilist   *ilist,
                 }
                 /* Copy this vsite to the thread data struct of thread th */
                 il_th = &vsite->tdata[th].ilist[ftype];
-                for (j = i; j < i+inc; j++)
+                for (j = i; j < i + inc; j++)
                 {
                     il_th->iatoms[il_th->nr++] = iat[j];
                 }
@@ -2185,6 +2208,7 @@ void set_vsite_top(gmx_vsite_t *vsite, gmx_localtop_t *top, t_mdatoms *md,
             gmx_incons("Can not use threads virtual sites combined with charge groups or particle decomposition");
         }
 
-        split_vsites_over_threads(top->idef.il, md, !DOMAINDECOMP(cr), vsite);
+        split_vsites_over_threads(top->idef.il, top->idef.iparams,
+                                  md, !DOMAINDECOMP(cr), vsite);
     }
 }
-- 
2.22.0