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