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