src/kernel/steep.c

   1 /*
   2  *       $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1997
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://rugmd0.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_steep_c = "$Id$";
  30
  31 #define FORCE_CRIT
  32
  33 #include <string.h>
  34 #include <time.h>
  35 #include <math.h>
  36 #include "sysstuff.h"
  37 #include "string2.h"
  38 #include "led.h"
  39 #include "network.h"
  40 #include "confio.h"
  41 #include "copyrite.h"
  42 #include "smalloc.h"
  43 #include "nrnb.h"
  44 #include "main.h"
  45 #include "pbc.h"
  46 #include "force.h"
  47 #include "macros.h"
  48 #include "random.h"
  49 #include "names.h"
  50 #include "fatal.h"
  51 #include "txtdump.h"
  52 #include "typedefs.h"
  53 #include "update.h"
  54 #include "random.h"
  55 #include "vec.h"
  56 #include "enxio.h"
  57 #include "tgroup.h"
  58 #include "mdebin.h"
  59 #include "mdrun.h"
  60 #include "pppm.h"
  61 #include "dummies.h"
  62
  63 #ifdef FORCE_CRIT
  64 static void sp_header(FILE *out,real epot,real fsqrt,real step,real ftol)
  65 {
  66   fprintf(out,"STEEPEST DESCENTS:\n");
  67   fprintf(out,"   Stepsize          = %12.5e\n",step);
  68   fprintf(out,"   Tolerance         = %12.5e\n",ftol);
  69   fprintf(out,"   Starting rmsF     = %30.20e\n",fsqrt);
  70   fprintf(out,"   Starting Energy   = %30.20e\n",epot);
  71 }
  72 #else
  73 static void sp_header(FILE *out,real epot,real step,real ftol)
  74 {
  75   fprintf(out,"STEEPEST DESCENTS:\n");
  76   fprintf(out,"   Stepsize          = %12.5e\n",step);
  77   fprintf(out,"   Tolerance         = %12.5e\n",ftol);
  78   fprintf(out,"   Starting Energy   = %30.20e\n",epot);
  79 }
  80 #endif
  81
  82 real f_max(FILE *log,
  83            int left,int right,int nprocs,
  84            int start,int end,rvec grad[])
  85 {
  86   real fmax,fmax_0,fam;
  87   int  i,m;
  88
  89   /* This routine finds the largest force component (abs value)
  90    * and returns it. On parallel machines the global max is taken.
  91    */
  92   fmax = fabs(grad[start][0]);
  93
  94   for(i=start; (i<end); i++)
  95     for(m=0; (m<DIM); m++) {
  96       fam=fabs(grad[i][m]);
  97       fmax=max(fmax,fam);
  98     }
  99   if (nprocs == 1)
 100     return fmax;
 101
 102   for(i=0; (i<nprocs-1); i++) {
 103     gmx_tx(left,(void *)&fmax,sizeof(fmax));
 104     gmx_rx(right,(void *)&fmax_0,sizeof(fmax_0));
 105     gmx_wait(left,right);
 106     if (fmax_0 > fmax)
 107       fmax=fmax_0;
 108   }
 109
 110   return fmax;
 111 }
 112
 113 real f_norm(FILE *log,
 114            int left,int right,int nprocs,
 115            int start,int end,rvec grad[])
 116 {
 117   real fnorm;
 118   int  i,m;
 119
 120   /* This routine finds the norm of the force
 121    * and returns it. Parallel machines not supported.
 122    */
 123   fnorm = 0;
 124
 125   for(i=start; (i<end); i++)
 126     for(m=0; (m<DIM); m++) {
 127       fnorm=fnorm+grad[i][m]*grad[i][m];
 128     }
 129   fnorm=sqrt(fnorm);
 130
 131   if (nprocs == 1)
 132     return fnorm;
 133
 134   fatal_error(0,"This version of Steepest Descents cannot be run in parallel");
 135   return fnorm;
 136 }
 137
 138 real f_sqrt(FILE *log,
 139             int left,int right,int nprocs,
 140             int start,int end,rvec f[])
 141 {
 142   int i,m;
 143   real fsqr;
 144
 145
 146   /* this routine calculates the route mean square force  */
 147
 148   fsqr=0;
 149
 150   /* calculate the sum of the square force of this processor  */
 151   for(i=start;(i<end);i++)
 152     for(m=0;(m<DIM);m++)
 153       fsqr += sqr(f[i][m]);
 154
 155
 156   if (nprocs == 1)
 157     return sqrt(fsqr)/(end - start);
 158
 159
 160   fatal_error(0,"This version of f_sqrt cannot run in parrallel");
 161
 162
 163   return ( sqrt(fsqr) );
 164 }
 165
 166 static void do_step(int start,int end,rvec x[],rvec f[],
 167                     bool bRandom,real step)
 168 {
 169   static int seed=1993;
 170   int    i,m;
 171   real   r;
 172
 173   if (bRandom) {
 174     fprintf(stderr,"\rRandom step\n");
 175     for(i=start; (i<end); i++) {
 176       for(m=0; (m<DIM); m++) {
 177         r=rando(&seed);
 178         x[i][m]=x[i][m]+step*r;
 179       }
 180     }
 181   }
 182   else {
 183     /* New positions to try  */
 184     for(i=start; (i<end);i++)
 185       for(m=0; (m<DIM); m++)
 186         x[i][m] = x[i][m]+step*f[i][m];
 187   }
 188 }
 189
 190 time_t do_steep(FILE *log,int nfile,t_filenm fnm[],
 191                 t_parm *parm,t_topology *top,
 192                 t_groups *grps,t_nsborder *nsb,
 193                 rvec x[],rvec grad[],rvec buf[],t_mdatoms *mdatoms,
 194                 tensor ekin,real ener[],t_nrnb nrnb[],
 195                 bool bVerbose,bool bDummies,t_commrec *cr,t_graph *graph,
 196                 t_forcerec *fr,rvec box_size)
 197 {
 198   static char *SD="STEEPEST DESCENTS";
 199   real   step,lambda,ftol,fmax;
 200   rvec   *pos[2],*force[2];
 201   rvec   *xx,*ff;
 202 #ifdef FORCE_CRIT
 203   real   Fsqrt[2];
 204 #endif
 205   real   Epot[2];
 206   real   vcm[4],fnorm,ustep;
 207   int        fp_ene;
 208   t_mdebin   *mdebin;
 209   t_nrnb mynrnb;
 210   bool   bNS=TRUE,bDone,bLR,bLJLR,bBHAM,b14;
 211   time_t start_t;
 212   tensor force_vir,shake_vir;
 213   rvec   mu_tot;
 214   int    number_steps;
 215   int    count=0;
 216   int    i,m,start,end;
 217   int    Min=0;
 218   int    steps_accepted=0;
 219 #define  TRY (1-Min)
 220
 221   /* Initiate some variables  */
 222   if (parm->ir.bPert)
 223     lambda       = parm->ir.init_lambda;
 224   else
 225     lambda = 0.0;
 226
 227   clear_rvec(mu_tot);
 228   calc_shifts(parm->box,box_size,fr->shift_vec,FALSE);
 229
 230   vcm[0]=vcm[1]=vcm[2]=vcm[3]=0.0;
 231
 232   /* Print to log file  */
 233   start_t=print_date_and_time(log,cr->pid,"Started Steepest Descents");
 234
 235   /* We need two coordinate arrays and two force arrays  */
 236   for(i=0; (i<2); i++) {
 237     snew(pos[i],nsb->natoms);
 238     snew(force[i],nsb->natoms);
 239   }
 240
 241   start=nsb->index[cr->pid];
 242   end=nsb->homenr[cr->pid]-start;
 243
 244   /* Open the enrgy file */
 245   if (MASTER(cr))
 246     fp_ene=open_enx(ftp2fn(efENX,nfile,fnm),"w");
 247   else
 248     fp_ene=-1;
 249
 250   /* Set some booleans for the epot routines  */
 251   set_pot_bools(&(parm->ir),top,&bLR,&bLJLR,&bBHAM,&b14);
 252
 253   /* Init bin for energy stuff  */
 254   mdebin=init_mdebin(fp_ene,grps,&(top->atoms),&(top->idef),bLR,bLJLR,
 255                      bBHAM,b14,parm->ir.bPert,parm->ir.epc,
 256                      parm->ir.bDispCorr);
 257
 258   /* Clear some matrix variables  */
 259   clear_mat(force_vir);
 260   clear_mat(shake_vir);
 261
 262   /* Copy coord vectors to our temp array  */
 263   for(i=0; (i<nsb->natoms); i++) {
 264     copy_rvec(x[i],pos[Min][i]);
 265     copy_rvec(x[i],pos[TRY][i]);
 266   }
 267
 268   /* Set variables for stepsize (in nm). This is the largest
 269    * step that we are going to make in any direction.
 270    */
 271   step=ustep=parm->ir.em_stepsize;
 272
 273   /* Tolerance for conversion  */
 274   ftol=parm->ir.em_tol;
 275
 276   /* Max number of steps  */
 277   number_steps=parm->ir.nsteps;
 278
 279   if (MASTER(cr)) {
 280     /* Print to the screen  */
 281     start_t=print_date_and_time(log,cr->pid,"Started EM");
 282     fprintf(stderr,"STEEPEST DESCENTS:\n");
 283     fprintf(log,"STEEPEST DESCENTS:\n");
 284     fprintf(stderr,"   Tolerance         = %12.5g\n",ftol);
 285     fprintf(log,"   Tolerance         = %12.5g\n",ftol);
 286   }
 287
 288   /* Here starts the loop, count is the counter for the number of steps
 289    * bDone is a BOOLEAN variable, which will be set TRUE when
 290    * the minimization has converged.
 291    */
 292   for(count=0,bDone=FALSE;
 293       !(bDone || ((number_steps > 0) && (count==number_steps))); count++) {
 294
 295     /* set new coordinates, except for first step */
 296     if (count>0)
 297       for(i=start; (i<end); i++)
 298         for(m=0;(m<DIM);m++)
 299           pos[TRY][i][m] = pos[Min][i][m] + step * force[Min][i][m];
 300
 301     if (bDummies) {
 302       /* Construct dummy particles */
 303       shift_self(graph,fr->shift_vec,pos[TRY]);
 304       construct_dummies(log,pos[TRY],&(nrnb[cr->pid]),1,NULL,&top->idef);
 305       unshift_self(graph,fr->shift_vec,pos[TRY]);
 306     }
 307
 308     /* Calc force & energy on new positions  */
 309     do_force(log,cr,parm,nsb,force_vir,
 310              count,&(nrnb[cr->pid]),top,grps,pos[TRY],buf,force[TRY],buf,
 311              mdatoms,ener,bVerbose && !(PAR(cr)),
 312              lambda,graph,bNS,FALSE,fr);
 313     unshift_self(graph,fr->shift_vec,pos[TRY]);
 314
 315     /* Spread the force on dummy particle to the other particles... */
 316     spread_dummy_f(log,pos[TRY],force[TRY],&(nrnb[cr->pid]),&top->idef);
 317
 318     /* Sum the potential energy terms from group contributions  */
 319     sum_epot(&(parm->ir.opts),grps,ener);
 320
 321     /* Clear stuff again  */
 322     clear_mat(force_vir);
 323     clear_mat(shake_vir);
 324
 325     /* Communicat stuff when parallel  */
 326     if (PAR(cr))
 327       global_stat(log,cr,ener,force_vir,shake_vir,
 328                   &(parm->ir.opts),grps,&mynrnb,nrnb,vcm,mu_tot);
 329
 330     /* This is the new energy  */
 331 #ifdef FORCE_CRIT
 332     Fsqrt[TRY]=f_sqrt(log,cr->left,cr->right,nsb->nprocs,start,end,force[TRY]);
 333 #endif
 334     Epot[TRY]=ener[F_EPOT];
 335
 336     /* Print it if necessary  */
 337 #ifdef FORCE_CRIT
 338     if ((bVerbose || (Fsqrt[TRY] < Fsqrt[Min])) && MASTER(cr)) {
 339       fprintf(stderr,
 340               "\rStep = %5d, Dx = %12.5e, Epot = %12.5e rmsF = %12.5e\n",
 341               count,step,Epot[TRY],Fsqrt[TRY]);
 342       /* Store the new (lower) energies  */
 343       upd_mdebin(mdebin,mdatoms->tmass,count,ener,parm->box,shake_vir,
 344                  force_vir,parm->vir,parm->pres,grps,mu_tot);
 345       /* Print the energies allways when we should be verbose  */
 346       if (MASTER(cr))
 347         print_ebin(fp_ene,TRUE,FALSE,log,count,count,lambda,
 348                    0.0,eprNORMAL,TRUE,mdebin,grps,&(top->atoms));
 349     }
 350 #else
 351     if ((bVerbose || (Epot[TRY] < Epot[Min])) && MASTER(cr)) {
 352       fprintf(stderr,"\rStep = %5d, Dx = %12.5e, E-Pot = %30.20e\n",
 353               count,step,Epot[TRY]);
 354
 355       /* Store the new (lower) energies  */
 356       upd_mdebin(mdebin,mdatoms->tmass,count,ener,parm->box,shake_vir,
 357                  force_vir,parm->vir,parm->pres,grps);
 358       /* Print the energies allways when we should be verbose  */
 359       if (MASTER(cr))
 360         print_ebin(fp_ene,log,count,count,lambda,0.0,eprNORMAL,TRUE,
 361                    mdebin,grps,&(top->atoms));
 362     }
 363 #endif
 364
 365     /* Now if the new energy is smaller than the previous...
 366      * or if this is the first step!
 367      * or if we did random steps!
 368      */
 369
 370 #ifdef FORCE_CRIT
 371     if ( (count==0) || (Fsqrt[TRY] < Fsqrt[Min]) ) {
 372 #else
 373     if ( (count==0) || (Epot[TRY] < Epot[Min]) ) {
 374 #endif
 375       steps_accepted++;
 376       if (do_per_step(steps_accepted,parm->ir.nstfout))
 377         ff=force[TRY];
 378       else
 379         ff=NULL;
 380       if (do_per_step(steps_accepted,parm->ir.nstxout)) {
 381         xx=pos[TRY];
 382         write_traj(log,cr,ftp2fn(efTRN,nfile,fnm),
 383                    nsb,count,(real) count,
 384                    lambda,nrnb,nsb->natoms,xx,NULL,ff,parm->box);
 385       } else
 386         xx=NULL;
 387
 388       /* Test whether the convergence criterion is met...  */
 389 #ifdef FORCE_CRIT
 390       bDone=(Fsqrt[TRY] < ftol);
 391 #else
 392       /* Stop when the difference between two tries is less than
 393        * the tolerance times the energy.
 394        */
 395       bDone=(fabs(Epot[Min]-Epot[TRY]) < ftol*fabs(Epot[Min]));
 396 #endif
 397
 398       /* Copy the arrays for force, positions and energy  */
 399       /* The 'Min' array always holds the coords and forces of the minimal
 400          sampled energy  */
 401       Min = TRY;
 402
 403       /* increase step size  */
 404       if (count>0)
 405         ustep *= 1.2;
 406
 407       /*if (MASTER(cr))
 408         fprintf(stderr,"\n"); */
 409     }
 410     else {
 411       /* If energy is not smaller make the step smaller...  */
 412       if (ustep > 1.0e-14)
 413         ustep *= 0.5;
 414       else
 415         ustep = 1.0e-14;
 416     }
 417
 418     /* Determine new step  */
 419     fnorm = f_norm(log,cr->left,cr->right,nsb->nprocs,start,end,force[Min]);
 420     step=ustep/fnorm;
 421
 422   } /* End of the loop  */
 423
 424   /* Print some shit...  */
 425   if (MASTER(cr)) {
 426     fprintf(stderr,"\nwriting lowest energy coordinates.\n");
 427     xx=pos[Min];
 428     ff=force[Min];
 429     write_traj(log,cr,ftp2fn(efTRN,nfile,fnm),
 430                nsb,count,(real) count,
 431                lambda,nrnb,nsb->natoms,xx,NULL,ff,parm->box);
 432     write_sto_conf(ftp2fn(efSTO,nfile,fnm),
 433                    *top->name, &(top->atoms),xx,NULL,parm->box);
 434
 435     fmax=f_max(log,cr->left,cr->right,nsb->nprocs,start,end,force[Min]);
 436     fprintf(stderr,"Maximum force: %12.5e\n",fmax);
 437     if (bDone) {
 438       fprintf(stderr,"\n%s converged to %8.6f \n",SD,ftol);
 439       fprintf(log,"%s converged to %8.6f \n",SD,ftol);
 440     }
 441     else {
 442       fprintf(stderr,"\n%s did not converge in %d steps\n",SD,number_steps);
 443       fprintf(log,"%s did not converge in %d steps\n",SD,number_steps);
 444     }
 445 #ifdef FORCE_CRIT
 446     fprintf(stderr,"  Minimum Root Mean Square Force  = %12.4e\n",Fsqrt[Min]);
 447     fprintf(log,"  Minimum Root Mean Square Force = %12.4e\n",Fsqrt[Min]);
 448 #endif
 449     fprintf(stderr,"  Function value at minimum = %12.4e\n",Epot[Min]);
 450     fprintf(log,"  Function value at minimum = %12.4e\n",Epot[Min]);
 451   }
 452   if (MASTER(cr))
 453     close_enx(fp_ene);
 454
 455   /* Put the coordinates back in the x array (otherwise the whole
 456    * minimization would be in vain)
 457    */
 458   for(i=0; (i<nsb->natoms); i++)
 459     copy_rvec(pos[Min][i],x[i]);
 460
 461   /* To print the actual number of steps we needed somewhere */
 462   parm->ir.nsteps=count;
 463
 464   return start_t;
 465 } /* That's all folks */
 466