Merge release-4-5-patches into release-4-6
[alexxy/gromacs.git] / src / tools / gmx_membed.c
index e3d4c671ff9bc2557e7c33c71250bba0e2ac80d9..cb73909dc4df5d4c1d4a309a138d340887ae12b9 100644 (file)
 #include <config.h>
 #endif
 
+#ifdef __linux
+#define _GNU_SOURCE
+#include <sched.h>
+#include <sys/syscall.h>
+#endif
 #include <signal.h>
 #include <stdlib.h>
+
 #include "typedefs.h"
 #include "smalloc.h"
 #include "sysstuff.h"
@@ -69,7 +75,6 @@
 #include "disre.h"
 #include "orires.h"
 #include "dihre.h"
-#include "pppm.h"
 #include "pme.h"
 #include "mdatoms.h"
 #include "qmmm.h"
 #ifdef GMX_LIB_MPI
 #include <mpi.h>
 #endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
 #include "tmpi.h"
 #endif
 
+#ifdef GMX_OPENMP
+#include <omp.h>
+#endif
+
 /* afm stuf */
 #include "pull.h"
 
@@ -244,23 +253,6 @@ int get_tpr_version(const char *infile)
        return fver;
 }
 
-void set_inbox(int natom, rvec *x)
-{
-       rvec tmp;
-       int  i;
-
-       tmp[XX]=tmp[YY]=tmp[ZZ]=0.0;
-       for(i=0;i<natom;i++)
-       {
-               if(x[i][XX]<tmp[XX])            tmp[XX]=x[i][XX];
-               if(x[i][YY]<tmp[YY])            tmp[YY]=x[i][YY];
-               if(x[i][ZZ]<tmp[ZZ])            tmp[ZZ]=x[i][ZZ];
-       }
-
-       for(i=0;i<natom;i++)
-                       rvec_inc(x[i],tmp);
-}
-
 int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
 {
        int i,j,nr,mol_id;
@@ -410,28 +402,6 @@ real est_prot_area(pos_ins_t *pos_ins,rvec *r,t_block *ins_at, mem_t *mem_p)
        return area;
 }
 
-void init_lip(matrix box, gmx_mtop_t *mtop, lip_t *lip)
-{
-       int i;
-       real mem_area;
-       int mol1=0;
-
-       mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
-       for(i=0;i<mtop->nmolblock;i++)
-       {
-               if(mtop->molblock[i].type == lip->id)
-               {
-                       lip->nr=mtop->molblock[i].nmol;
-                       lip->natoms=mtop->molblock[i].natoms_mol;
-               }
-       }
-       lip->area=2.0*mem_area/(double)lip->nr;
-
-       for (i=0;i<lip->id;i++)
-               mol1+=mtop->molblock[i].nmol;
-       lip->mol1=mol1;
-}
-
 int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
 {
        int i,j,at,mol,nmol,nmolbox,count;
@@ -936,10 +906,11 @@ void top_update(const char *topfile, char *ins, rmm_t *rm_p, gmx_mtop_t *mtop)
                }
        }
 
-       fclose(fpout);
+       ffclose(fpout);
+       ffclose(fpin);
        /* use ffopen to generate backup of topinout */
        fpout=ffopen(topfile,"w");
-       fclose(fpout);
+       ffclose(fpout);
        rename(TEMP_FILENM,topfile);
 #undef TEMP_FILENM
 }
@@ -1138,13 +1109,13 @@ double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 
 /*    if (DEFORM(*ir))
     {
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
 #endif
         set_deform_reference_box(upd,
                                  deform_init_init_step_tpx,
                                  deform_init_box_tpx);
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
 #endif
     }*/
@@ -2118,7 +2089,7 @@ double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 
         /* Check whether everything is still allright */
         if (((int)gmx_get_stop_condition() > handled_stop_condition)
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
            && MASTER(cr)
 #endif
            )
@@ -2747,7 +2718,7 @@ int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
     gmx_edsam_t ed=NULL;
     t_commrec   *cr_old=cr;
     int        nthreads=1,nthreads_requested=1;
-
+    int         omp_nthreads = 1;
 
        char                    *ins;
        int                     rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
@@ -2787,7 +2758,7 @@ int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
         read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
 
         /* NOW the threads will be started: */
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
 #endif
     }
     /* END OF CAUTION: cr is now reliable */
@@ -3160,7 +3131,32 @@ int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
         gmx_setup_nodecomm(fplog,cr);
     }
 
-    wcycle = wallcycle_init(fplog,resetstep,cr);
+    /* get number of OpenMP/PME threads
+    * env variable should be read only on one node to make sure it is identical everywhere */
+#ifdef GMX_OPENMP
+   if (EEL_PME(inputrec->coulombtype))
+   {
+       if (MASTER(cr))
+       {
+           char *ptr;
+           omp_nthreads = omp_get_max_threads();
+           if ((ptr=getenv("GMX_PME_NTHREADS")) != NULL)
+           {
+               sscanf(ptr,"%d",&omp_nthreads);
+           }
+           if (fplog!=NULL)
+           {
+               fprintf(fplog,"Using %d threads for PME\n",omp_nthreads);
+           }
+       }
+       if (PAR(cr))
+       {
+           gmx_bcast_sim(sizeof(omp_nthreads),&omp_nthreads,cr);
+       }
+   }
+#endif
+
+    wcycle = wallcycle_init(fplog,resetstep,cr, omp_nthreads);
     if (PAR(cr))
     {
         /* Master synchronizes its value of reset_counters with all nodes
@@ -3203,12 +3199,13 @@ int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
         fr = mk_forcerec();
         init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
                       opt2fn("-table",nfile,fnm),
+                      opt2fn("-tabletf",nfile,fnm),
                       opt2fn("-tablep",nfile,fnm),
                       opt2fn("-tableb",nfile,fnm),FALSE,pforce);
 
         /* version for PCA_NOT_READ_NODE (see md.c) */
         /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
-          "nofile","nofile","nofile",FALSE,pforce);
+          "nofile","nofile","nofile","nofile",FALSE,pforce);
           */
         fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
 
@@ -3250,22 +3247,6 @@ int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
             }
         }
 
-        /* Initiate PPPM if necessary */
-        if (fr->eeltype == eelPPPM)
-        {
-            if (mdatoms->nChargePerturbed)
-            {
-                gmx_fatal(FARGS,"Free energy with %s is not implemented",
-                          eel_names[fr->eeltype]);
-            }
-            status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
-                                   getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
-            if (status != 0)
-            {
-                gmx_fatal(FARGS,"Error %d initializing PPPM",status);
-            }
-        }
-
         if (EEL_PME(fr->eeltype))
         {
             ewaldcoeff = fr->ewaldcoeff;
@@ -3300,11 +3281,37 @@ int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
             /* The PME only nodes need to know nChargePerturbed */
             gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
         }
+
+
+        /*set CPU affinity*/
+#ifdef GMX_OPENMP
+#ifdef __linux
+#ifdef GMX_LIB_MPI
+        {
+            int core;
+            MPI_Comm comm_intra; /*intra communicator (but different to nc.comm_intra includes PME nodes)*/
+            MPI_Comm_split(MPI_COMM_WORLD,gmx_hostname_num(),gmx_node_rank(),&comm_intra);
+            int local_omp_nthreads = (cr->duty & DUTY_PME) ? omp_nthreads : 1; /*threads on this node*/
+            MPI_Scan(&local_omp_nthreads,&core, 1, MPI_INT, MPI_SUM, comm_intra);
+            core-=local_omp_nthreads; /*make exclusive scan*/
+    #pragma omp parallel firstprivate(core) num_threads(local_omp_nthreads)
+            {
+                cpu_set_t mask;
+                CPU_ZERO(&mask);
+                core+=omp_get_thread_num();
+                CPU_SET(core,&mask);
+                sched_setaffinity((pid_t) syscall (SYS_gettid),sizeof(cpu_set_t),&mask);
+            }
+        }
+#endif /*GMX_MPI*/
+#endif /*__linux*/
+#endif /*GMX_OPENMP*/
+
         if (cr->duty & DUTY_PME)
         {
             status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
                                   mtop ? mtop->natoms : 0,nChargePerturbed,
-                                  (Flags & MD_REPRODUCIBLE));
+                                  (Flags & MD_REPRODUCIBLE),omp_nthreads);
             if (status != 0)
             {
                 gmx_fatal(FARGS,"Error %d initializing PME",status);
@@ -3675,7 +3682,7 @@ int gmx_membed(int argc,char *argv[])
        dd_node_order = nenum(ddno_opt);
        cr->npmenodes = npme;
 
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
        /* now determine the number of threads automatically. The threads are
    only started at mdrunner_threads, though. */
        if (nthreads<1)
@@ -3702,7 +3709,7 @@ int gmx_membed(int argc,char *argv[])
                gmx_fatal(FARGS,"Need at least two replicas for replica exchange (option -multi)");
 
        if (nmultisim > 1) {
-#ifndef GMX_THREADS
+#ifndef GMX_THREAD_MPI
                 gmx_bool bParFn = (multidir == NULL);
                init_multisystem(cr,nmultisim,multidir,NFILE,fnm,TRUE);
 #else