src/gromacs/mdrun/minimize.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  *
  39  * \brief This file defines integrators for energy minimization
  40  *
  41  * \author Berk Hess <hess@kth.se>
  42  * \author Erik Lindahl <erik@kth.se>
  43  * \ingroup module_mdrun
  44  */
  45 #include "gmxpre.h"
  46
  47 #include "config.h"
  48
  49 #include <cmath>
  50 #include <cstring>
  51 #include <ctime>
  52
  53 #include <algorithm>
  54 #include <vector>
  55
  56 #include "gromacs/commandline/filenm.h"
  57 #include "gromacs/domdec/collect.h"
  58 #include "gromacs/domdec/dlbtiming.h"
  59 #include "gromacs/domdec/domdec.h"
  60 #include "gromacs/domdec/domdec_struct.h"
  61 #include "gromacs/domdec/partition.h"
  62 #include "gromacs/ewald/pme.h"
  63 #include "gromacs/fileio/confio.h"
  64 #include "gromacs/fileio/mtxio.h"
  65 #include "gromacs/gmxlib/network.h"
  66 #include "gromacs/gmxlib/nrnb.h"
  67 #include "gromacs/imd/imd.h"
  68 #include "gromacs/linearalgebra/sparsematrix.h"
  69 #include "gromacs/listed_forces/manage_threading.h"
  70 #include "gromacs/math/functions.h"
  71 #include "gromacs/math/vec.h"
  72 #include "gromacs/mdlib/constr.h"
  73 #include "gromacs/mdlib/dispersioncorrection.h"
  74 #include "gromacs/mdlib/ebin.h"
  75 #include "gromacs/mdlib/energyoutput.h"
  76 #include "gromacs/mdlib/force.h"
  77 #include "gromacs/mdlib/forcerec.h"
  78 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  79 #include "gromacs/mdlib/md_support.h"
  80 #include "gromacs/mdlib/mdatoms.h"
  81 #include "gromacs/mdlib/mdsetup.h"
  82 #include "gromacs/mdlib/ns.h"
  83 #include "gromacs/mdlib/shellfc.h"
  84 #include "gromacs/mdlib/stat.h"
  85 #include "gromacs/mdlib/tgroup.h"
  86 #include "gromacs/mdlib/trajectory_writing.h"
  87 #include "gromacs/mdlib/update.h"
  88 #include "gromacs/mdlib/vsite.h"
  89 #include "gromacs/mdrunutility/printtime.h"
  90 #include "gromacs/mdtypes/commrec.h"
  91 #include "gromacs/mdtypes/inputrec.h"
  92 #include "gromacs/mdtypes/md_enums.h"
  93 #include "gromacs/mdtypes/mdrunoptions.h"
  94 #include "gromacs/mdtypes/state.h"
  95 #include "gromacs/pbcutil/mshift.h"
  96 #include "gromacs/pbcutil/pbc.h"
  97 #include "gromacs/timing/wallcycle.h"
  98 #include "gromacs/timing/walltime_accounting.h"
  99 #include "gromacs/topology/mtop_util.h"
 100 #include "gromacs/topology/topology.h"
 101 #include "gromacs/utility/cstringutil.h"
 102 #include "gromacs/utility/exceptions.h"
 103 #include "gromacs/utility/fatalerror.h"
 104 #include "gromacs/utility/logger.h"
 105 #include "gromacs/utility/smalloc.h"
 106
 107 #include "integrator.h"
 108
 109 //! Utility structure for manipulating states during EM
 110 typedef struct {
 111     //! Copy of the global state
 112     t_state                 s;
 113     //! Force array
 114     PaddedVector<gmx::RVec> f;
 115     //! Potential energy
 116     real                    epot;
 117     //! Norm of the force
 118     real                    fnorm;
 119     //! Maximum force
 120     real                    fmax;
 121     //! Direction
 122     int                     a_fmax;
 123 } em_state_t;
 124
 125 //! Print the EM starting conditions
 126 static void print_em_start(FILE                     *fplog,
 127                            const t_commrec          *cr,
 128                            gmx_walltime_accounting_t walltime_accounting,
 129                            gmx_wallcycle_t           wcycle,
 130                            const char               *name)
 131 {
 132     walltime_accounting_start_time(walltime_accounting);
 133     wallcycle_start(wcycle, ewcRUN);
 134     print_start(fplog, cr, walltime_accounting, name);
 135 }
 136
 137 //! Stop counting time for EM
 138 static void em_time_end(gmx_walltime_accounting_t walltime_accounting,
 139                         gmx_wallcycle_t           wcycle)
 140 {
 141     wallcycle_stop(wcycle, ewcRUN);
 142
 143     walltime_accounting_end_time(walltime_accounting);
 144 }
 145
 146 //! Printing a log file and console header
 147 static void sp_header(FILE *out, const char *minimizer, real ftol, int nsteps)
 148 {
 149     fprintf(out, "\n");
 150     fprintf(out, "%s:\n", minimizer);
 151     fprintf(out, "   Tolerance (Fmax)   = %12.5e\n", ftol);
 152     fprintf(out, "   Number of steps    = %12d\n", nsteps);
 153 }
 154
 155 //! Print warning message
 156 static void warn_step(FILE     *fp,
 157                       real      ftol,
 158                       real      fmax,
 159                       gmx_bool  bLastStep,
 160                       gmx_bool  bConstrain)
 161 {
 162     constexpr bool realIsDouble = GMX_DOUBLE;
 163     char           buffer[2048];
 164
 165     if (!std::isfinite(fmax))
 166     {
 167         sprintf(buffer,
 168                 "\nEnergy minimization has stopped because the force "
 169                 "on at least one atom is not finite. This usually means "
 170                 "atoms are overlapping. Modify the input coordinates to "
 171                 "remove atom overlap or use soft-core potentials with "
 172                 "the free energy code to avoid infinite forces.\n%s",
 173                 !realIsDouble ?
 174                 "You could also be lucky that switching to double precision "
 175                 "is sufficient to obtain finite forces.\n" :
 176                 "");
 177     }
 178     else if (bLastStep)
 179     {
 180         sprintf(buffer,
 181                 "\nEnergy minimization reached the maximum number "
 182                 "of steps before the forces reached the requested "
 183                 "precision Fmax < %g.\n", ftol);
 184     }
 185     else
 186     {
 187         sprintf(buffer,
 188                 "\nEnergy minimization has stopped, but the forces have "
 189                 "not converged to the requested precision Fmax < %g (which "
 190                 "may not be possible for your system). It stopped "
 191                 "because the algorithm tried to make a new step whose size "
 192                 "was too small, or there was no change in the energy since "
 193                 "last step. Either way, we regard the minimization as "
 194                 "converged to within the available machine precision, "
 195                 "given your starting configuration and EM parameters.\n%s%s",
 196                 ftol,
 197                 !realIsDouble ?
 198                 "\nDouble precision normally gives you higher accuracy, but "
 199                 "this is often not needed for preparing to run molecular "
 200                 "dynamics.\n" :
 201                 "",
 202                 bConstrain ?
 203                 "You might need to increase your constraint accuracy, or turn\n"
 204                 "off constraints altogether (set constraints = none in mdp file)\n" :
 205                 "");
 206     }
 207
 208     fputs(wrap_lines(buffer, 78, 0, FALSE), stderr);
 209     fputs(wrap_lines(buffer, 78, 0, FALSE), fp);
 210 }
 211
 212 //! Print message about convergence of the EM
 213 static void print_converged(FILE *fp, const char *alg, real ftol,
 214                             int64_t count, gmx_bool bDone, int64_t nsteps,
 215                             const em_state_t *ems, double sqrtNumAtoms)
 216 {
 217     char buf[STEPSTRSIZE];
 218
 219     if (bDone)
 220     {
 221         fprintf(fp, "\n%s converged to Fmax < %g in %s steps\n",
 222                 alg, ftol, gmx_step_str(count, buf));
 223     }
 224     else if (count < nsteps)
 225     {
 226         fprintf(fp, "\n%s converged to machine precision in %s steps,\n"
 227                 "but did not reach the requested Fmax < %g.\n",
 228                 alg, gmx_step_str(count, buf), ftol);
 229     }
 230     else
 231     {
 232         fprintf(fp, "\n%s did not converge to Fmax < %g in %s steps.\n",
 233                 alg, ftol, gmx_step_str(count, buf));
 234     }
 235
 236 #if GMX_DOUBLE
 237     fprintf(fp, "Potential Energy  = %21.14e\n", ems->epot);
 238     fprintf(fp, "Maximum force     = %21.14e on atom %d\n", ems->fmax, ems->a_fmax + 1);
 239     fprintf(fp, "Norm of force     = %21.14e\n", ems->fnorm/sqrtNumAtoms);
 240 #else
 241     fprintf(fp, "Potential Energy  = %14.7e\n", ems->epot);
 242     fprintf(fp, "Maximum force     = %14.7e on atom %d\n", ems->fmax, ems->a_fmax + 1);
 243     fprintf(fp, "Norm of force     = %14.7e\n", ems->fnorm/sqrtNumAtoms);
 244 #endif
 245 }
 246
 247 //! Compute the norm and max of the force array in parallel
 248 static void get_f_norm_max(const t_commrec *cr,
 249                            t_grpopts *opts, t_mdatoms *mdatoms, const rvec *f,
 250                            real *fnorm, real *fmax, int *a_fmax)
 251 {
 252     double fnorm2, *sum;
 253     real   fmax2, fam;
 254     int    la_max, a_max, start, end, i, m, gf;
 255
 256     /* This routine finds the largest force and returns it.
 257      * On parallel machines the global max is taken.
 258      */
 259     fnorm2 = 0;
 260     fmax2  = 0;
 261     la_max = -1;
 262     start  = 0;
 263     end    = mdatoms->homenr;
 264     if (mdatoms->cFREEZE)
 265     {
 266         for (i = start; i < end; i++)
 267         {
 268             gf  = mdatoms->cFREEZE[i];
 269             fam = 0;
 270             for (m = 0; m < DIM; m++)
 271             {
 272                 if (!opts->nFreeze[gf][m])
 273                 {
 274                     fam += gmx::square(f[i][m]);
 275                 }
 276             }
 277             fnorm2 += fam;
 278             if (fam > fmax2)
 279             {
 280                 fmax2  = fam;
 281                 la_max = i;
 282             }
 283         }
 284     }
 285     else
 286     {
 287         for (i = start; i < end; i++)
 288         {
 289             fam     = norm2(f[i]);
 290             fnorm2 += fam;
 291             if (fam > fmax2)
 292             {
 293                 fmax2  = fam;
 294                 la_max = i;
 295             }
 296         }
 297     }
 298
 299     if (la_max >= 0 && DOMAINDECOMP(cr))
 300     {
 301         a_max = cr->dd->globalAtomIndices[la_max];
 302     }
 303     else
 304     {
 305         a_max = la_max;
 306     }
 307     if (PAR(cr))
 308     {
 309         snew(sum, 2*cr->nnodes+1);
 310         sum[2*cr->nodeid]   = fmax2;
 311         sum[2*cr->nodeid+1] = a_max;
 312         sum[2*cr->nnodes]   = fnorm2;
 313         gmx_sumd(2*cr->nnodes+1, sum, cr);
 314         fnorm2 = sum[2*cr->nnodes];
 315         /* Determine the global maximum */
 316         for (i = 0; i < cr->nnodes; i++)
 317         {
 318             if (sum[2*i] > fmax2)
 319             {
 320                 fmax2 = sum[2*i];
 321                 a_max = gmx::roundToInt(sum[2*i+1]);
 322             }
 323         }
 324         sfree(sum);
 325     }
 326
 327     if (fnorm)
 328     {
 329         *fnorm = sqrt(fnorm2);
 330     }
 331     if (fmax)
 332     {
 333         *fmax  = sqrt(fmax2);
 334     }
 335     if (a_fmax)
 336     {
 337         *a_fmax = a_max;
 338     }
 339 }
 340
 341 //! Compute the norm of the force
 342 static void get_state_f_norm_max(const t_commrec *cr,
 343                                  t_grpopts *opts, t_mdatoms *mdatoms,
 344                                  em_state_t *ems)
 345 {
 346     get_f_norm_max(cr, opts, mdatoms, ems->f.rvec_array(),
 347                    &ems->fnorm, &ems->fmax, &ems->a_fmax);
 348 }
 349
 350 //! Initialize the energy minimization
 351 static void init_em(FILE *fplog,
 352                     const gmx::MDLogger &mdlog,
 353                     const char *title,
 354                     const t_commrec *cr,
 355                     t_inputrec *ir,
 356                     gmx::ImdSession *imdSession,
 357                     t_state *state_global, gmx_mtop_t *top_global,
 358                     em_state_t *ems, gmx_localtop_t *top,
 359                     t_nrnb *nrnb,
 360                     t_forcerec *fr,
 361                     t_graph **graph, gmx::MDAtoms *mdAtoms, gmx_global_stat_t *gstat,
 362                     gmx_vsite_t *vsite, gmx::Constraints *constr, gmx_shellfc_t **shellfc)
 363 {
 364     real dvdl_constr;
 365
 366     if (fplog)
 367     {
 368         fprintf(fplog, "Initiating %s\n", title);
 369     }
 370
 371     if (MASTER(cr))
 372     {
 373         state_global->ngtc = 0;
 374     }
 375     initialize_lambdas(fplog, *ir, MASTER(cr), &(state_global->fep_state), state_global->lambda, nullptr);
 376
 377     init_nrnb(nrnb);
 378
 379     if (ir->eI == eiNM)
 380     {
 381         GMX_ASSERT(shellfc != nullptr, "With NM we always support shells");
 382
 383         *shellfc = init_shell_flexcon(stdout,
 384                                       top_global,
 385                                       constr ? constr->numFlexibleConstraints() : 0,
 386                                       ir->nstcalcenergy,
 387                                       DOMAINDECOMP(cr));
 388     }
 389     else
 390     {
 391         GMX_ASSERT(EI_ENERGY_MINIMIZATION(ir->eI), "This else currently only handles energy minimizers, consider if your algorithm needs shell/flexible-constraint support");
 392
 393         /* With energy minimization, shells and flexible constraints are
 394          * automatically minimized when treated like normal DOFS.
 395          */
 396         if (shellfc != nullptr)
 397         {
 398             *shellfc = nullptr;
 399         }
 400     }
 401
 402     auto mdatoms = mdAtoms->mdatoms();
 403     if (DOMAINDECOMP(cr))
 404     {
 405         top->useInDomainDecomp_ = true;
 406         dd_init_local_top(*top_global, top);
 407
 408         dd_init_local_state(cr->dd, state_global, &ems->s);
 409
 410         /* Distribute the charge groups over the nodes from the master node */
 411         dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
 412                             state_global, *top_global, ir, imdSession,
 413                             &ems->s, &ems->f, mdAtoms, top,
 414                             fr, vsite, constr,
 415                             nrnb, nullptr, FALSE);
 416         dd_store_state(cr->dd, &ems->s);
 417
 418         *graph = nullptr;
 419     }
 420     else
 421     {
 422         state_change_natoms(state_global, state_global->natoms);
 423         /* Just copy the state */
 424         ems->s = *state_global;
 425         state_change_natoms(&ems->s, ems->s.natoms);
 426         ems->f.resizeWithPadding(ems->s.natoms);
 427
 428         mdAlgorithmsSetupAtomData(cr, ir, *top_global, top, fr,
 429                                   graph, mdAtoms,
 430                                   constr, vsite, shellfc ? *shellfc : nullptr);
 431
 432         if (vsite)
 433         {
 434             set_vsite_top(vsite, top, mdatoms);
 435         }
 436     }
 437
 438     update_mdatoms(mdAtoms->mdatoms(), ems->s.lambda[efptMASS]);
 439
 440     if (constr)
 441     {
 442         // TODO how should this cross-module support dependency be managed?
 443         if (ir->eConstrAlg == econtSHAKE &&
 444             gmx_mtop_ftype_count(top_global, F_CONSTR) > 0)
 445         {
 446             gmx_fatal(FARGS, "Can not do energy minimization with %s, use %s\n",
 447                       econstr_names[econtSHAKE], econstr_names[econtLINCS]);
 448         }
 449
 450         if (!ir->bContinuation)
 451         {
 452             /* Constrain the starting coordinates */
 453             dvdl_constr = 0;
 454             constr->apply(TRUE, TRUE,
 455                           -1, 0, 1.0,
 456                           ems->s.x.rvec_array(),
 457                           ems->s.x.rvec_array(),
 458                           nullptr,
 459                           ems->s.box,
 460                           ems->s.lambda[efptFEP], &dvdl_constr,
 461                           nullptr, nullptr, gmx::ConstraintVariable::Positions);
 462         }
 463     }
 464
 465     if (PAR(cr))
 466     {
 467         *gstat = global_stat_init(ir);
 468     }
 469     else
 470     {
 471         *gstat = nullptr;
 472     }
 473
 474     calc_shifts(ems->s.box, fr->shift_vec);
 475 }
 476
 477 //! Finalize the minimization
 478 static void finish_em(const t_commrec *cr, gmx_mdoutf_t outf,
 479                       gmx_walltime_accounting_t walltime_accounting,
 480                       gmx_wallcycle_t wcycle)
 481 {
 482     if (!thisRankHasDuty(cr, DUTY_PME))
 483     {
 484         /* Tell the PME only node to finish */
 485         gmx_pme_send_finish(cr);
 486     }
 487
 488     done_mdoutf(outf);
 489
 490     em_time_end(walltime_accounting, wcycle);
 491 }
 492
 493 //! Swap two different EM states during minimization
 494 static void swap_em_state(em_state_t **ems1, em_state_t **ems2)
 495 {
 496     em_state_t *tmp;
 497
 498     tmp   = *ems1;
 499     *ems1 = *ems2;
 500     *ems2 = tmp;
 501 }
 502
 503 //! Save the EM trajectory
 504 static void write_em_traj(FILE *fplog, const t_commrec *cr,
 505                           gmx_mdoutf_t outf,
 506                           gmx_bool bX, gmx_bool bF, const char *confout,
 507                           gmx_mtop_t *top_global,
 508                           t_inputrec *ir, int64_t step,
 509                           em_state_t *state,
 510                           t_state *state_global,
 511                           ObservablesHistory *observablesHistory)
 512 {
 513     int mdof_flags = 0;
 514
 515     if (bX)
 516     {
 517         mdof_flags |= MDOF_X;
 518     }
 519     if (bF)
 520     {
 521         mdof_flags |= MDOF_F;
 522     }
 523
 524     /* If we want IMD output, set appropriate MDOF flag */
 525     if (ir->bIMD)
 526     {
 527         mdof_flags |= MDOF_IMD;
 528     }
 529
 530     mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
 531                                      top_global, step, static_cast<double>(step),
 532                                      &state->s, state_global, observablesHistory,
 533                                      state->f);
 534
 535     if (confout != nullptr)
 536     {
 537         if (DOMAINDECOMP(cr))
 538         {
 539             /* If bX=true, x was collected to state_global in the call above */
 540             if (!bX)
 541             {
 542                 auto globalXRef = MASTER(cr) ? state_global->x : gmx::ArrayRef<gmx::RVec>();
 543                 dd_collect_vec(cr->dd, &state->s, state->s.x, globalXRef);
 544             }
 545         }
 546         else
 547         {
 548             /* Copy the local state pointer */
 549             state_global = &state->s;
 550         }
 551
 552         if (MASTER(cr))
 553         {
 554             if (ir->ePBC != epbcNONE && !ir->bPeriodicMols && DOMAINDECOMP(cr))
 555             {
 556                 /* Make molecules whole only for confout writing */
 557                 do_pbc_mtop(ir->ePBC, state->s.box, top_global,
 558                             state_global->x.rvec_array());
 559             }
 560
 561             write_sto_conf_mtop(confout,
 562                                 *top_global->name, top_global,
 563                                 state_global->x.rvec_array(), nullptr, ir->ePBC, state->s.box);
 564         }
 565     }
 566 }
 567
 568 //! \brief Do one minimization step
 569 //
 570 // \returns true when the step succeeded, false when a constraint error occurred
 571 static bool do_em_step(const t_commrec *cr,
 572                        t_inputrec *ir, t_mdatoms *md,
 573                        em_state_t *ems1, real a, const PaddedVector<gmx::RVec> *force,
 574                        em_state_t *ems2,
 575                        gmx::Constraints *constr,
 576                        int64_t count)
 577
 578 {
 579     t_state *s1, *s2;
 580     int      start, end;
 581     real     dvdl_constr;
 582     int      nthreads gmx_unused;
 583
 584     bool     validStep = true;
 585
 586     s1 = &ems1->s;
 587     s2 = &ems2->s;
 588
 589     if (DOMAINDECOMP(cr) && s1->ddp_count != cr->dd->ddp_count)
 590     {
 591         gmx_incons("state mismatch in do_em_step");
 592     }
 593
 594     s2->flags = s1->flags;
 595
 596     if (s2->natoms != s1->natoms)
 597     {
 598         state_change_natoms(s2, s1->natoms);
 599         ems2->f.resizeWithPadding(s2->natoms);
 600     }
 601     if (DOMAINDECOMP(cr) && s2->cg_gl.size() != s1->cg_gl.size())
 602     {
 603         s2->cg_gl.resize(s1->cg_gl.size());
 604     }
 605
 606     copy_mat(s1->box, s2->box);
 607     /* Copy free energy state */
 608     s2->lambda = s1->lambda;
 609     copy_mat(s1->box, s2->box);
 610
 611     start = 0;
 612     end   = md->homenr;
 613
 614     nthreads = gmx_omp_nthreads_get(emntUpdate);
 615 #pragma omp parallel num_threads(nthreads)
 616     {
 617         const rvec *x1 = s1->x.rvec_array();
 618         rvec       *x2 = s2->x.rvec_array();
 619         const rvec *f  = force->rvec_array();
 620
 621         int         gf = 0;
 622 #pragma omp for schedule(static) nowait
 623         for (int i = start; i < end; i++)
 624         {
 625             try
 626             {
 627                 if (md->cFREEZE)
 628                 {
 629                     gf = md->cFREEZE[i];
 630                 }
 631                 for (int m = 0; m < DIM; m++)
 632                 {
 633                     if (ir->opts.nFreeze[gf][m])
 634                     {
 635                         x2[i][m] = x1[i][m];
 636                     }
 637                     else
 638                     {
 639                         x2[i][m] = x1[i][m] + a*f[i][m];
 640                     }
 641                 }
 642             }
 643             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 644         }
 645
 646         if (s2->flags & (1<<estCGP))
 647         {
 648             /* Copy the CG p vector */
 649             const rvec *p1 = s1->cg_p.rvec_array();
 650             rvec       *p2 = s2->cg_p.rvec_array();
 651 #pragma omp for schedule(static) nowait
 652             for (int i = start; i < end; i++)
 653             {
 654                 // Trivial OpenMP block that does not throw
 655                 copy_rvec(p1[i], p2[i]);
 656             }
 657         }
 658
 659         if (DOMAINDECOMP(cr))
 660         {
 661             s2->ddp_count = s1->ddp_count;
 662
 663             /* OpenMP does not supported unsigned loop variables */
 664 #pragma omp for schedule(static) nowait
 665             for (int i = 0; i < gmx::ssize(s2->cg_gl); i++)
 666             {
 667                 s2->cg_gl[i] = s1->cg_gl[i];
 668             }
 669             s2->ddp_count_cg_gl = s1->ddp_count_cg_gl;
 670         }
 671     }
 672
 673     if (constr)
 674     {
 675         dvdl_constr = 0;
 676         validStep   =
 677             constr->apply(TRUE, TRUE,
 678                           count, 0, 1.0,
 679                           s1->x.rvec_array(), s2->x.rvec_array(),
 680                           nullptr, s2->box,
 681                           s2->lambda[efptBONDED], &dvdl_constr,
 682                           nullptr, nullptr, gmx::ConstraintVariable::Positions);
 683
 684         if (cr->nnodes > 1)
 685         {
 686             /* This global reduction will affect performance at high
 687              * parallelization, but we can not really avoid it.
 688              * But usually EM is not run at high parallelization.
 689              */
 690             int reductionBuffer = static_cast<int>(!validStep);
 691             gmx_sumi(1, &reductionBuffer, cr);
 692             validStep           = (reductionBuffer == 0);
 693         }
 694
 695         // We should move this check to the different minimizers
 696         if (!validStep && ir->eI != eiSteep)
 697         {
 698             gmx_fatal(FARGS, "The coordinates could not be constrained. Minimizer '%s' can not handle constraint failures, use minimizer '%s' before using '%s'.",
 699                       EI(ir->eI), EI(eiSteep), EI(ir->eI));
 700         }
 701     }
 702
 703     return validStep;
 704 }
 705
 706 //! Prepare EM for using domain decomposition parallellization
 707 static void em_dd_partition_system(FILE *fplog,
 708                                    const gmx::MDLogger &mdlog,
 709                                    int step, const t_commrec *cr,
 710                                    gmx_mtop_t *top_global, t_inputrec *ir,
 711                                    gmx::ImdSession *imdSession,
 712                                    em_state_t *ems, gmx_localtop_t *top,
 713                                    gmx::MDAtoms *mdAtoms, t_forcerec *fr,
 714                                    gmx_vsite_t *vsite, gmx::Constraints *constr,
 715                                    t_nrnb *nrnb, gmx_wallcycle_t wcycle)
 716 {
 717     /* Repartition the domain decomposition */
 718     dd_partition_system(fplog, mdlog, step, cr, FALSE, 1,
 719                         nullptr, *top_global, ir, imdSession,
 720                         &ems->s, &ems->f,
 721                         mdAtoms, top, fr, vsite, constr,
 722                         nrnb, wcycle, FALSE);
 723     dd_store_state(cr->dd, &ems->s);
 724 }
 725
 726 namespace
 727 {
 728
 729 /*! \brief Class to handle the work of setting and doing an energy evaluation.
 730  *
 731  * This class is a mere aggregate of parameters to pass to evaluate an
 732  * energy, so that future changes to names and types of them consume
 733  * less time when refactoring other code.
 734  *
 735  * Aggregate initialization is used, for which the chief risk is that
 736  * if a member is added at the end and not all initializer lists are
 737  * updated, then the member will be value initialized, which will
 738  * typically mean initialization to zero.
 739  *
 740  * We only want to construct one of these with an initializer list, so
 741  * we explicitly delete the default constructor. */
 742 class EnergyEvaluator
 743 {
 744     public:
 745         //! We only intend to construct such objects with an initializer list.
 746         EnergyEvaluator() = delete;
 747         /*! \brief Evaluates an energy on the state in \c ems.
 748          *
 749          * \todo In practice, the same objects mu_tot, vir, and pres
 750          * are always passed to this function, so we would rather have
 751          * them as data members. However, their C-array types are
 752          * unsuited for aggregate initialization. When the types
 753          * improve, the call signature of this method can be reduced.
 754          */
 755         void run(em_state_t *ems, rvec mu_tot,
 756                  tensor vir, tensor pres,
 757                  int64_t count, gmx_bool bFirst);
 758         //! Handles logging (deprecated).
 759         FILE                 *fplog;
 760         //! Handles logging.
 761         const gmx::MDLogger  &mdlog;
 762         //! Handles communication.
 763         const t_commrec      *cr;
 764         //! Coordinates multi-simulations.
 765         const gmx_multisim_t *ms;
 766         //! Holds the simulation topology.
 767         gmx_mtop_t           *top_global;
 768         //! Holds the domain topology.
 769         gmx_localtop_t       *top;
 770         //! User input options.
 771         t_inputrec           *inputrec;
 772         //! The Interactive Molecular Dynamics session.
 773         gmx::ImdSession      *imdSession;
 774         //! Manages flop accounting.
 775         t_nrnb               *nrnb;
 776         //! Manages wall cycle accounting.
 777         gmx_wallcycle_t       wcycle;
 778         //! Coordinates global reduction.
 779         gmx_global_stat_t     gstat;
 780         //! Handles virtual sites.
 781         gmx_vsite_t          *vsite;
 782         //! Handles constraints.
 783         gmx::Constraints     *constr;
 784         //! Handles strange things.
 785         t_fcdata             *fcd;
 786         //! Molecular graph for SHAKE.
 787         t_graph              *graph;
 788         //! Per-atom data for this domain.
 789         gmx::MDAtoms         *mdAtoms;
 790         //! Handles how to calculate the forces.
 791         t_forcerec           *fr;
 792         //! Schedule of force-calculation work each step for this task.
 793         gmx::PpForceWorkload *ppForceWorkload;
 794         //! Stores the computed energies.
 795         gmx_enerdata_t       *enerd;
 796 };
 797
 798 void
 799 EnergyEvaluator::run(em_state_t *ems, rvec mu_tot,
 800                      tensor vir, tensor pres,
 801                      int64_t count, gmx_bool bFirst)
 802 {
 803     real     t;
 804     gmx_bool bNS;
 805     tensor   force_vir, shake_vir, ekin;
 806     real     dvdl_constr, prescorr, enercorr, dvdlcorr;
 807     real     terminate = 0;
 808
 809     /* Set the time to the initial time, the time does not change during EM */
 810     t = inputrec->init_t;
 811
 812     if (bFirst ||
 813         (DOMAINDECOMP(cr) && ems->s.ddp_count < cr->dd->ddp_count))
 814     {
 815         /* This is the first state or an old state used before the last ns */
 816         bNS = TRUE;
 817     }
 818     else
 819     {
 820         bNS = FALSE;
 821         if (inputrec->nstlist > 0)
 822         {
 823             bNS = TRUE;
 824         }
 825     }
 826
 827     if (vsite)
 828     {
 829         construct_vsites(vsite, ems->s.x.rvec_array(), 1, nullptr,
 830                          top->idef.iparams, top->idef.il,
 831                          fr->ePBC, fr->bMolPBC, cr, ems->s.box);
 832     }
 833
 834     if (DOMAINDECOMP(cr) && bNS)
 835     {
 836         /* Repartition the domain decomposition */
 837         em_dd_partition_system(fplog, mdlog, count, cr, top_global, inputrec, imdSession,
 838                                ems, top, mdAtoms, fr, vsite, constr,
 839                                nrnb, wcycle);
 840     }
 841
 842     /* Calc force & energy on new trial position  */
 843     /* do_force always puts the charge groups in the box and shifts again
 844      * We do not unshift, so molecules are always whole in congrad.c
 845      */
 846     do_force(fplog, cr, ms, inputrec, nullptr, nullptr, imdSession,
 847              count, nrnb, wcycle, top,
 848              ems->s.box, ems->s.x.arrayRefWithPadding(), &ems->s.hist,
 849              ems->f.arrayRefWithPadding(), force_vir, mdAtoms->mdatoms(), enerd, fcd,
 850              ems->s.lambda, graph, fr, ppForceWorkload, vsite, mu_tot, t, nullptr,
 851              GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES |
 852              GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY |
 853              (bNS ? GMX_FORCE_NS : 0),
 854              DDBalanceRegionHandler(cr));
 855
 856     /* Clear the unused shake virial and pressure */
 857     clear_mat(shake_vir);
 858     clear_mat(pres);
 859
 860     /* Communicate stuff when parallel */
 861     if (PAR(cr) && inputrec->eI != eiNM)
 862     {
 863         wallcycle_start(wcycle, ewcMoveE);
 864
 865         global_stat(gstat, cr, enerd, force_vir, shake_vir, mu_tot,
 866                     inputrec, nullptr, nullptr, nullptr, 1, &terminate,
 867                     nullptr, FALSE,
 868                     CGLO_ENERGY |
 869                     CGLO_PRESSURE |
 870                     CGLO_CONSTRAINT);
 871
 872         wallcycle_stop(wcycle, ewcMoveE);
 873     }
 874
 875     /* Calculate long range corrections to pressure and energy */
 876     calc_dispcorr(inputrec, fr, ems->s.box, ems->s.lambda[efptVDW],
 877                   pres, force_vir, &prescorr, &enercorr, &dvdlcorr);
 878     enerd->term[F_DISPCORR] = enercorr;
 879     enerd->term[F_EPOT]    += enercorr;
 880     enerd->term[F_PRES]    += prescorr;
 881     enerd->term[F_DVDL]    += dvdlcorr;
 882
 883     ems->epot = enerd->term[F_EPOT];
 884
 885     if (constr)
 886     {
 887         /* Project out the constraint components of the force */
 888         dvdl_constr = 0;
 889         rvec *f_rvec = ems->f.rvec_array();
 890         constr->apply(FALSE, FALSE,
 891                       count, 0, 1.0,
 892                       ems->s.x.rvec_array(), f_rvec, f_rvec,
 893                       ems->s.box,
 894                       ems->s.lambda[efptBONDED], &dvdl_constr,
 895                       nullptr, &shake_vir, gmx::ConstraintVariable::ForceDispl);
 896         enerd->term[F_DVDL_CONSTR] += dvdl_constr;
 897         m_add(force_vir, shake_vir, vir);
 898     }
 899     else
 900     {
 901         copy_mat(force_vir, vir);
 902     }
 903
 904     clear_mat(ekin);
 905     enerd->term[F_PRES] =
 906         calc_pres(fr->ePBC, inputrec->nwall, ems->s.box, ekin, vir, pres);
 907
 908     sum_dhdl(enerd, ems->s.lambda, inputrec->fepvals);
 909
 910     if (EI_ENERGY_MINIMIZATION(inputrec->eI))
 911     {
 912         get_state_f_norm_max(cr, &(inputrec->opts), mdAtoms->mdatoms(), ems);
 913     }
 914 }
 915
 916 } // namespace
 917
 918 //! Parallel utility summing energies and forces
 919 static double reorder_partsum(const t_commrec *cr, t_grpopts *opts,
 920                               gmx_mtop_t *top_global,
 921                               em_state_t *s_min, em_state_t *s_b)
 922 {
 923     t_block       *cgs_gl;
 924     int            ncg, *cg_gl, *index, c, cg, i, a0, a1, a, gf, m;
 925     double         partsum;
 926
 927     if (debug)
 928     {
 929         fprintf(debug, "Doing reorder_partsum\n");
 930     }
 931
 932     const rvec *fm = s_min->f.rvec_array();
 933     const rvec *fb = s_b->f.rvec_array();
 934
 935     cgs_gl = dd_charge_groups_global(cr->dd);
 936     index  = cgs_gl->index;
 937
 938     /* Collect fm in a global vector fmg.
 939      * This conflicts with the spirit of domain decomposition,
 940      * but to fully optimize this a much more complicated algorithm is required.
 941      */
 942     rvec *fmg;
 943     snew(fmg, top_global->natoms);
 944
 945     ncg   = s_min->s.cg_gl.size();
 946     cg_gl = s_min->s.cg_gl.data();
 947     i     = 0;
 948     for (c = 0; c < ncg; c++)
 949     {
 950         cg = cg_gl[c];
 951         a0 = index[cg];
 952         a1 = index[cg+1];
 953         for (a = a0; a < a1; a++)
 954         {
 955             copy_rvec(fm[i], fmg[a]);
 956             i++;
 957         }
 958     }
 959     gmx_sum(top_global->natoms*3, fmg[0], cr);
 960
 961     /* Now we will determine the part of the sum for the cgs in state s_b */
 962     ncg         = s_b->s.cg_gl.size();
 963     cg_gl       = s_b->s.cg_gl.data();
 964     partsum     = 0;
 965     i           = 0;
 966     gf          = 0;
 967     gmx::ArrayRef<unsigned char> grpnrFREEZE = top_global->groups.groupNumbers[SimulationAtomGroupType::Freeze];
 968     for (c = 0; c < ncg; c++)
 969     {
 970         cg = cg_gl[c];
 971         a0 = index[cg];
 972         a1 = index[cg+1];
 973         for (a = a0; a < a1; a++)
 974         {
 975             if (!grpnrFREEZE.empty())
 976             {
 977                 gf = grpnrFREEZE[i];
 978             }
 979             for (m = 0; m < DIM; m++)
 980             {
 981                 if (!opts->nFreeze[gf][m])
 982                 {
 983                     partsum += (fb[i][m] - fmg[a][m])*fb[i][m];
 984                 }
 985             }
 986             i++;
 987         }
 988     }
 989
 990     sfree(fmg);
 991
 992     return partsum;
 993 }
 994
 995 //! Print some stuff, like beta, whatever that means.
 996 static real pr_beta(const t_commrec *cr, t_grpopts *opts, t_mdatoms *mdatoms,
 997                     gmx_mtop_t *top_global,
 998                     em_state_t *s_min, em_state_t *s_b)
 999 {
1000     double sum;
1001
1002     /* This is just the classical Polak-Ribiere calculation of beta;
1003      * it looks a bit complicated since we take freeze groups into account,
1004      * and might have to sum it in parallel runs.
1005      */
1006
1007     if (!DOMAINDECOMP(cr) ||
1008         (s_min->s.ddp_count == cr->dd->ddp_count &&
1009          s_b->s.ddp_count   == cr->dd->ddp_count))
1010     {
1011         const rvec *fm  = s_min->f.rvec_array();
1012         const rvec *fb  = s_b->f.rvec_array();
1013         sum             = 0;
1014         int         gf  = 0;
1015         /* This part of code can be incorrect with DD,
1016          * since the atom ordering in s_b and s_min might differ.
1017          */
1018         for (int i = 0; i < mdatoms->homenr; i++)
1019         {
1020             if (mdatoms->cFREEZE)
1021             {
1022                 gf = mdatoms->cFREEZE[i];
1023             }
1024             for (int m = 0; m < DIM; m++)
1025             {
1026                 if (!opts->nFreeze[gf][m])
1027                 {
1028                     sum += (fb[i][m] - fm[i][m])*fb[i][m];
1029                 }
1030             }
1031         }
1032     }
1033     else
1034     {
1035         /* We need to reorder cgs while summing */
1036         sum = reorder_partsum(cr, opts, top_global, s_min, s_b);
1037     }
1038     if (PAR(cr))
1039     {
1040         gmx_sumd(1, &sum, cr);
1041     }
1042
1043     return sum/gmx::square(s_min->fnorm);
1044 }
1045
1046 namespace gmx
1047 {
1048
1049 void
1050 Integrator::do_cg()
1051 {
1052     const char        *CG = "Polak-Ribiere Conjugate Gradients";
1053
1054     gmx_localtop_t     top;
1055     gmx_global_stat_t  gstat;
1056     t_graph           *graph;
1057     double             tmp, minstep;
1058     real               stepsize;
1059     real               a, b, c, beta = 0.0;
1060     real               epot_repl = 0;
1061     real               pnorm;
1062     gmx_bool           converged, foundlower;
1063     rvec               mu_tot = {0};
1064     gmx_bool           do_log = FALSE, do_ene = FALSE, do_x, do_f;
1065     tensor             vir, pres;
1066     int                number_steps, neval = 0, nstcg = inputrec->nstcgsteep;
1067     int                m, step, nminstep;
1068     auto               mdatoms = mdAtoms->mdatoms();
1069
1070     GMX_LOG(mdlog.info).asParagraph().
1071         appendText("Note that activating conjugate gradient energy minimization via the "
1072                    "integrator .mdp option and the command gmx mdrun may "
1073                    "be available in a different form in a future version of GROMACS, "
1074                    "e.g. gmx minimize and an .mdp option.");
1075
1076     step = 0;
1077
1078     if (MASTER(cr))
1079     {
1080         // In CG, the state is extended with a search direction
1081         state_global->flags |= (1<<estCGP);
1082
1083         // Ensure the extra per-atom state array gets allocated
1084         state_change_natoms(state_global, state_global->natoms);
1085
1086         // Initialize the search direction to zero
1087         for (RVec &cg_p : state_global->cg_p)
1088         {
1089             cg_p = { 0, 0, 0 };
1090         }
1091     }
1092
1093     /* Create 4 states on the stack and extract pointers that we will swap */
1094     em_state_t  s0 {}, s1 {}, s2 {}, s3 {};
1095     em_state_t *s_min = &s0;
1096     em_state_t *s_a   = &s1;
1097     em_state_t *s_b   = &s2;
1098     em_state_t *s_c   = &s3;
1099
1100     /* Init em and store the local state in s_min */
1101     init_em(fplog, mdlog, CG, cr, inputrec, imdSession,
1102             state_global, top_global, s_min, &top,
1103             nrnb, fr, &graph, mdAtoms, &gstat,
1104             vsite, constr, nullptr);
1105     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, inputrec, top_global, nullptr, wcycle);
1106     gmx::EnergyOutput energyOutput;
1107     energyOutput.prepare(mdoutf_get_fp_ene(outf), top_global, inputrec, nullptr);
1108
1109     /* Print to log file */
1110     print_em_start(fplog, cr, walltime_accounting, wcycle, CG);
1111
1112     /* Max number of steps */
1113     number_steps = inputrec->nsteps;
1114
1115     if (MASTER(cr))
1116     {
1117         sp_header(stderr, CG, inputrec->em_tol, number_steps);
1118     }
1119     if (fplog)
1120     {
1121         sp_header(fplog, CG, inputrec->em_tol, number_steps);
1122     }
1123
1124     EnergyEvaluator energyEvaluator {
1125         fplog, mdlog, cr, ms,
1126         top_global, &top,
1127         inputrec, imdSession, nrnb, wcycle, gstat,
1128         vsite, constr, fcd, graph,
1129         mdAtoms, fr, ppForceWorkload, enerd
1130     };
1131     /* Call the force routine and some auxiliary (neighboursearching etc.) */
1132     /* do_force always puts the charge groups in the box and shifts again
1133      * We do not unshift, so molecules are always whole in congrad.c
1134      */
1135     energyEvaluator.run(s_min, mu_tot, vir, pres, -1, TRUE);
1136
1137     if (MASTER(cr))
1138     {
1139         /* Copy stuff to the energy bin for easy printing etc. */
1140         matrix nullBox = {};
1141         energyOutput.addDataAtEnergyStep(false, false, static_cast<double>(step),
1142                                          mdatoms->tmass, enerd, nullptr, nullptr, nullptr, nullBox,
1143                                          nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
1144
1145         print_ebin_header(fplog, step, step);
1146         energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE,
1147                                            fplog, step, step, eprNORMAL,
1148                                            fcd, &(top_global->groups), &(inputrec->opts), nullptr);
1149     }
1150
1151     /* Estimate/guess the initial stepsize */
1152     stepsize = inputrec->em_stepsize/s_min->fnorm;
1153
1154     if (MASTER(cr))
1155     {
1156         double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
1157         fprintf(stderr, "   F-max             = %12.5e on atom %d\n",
1158                 s_min->fmax, s_min->a_fmax+1);
1159         fprintf(stderr, "   F-Norm            = %12.5e\n",
1160                 s_min->fnorm/sqrtNumAtoms);
1161         fprintf(stderr, "\n");
1162         /* and copy to the log file too... */
1163         fprintf(fplog, "   F-max             = %12.5e on atom %d\n",
1164                 s_min->fmax, s_min->a_fmax+1);
1165         fprintf(fplog, "   F-Norm            = %12.5e\n",
1166                 s_min->fnorm/sqrtNumAtoms);
1167         fprintf(fplog, "\n");
1168     }
1169     /* Start the loop over CG steps.
1170      * Each successful step is counted, and we continue until
1171      * we either converge or reach the max number of steps.
1172      */
1173     converged = FALSE;
1174     for (step = 0; (number_steps < 0 || step <= number_steps) && !converged; step++)
1175     {
1176
1177         /* start taking steps in a new direction
1178          * First time we enter the routine, beta=0, and the direction is
1179          * simply the negative gradient.
1180          */
1181
1182         /* Calculate the new direction in p, and the gradient in this direction, gpa */
1183         rvec       *pm  = s_min->s.cg_p.rvec_array();
1184         const rvec *sfm = s_min->f.rvec_array();
1185         double      gpa = 0;
1186         int         gf  = 0;
1187         for (int i = 0; i < mdatoms->homenr; i++)
1188         {
1189             if (mdatoms->cFREEZE)
1190             {
1191                 gf = mdatoms->cFREEZE[i];
1192             }
1193             for (m = 0; m < DIM; m++)
1194             {
1195                 if (!inputrec->opts.nFreeze[gf][m])
1196                 {
1197                     pm[i][m] = sfm[i][m] + beta*pm[i][m];
1198                     gpa     -= pm[i][m]*sfm[i][m];
1199                     /* f is negative gradient, thus the sign */
1200                 }
1201                 else
1202                 {
1203                     pm[i][m] = 0;
1204                 }
1205             }
1206         }
1207
1208         /* Sum the gradient along the line across CPUs */
1209         if (PAR(cr))
1210         {
1211             gmx_sumd(1, &gpa, cr);
1212         }
1213
1214         /* Calculate the norm of the search vector */
1215         get_f_norm_max(cr, &(inputrec->opts), mdatoms, pm, &pnorm, nullptr, nullptr);
1216
1217         /* Just in case stepsize reaches zero due to numerical precision... */
1218         if (stepsize <= 0)
1219         {
1220             stepsize = inputrec->em_stepsize/pnorm;
1221         }
1222
1223         /*
1224          * Double check the value of the derivative in the search direction.
1225          * If it is positive it must be due to the old information in the
1226          * CG formula, so just remove that and start over with beta=0.
1227          * This corresponds to a steepest descent step.
1228          */
1229         if (gpa > 0)
1230         {
1231             beta = 0;
1232             step--;   /* Don't count this step since we are restarting */
1233             continue; /* Go back to the beginning of the big for-loop */
1234         }
1235
1236         /* Calculate minimum allowed stepsize, before the average (norm)
1237          * relative change in coordinate is smaller than precision
1238          */
1239         minstep = 0;
1240         auto s_min_x = makeArrayRef(s_min->s.x);
1241         for (int i = 0; i < mdatoms->homenr; i++)
1242         {
1243             for (m = 0; m < DIM; m++)
1244             {
1245                 tmp = fabs(s_min_x[i][m]);
1246                 if (tmp < 1.0)
1247                 {
1248                     tmp = 1.0;
1249                 }
1250                 tmp      = pm[i][m]/tmp;
1251                 minstep += tmp*tmp;
1252             }
1253         }
1254         /* Add up from all CPUs */
1255         if (PAR(cr))
1256         {
1257             gmx_sumd(1, &minstep, cr);
1258         }
1259
1260         minstep = GMX_REAL_EPS/sqrt(minstep/(3*top_global->natoms));
1261
1262         if (stepsize < minstep)
1263         {
1264             converged = TRUE;
1265             break;
1266         }
1267
1268         /* Write coordinates if necessary */
1269         do_x = do_per_step(step, inputrec->nstxout);
1270         do_f = do_per_step(step, inputrec->nstfout);
1271
1272         write_em_traj(fplog, cr, outf, do_x, do_f, nullptr,
1273                       top_global, inputrec, step,
1274                       s_min, state_global, observablesHistory);
1275
1276         /* Take a step downhill.
1277          * In theory, we should minimize the function along this direction.
1278          * That is quite possible, but it turns out to take 5-10 function evaluations
1279          * for each line. However, we dont really need to find the exact minimum -
1280          * it is much better to start a new CG step in a modified direction as soon
1281          * as we are close to it. This will save a lot of energy evaluations.
1282          *
1283          * In practice, we just try to take a single step.
1284          * If it worked (i.e. lowered the energy), we increase the stepsize but
1285          * the continue straight to the next CG step without trying to find any minimum.
1286          * If it didn't work (higher energy), there must be a minimum somewhere between
1287          * the old position and the new one.
1288          *
1289          * Due to the finite numerical accuracy, it turns out that it is a good idea
1290          * to even accept a SMALL increase in energy, if the derivative is still downhill.
1291          * This leads to lower final energies in the tests I've done. / Erik
1292          */
1293         s_a->epot = s_min->epot;
1294         a         = 0.0;
1295         c         = a + stepsize; /* reference position along line is zero */
1296
1297         if (DOMAINDECOMP(cr) && s_min->s.ddp_count < cr->dd->ddp_count)
1298         {
1299             em_dd_partition_system(fplog, mdlog, step, cr, top_global, inputrec, imdSession,
1300                                    s_min, &top, mdAtoms, fr, vsite, constr,
1301                                    nrnb, wcycle);
1302         }
1303
1304         /* Take a trial step (new coords in s_c) */
1305         do_em_step(cr, inputrec, mdatoms, s_min, c, &s_min->s.cg_p, s_c,
1306                    constr, -1);
1307
1308         neval++;
1309         /* Calculate energy for the trial step */
1310         energyEvaluator.run(s_c, mu_tot, vir, pres, -1, FALSE);
1311
1312         /* Calc derivative along line */
1313         const rvec *pc  = s_c->s.cg_p.rvec_array();
1314         const rvec *sfc = s_c->f.rvec_array();
1315         double      gpc = 0;
1316         for (int i = 0; i < mdatoms->homenr; i++)
1317         {
1318             for (m = 0; m < DIM; m++)
1319             {
1320                 gpc -= pc[i][m]*sfc[i][m]; /* f is negative gradient, thus the sign */
1321             }
1322         }
1323         /* Sum the gradient along the line across CPUs */
1324         if (PAR(cr))
1325         {
1326             gmx_sumd(1, &gpc, cr);
1327         }
1328
1329         /* This is the max amount of increase in energy we tolerate */
1330         tmp = std::sqrt(GMX_REAL_EPS)*fabs(s_a->epot);
1331
1332         /* Accept the step if the energy is lower, or if it is not significantly higher
1333          * and the line derivative is still negative.
1334          */
1335         if (s_c->epot < s_a->epot || (gpc < 0 && s_c->epot < (s_a->epot + tmp)))
1336         {
1337             foundlower = TRUE;
1338             /* Great, we found a better energy. Increase step for next iteration
1339              * if we are still going down, decrease it otherwise
1340              */
1341             if (gpc < 0)
1342             {
1343                 stepsize *= 1.618034; /* The golden section */
1344             }
1345             else
1346             {
1347                 stepsize *= 0.618034; /* 1/golden section */
1348             }
1349         }
1350         else
1351         {
1352             /* New energy is the same or higher. We will have to do some work
1353              * to find a smaller value in the interval. Take smaller step next time!
1354              */
1355             foundlower = FALSE;
1356             stepsize  *= 0.618034;
1357         }
1358
1359
1360
1361
1362         /* OK, if we didn't find a lower value we will have to locate one now - there must
1363          * be one in the interval [a=0,c].
1364          * The same thing is valid here, though: Don't spend dozens of iterations to find
1365          * the line minimum. We try to interpolate based on the derivative at the endpoints,
1366          * and only continue until we find a lower value. In most cases this means 1-2 iterations.
1367          *
1368          * I also have a safeguard for potentially really pathological functions so we never
1369          * take more than 20 steps before we give up ...
1370          *
1371          * If we already found a lower value we just skip this step and continue to the update.
1372          */
1373         double gpb;
1374         if (!foundlower)
1375         {
1376             nminstep = 0;
1377
1378             do
1379             {
1380                 /* Select a new trial point.
1381                  * If the derivatives at points a & c have different sign we interpolate to zero,
1382                  * otherwise just do a bisection.
1383                  */
1384                 if (gpa < 0 && gpc > 0)
1385                 {
1386                     b = a + gpa*(a-c)/(gpc-gpa);
1387                 }
1388                 else
1389                 {
1390                     b = 0.5*(a+c);
1391                 }
1392
1393                 /* safeguard if interpolation close to machine accuracy causes errors:
1394                  * never go outside the interval
1395                  */
1396                 if (b <= a || b >= c)
1397                 {
1398                     b = 0.5*(a+c);
1399                 }
1400
1401                 if (DOMAINDECOMP(cr) && s_min->s.ddp_count != cr->dd->ddp_count)
1402                 {
1403                     /* Reload the old state */
1404                     em_dd_partition_system(fplog, mdlog, -1, cr, top_global, inputrec, imdSession,
1405                                            s_min, &top, mdAtoms, fr, vsite, constr,
1406                                            nrnb, wcycle);
1407                 }
1408
1409                 /* Take a trial step to this new point - new coords in s_b */
1410                 do_em_step(cr, inputrec, mdatoms, s_min, b, &s_min->s.cg_p, s_b,
1411                            constr, -1);
1412
1413                 neval++;
1414                 /* Calculate energy for the trial step */
1415                 energyEvaluator.run(s_b, mu_tot, vir, pres, -1, FALSE);
1416
1417                 /* p does not change within a step, but since the domain decomposition
1418                  * might change, we have to use cg_p of s_b here.
1419                  */
1420                 const rvec *pb  = s_b->s.cg_p.rvec_array();
1421                 const rvec *sfb = s_b->f.rvec_array();
1422                 gpb             = 0;
1423                 for (int i = 0; i < mdatoms->homenr; i++)
1424                 {
1425                     for (m = 0; m < DIM; m++)
1426                     {
1427                         gpb -= pb[i][m]*sfb[i][m]; /* f is negative gradient, thus the sign */
1428                     }
1429                 }
1430                 /* Sum the gradient along the line across CPUs */
1431                 if (PAR(cr))
1432                 {
1433                     gmx_sumd(1, &gpb, cr);
1434                 }
1435
1436                 if (debug)
1437                 {
1438                     fprintf(debug, "CGE: EpotA %f EpotB %f EpotC %f gpb %f\n",
1439                             s_a->epot, s_b->epot, s_c->epot, gpb);
1440                 }
1441
1442                 epot_repl = s_b->epot;
1443
1444                 /* Keep one of the intervals based on the value of the derivative at the new point */
1445                 if (gpb > 0)
1446                 {
1447                     /* Replace c endpoint with b */
1448                     swap_em_state(&s_b, &s_c);
1449                     c   = b;
1450                     gpc = gpb;
1451                 }
1452                 else
1453                 {
1454                     /* Replace a endpoint with b */
1455                     swap_em_state(&s_b, &s_a);
1456                     a   = b;
1457                     gpa = gpb;
1458                 }
1459
1460                 /*
1461                  * Stop search as soon as we find a value smaller than the endpoints.
1462                  * Never run more than 20 steps, no matter what.
1463                  */
1464                 nminstep++;
1465             }
1466             while ((epot_repl > s_a->epot || epot_repl > s_c->epot) &&
1467                    (nminstep < 20));
1468
1469             if (std::fabs(epot_repl - s_min->epot) < fabs(s_min->epot)*GMX_REAL_EPS ||
1470                 nminstep >= 20)
1471             {
1472                 /* OK. We couldn't find a significantly lower energy.
1473                  * If beta==0 this was steepest descent, and then we give up.
1474                  * If not, set beta=0 and restart with steepest descent before quitting.
1475                  */
1476                 if (beta == 0.0)
1477                 {
1478                     /* Converged */
1479                     converged = TRUE;
1480                     break;
1481                 }
1482                 else
1483                 {
1484                     /* Reset memory before giving up */
1485                     beta = 0.0;
1486                     continue;
1487                 }
1488             }
1489
1490             /* Select min energy state of A & C, put the best in B.
1491              */
1492             if (s_c->epot < s_a->epot)
1493             {
1494                 if (debug)
1495                 {
1496                     fprintf(debug, "CGE: C (%f) is lower than A (%f), moving C to B\n",
1497                             s_c->epot, s_a->epot);
1498                 }
1499                 swap_em_state(&s_b, &s_c);
1500                 gpb = gpc;
1501             }
1502             else
1503             {
1504                 if (debug)
1505                 {
1506                     fprintf(debug, "CGE: A (%f) is lower than C (%f), moving A to B\n",
1507                             s_a->epot, s_c->epot);
1508                 }
1509                 swap_em_state(&s_b, &s_a);
1510                 gpb = gpa;
1511             }
1512
1513         }
1514         else
1515         {
1516             if (debug)
1517             {
1518                 fprintf(debug, "CGE: Found a lower energy %f, moving C to B\n",
1519                         s_c->epot);
1520             }
1521             swap_em_state(&s_b, &s_c);
1522             gpb = gpc;
1523         }
1524
1525         /* new search direction */
1526         /* beta = 0 means forget all memory and restart with steepest descents. */
1527         if (nstcg && ((step % nstcg) == 0))
1528         {
1529             beta = 0.0;
1530         }
1531         else
1532         {
1533             /* s_min->fnorm cannot be zero, because then we would have converged
1534              * and broken out.
1535              */
1536
1537             /* Polak-Ribiere update.
1538              * Change to fnorm2/fnorm2_old for Fletcher-Reeves
1539              */
1540             beta = pr_beta(cr, &inputrec->opts, mdatoms, top_global, s_min, s_b);
1541         }
1542         /* Limit beta to prevent oscillations */
1543         if (fabs(beta) > 5.0)
1544         {
1545             beta = 0.0;
1546         }
1547
1548
1549         /* update positions */
1550         swap_em_state(&s_min, &s_b);
1551         gpa = gpb;
1552
1553         /* Print it if necessary */
1554         if (MASTER(cr))
1555         {
1556             if (mdrunOptions.verbose)
1557             {
1558                 double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
1559                 fprintf(stderr, "\rStep %d, Epot=%12.6e, Fnorm=%9.3e, Fmax=%9.3e (atom %d)\n",
1560                         step, s_min->epot, s_min->fnorm/sqrtNumAtoms,
1561                         s_min->fmax, s_min->a_fmax+1);
1562                 fflush(stderr);
1563             }
1564             /* Store the new (lower) energies */
1565             matrix nullBox = {};
1566             energyOutput.addDataAtEnergyStep(false, false, static_cast<double>(step),
1567                                              mdatoms->tmass, enerd, nullptr, nullptr, nullptr, nullBox,
1568                                              nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
1569
1570             do_log = do_per_step(step, inputrec->nstlog);
1571             do_ene = do_per_step(step, inputrec->nstenergy);
1572
1573             imdSession->fillEnergyRecord(step, TRUE);
1574
1575             if (do_log)
1576             {
1577                 print_ebin_header(fplog, step, step);
1578             }
1579             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
1580                                                do_log ? fplog : nullptr, step, step, eprNORMAL,
1581                                                fcd, &(top_global->groups), &(inputrec->opts), nullptr);
1582         }
1583
1584         /* Send energies and positions to the IMD client if bIMD is TRUE. */
1585         if (MASTER(cr) && imdSession->run(step, TRUE, state_global->box, state_global->x.rvec_array(), 0))
1586         {
1587             imdSession->sendPositionsAndEnergies();
1588         }
1589
1590         /* Stop when the maximum force lies below tolerance.
1591          * If we have reached machine precision, converged is already set to true.
1592          */
1593         converged = converged || (s_min->fmax < inputrec->em_tol);
1594
1595     }   /* End of the loop */
1596
1597     if (converged)
1598     {
1599         step--; /* we never took that last step in this case */
1600
1601     }
1602     if (s_min->fmax > inputrec->em_tol)
1603     {
1604         if (MASTER(cr))
1605         {
1606             warn_step(fplog, inputrec->em_tol, s_min->fmax,
1607                       step-1 == number_steps, FALSE);
1608         }
1609         converged = FALSE;
1610     }
1611
1612     if (MASTER(cr))
1613     {
1614         /* If we printed energy and/or logfile last step (which was the last step)
1615          * we don't have to do it again, but otherwise print the final values.
1616          */
1617         if (!do_log)
1618         {
1619             /* Write final value to log since we didn't do anything the last step */
1620             print_ebin_header(fplog, step, step);
1621         }
1622         if (!do_ene || !do_log)
1623         {
1624             /* Write final energy file entries */
1625             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
1626                                                !do_log ? fplog : nullptr, step, step, eprNORMAL,
1627                                                fcd, &(top_global->groups), &(inputrec->opts), nullptr);
1628         }
1629     }
1630
1631     /* Print some stuff... */
1632     if (MASTER(cr))
1633     {
1634         fprintf(stderr, "\nwriting lowest energy coordinates.\n");
1635     }
1636
1637     /* IMPORTANT!
1638      * For accurate normal mode calculation it is imperative that we
1639      * store the last conformation into the full precision binary trajectory.
1640      *
1641      * However, we should only do it if we did NOT already write this step
1642      * above (which we did if do_x or do_f was true).
1643      */
1644     /* Note that with 0 < nstfout != nstxout we can end up with two frames
1645      * in the trajectory with the same step number.
1646      */
1647     do_x = !do_per_step(step, inputrec->nstxout);
1648     do_f = (inputrec->nstfout > 0 && !do_per_step(step, inputrec->nstfout));
1649
1650     write_em_traj(fplog, cr, outf, do_x, do_f, ftp2fn(efSTO, nfile, fnm),
1651                   top_global, inputrec, step,
1652                   s_min, state_global, observablesHistory);
1653
1654
1655     if (MASTER(cr))
1656     {
1657         double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
1658         print_converged(stderr, CG, inputrec->em_tol, step, converged, number_steps,
1659                         s_min, sqrtNumAtoms);
1660         print_converged(fplog, CG, inputrec->em_tol, step, converged, number_steps,
1661                         s_min, sqrtNumAtoms);
1662
1663         fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
1664     }
1665
1666     finish_em(cr, outf, walltime_accounting, wcycle);
1667
1668     /* To print the actual number of steps we needed somewhere */
1669     walltime_accounting_set_nsteps_done(walltime_accounting, step);
1670 }
1671
1672
1673 void
1674 Integrator::do_lbfgs()
1675 {
1676     static const char *LBFGS = "Low-Memory BFGS Minimizer";
1677     em_state_t         ems;
1678     gmx_localtop_t     top;
1679     gmx_global_stat_t  gstat;
1680     t_graph           *graph;
1681     int                ncorr, nmaxcorr, point, cp, neval, nminstep;
1682     double             stepsize, step_taken, gpa, gpb, gpc, tmp, minstep;
1683     real              *rho, *alpha, *p, *s, **dx, **dg;
1684     real               a, b, c, maxdelta, delta;
1685     real               diag, Epot0;
1686     real               dgdx, dgdg, sq, yr, beta;
1687     gmx_bool           converged;
1688     rvec               mu_tot = {0};
1689     gmx_bool           do_log, do_ene, do_x, do_f, foundlower, *frozen;
1690     tensor             vir, pres;
1691     int                start, end, number_steps;
1692     int                i, k, m, n, gf, step;
1693     int                mdof_flags;
1694     auto               mdatoms = mdAtoms->mdatoms();
1695
1696     GMX_LOG(mdlog.info).asParagraph().
1697         appendText("Note that activating L-BFGS energy minimization via the "
1698                    "integrator .mdp option and the command gmx mdrun may "
1699                    "be available in a different form in a future version of GROMACS, "
1700                    "e.g. gmx minimize and an .mdp option.");
1701
1702     if (PAR(cr))
1703     {
1704         gmx_fatal(FARGS, "L-BFGS minimization only supports a single rank");
1705     }
1706
1707     if (nullptr != constr)
1708     {
1709         gmx_fatal(FARGS, "The combination of constraints and L-BFGS minimization is not implemented. Either do not use constraints, or use another minimizer (e.g. steepest descent).");
1710     }
1711
1712     n        = 3*state_global->natoms;
1713     nmaxcorr = inputrec->nbfgscorr;
1714
1715     snew(frozen, n);
1716
1717     snew(p, n);
1718     snew(rho, nmaxcorr);
1719     snew(alpha, nmaxcorr);
1720
1721     snew(dx, nmaxcorr);
1722     for (i = 0; i < nmaxcorr; i++)
1723     {
1724         snew(dx[i], n);
1725     }
1726
1727     snew(dg, nmaxcorr);
1728     for (i = 0; i < nmaxcorr; i++)
1729     {
1730         snew(dg[i], n);
1731     }
1732
1733     step  = 0;
1734     neval = 0;
1735
1736     /* Init em */
1737     init_em(fplog, mdlog, LBFGS, cr, inputrec, imdSession,
1738             state_global, top_global, &ems, &top,
1739             nrnb, fr, &graph, mdAtoms, &gstat,
1740             vsite, constr, nullptr);
1741     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, inputrec, top_global, nullptr, wcycle);
1742     gmx::EnergyOutput energyOutput;
1743     energyOutput.prepare(mdoutf_get_fp_ene(outf), top_global, inputrec, nullptr);
1744
1745     start = 0;
1746     end   = mdatoms->homenr;
1747
1748     /* We need 4 working states */
1749     em_state_t  s0 {}, s1 {}, s2 {}, s3 {};
1750     em_state_t *sa   = &s0;
1751     em_state_t *sb   = &s1;
1752     em_state_t *sc   = &s2;
1753     em_state_t *last = &s3;
1754     /* Initialize by copying the state from ems (we could skip x and f here) */
1755     *sa              = ems;
1756     *sb              = ems;
1757     *sc              = ems;
1758
1759     /* Print to log file */
1760     print_em_start(fplog, cr, walltime_accounting, wcycle, LBFGS);
1761
1762     do_log = do_ene = do_x = do_f = TRUE;
1763
1764     /* Max number of steps */
1765     number_steps = inputrec->nsteps;
1766
1767     /* Create a 3*natoms index to tell whether each degree of freedom is frozen */
1768     gf = 0;
1769     for (i = start; i < end; i++)
1770     {
1771         if (mdatoms->cFREEZE)
1772         {
1773             gf = mdatoms->cFREEZE[i];
1774         }
1775         for (m = 0; m < DIM; m++)
1776         {
1777             frozen[3*i+m] = (inputrec->opts.nFreeze[gf][m] != 0);
1778         }
1779     }
1780     if (MASTER(cr))
1781     {
1782         sp_header(stderr, LBFGS, inputrec->em_tol, number_steps);
1783     }
1784     if (fplog)
1785     {
1786         sp_header(fplog, LBFGS, inputrec->em_tol, number_steps);
1787     }
1788
1789     if (vsite)
1790     {
1791         construct_vsites(vsite, state_global->x.rvec_array(), 1, nullptr,
1792                          top.idef.iparams, top.idef.il,
1793                          fr->ePBC, fr->bMolPBC, cr, state_global->box);
1794     }
1795
1796     /* Call the force routine and some auxiliary (neighboursearching etc.) */
1797     /* do_force always puts the charge groups in the box and shifts again
1798      * We do not unshift, so molecules are always whole
1799      */
1800     neval++;
1801     EnergyEvaluator energyEvaluator {
1802         fplog, mdlog, cr, ms,
1803         top_global, &top,
1804         inputrec, imdSession, nrnb, wcycle, gstat,
1805         vsite, constr, fcd, graph,
1806         mdAtoms, fr, ppForceWorkload, enerd
1807     };
1808     energyEvaluator.run(&ems, mu_tot, vir, pres, -1, TRUE);
1809
1810     if (MASTER(cr))
1811     {
1812         /* Copy stuff to the energy bin for easy printing etc. */
1813         matrix nullBox = {};
1814         energyOutput.addDataAtEnergyStep(false, false, static_cast<double>(step),
1815                                          mdatoms->tmass, enerd, nullptr, nullptr, nullptr, nullBox,
1816                                          nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
1817
1818         print_ebin_header(fplog, step, step);
1819         energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE,
1820                                            fplog, step, step, eprNORMAL,
1821                                            fcd, &(top_global->groups), &(inputrec->opts), nullptr);
1822     }
1823
1824     /* Set the initial step.
1825      * since it will be multiplied by the non-normalized search direction
1826      * vector (force vector the first time), we scale it by the
1827      * norm of the force.
1828      */
1829
1830     if (MASTER(cr))
1831     {
1832         double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
1833         fprintf(stderr, "Using %d BFGS correction steps.\n\n", nmaxcorr);
1834         fprintf(stderr, "   F-max             = %12.5e on atom %d\n", ems.fmax, ems.a_fmax + 1);
1835         fprintf(stderr, "   F-Norm            = %12.5e\n", ems.fnorm/sqrtNumAtoms);
1836         fprintf(stderr, "\n");
1837         /* and copy to the log file too... */
1838         fprintf(fplog, "Using %d BFGS correction steps.\n\n", nmaxcorr);
1839         fprintf(fplog, "   F-max             = %12.5e on atom %d\n", ems.fmax, ems.a_fmax + 1);
1840         fprintf(fplog, "   F-Norm            = %12.5e\n", ems.fnorm/sqrtNumAtoms);
1841         fprintf(fplog, "\n");
1842     }
1843
1844     // Point is an index to the memory of search directions, where 0 is the first one.
1845     point = 0;
1846
1847     // Set initial search direction to the force (-gradient), or 0 for frozen particles.
1848     real *fInit = static_cast<real *>(ems.f.rvec_array()[0]);
1849     for (i = 0; i < n; i++)
1850     {
1851         if (!frozen[i])
1852         {
1853             dx[point][i] = fInit[i]; /* Initial search direction */
1854         }
1855         else
1856         {
1857             dx[point][i] = 0;
1858         }
1859     }
1860
1861     // Stepsize will be modified during the search, and actually it is not critical
1862     // (the main efficiency in the algorithm comes from changing directions), but
1863     // we still need an initial value, so estimate it as the inverse of the norm
1864     // so we take small steps where the potential fluctuates a lot.
1865     stepsize  = 1.0/ems.fnorm;
1866
1867     /* Start the loop over BFGS steps.
1868      * Each successful step is counted, and we continue until
1869      * we either converge or reach the max number of steps.
1870      */
1871
1872     ncorr = 0;
1873
1874     /* Set the gradient from the force */
1875     converged = FALSE;
1876     for (step = 0; (number_steps < 0 || step <= number_steps) && !converged; step++)
1877     {
1878
1879         /* Write coordinates if necessary */
1880         do_x = do_per_step(step, inputrec->nstxout);
1881         do_f = do_per_step(step, inputrec->nstfout);
1882
1883         mdof_flags = 0;
1884         if (do_x)
1885         {
1886             mdof_flags |= MDOF_X;
1887         }
1888
1889         if (do_f)
1890         {
1891             mdof_flags |= MDOF_F;
1892         }
1893
1894         if (inputrec->bIMD)
1895         {
1896             mdof_flags |= MDOF_IMD;
1897         }
1898
1899         mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
1900                                          top_global, step, static_cast<real>(step), &ems.s, state_global, observablesHistory, ems.f);
1901
1902         /* Do the linesearching in the direction dx[point][0..(n-1)] */
1903
1904         /* make s a pointer to current search direction - point=0 first time we get here */
1905         s = dx[point];
1906
1907         real *xx = static_cast<real *>(ems.s.x.rvec_array()[0]);
1908         real *ff = static_cast<real *>(ems.f.rvec_array()[0]);
1909
1910         // calculate line gradient in position A
1911         for (gpa = 0, i = 0; i < n; i++)
1912         {
1913             gpa -= s[i]*ff[i];
1914         }
1915
1916         /* Calculate minimum allowed stepsize along the line, before the average (norm)
1917          * relative change in coordinate is smaller than precision
1918          */
1919         for (minstep = 0, i = 0; i < n; i++)
1920         {
1921             tmp = fabs(xx[i]);
1922             if (tmp < 1.0)
1923             {
1924                 tmp = 1.0;
1925             }
1926             tmp      = s[i]/tmp;
1927             minstep += tmp*tmp;
1928         }
1929         minstep = GMX_REAL_EPS/sqrt(minstep/n);
1930
1931         if (stepsize < minstep)
1932         {
1933             converged = TRUE;
1934             break;
1935         }
1936
1937         // Before taking any steps along the line, store the old position
1938         *last       = ems;
1939         real *lastx = static_cast<real *>(last->s.x.data()[0]);
1940         real *lastf = static_cast<real *>(last->f.data()[0]);
1941         Epot0       = ems.epot;
1942
1943         *sa         = ems;
1944
1945         /* Take a step downhill.
1946          * In theory, we should find the actual minimum of the function in this
1947          * direction, somewhere along the line.
1948          * That is quite possible, but it turns out to take 5-10 function evaluations
1949          * for each line. However, we dont really need to find the exact minimum -
1950          * it is much better to start a new BFGS step in a modified direction as soon
1951          * as we are close to it. This will save a lot of energy evaluations.
1952          *
1953          * In practice, we just try to take a single step.
1954          * If it worked (i.e. lowered the energy), we increase the stepsize but
1955          * continue straight to the next BFGS step without trying to find any minimum,
1956          * i.e. we change the search direction too. If the line was smooth, it is
1957          * likely we are in a smooth region, and then it makes sense to take longer
1958          * steps in the modified search direction too.
1959          *
1960          * If it didn't work (higher energy), there must be a minimum somewhere between
1961          * the old position and the new one. Then we need to start by finding a lower
1962          * value before we change search direction. Since the energy was apparently
1963          * quite rough, we need to decrease the step size.
1964          *
1965          * Due to the finite numerical accuracy, it turns out that it is a good idea
1966          * to accept a SMALL increase in energy, if the derivative is still downhill.
1967          * This leads to lower final energies in the tests I've done. / Erik
1968          */
1969
1970         // State "A" is the first position along the line.
1971         // reference position along line is initially zero
1972         a          = 0.0;
1973
1974         // Check stepsize first. We do not allow displacements
1975         // larger than emstep.
1976         //
1977         do
1978         {
1979             // Pick a new position C by adding stepsize to A.
1980             c        = a + stepsize;
1981
1982             // Calculate what the largest change in any individual coordinate
1983             // would be (translation along line * gradient along line)
1984             maxdelta = 0;
1985             for (i = 0; i < n; i++)
1986             {
1987                 delta = c*s[i];
1988                 if (delta > maxdelta)
1989                 {
1990                     maxdelta = delta;
1991                 }
1992             }
1993             // If any displacement is larger than the stepsize limit, reduce the step
1994             if (maxdelta > inputrec->em_stepsize)
1995             {
1996                 stepsize *= 0.1;
1997             }
1998         }
1999         while (maxdelta > inputrec->em_stepsize);
2000
2001         // Take a trial step and move the coordinate array xc[] to position C
2002         real *xc = static_cast<real *>(sc->s.x.rvec_array()[0]);
2003         for (i = 0; i < n; i++)
2004         {
2005             xc[i] = lastx[i] + c*s[i];
2006         }
2007
2008         neval++;
2009         // Calculate energy for the trial step in position C
2010         energyEvaluator.run(sc, mu_tot, vir, pres, step, FALSE);
2011
2012         // Calc line gradient in position C
2013         real *fc = static_cast<real *>(sc->f.rvec_array()[0]);
2014         for (gpc = 0, i = 0; i < n; i++)
2015         {
2016             gpc -= s[i]*fc[i]; /* f is negative gradient, thus the sign */
2017         }
2018         /* Sum the gradient along the line across CPUs */
2019         if (PAR(cr))
2020         {
2021             gmx_sumd(1, &gpc, cr);
2022         }
2023
2024         // This is the max amount of increase in energy we tolerate.
2025         // By allowing VERY small changes (close to numerical precision) we
2026         // frequently find even better (lower) final energies.
2027         tmp = std::sqrt(GMX_REAL_EPS)*fabs(sa->epot);
2028
2029         // Accept the step if the energy is lower in the new position C (compared to A),
2030         // or if it is not significantly higher and the line derivative is still negative.
2031         foundlower = sc->epot < sa->epot || (gpc < 0 && sc->epot < (sa->epot + tmp));
2032         // If true, great, we found a better energy. We no longer try to alter the
2033         // stepsize, but simply accept this new better position. The we select a new
2034         // search direction instead, which will be much more efficient than continuing
2035         // to take smaller steps along a line. Set fnorm based on the new C position,
2036         // which will be used to update the stepsize to 1/fnorm further down.
2037
2038         // If false, the energy is NOT lower in point C, i.e. it will be the same
2039         // or higher than in point A. In this case it is pointless to move to point C,
2040         // so we will have to do more iterations along the same line to find a smaller
2041         // value in the interval [A=0.0,C].
2042         // Here, A is still 0.0, but that will change when we do a search in the interval
2043         // [0.0,C] below. That search we will do by interpolation or bisection rather
2044         // than with the stepsize, so no need to modify it. For the next search direction
2045         // it will be reset to 1/fnorm anyway.
2046
2047         if (!foundlower)
2048         {
2049             // OK, if we didn't find a lower value we will have to locate one now - there must
2050             // be one in the interval [a,c].
2051             // The same thing is valid here, though: Don't spend dozens of iterations to find
2052             // the line minimum. We try to interpolate based on the derivative at the endpoints,
2053             // and only continue until we find a lower value. In most cases this means 1-2 iterations.
2054             // I also have a safeguard for potentially really pathological functions so we never
2055             // take more than 20 steps before we give up.
2056             // If we already found a lower value we just skip this step and continue to the update.
2057             real fnorm = 0;
2058             nminstep   = 0;
2059             do
2060             {
2061                 // Select a new trial point B in the interval [A,C].
2062                 // If the derivatives at points a & c have different sign we interpolate to zero,
2063                 // otherwise just do a bisection since there might be multiple minima/maxima
2064                 // inside the interval.
2065                 if (gpa < 0 && gpc > 0)
2066                 {
2067                     b = a + gpa*(a-c)/(gpc-gpa);
2068                 }
2069                 else
2070                 {
2071                     b = 0.5*(a+c);
2072                 }
2073
2074                 /* safeguard if interpolation close to machine accuracy causes errors:
2075                  * never go outside the interval
2076                  */
2077                 if (b <= a || b >= c)
2078                 {
2079                     b = 0.5*(a+c);
2080                 }
2081
2082                 // Take a trial step to point B
2083                 real *xb = static_cast<real *>(sb->s.x.rvec_array()[0]);
2084                 for (i = 0; i < n; i++)
2085                 {
2086                     xb[i] = lastx[i] + b*s[i];
2087                 }
2088
2089                 neval++;
2090                 // Calculate energy for the trial step in point B
2091                 energyEvaluator.run(sb, mu_tot, vir, pres, step, FALSE);
2092                 fnorm = sb->fnorm;
2093
2094                 // Calculate gradient in point B
2095                 real *fb = static_cast<real *>(sb->f.rvec_array()[0]);
2096                 for (gpb = 0, i = 0; i < n; i++)
2097                 {
2098                     gpb -= s[i]*fb[i]; /* f is negative gradient, thus the sign */
2099
2100                 }
2101                 /* Sum the gradient along the line across CPUs */
2102                 if (PAR(cr))
2103                 {
2104                     gmx_sumd(1, &gpb, cr);
2105                 }
2106
2107                 // Keep one of the intervals [A,B] or [B,C] based on the value of the derivative
2108                 // at the new point B, and rename the endpoints of this new interval A and C.
2109                 if (gpb > 0)
2110                 {
2111                     /* Replace c endpoint with b */
2112                     c   = b;
2113                     /* swap states b and c */
2114                     swap_em_state(&sb, &sc);
2115                 }
2116                 else
2117                 {
2118                     /* Replace a endpoint with b */
2119                     a   = b;
2120                     /* swap states a and b */
2121                     swap_em_state(&sa, &sb);
2122                 }
2123
2124                 /*
2125                  * Stop search as soon as we find a value smaller than the endpoints,
2126                  * or if the tolerance is below machine precision.
2127                  * Never run more than 20 steps, no matter what.
2128                  */
2129                 nminstep++;
2130             }
2131             while ((sb->epot > sa->epot || sb->epot > sc->epot) && (nminstep < 20));
2132
2133             if (std::fabs(sb->epot - Epot0) < GMX_REAL_EPS || nminstep >= 20)
2134             {
2135                 /* OK. We couldn't find a significantly lower energy.
2136                  * If ncorr==0 this was steepest descent, and then we give up.
2137                  * If not, reset memory to restart as steepest descent before quitting.
2138                  */
2139                 if (ncorr == 0)
2140                 {
2141                     /* Converged */
2142                     converged = TRUE;
2143                     break;
2144                 }
2145                 else
2146                 {
2147                     /* Reset memory */
2148                     ncorr = 0;
2149                     /* Search in gradient direction */
2150                     for (i = 0; i < n; i++)
2151                     {
2152                         dx[point][i] = ff[i];
2153                     }
2154                     /* Reset stepsize */
2155                     stepsize = 1.0/fnorm;
2156                     continue;
2157                 }
2158             }
2159
2160             /* Select min energy state of A & C, put the best in xx/ff/Epot
2161              */
2162             if (sc->epot < sa->epot)
2163             {
2164                 /* Use state C */
2165                 ems        = *sc;
2166                 step_taken = c;
2167             }
2168             else
2169             {
2170                 /* Use state A */
2171                 ems        = *sa;
2172                 step_taken = a;
2173             }
2174
2175         }
2176         else
2177         {
2178             /* found lower */
2179             /* Use state C */
2180             ems        = *sc;
2181             step_taken = c;
2182         }
2183
2184         /* Update the memory information, and calculate a new
2185          * approximation of the inverse hessian
2186          */
2187
2188         /* Have new data in Epot, xx, ff */
2189         if (ncorr < nmaxcorr)
2190         {
2191             ncorr++;
2192         }
2193
2194         for (i = 0; i < n; i++)
2195         {
2196             dg[point][i]  = lastf[i]-ff[i];
2197             dx[point][i] *= step_taken;
2198         }
2199
2200         dgdg = 0;
2201         dgdx = 0;
2202         for (i = 0; i < n; i++)
2203         {
2204             dgdg += dg[point][i]*dg[point][i];
2205             dgdx += dg[point][i]*dx[point][i];
2206         }
2207
2208         diag = dgdx/dgdg;
2209
2210         rho[point] = 1.0/dgdx;
2211         point++;
2212
2213         if (point >= nmaxcorr)
2214         {
2215             point = 0;
2216         }
2217
2218         /* Update */
2219         for (i = 0; i < n; i++)
2220         {
2221             p[i] = ff[i];
2222         }
2223
2224         cp = point;
2225
2226         /* Recursive update. First go back over the memory points */
2227         for (k = 0; k < ncorr; k++)
2228         {
2229             cp--;
2230             if (cp < 0)
2231             {
2232                 cp = ncorr-1;
2233             }
2234
2235             sq = 0;
2236             for (i = 0; i < n; i++)
2237             {
2238                 sq += dx[cp][i]*p[i];
2239             }
2240
2241             alpha[cp] = rho[cp]*sq;
2242
2243             for (i = 0; i < n; i++)
2244             {
2245                 p[i] -= alpha[cp]*dg[cp][i];
2246             }
2247         }
2248
2249         for (i = 0; i < n; i++)
2250         {
2251             p[i] *= diag;
2252         }
2253
2254         /* And then go forward again */
2255         for (k = 0; k < ncorr; k++)
2256         {
2257             yr = 0;
2258             for (i = 0; i < n; i++)
2259             {
2260                 yr += p[i]*dg[cp][i];
2261             }
2262
2263             beta = rho[cp]*yr;
2264             beta = alpha[cp]-beta;
2265
2266             for (i = 0; i < n; i++)
2267             {
2268                 p[i] += beta*dx[cp][i];
2269             }
2270
2271             cp++;
2272             if (cp >= ncorr)
2273             {
2274                 cp = 0;
2275             }
2276         }
2277
2278         for (i = 0; i < n; i++)
2279         {
2280             if (!frozen[i])
2281             {
2282                 dx[point][i] = p[i];
2283             }
2284             else
2285             {
2286                 dx[point][i] = 0;
2287             }
2288         }
2289
2290         /* Print it if necessary */
2291         if (MASTER(cr))
2292         {
2293             if (mdrunOptions.verbose)
2294             {
2295                 double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
2296                 fprintf(stderr, "\rStep %d, Epot=%12.6e, Fnorm=%9.3e, Fmax=%9.3e (atom %d)\n",
2297                         step, ems.epot, ems.fnorm/sqrtNumAtoms, ems.fmax, ems.a_fmax + 1);
2298                 fflush(stderr);
2299             }
2300             /* Store the new (lower) energies */
2301             matrix nullBox = {};
2302             energyOutput.addDataAtEnergyStep(false, false, static_cast<double>(step),
2303                                              mdatoms->tmass, enerd, nullptr, nullptr, nullptr, nullBox,
2304                                              nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
2305
2306             do_log = do_per_step(step, inputrec->nstlog);
2307             do_ene = do_per_step(step, inputrec->nstenergy);
2308
2309             imdSession->fillEnergyRecord(step, TRUE);
2310
2311             if (do_log)
2312             {
2313                 print_ebin_header(fplog, step, step);
2314             }
2315             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
2316                                                do_log ? fplog : nullptr, step, step, eprNORMAL,
2317                                                fcd, &(top_global->groups), &(inputrec->opts), nullptr);
2318         }
2319
2320         /* Send x and E to IMD client, if bIMD is TRUE. */
2321         if (imdSession->run(step, TRUE, state_global->box, state_global->x.rvec_array(), 0) && MASTER(cr))
2322         {
2323             imdSession->sendPositionsAndEnergies();
2324         }
2325
2326         // Reset stepsize in we are doing more iterations
2327         stepsize = 1.0/ems.fnorm;
2328
2329         /* Stop when the maximum force lies below tolerance.
2330          * If we have reached machine precision, converged is already set to true.
2331          */
2332         converged = converged || (ems.fmax < inputrec->em_tol);
2333
2334     }   /* End of the loop */
2335
2336     if (converged)
2337     {
2338         step--; /* we never took that last step in this case */
2339
2340     }
2341     if (ems.fmax > inputrec->em_tol)
2342     {
2343         if (MASTER(cr))
2344         {
2345             warn_step(fplog, inputrec->em_tol, ems.fmax,
2346                       step-1 == number_steps, FALSE);
2347         }
2348         converged = FALSE;
2349     }
2350
2351     /* If we printed energy and/or logfile last step (which was the last step)
2352      * we don't have to do it again, but otherwise print the final values.
2353      */
2354     if (!do_log) /* Write final value to log since we didn't do anythin last step */
2355     {
2356         print_ebin_header(fplog, step, step);
2357     }
2358     if (!do_ene || !do_log) /* Write final energy file entries */
2359     {
2360         energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
2361                                            !do_log ? fplog : nullptr, step, step, eprNORMAL,
2362                                            fcd, &(top_global->groups), &(inputrec->opts), nullptr);
2363     }
2364
2365     /* Print some stuff... */
2366     if (MASTER(cr))
2367     {
2368         fprintf(stderr, "\nwriting lowest energy coordinates.\n");
2369     }
2370
2371     /* IMPORTANT!
2372      * For accurate normal mode calculation it is imperative that we
2373      * store the last conformation into the full precision binary trajectory.
2374      *
2375      * However, we should only do it if we did NOT already write this step
2376      * above (which we did if do_x or do_f was true).
2377      */
2378     do_x = !do_per_step(step, inputrec->nstxout);
2379     do_f = !do_per_step(step, inputrec->nstfout);
2380     write_em_traj(fplog, cr, outf, do_x, do_f, ftp2fn(efSTO, nfile, fnm),
2381                   top_global, inputrec, step,
2382                   &ems, state_global, observablesHistory);
2383
2384     if (MASTER(cr))
2385     {
2386         double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
2387         print_converged(stderr, LBFGS, inputrec->em_tol, step, converged,
2388                         number_steps, &ems, sqrtNumAtoms);
2389         print_converged(fplog, LBFGS, inputrec->em_tol, step, converged,
2390                         number_steps, &ems, sqrtNumAtoms);
2391
2392         fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
2393     }
2394
2395     finish_em(cr, outf, walltime_accounting, wcycle);
2396
2397     /* To print the actual number of steps we needed somewhere */
2398     walltime_accounting_set_nsteps_done(walltime_accounting, step);
2399 }
2400
2401 void
2402 Integrator::do_steep()
2403 {
2404     const char       *SD  = "Steepest Descents";
2405     gmx_localtop_t    top;
2406     gmx_global_stat_t gstat;
2407     t_graph          *graph;
2408     real              stepsize;
2409     real              ustep;
2410     gmx_bool          bDone, bAbort, do_x, do_f;
2411     tensor            vir, pres;
2412     rvec              mu_tot = {0};
2413     int               nsteps;
2414     int               count          = 0;
2415     int               steps_accepted = 0;
2416     auto              mdatoms        = mdAtoms->mdatoms();
2417
2418     GMX_LOG(mdlog.info).asParagraph().
2419         appendText("Note that activating steepest-descent energy minimization via the "
2420                    "integrator .mdp option and the command gmx mdrun may "
2421                    "be available in a different form in a future version of GROMACS, "
2422                    "e.g. gmx minimize and an .mdp option.");
2423
2424     /* Create 2 states on the stack and extract pointers that we will swap */
2425     em_state_t  s0 {}, s1 {};
2426     em_state_t *s_min = &s0;
2427     em_state_t *s_try = &s1;
2428
2429     /* Init em and store the local state in s_try */
2430     init_em(fplog, mdlog, SD, cr, inputrec, imdSession,
2431             state_global, top_global, s_try, &top,
2432             nrnb, fr, &graph, mdAtoms, &gstat,
2433             vsite, constr, nullptr);
2434     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, inputrec, top_global, nullptr, wcycle);
2435     gmx::EnergyOutput energyOutput;
2436     energyOutput.prepare(mdoutf_get_fp_ene(outf), top_global, inputrec, nullptr);
2437
2438     /* Print to log file  */
2439     print_em_start(fplog, cr, walltime_accounting, wcycle, SD);
2440
2441     /* Set variables for stepsize (in nm). This is the largest
2442      * step that we are going to make in any direction.
2443      */
2444     ustep    = inputrec->em_stepsize;
2445     stepsize = 0;
2446
2447     /* Max number of steps  */
2448     nsteps = inputrec->nsteps;
2449
2450     if (MASTER(cr))
2451     {
2452         /* Print to the screen  */
2453         sp_header(stderr, SD, inputrec->em_tol, nsteps);
2454     }
2455     if (fplog)
2456     {
2457         sp_header(fplog, SD, inputrec->em_tol, nsteps);
2458     }
2459     EnergyEvaluator energyEvaluator {
2460         fplog, mdlog, cr, ms,
2461         top_global, &top,
2462         inputrec, imdSession, nrnb, wcycle, gstat,
2463         vsite, constr, fcd, graph,
2464         mdAtoms, fr, ppForceWorkload, enerd
2465     };
2466
2467     /**** HERE STARTS THE LOOP ****
2468      * count is the counter for the number of steps
2469      * bDone will be TRUE when the minimization has converged
2470      * bAbort will be TRUE when nsteps steps have been performed or when
2471      * the stepsize becomes smaller than is reasonable for machine precision
2472      */
2473     count  = 0;
2474     bDone  = FALSE;
2475     bAbort = FALSE;
2476     while (!bDone && !bAbort)
2477     {
2478         bAbort = (nsteps >= 0) && (count == nsteps);
2479
2480         /* set new coordinates, except for first step */
2481         bool validStep = true;
2482         if (count > 0)
2483         {
2484             validStep =
2485                 do_em_step(cr, inputrec, mdatoms,
2486                            s_min, stepsize, &s_min->f, s_try,
2487                            constr, count);
2488         }
2489
2490         if (validStep)
2491         {
2492             energyEvaluator.run(s_try, mu_tot, vir, pres, count, count == 0);
2493         }
2494         else
2495         {
2496             // Signal constraint error during stepping with energy=inf
2497             s_try->epot = std::numeric_limits<real>::infinity();
2498         }
2499
2500         if (MASTER(cr))
2501         {
2502             print_ebin_header(fplog, count, count);
2503         }
2504
2505         if (count == 0)
2506         {
2507             s_min->epot = s_try->epot;
2508         }
2509
2510         /* Print it if necessary  */
2511         if (MASTER(cr))
2512         {
2513             if (mdrunOptions.verbose)
2514             {
2515                 fprintf(stderr, "Step=%5d, Dmax= %6.1e nm, Epot= %12.5e Fmax= %11.5e, atom= %d%c",
2516                         count, ustep, s_try->epot, s_try->fmax, s_try->a_fmax+1,
2517                         ( (count == 0) || (s_try->epot < s_min->epot) ) ? '\n' : '\r');
2518                 fflush(stderr);
2519             }
2520
2521             if ( (count == 0) || (s_try->epot < s_min->epot) )
2522             {
2523                 /* Store the new (lower) energies  */
2524                 matrix nullBox = {};
2525                 energyOutput.addDataAtEnergyStep(false, false, static_cast<double>(count),
2526                                                  mdatoms->tmass, enerd, nullptr, nullptr, nullptr, nullBox,
2527                                                  nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
2528
2529                 imdSession->fillEnergyRecord(count, TRUE);
2530
2531                 const bool do_dr = do_per_step(steps_accepted, inputrec->nstdisreout);
2532                 const bool do_or = do_per_step(steps_accepted, inputrec->nstorireout);
2533                 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE,
2534                                                    do_dr, do_or,
2535                                                    fplog, count, count, eprNORMAL,
2536                                                    fcd, &(top_global->groups),
2537                                                    &(inputrec->opts), nullptr);
2538                 fflush(fplog);
2539             }
2540         }
2541
2542         /* Now if the new energy is smaller than the previous...
2543          * or if this is the first step!
2544          * or if we did random steps!
2545          */
2546
2547         if ( (count == 0) || (s_try->epot < s_min->epot) )
2548         {
2549             steps_accepted++;
2550
2551             /* Test whether the convergence criterion is met...  */
2552             bDone = (s_try->fmax < inputrec->em_tol);
2553
2554             /* Copy the arrays for force, positions and energy  */
2555             /* The 'Min' array always holds the coords and forces of the minimal
2556                sampled energy  */
2557             swap_em_state(&s_min, &s_try);
2558             if (count > 0)
2559             {
2560                 ustep *= 1.2;
2561             }
2562
2563             /* Write to trn, if necessary */
2564             do_x = do_per_step(steps_accepted, inputrec->nstxout);
2565             do_f = do_per_step(steps_accepted, inputrec->nstfout);
2566             write_em_traj(fplog, cr, outf, do_x, do_f, nullptr,
2567                           top_global, inputrec, count,
2568                           s_min, state_global, observablesHistory);
2569         }
2570         else
2571         {
2572             /* If energy is not smaller make the step smaller...  */
2573             ustep *= 0.5;
2574
2575             if (DOMAINDECOMP(cr) && s_min->s.ddp_count != cr->dd->ddp_count)
2576             {
2577                 /* Reload the old state */
2578                 em_dd_partition_system(fplog, mdlog, count, cr, top_global, inputrec, imdSession,
2579                                        s_min, &top, mdAtoms, fr, vsite, constr,
2580                                        nrnb, wcycle);
2581             }
2582         }
2583
2584         /* Determine new step  */
2585         stepsize = ustep/s_min->fmax;
2586
2587         /* Check if stepsize is too small, with 1 nm as a characteristic length */
2588 #if GMX_DOUBLE
2589         if (count == nsteps || ustep < 1e-12)
2590 #else
2591         if (count == nsteps || ustep < 1e-6)
2592 #endif
2593         {
2594             if (MASTER(cr))
2595             {
2596                 warn_step(fplog, inputrec->em_tol, s_min->fmax,
2597                           count == nsteps, constr != nullptr);
2598             }
2599             bAbort = TRUE;
2600         }
2601
2602         /* Send IMD energies and positions, if bIMD is TRUE. */
2603         if (imdSession->run(count, TRUE, state_global->box,
2604                             MASTER(cr) ? state_global->x.rvec_array() : nullptr,
2605                             0) &&
2606             MASTER(cr))
2607         {
2608             imdSession->sendPositionsAndEnergies();
2609         }
2610
2611         count++;
2612     }   /* End of the loop  */
2613
2614     /* Print some data...  */
2615     if (MASTER(cr))
2616     {
2617         fprintf(stderr, "\nwriting lowest energy coordinates.\n");
2618     }
2619     write_em_traj(fplog, cr, outf, TRUE, inputrec->nstfout != 0, ftp2fn(efSTO, nfile, fnm),
2620                   top_global, inputrec, count,
2621                   s_min, state_global, observablesHistory);
2622
2623     if (MASTER(cr))
2624     {
2625         double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
2626
2627         print_converged(stderr, SD, inputrec->em_tol, count, bDone, nsteps,
2628                         s_min, sqrtNumAtoms);
2629         print_converged(fplog, SD, inputrec->em_tol, count, bDone, nsteps,
2630                         s_min, sqrtNumAtoms);
2631     }
2632
2633     finish_em(cr, outf, walltime_accounting, wcycle);
2634
2635     /* To print the actual number of steps we needed somewhere */
2636     inputrec->nsteps = count;
2637
2638     walltime_accounting_set_nsteps_done(walltime_accounting, count);
2639 }
2640
2641 void
2642 Integrator::do_nm()
2643 {
2644     const char          *NM = "Normal Mode Analysis";
2645     int                  nnodes, node;
2646     gmx_localtop_t       top;
2647     gmx_global_stat_t    gstat;
2648     t_graph             *graph;
2649     tensor               vir, pres;
2650     rvec                 mu_tot = {0};
2651     rvec                *dfdx;
2652     gmx_bool             bSparse; /* use sparse matrix storage format */
2653     size_t               sz;
2654     gmx_sparsematrix_t * sparse_matrix           = nullptr;
2655     real           *     full_matrix             = nullptr;
2656
2657     /* added with respect to mdrun */
2658     int                       row, col;
2659     real                      der_range = 10.0*std::sqrt(GMX_REAL_EPS);
2660     real                      x_min;
2661     bool                      bIsMaster = MASTER(cr);
2662     auto                      mdatoms   = mdAtoms->mdatoms();
2663
2664     GMX_LOG(mdlog.info).asParagraph().
2665         appendText("Note that activating normal-mode analysis via the integrator "
2666                    ".mdp option and the command gmx mdrun may "
2667                    "be available in a different form in a future version of GROMACS, "
2668                    "e.g. gmx normal-modes.");
2669
2670     if (constr != nullptr)
2671     {
2672         gmx_fatal(FARGS, "Constraints present with Normal Mode Analysis, this combination is not supported");
2673     }
2674
2675     gmx_shellfc_t *shellfc;
2676
2677     em_state_t     state_work {};
2678
2679     /* Init em and store the local state in state_minimum */
2680     init_em(fplog, mdlog, NM, cr, inputrec, imdSession,
2681             state_global, top_global, &state_work, &top,
2682             nrnb, fr, &graph, mdAtoms, &gstat,
2683             vsite, constr, &shellfc);
2684     gmx_mdoutf            *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, inputrec, top_global, nullptr, wcycle);
2685
2686     std::vector<int>       atom_index = get_atom_index(top_global);
2687     std::vector<gmx::RVec> fneg(atom_index.size(), {0, 0, 0});
2688     snew(dfdx, atom_index.size());
2689
2690 #if !GMX_DOUBLE
2691     if (bIsMaster)
2692     {
2693         fprintf(stderr,
2694                 "NOTE: This version of GROMACS has been compiled in single precision,\n"
2695                 "      which MIGHT not be accurate enough for normal mode analysis.\n"
2696                 "      GROMACS now uses sparse matrix storage, so the memory requirements\n"
2697                 "      are fairly modest even if you recompile in double precision.\n\n");
2698     }
2699 #endif
2700
2701     /* Check if we can/should use sparse storage format.
2702      *
2703      * Sparse format is only useful when the Hessian itself is sparse, which it
2704      * will be when we use a cutoff.
2705      * For small systems (n<1000) it is easier to always use full matrix format, though.
2706      */
2707     if (EEL_FULL(fr->ic->eeltype) || fr->rlist == 0.0)
2708     {
2709         GMX_LOG(mdlog.warning).appendText("Non-cutoff electrostatics used, forcing full Hessian format.");
2710         bSparse = FALSE;
2711     }
2712     else if (atom_index.size() < 1000)
2713     {
2714         GMX_LOG(mdlog.warning).appendTextFormatted("Small system size (N=%zu), using full Hessian format.",
2715                                                    atom_index.size());
2716         bSparse = FALSE;
2717     }
2718     else
2719     {
2720         GMX_LOG(mdlog.warning).appendText("Using compressed symmetric sparse Hessian format.");
2721         bSparse = TRUE;
2722     }
2723
2724     /* Number of dimensions, based on real atoms, that is not vsites or shell */
2725     sz = DIM*atom_index.size();
2726
2727     fprintf(stderr, "Allocating Hessian memory...\n\n");
2728
2729     if (bSparse)
2730     {
2731         sparse_matrix = gmx_sparsematrix_init(sz);
2732         sparse_matrix->compressed_symmetric = TRUE;
2733     }
2734     else
2735     {
2736         snew(full_matrix, sz*sz);
2737     }
2738
2739     init_nrnb(nrnb);
2740
2741
2742     /* Write start time and temperature */
2743     print_em_start(fplog, cr, walltime_accounting, wcycle, NM);
2744
2745     /* fudge nr of steps to nr of atoms */
2746     inputrec->nsteps = atom_index.size()*2;
2747
2748     if (bIsMaster)
2749     {
2750         fprintf(stderr, "starting normal mode calculation '%s'\n%" PRId64 " steps.\n\n",
2751                 *(top_global->name), inputrec->nsteps);
2752     }
2753
2754     nnodes = cr->nnodes;
2755
2756     /* Make evaluate_energy do a single node force calculation */
2757     cr->nnodes = 1;
2758     EnergyEvaluator energyEvaluator {
2759         fplog, mdlog, cr, ms,
2760         top_global, &top,
2761         inputrec, imdSession, nrnb, wcycle, gstat,
2762         vsite, constr, fcd, graph,
2763         mdAtoms, fr, ppForceWorkload, enerd
2764     };
2765     energyEvaluator.run(&state_work, mu_tot, vir, pres, -1, TRUE);
2766     cr->nnodes = nnodes;
2767
2768     /* if forces are not small, warn user */
2769     get_state_f_norm_max(cr, &(inputrec->opts), mdatoms, &state_work);
2770
2771     GMX_LOG(mdlog.warning).appendTextFormatted("Maximum force:%12.5e", state_work.fmax);
2772     if (state_work.fmax > 1.0e-3)
2773     {
2774         GMX_LOG(mdlog.warning).appendText(
2775                 "The force is probably not small enough to "
2776                 "ensure that you are at a minimum.\n"
2777                 "Be aware that negative eigenvalues may occur\n"
2778                 "when the resulting matrix is diagonalized.");
2779     }
2780
2781     /***********************************************************
2782      *
2783      *      Loop over all pairs in matrix
2784      *
2785      *      do_force called twice. Once with positive and
2786      *      once with negative displacement
2787      *
2788      ************************************************************/
2789
2790     /* Steps are divided one by one over the nodes */
2791     bool bNS          = true;
2792     auto state_work_x = makeArrayRef(state_work.s.x);
2793     auto state_work_f = makeArrayRef(state_work.f);
2794     for (unsigned int aid = cr->nodeid; aid < atom_index.size(); aid += nnodes)
2795     {
2796         size_t atom = atom_index[aid];
2797         for (size_t d = 0; d < DIM; d++)
2798         {
2799             int64_t     step        = 0;
2800             int         force_flags = GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES;
2801             double      t           = 0;
2802
2803             x_min = state_work_x[atom][d];
2804
2805             for (unsigned int dx = 0; (dx < 2); dx++)
2806             {
2807                 if (dx == 0)
2808                 {
2809                     state_work_x[atom][d] = x_min - der_range;
2810                 }
2811                 else
2812                 {
2813                     state_work_x[atom][d] = x_min + der_range;
2814                 }
2815
2816                 /* Make evaluate_energy do a single node force calculation */
2817                 cr->nnodes = 1;
2818                 if (shellfc)
2819                 {
2820                     /* Now is the time to relax the shells */
2821                     relax_shell_flexcon(fplog,
2822                                         cr,
2823                                         ms,
2824                                         mdrunOptions.verbose,
2825                                         nullptr,
2826                                         step,
2827                                         inputrec,
2828                                         imdSession,
2829                                         bNS,
2830                                         force_flags,
2831                                         &top,
2832                                         constr,
2833                                         enerd,
2834                                         fcd,
2835                                         &state_work.s,
2836                                         state_work.f.arrayRefWithPadding(),
2837                                         vir,
2838                                         mdatoms,
2839                                         nrnb,
2840                                         wcycle,
2841                                         graph,
2842                                         shellfc,
2843                                         fr,
2844                                         ppForceWorkload,
2845                                         t,
2846                                         mu_tot,
2847                                         vsite,
2848                                         DDBalanceRegionHandler(nullptr));
2849                     bNS = false;
2850                     step++;
2851                 }
2852                 else
2853                 {
2854                     energyEvaluator.run(&state_work, mu_tot, vir, pres, aid*2+dx, FALSE);
2855                 }
2856
2857                 cr->nnodes = nnodes;
2858
2859                 if (dx == 0)
2860                 {
2861                     std::copy(state_work_f.begin(), state_work_f.begin()+atom_index.size(), fneg.begin());
2862                 }
2863             }
2864
2865             /* x is restored to original */
2866             state_work_x[atom][d] = x_min;
2867
2868             for (size_t j = 0; j < atom_index.size(); j++)
2869             {
2870                 for (size_t k = 0; (k < DIM); k++)
2871                 {
2872                     dfdx[j][k] =
2873                         -(state_work_f[atom_index[j]][k] - fneg[j][k])/(2*der_range);
2874                 }
2875             }
2876
2877             if (!bIsMaster)
2878             {
2879 #if GMX_MPI
2880 #define mpi_type GMX_MPI_REAL
2881                 MPI_Send(dfdx[0], atom_index.size()*DIM, mpi_type, MASTER(cr),
2882                          cr->nodeid, cr->mpi_comm_mygroup);
2883 #endif
2884             }
2885             else
2886             {
2887                 for (node = 0; (node < nnodes && aid+node < atom_index.size()); node++)
2888                 {
2889                     if (node > 0)
2890                     {
2891 #if GMX_MPI
2892                         MPI_Status stat;
2893                         MPI_Recv(dfdx[0], atom_index.size()*DIM, mpi_type, node, node,
2894                                  cr->mpi_comm_mygroup, &stat);
2895 #undef mpi_type
2896 #endif
2897                     }
2898
2899                     row = (aid + node)*DIM + d;
2900
2901                     for (size_t j = 0; j < atom_index.size(); j++)
2902                     {
2903                         for (size_t k = 0; k < DIM; k++)
2904                         {
2905                             col = j*DIM + k;
2906
2907                             if (bSparse)
2908                             {
2909                                 if (col >= row && dfdx[j][k] != 0.0)
2910                                 {
2911                                     gmx_sparsematrix_increment_value(sparse_matrix,
2912                                                                      row, col, dfdx[j][k]);
2913                                 }
2914                             }
2915                             else
2916                             {
2917                                 full_matrix[row*sz+col] = dfdx[j][k];
2918                             }
2919                         }
2920                     }
2921                 }
2922             }
2923
2924             if (mdrunOptions.verbose && fplog)
2925             {
2926                 fflush(fplog);
2927             }
2928         }
2929         /* write progress */
2930         if (bIsMaster && mdrunOptions.verbose)
2931         {
2932             fprintf(stderr, "\rFinished step %d out of %td",
2933                     std::min<int>(atom+nnodes, atom_index.size()),
2934                     ssize(atom_index));
2935             fflush(stderr);
2936         }
2937     }
2938
2939     if (bIsMaster)
2940     {
2941         fprintf(stderr, "\n\nWriting Hessian...\n");
2942         gmx_mtxio_write(ftp2fn(efMTX, nfile, fnm), sz, sz, full_matrix, sparse_matrix);
2943     }
2944
2945     finish_em(cr, outf, walltime_accounting, wcycle);
2946
2947     walltime_accounting_set_nsteps_done(walltime_accounting, atom_index.size()*2);
2948 }
2949
2950 } // namespace gmx