src/gromacs/mdlib/minimize.c

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