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