Avoid using function calls in OpenMP directives
[alexxy/gromacs.git] / src / gromacs / mdlib / minimize.c
index 26122e26d4c6f0ce4ed295b0ab80d246e1b43747..69008f53fabcb3388951c01d7554d479c67ccdca 100644 (file)
@@ -42,9 +42,9 @@
 #include <time.h>
 #include <math.h>
 #include "sysstuff.h"
-#include "string2.h"
+#include "gromacs/utility/cstringutil.h"
 #include "network.h"
-#include "smalloc.h"
+#include "gromacs/utility/smalloc.h"
 #include "nrnb.h"
 #include "main.h"
 #include "force.h"
@@ -78,6 +78,7 @@
 #include "gromacs/linearalgebra/sparsematrix.h"
 #include "gromacs/timing/wallcycle.h"
 #include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/imd/imd.h"
 
 typedef struct {
     t_state  s;
@@ -308,7 +309,9 @@ void init_em(FILE *fplog, const char *title,
              t_graph **graph, t_mdatoms *mdatoms, gmx_global_stat_t *gstat,
              gmx_vsite_t *vsite, gmx_constr_t constr,
              int nfile, const t_filenm fnm[],
-             gmx_mdoutf_t *outf, t_mdebin **mdebin)
+             gmx_mdoutf_t *outf, t_mdebin **mdebin,
+             int imdport, unsigned long gmx_unused Flags,
+             gmx_wallcycle_t wcycle)
 {
     int  i;
     real dvdl_constr;
@@ -325,6 +328,10 @@ void init_em(FILE *fplog, const char *title,
 
     init_nrnb(nrnb);
 
+    /* Interactive molecular dynamics */
+    init_IMD(ir, cr, top_global, fplog, 1, state_global->x,
+             nfile, fnm, NULL, imdport, Flags);
+
     if (DOMAINDECOMP(cr))
     {
         *top = dd_init_local_top(top_global);
@@ -409,7 +416,7 @@ void init_em(FILE *fplog, const char *title,
             /* Constrain the starting coordinates */
             dvdl_constr = 0;
             constrain(PAR(cr) ? NULL : fplog, TRUE, TRUE, constr, &(*top)->idef,
-                      ir, NULL, cr, -1, 0, mdatoms,
+                      ir, NULL, cr, -1, 0, 1.0, mdatoms,
                       ems->s.x, ems->s.x, NULL, fr->bMolPBC, ems->s.box,
                       ems->s.lambda[efptFEP], &dvdl_constr,
                       NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
@@ -421,7 +428,7 @@ void init_em(FILE *fplog, const char *title,
         *gstat = global_stat_init(ir);
     }
 
-    *outf = init_mdoutf(fplog, nfile, fnm, 0, cr, ir, top_global, NULL);
+    *outf = init_mdoutf(fplog, nfile, fnm, 0, cr, ir, top_global, NULL, wcycle);
 
     snew(*enerd, 1);
     init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
@@ -479,9 +486,17 @@ static void write_em_traj(FILE *fplog, t_commrec *cr,
                           em_state_t *state,
                           t_state *state_global, rvec *f_global)
 {
-    int mdof_flags;
+    int      mdof_flags;
+    gmx_bool bIMDout = FALSE;
+
 
-    if ((bX || bF || confout != NULL) && !DOMAINDECOMP(cr))
+    /* Shall we do IMD output? */
+    if (ir->bIMD)
+    {
+        bIMDout = do_per_step(step, IMD_get_step(ir->imd->setup));
+    }
+
+    if ((bX || bF || bIMDout || confout != NULL) && !DOMAINDECOMP(cr))
     {
         copy_em_coords(state, state_global);
         f_global = state->f;
@@ -496,6 +511,13 @@ static void write_em_traj(FILE *fplog, t_commrec *cr,
     {
         mdof_flags |= MDOF_F;
     }
+
+    /* If we want IMD output, set appropriate MDOF flag */
+    if (ir->bIMD)
+    {
+        mdof_flags |= MDOF_IMD;
+    }
+
     mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
                                      top_global, step, (double)step,
                                      &state->s, state_global, state->f, f_global);
@@ -528,6 +550,7 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
     int      start, end;
     rvec    *x1, *x2;
     real     dvdl_constr;
+    int      nthreads gmx_unused;
 
     s1 = &ems1->s;
     s2 = &ems2->s;
@@ -565,7 +588,8 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
     x1 = s1->x;
     x2 = s2->x;
 
-#pragma omp parallel num_threads(gmx_omp_nthreads_get(emntUpdate))
+    nthreads = gmx_omp_nthreads_get(emntUpdate);
+#pragma omp parallel num_threads(nthreads)
     {
         int gf, i, m;
 
@@ -627,7 +651,7 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
         wallcycle_start(wcycle, ewcCONSTR);
         dvdl_constr = 0;
         constrain(NULL, TRUE, TRUE, constr, &top->idef,
-                  ir, NULL, cr, count, 0, md,
+                  ir, NULL, cr, count, 0, 1.0, md,
                   s1->x, s2->x, NULL, bMolPBC, s2->box,
                   s2->lambda[efptBONDED], &dvdl_constr,
                   NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
@@ -764,7 +788,7 @@ static void evaluate_energy(FILE *fplog, t_commrec *cr,
         wallcycle_start(wcycle, ewcCONSTR);
         dvdl_constr = 0;
         constrain(NULL, FALSE, FALSE, constr, &top->idef,
-                  inputrec, NULL, cr, count, 0, mdatoms,
+                  inputrec, NULL, cr, count, 0, 1.0, mdatoms,
                   ems->s.x, ems->f, ems->f, fr->bMolPBC, ems->s.box,
                   ems->s.lambda[efptBONDED], &dvdl_constr,
                   NULL, &shake_vir, nrnb, econqForceDispl, FALSE, 0, 0);
@@ -939,6 +963,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
              gmx_membed_t gmx_unused membed,
              real gmx_unused cpt_period, real gmx_unused max_hours,
              const char gmx_unused *deviceOptions,
+             int imdport,
              unsigned long gmx_unused Flags,
              gmx_walltime_accounting_t walltime_accounting)
 {
@@ -978,7 +1003,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
     init_em(fplog, CG, cr, inputrec,
             state_global, top_global, s_min, &top, &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, &mdebin);
+            nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
 
     /* Print to log file */
     print_em_start(fplog, cr, walltime_accounting, wcycle, CG);
@@ -1445,6 +1470,10 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
             do_log = do_per_step(step, inputrec->nstlog);
             do_ene = do_per_step(step, inputrec->nstenergy);
+
+            /* Prepare IMD energy record, if bIMD is TRUE. */
+            IMD_fill_energy_record(inputrec->bIMD, inputrec->imd, enerd, step, TRUE);
+
             if (do_log)
             {
                 print_ebin_header(fplog, step, step, s_min->s.lambda[efptFEP]);
@@ -1454,6 +1483,12 @@ double do_cg(FILE *fplog, t_commrec *cr,
                        TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
         }
 
+        /* Send energies and positions to the IMD client if bIMD is TRUE. */
+        if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+        {
+            IMD_send_positions(inputrec->imd);
+        }
+
         /* Stop when the maximum force lies below tolerance.
          * If we have reached machine precision, converged is already set to true.
          */
@@ -1461,6 +1496,9 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
     } /* End of the loop */
 
+    /* IMD cleanup, if bIMD is TRUE. */
+    IMD_finalize(inputrec->bIMD, inputrec->imd);
+
     if (converged)
     {
         step--; /* we never took that last step in this case */
@@ -1553,6 +1591,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
                 gmx_membed_t gmx_unused membed,
                 real gmx_unused cpt_period, real gmx_unused max_hours,
                 const char gmx_unused *deviceOptions,
+                int imdport,
                 unsigned long gmx_unused Flags,
                 gmx_walltime_accounting_t walltime_accounting)
 {
@@ -1636,7 +1675,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     init_em(fplog, LBFGS, cr, inputrec,
             state, top_global, &ems, &top, &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, &mdebin);
+            nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
     /* Do_lbfgs is not completely updated like do_steep and do_cg,
      * so we free some memory again.
      */
@@ -1782,6 +1821,11 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
             mdof_flags |= MDOF_F;
         }
 
+        if (inputrec->bIMD)
+        {
+            mdof_flags |= MDOF_IMD;
+        }
+
         mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
                                          top_global, step, (real)step, state, state, f, f);
 
@@ -2243,6 +2287,12 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
                        TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
         }
 
+        /* Send x and E to IMD client, if bIMD is TRUE. */
+        if (do_IMD(inputrec->bIMD, step, cr, TRUE, state->box, state->x, inputrec, 0, wcycle) && MASTER(cr))
+        {
+            IMD_send_positions(inputrec->imd);
+        }
+
         /* Stop when the maximum force lies below tolerance.
          * If we have reached machine precision, converged is already set to true.
          */
@@ -2251,6 +2301,9 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
 
     } /* End of the loop */
 
+    /* IMD cleanup, if bIMD is TRUE. */
+    IMD_finalize(inputrec->bIMD, inputrec->imd);
+
     if (converged)
     {
         step--; /* we never took that last step in this case */
@@ -2335,6 +2388,7 @@ double do_steep(FILE *fplog, t_commrec *cr,
                 gmx_membed_t gmx_unused membed,
                 real gmx_unused cpt_period, real gmx_unused max_hours,
                 const char  gmx_unused *deviceOptions,
+                int imdport,
                 unsigned long gmx_unused Flags,
                 gmx_walltime_accounting_t walltime_accounting)
 {
@@ -2366,7 +2420,7 @@ double do_steep(FILE *fplog, t_commrec *cr,
     init_em(fplog, SD, cr, inputrec,
             state_global, top_global, s_try, &top, &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, &mdebin);
+            nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
 
     /* Print to log file  */
     print_em_start(fplog, cr, walltime_accounting, wcycle, SD);
@@ -2443,6 +2497,10 @@ double do_steep(FILE *fplog, t_commrec *cr,
                 upd_mdebin(mdebin, FALSE, FALSE, (double)count,
                            mdatoms->tmass, enerd, &s_try->s, inputrec->fepvals, inputrec->expandedvals,
                            s_try->s.box, NULL, NULL, vir, pres, NULL, mu_tot, constr);
+
+                /* Prepare IMD energy record, if bIMD is TRUE. */
+                IMD_fill_energy_record(inputrec->bIMD, inputrec->imd, enerd, count, TRUE);
+
                 print_ebin(mdoutf_get_fp_ene(outf), TRUE,
                            do_per_step(steps_accepted, inputrec->nstdisreout),
                            do_per_step(steps_accepted, inputrec->nstorireout),
@@ -2512,10 +2570,19 @@ double do_steep(FILE *fplog, t_commrec *cr,
             bAbort = TRUE;
         }
 
+        /* Send IMD energies and positions, if bIMD is TRUE. */
+        if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+        {
+            IMD_send_positions(inputrec->imd);
+        }
+
         count++;
     } /* End of the loop  */
 
-    /* Print some shit...  */
+    /* IMD cleanup, if bIMD is TRUE. */
+    IMD_finalize(inputrec->bIMD, inputrec->imd);
+
+    /* Print some data...  */
     if (MASTER(cr))
     {
         fprintf(stderr, "\nwriting lowest energy coordinates.\n");
@@ -2562,6 +2629,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
              gmx_membed_t gmx_unused membed,
              real gmx_unused cpt_period, real gmx_unused max_hours,
              const char gmx_unused *deviceOptions,
+             int imdport,
              unsigned long gmx_unused Flags,
              gmx_walltime_accounting_t walltime_accounting)
 {
@@ -2604,7 +2672,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
             state_global, top_global, state_work, &top,
             &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, NULL);
+            nfile, fnm, &outf, NULL, imdport, Flags, wcycle);
 
     natoms = top_global->natoms;
     snew(fneg, natoms);