Merge branch 'release-4-5-patches' into release-4-6

[alexxy/gromacs.git] / src / mdlib / sim_util.c

diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c
index 590aaaf5672ac4d7507a00e1ddffb8d57784ecdf..f5cbedd5f9380961b0fae4b40dc588d30b5ca4cc 100644 (file)
--- a/src/mdlib/sim_util.c
+++ b/src/mdlib/sim_util.c
@@ -59,7 +59,7 @@
 #include "pbc.h"
 #include "chargegroup.h"
 #include "vec.h"
-#include "time.h"
+#include <time.h>
 #include "nrnb.h"
 #include "mshift.h"
 #include "mdrun.h"
@@ -80,7 +80,7 @@
 #include "trnio.h"
 #include "xtcio.h"
 #include "copyrite.h"
-
+#include "pull_rotation.h"
 #include "mpelogging.h"
 #include "domdec.h"
 #include "partdec.h"
@@ -90,10 +90,11 @@
 #ifdef GMX_LIB_MPI
 #include <mpi.h>
 #endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
 #include "tmpi.h"
 #endif
 
+#include "adress.h"
 #include "qmmm.h"
 
 #if 0
@@ -139,7 +140,7 @@ void print_time(FILE *out,gmx_runtime_t *runtime,gmx_large_int_t step,
     double dt;
     char buf[48];
     
-#ifndef GMX_THREADS
+#ifndef GMX_THREAD_MPI
     if (!PAR(cr))
 #endif
     {
@@ -176,7 +177,7 @@ void print_time(FILE *out,gmx_runtime_t *runtime,gmx_large_int_t step,
                     ir->delta_t/1000*24*60*60/runtime->time_per_step);
         }
     }
-#ifndef GMX_THREADS
+#ifndef GMX_THREAD_MPI
     if (PAR(cr))
     {
         fprintf(out,"\n");
@@ -427,6 +428,7 @@ void do_force(FILE *fplog,t_commrec *cr,
     double mu[2*DIM]; 
     gmx_bool   bSepDVDL,bStateChanged,bNS,bFillGrid,bCalcCGCM,bBS;
     gmx_bool   bDoLongRange,bDoForces,bSepLRF;
+    gmx_bool   bDoAdressWF;
     matrix boxs;
     real   e,v,dvdl;
     t_pbc  pbc;
@@ -467,6 +469,8 @@ void do_force(FILE *fplog,t_commrec *cr,
     bDoLongRange  = (fr->bTwinRange && bNS && (flags & GMX_FORCE_DOLR));
     bDoForces     = (flags & GMX_FORCE_FORCES);
     bSepLRF       = (bDoLongRange && bDoForces && (flags & GMX_FORCE_SEPLRF));
+    /* should probably move this to the forcerec since it doesn't change */
+    bDoAdressWF   = ((fr->adress_type!=eAdressOff));
 
     if (bStateChanged)
     {
@@ -557,6 +561,33 @@ void do_force(FILE *fplog,t_commrec *cr,
     }
     if (bStateChanged)
     {
+
+        /* update adress weight beforehand */
+        if(bDoAdressWF)
+        {
+            /* need pbc for adress weight calculation with pbc_dx */
+            set_pbc(&pbc,inputrec->ePBC,box);
+            if(fr->adress_site == eAdressSITEcog)
+            {
+                update_adress_weights_cog(top->idef.iparams,top->idef.il,x,fr,mdatoms,
+                                          inputrec->ePBC==epbcNONE ? NULL : &pbc);
+            }
+            else if (fr->adress_site == eAdressSITEcom)
+            {
+                update_adress_weights_com(fplog,cg0,cg1,&(top->cgs),x,fr,mdatoms,
+                                          inputrec->ePBC==epbcNONE ? NULL : &pbc);
+            }
+            else if (fr->adress_site == eAdressSITEatomatom){
+                update_adress_weights_atom_per_atom(cg0,cg1,&(top->cgs),x,fr,mdatoms,
+                                          inputrec->ePBC==epbcNONE ? NULL : &pbc);
+            }
+            else
+            {
+                update_adress_weights_atom(cg0,cg1,&(top->cgs),x,fr,mdatoms,
+                                           inputrec->ePBC==epbcNONE ? NULL : &pbc);
+            }
+        }
+
         for(i=0; i<2; i++)
         {
             for(j=0;j<DIM;j++)
@@ -630,14 +661,24 @@ void do_force(FILE *fplog,t_commrec *cr,
             dd_force_flop_start(cr->dd,nrnb);
         }
     }
-       
+    
+    if (inputrec->bRot)
+    {
+        /* Enforced rotation has its own cycle counter that starts after the collective
+         * coordinates have been communicated. It is added to ddCyclF to allow
+         * for proper load-balancing */
+        wallcycle_start(wcycle,ewcROT);
+        do_rotation(cr,inputrec,box,x,t,step,wcycle,bNS);
+        wallcycle_stop(wcycle,ewcROT);
+    }
+
     /* Start the force cycle counter.
      * This counter is stopped in do_forcelow_level.
      * No parallel communication should occur while this counter is running,
      * since that will interfere with the dynamic load balancing.
      */
     wallcycle_start(wcycle,ewcFORCE);
-    
+
     if (bDoForces)
     {
         /* Reset forces for which the virial is calculated separately:
@@ -699,8 +740,11 @@ void do_force(FILE *fplog,t_commrec *cr,
 
     if ((flags & GMX_FORCE_BONDED) && top->idef.il[F_POSRES].nr > 0)
     {
-        /* Position restraints always require full pbc */
-        set_pbc(&pbc,inputrec->ePBC,box);
+        /* Position restraints always require full pbc. Check if we already did it for Adress */
+        if(!(bStateChanged && bDoAdressWF))
+        {
+            set_pbc(&pbc,inputrec->ePBC,box);
+        }
         v = posres(top->idef.il[F_POSRES].nr,top->idef.il[F_POSRES].iatoms,
                    top->idef.iparams_posres,
                    (const rvec*)x,fr->f_novirsum,fr->vir_diag_posres,
@@ -755,6 +799,13 @@ void do_force(FILE *fplog,t_commrec *cr,
                       start,homenr,mdatoms->chargeA,x,fr->f_novirsum,
                       inputrec->ex,inputrec->et,t);
         }
+
+        if (bDoAdressWF && fr->adress_icor == eAdressICThermoForce)
+        {
+            /* Compute thermodynamic force in hybrid AdResS region */
+            adress_thermo_force(start,homenr,&(top->cgs),x,fr->f_novirsum,fr,mdatoms,
+                                inputrec->ePBC==epbcNONE ? NULL : &pbc);
+        }
         
         /* Communicate the forces */
         if (PAR(cr))
@@ -822,6 +873,7 @@ void do_force(FILE *fplog,t_commrec *cr,
         }
     }
 
+    enerd->term[F_COM_PULL] = 0;
     if (inputrec->ePull == epullUMBRELLA || inputrec->ePull == epullCONST_F)
     {
         /* Calculate the center of mass forces, this requires communication,
@@ -831,7 +883,7 @@ void do_force(FILE *fplog,t_commrec *cr,
          */
         set_pbc(&pbc,inputrec->ePBC,box);
         dvdl = 0; 
-        enerd->term[F_COM_PULL] =
+        enerd->term[F_COM_PULL] +=
             pull_potential(inputrec->ePull,inputrec->pull,mdatoms,&pbc,
                            cr,t,lambda,x,f,vir_force,&dvdl);
         if (bSepDVDL)
@@ -840,6 +892,14 @@ void do_force(FILE *fplog,t_commrec *cr,
         }
         enerd->dvdl_lin += dvdl;
     }
+    
+    /* Add the forces from enforced rotation potentials (if any) */
+    if (inputrec->bRot)
+    {
+        wallcycle_start(wcycle,ewcROTadd);
+        enerd->term[F_COM_PULL] += add_rot_forces(inputrec->rot, f, cr,step,t);
+        wallcycle_stop(wcycle,ewcROTadd);
+    }
 
     if (PAR(cr) && !(cr->duty & DUTY_PME))
     {
@@ -1526,6 +1586,11 @@ void init_md(FILE *fplog,
                               mtop,ir, (*outf)->fp_dhdl);
     }
     
+    if (ir->bAdress)
+    {
+      please_cite(fplog,"Fritsch12");
+      please_cite(fplog,"Junghans10");
+    }
     /* Initiate variables */  
     clear_mat(force_vir);
     clear_mat(shake_vir);