Choose the PBC image of each rotation group atom depending

on the current position of it reference atom. Therefore,
1) in choose_pbc_image() the correctly rotated reference position
has to be used (the addition of the pivot vector was missing)
2) in init_rot_group(), the reference positions have to be
correctly determined for restarts. In these cases, the initial
angle of the reference group is different from zero.

Change-Id: I57f886c507317385619e42bf387e407479b0c08d

diff --git a/src/mdlib/pull_rotation.c b/src/mdlib/pull_rotation.c
index a490ccfff80532494d321841f271df358040c757..5b600b908efe27470d8e45b83ec33700b6602435 100644 (file)
--- a/src/mdlib/pull_rotation.c
+++ b/src/mdlib/pull_rotation.c
@@ -3386,10 +3386,10 @@ static inline void copy_correct_pbc_image(
 
 static void init_rot_group(FILE *fplog, t_commrec *cr, int g, t_rotgrp *rotg,
                            rvec *x, gmx_mtop_t *mtop, gmx_bool bVerbose, FILE *out_slabs, matrix box,
-                           gmx_bool bOutputCenters)
+                           t_inputrec *ir, gmx_bool bOutputCenters)
 {
     int                   i, ii;
-    rvec                  coord, *xdum;
+    rvec                  coord, xref, *xdum;
     gmx_bool              bFlex, bColl;
     t_atom               *atom;
     gmx_enfrotgrp_t       erg; /* Pointer to enforced rotation group data */
@@ -3397,6 +3397,7 @@ static void init_rot_group(FILE *fplog, t_commrec *cr, int g, t_rotgrp *rotg,
     gmx_mtop_atomlookup_t alook = NULL;
     real                  mass, totalmass;
     real                  start = 0.0;
+    double                t_start;
 
 
     /* Do we have a flexible axis? */
@@ -3412,24 +3413,7 @@ static void init_rot_group(FILE *fplog, t_commrec *cr, int g, t_rotgrp *rotg,
         snew(erg->xc, rotg->nat);
         snew(erg->xc_shifts, rotg->nat);
         snew(erg->xc_eshifts, rotg->nat);
-
-        /* Save the original (whole) set of positions such that later the
-         * molecule can always be made whole again */
         snew(erg->xc_old, rotg->nat);
-        if (MASTER(cr))
-        {
-            for (i = 0; i < rotg->nat; i++)
-            {
-                ii = rotg->ind[i];
-                copy_correct_pbc_image(x[ii], erg->xc_old[i], rotg->x_ref[i], box, 3);
-            }
-        }
-#ifdef GMX_MPI
-        if (PAR(cr))
-        {
-            gmx_bcast(rotg->nat*sizeof(erg->xc_old[0]), erg->xc_old, cr);
-        }
-#endif
 
         if (rotg->eFittype == erotgFitNORM)
         {
@@ -3545,6 +3529,41 @@ static void init_rot_group(FILE *fplog, t_commrec *cr, int g, t_rotgrp *rotg,
         }
 #endif
     }
+    
+    if (bColl)
+    {
+        /* Save the original (whole) set of positions in xc_old such that at later 
+         * steps the rotation group can always be made whole again. If the simulation is
+         * restarted, we compute the starting reference positions (given the time)
+         * and assume that the correct PBC image of each position is the one nearest
+         * to the current reference */
+        if (MASTER(cr))
+        {
+            /* Calculate the rotation matrix for this angle: */
+            t_start       = ir->init_t + ir->init_step*ir->delta_t;
+            erg->degangle = rotg->rate * t_start;
+            calc_rotmat(rotg->vec, erg->degangle, erg->rotmat);
+
+            for (i = 0; i < rotg->nat; i++)
+            {
+                ii = rotg->ind[i];
+
+                /* Subtract pivot, rotate, and add pivot again. This will yield the 
+                 * reference position for time t */
+                rvec_sub(rotg->x_ref[i], erg->xc_ref_center, coord);
+                mvmul(erg->rotmat, coord, xref);
+                rvec_inc(xref, erg->xc_ref_center);
+
+                copy_correct_pbc_image(x[ii], erg->xc_old[i], xref, box, 3);
+            }
+        }
+#ifdef GMX_MPI
+        if (PAR(cr))
+        {
+            gmx_bcast(rotg->nat*sizeof(erg->xc_old[0]), erg->xc_old, cr);
+        }
+#endif
+    }
 
     if ( (rotg->eType != erotgFLEX) && (rotg->eType != erotgFLEX2) )
     {
@@ -3749,7 +3768,7 @@ extern void init_rot(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[
                 erg->nat_loc = rotg->nat;
                 erg->ind_loc = rotg->ind;
             }
-            init_rot_group(fplog, cr, g, rotg, x_pbc, mtop, bVerbose, er->out_slabs, box,
+            init_rot_group(fplog, cr, g, rotg, x_pbc, mtop, bVerbose, er->out_slabs, box, ir,
                            !(er->Flags & MD_APPENDFILES) ); /* Do not output the reference centers
                                                              * again if we are appending */
         }
@@ -3866,11 +3885,13 @@ static void choose_pbc_image(rvec x[], t_rotgrp *rotg, matrix box, int npbcdim)
         /* Index of a rotation group atom  */
         ii = erg->ind_loc[i];
 
-        /* Get the reference position. The pivot was already
+        /* Get the correctly rotated reference position. The pivot was already
          * subtracted in init_rot_group() from the reference positions. Also,
          * the reference positions have already been rotated in
-         * rotate_local_reference() */
+         * rotate_local_reference(). For the current reference position we thus
+         * only need to add the pivot again. */
         copy_rvec(erg->xr_loc[i], xref);
+        rvec_inc(xref, erg->xc_ref_center);
 
         copy_correct_pbc_image(x[ii], erg->x_loc_pbc[i], xref, box, npbcdim);
     }