fixed continuation runs with essential dynamics

diff --git a/include/edsam.h b/include/edsam.h
index 1b8be2fc71a8fddc6c2ef12af06c64d7165551bd..c2092c76f60f0094eba43ce79295337ae93c3ebc 100644 (file)
--- a/include/edsam.h
+++ b/include/edsam.h
@@ -44,7 +44,7 @@
 extern "C" {
 #endif
 
-extern void do_edsam(t_inputrec *ir,int step,t_mdatoms *md,
+extern void do_edsam(t_inputrec *ir,gmx_large_int_t step,t_mdatoms *md,
                      t_commrec *cr,rvec xs[],rvec v[],matrix box,gmx_edsam_t ed);
 /* Essential dynamics constraints, called from constrain() */
 

diff --git a/src/mdlib/edsam.c b/src/mdlib/edsam.c
index c50d257509a91d0ee661219f0beb83f0cafac9c3..38783fba33d783fc4930bee21583c2913b67ae9b 100644 (file)
--- a/src/mdlib/edsam.c
+++ b/src/mdlib/edsam.c
@@ -289,8 +289,8 @@ static void project_to_eigvectors(rvec       *x,    /* The positions to project
 static void project(rvec      *x,     /* positions to project */ 
                     t_edpar   *edi)   /* edi data set */
 {
-    /* It is not more work to subtract the average position in every s
-     * ubroutine again, because these routines are rarely used simultanely */
+    /* It is not more work to subtract the average position in every
+     * subroutine again, because these routines are rarely used simultanely */
     project_to_eigvectors(x, &edi->vecs.mon   , edi);
     project_to_eigvectors(x, &edi->vecs.linfix, edi);
     project_to_eigvectors(x, &edi->vecs.linacc, edi);
@@ -1862,7 +1862,7 @@ static void do_radcon(rvec *xcoll, t_edpar *edi, t_commrec *cr)
 }
 
 
-static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
+static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step, t_commrec *cr)
 {
     int i;
 
@@ -2236,7 +2236,7 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
 
 
 void do_edsam(t_inputrec  *ir,
-              int         step,
+              gmx_large_int_t step,
               t_mdatoms   *md,
               t_commrec   *cr,
               rvec        xs[],   /* The local current positions on this processor */
@@ -2334,7 +2334,7 @@ void do_edsam(t_inputrec  *ir,
             }
 
             /* update radsam references, when required */
-            if (do_per_step(step,edi->maxedsteps) && step > edi->presteps)
+            if (do_per_step(step,edi->maxedsteps) && step >= edi->presteps)
             {
                 project(buf->xcoll, edi);
                 rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
@@ -2343,7 +2343,7 @@ void do_edsam(t_inputrec  *ir,
             }
 
             /* update radacc references, when required */
-            if (do_per_step(step,iupdate) && step > edi->presteps)
+            if (do_per_step(step,iupdate) && step >= edi->presteps)
             {
                 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
                 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
@@ -2356,8 +2356,12 @@ void do_edsam(t_inputrec  *ir,
             }
 
             /* apply the constraints */
-            if (step > edi->presteps && ed_constraints(ed->eEDtype, edi))
-                ed_apply_constraints(buf->xcoll, edi, step, cr);
+            if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
+            {
+                /* ED constraints should be applied already in the first MD step
+                 * (which is step 0), therfore we pass step+1 to the routine */
+                ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step, cr);
+            }
 
             /* write to edo, when required */
             if (do_per_step(step,edi->outfrq))