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