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