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