src/gromacs/mdlib/minimize.c

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