src/contrib/relax_sh.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-1999
  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://md.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_relax_sh_c = "$Id$";
  30
  31 #include <string.h>
  32 #include "assert.h"
  33 #include "typedefs.h"
  34 #include "smalloc.h"
  35 #include "fatal.h"
  36 #include "vec.h"
  37 #include "txtdump.h"
  38 #include "mdrun.h"
  39 #include "init_sh.h"
  40 #include "mdatoms.h"
  41 #include "network.h"
  42 #include "do_gct.h"
  43 #include "names.h"
  44 #include "constr.h"
  45
  46 static void do_1pos(rvec xnew,rvec xold,rvec f,real k_1,real step)
  47 {
  48   real xo,yo,zo;
  49   real dx,dy,dz,dx2;
  50
  51   xo=xold[XX];
  52   yo=xold[YY];
  53   zo=xold[ZZ];
  54
  55   dx=f[XX]*k_1;
  56   dy=f[YY]*k_1;
  57   dz=f[ZZ]*k_1;
  58
  59   xnew[XX]=xo+dx*step;
  60   xnew[YY]=yo+dy*step;
  61   xnew[ZZ]=zo+dz*step;
  62 }
  63
  64 static void directional_sd(FILE *log,real step,rvec xold[],rvec xnew[],
  65                            rvec min_f[],int start,int homenr,real k)
  66 {
  67   int  i;
  68
  69   for(i=start; i<homenr; i++)
  70     do_1pos(xnew[i],xold[i],min_f[i],-k,step);
  71 }
  72
  73 static void shell_pos_sd(FILE *log,real step,rvec xold[],rvec xnew[],rvec f[],
  74                          int ns,t_shell s[])
  75 {
  76   int  i,shell;
  77   real k_1;
  78
  79   for(i=0; (i<ns); i++) {
  80     shell = s[i].shell;
  81     k_1   = s[i].k_1;
  82     do_1pos(xnew[shell],xold[shell],f[shell],k_1,step);
  83     if (debug && 0) {
  84       pr_rvec(debug,0,"fshell",f[shell],DIM);
  85       pr_rvec(debug,0,"xold",xold[shell],DIM);
  86       pr_rvec(debug,0,"xnew",xnew[shell],DIM);
  87     }
  88   }
  89 }
  90
  91 static void zero_shell_forces(FILE *log,rvec f[],int ns,t_shell s[])
  92 {
  93   int i,s1,n1,n2,n3;
  94   real m1,m2,m3,tm;
  95
  96   /* Subtract remaining forces on shells from the attaching atoms.
  97    * This way energy conservation may get better, and it is not unreasonable:
  98    * It corresponds to an ad-hoc change in the force constants, which is exactly
  99    * large enough to compensate for the remaining force on the shell. The
 100    * shell forces then becomes zero.
 101    */
 102   for(i=0; (i<ns); i++) {
 103     s1 = s[i].shell;
 104     switch (s[i].nnucl) {
 105     case 1:
 106       n1 = s[i].nucl1;
 107       rvec_dec(f[n1],f[s1]);
 108       clear_rvec(f[s1]);
 109       break;
 110     case 2:
 111       n1 = s[i].nucl1;
 112       n2 = s[i].nucl2;
 113       /*m1 = mass[n1];
 114         m2 = mass[n2];
 115         tm = dt_1/(m1+m2);*/
 116       break;
 117     case 3:
 118       n1 = s[i].nucl1;
 119       n2 = s[i].nucl2;
 120       n3 = s[i].nucl3;
 121       /*m1 = mass[n1];
 122       m2 = mass[n2];
 123       m3 = mass[n3];
 124       tm = dt_1/(m1+m2+m3);*/
 125       break;
 126     default:
 127       fatal_error(0,"Shell %d has %d nuclei!",i,s[i].nnucl);
 128     }
 129   }
 130 }
 131
 132 static void predict_shells(FILE *log,rvec x[],rvec v[],real dt,
 133                            int ns,t_shell s[],
 134                            real mass[],bool bInit)
 135 {
 136   int  i,m,s1,n1,n2,n3;
 137   real dt_1,dt_2,dt_3,fudge,tm,m1,m2,m3;
 138   rvec *ptr;
 139
 140   /* We introduce a fudge factor for performance reasons: with this choice
 141    * the initial force on the shells is about a factor of two lower than
 142    * without
 143    */
 144   fudge = 1.0;
 145
 146   if (bInit) {
 147     fprintf(log,"RELAX: Using prediction for initial shell placement\n");
 148     ptr  = x;
 149     dt_1 = 1;
 150   }
 151   else {
 152     ptr  = v;
 153     dt_1 = fudge*dt;
 154   }
 155
 156   for(i=0; (i<ns); i++) {
 157     s1 = s[i].shell;
 158     if (bInit)
 159       clear_rvec(x[s1]);
 160     switch (s[i].nnucl) {
 161     case 1:
 162       n1 = s[i].nucl1;
 163       for(m=0; (m<DIM); m++)
 164         x[s1][m]+=ptr[n1][m]*dt_1;
 165       break;
 166     case 2:
 167       n1 = s[i].nucl1;
 168       n2 = s[i].nucl2;
 169       m1 = mass[n1];
 170       m2 = mass[n2];
 171       tm = dt_1/(m1+m2);
 172       for(m=0; (m<DIM); m++)
 173         x[s1][m]+=(m1*ptr[n1][m]+m2*ptr[n2][m])*tm;
 174       break;
 175     case 3:
 176       n1 = s[i].nucl1;
 177       n2 = s[i].nucl2;
 178       n3 = s[i].nucl3;
 179       m1 = mass[n1];
 180       m2 = mass[n2];
 181       m3 = mass[n3];
 182       tm = dt_1/(m1+m2+m3);
 183       for(m=0; (m<DIM); m++)
 184         x[s1][m]+=(m1*ptr[n1][m]+m2*ptr[n2][m]+m3*ptr[n3][m])*tm;
 185       break;
 186     default:
 187       fatal_error(0,"Shell %d has %d nuclei!",i,s[i].nnucl);
 188     }
 189   }
 190 }
 191
 192 static void print_epot(FILE *fp,int mdstep,int count,real step,real epot,
 193                        real df,bool bLast)
 194 {
 195   fprintf(fp,"MDStep=%5d/%2d lamb: %6g, EPot: %12.8e",
 196           mdstep,count,step,epot);
 197
 198   if (count != 0)
 199     fprintf(fp,", rmsF: %12.8e\n",df);
 200   else
 201     fprintf(fp,"\n");
 202 }
 203
 204
 205 static real rms_force(t_commrec *cr,rvec f[],int ns,t_shell s[])
 206 {
 207   int  i,shell;
 208   real df2;
 209
 210   if (!ns)
 211     return 0;
 212   df2=0.0;
 213   for(i=0; (i<ns); i++) {
 214     shell = s[i].shell;
 215     df2  += iprod(f[shell],f[shell]);
 216   }
 217   if (PAR(cr)) {
 218     gmx_sum(1,&df2,cr);
 219     gmx_sumi(1,&ns,cr);
 220   }
 221   return sqrt(df2/ns);
 222 }
 223
 224 static void check_pbc(FILE *fp,rvec x[],int shell)
 225 {
 226   int m,now;
 227
 228   now = shell-4;
 229   for(m=0; (m<DIM); m++)
 230     if (fabs(x[shell][m]-x[now][m]) > 0.3) {
 231       pr_rvecs(fp,0,"SHELL-X",x+now,5);
 232       break;
 233     }
 234 }
 235
 236 static void dump_shells(FILE *fp,rvec x[],rvec f[],real ftol,int ns,t_shell s[])
 237 {
 238   int  i,shell;
 239   real ft2,ff2;
 240
 241   ft2 = sqr(ftol);
 242
 243   for(i=0; (i<ns); i++) {
 244     shell = s[i].shell;
 245     ff2   = iprod(f[shell],f[shell]);
 246     if (ff2 > ft2)
 247       fprintf(fp,"SHELL %5d, force %10.5f  %10.5f  %10.5f, |f| %10.5f\n",
 248               shell,f[shell][XX],f[shell][YY],f[shell][ZZ],sqrt(ff2));
 249     check_pbc(fp,x,shell);
 250   }
 251 }
 252
 253 int relax_shells(FILE *log,t_commrec *cr,bool bVerbose,
 254                  int mdstep,t_parm *parm,bool bDoNS,bool bStopCM,
 255                  t_topology *top,real ener[],
 256                  rvec x[],rvec vold[],rvec v[],rvec vt[],rvec f[],
 257                  rvec buf[],t_mdatoms *md,t_nsborder *nsb,t_nrnb *nrnb,
 258                  t_graph *graph,t_groups *grps,tensor vir_part,
 259                  tensor pme_vir_part,
 260                  int nshell,t_shell shells[],t_forcerec *fr,
 261                  char *traj,real t,real lambda,rvec mu_tot,
 262                  int natoms,matrix box,bool *bConverged)
 263 {
 264   static bool bFirst=TRUE,bInit;
 265   static rvec *pos[2],*force[2];
 266   real   Epot[2],df[2],Estore[F_NRE];
 267   tensor my_vir[2],vir_last,pme_vir[2];
 268   rvec   *min_f_dir=NULL;
 269 #define NEPOT asize(Epot)
 270   real   ftol,step,step0,xiH,xiS,dum=0;
 271   char   cbuf[56];
 272   bool   bDone;
 273   int    g;
 274   int    number_steps;
 275   int    count=0;
 276   int    i,start=START(nsb),homenr=HOMENR(nsb),end=START(nsb)+HOMENR(nsb);
 277   int    Min=0;
 278 #define  Try (1-Min)             /* At start Try = 1 */
 279
 280   if (bFirst) {
 281     /* Allocate local arrays */
 282     for(i=0; (i<2); i++) {
 283       if (nshell > 0)
 284         snew(pos[i],nsb->natoms);
 285       snew(force[i],nsb->natoms);
 286     }
 287     bInit  = (getenv("FORCEINIT") != NULL);
 288     bFirst = FALSE;
 289   }
 290
 291   ftol         = parm->ir.em_tol;
 292   number_steps = parm->ir.niter;
 293   step0        = 1.0;
 294
 295   /* Do a prediction of the shell positions */
 296   predict_shells(log,x,v,parm->ir.delta_t,nshell,shells,
 297                  md->massT,bInit || (mdstep == 0));
 298
 299   /* Calculate the forces first time around */
 300   clear_mat(my_vir[Min]);
 301   clear_mat(pme_vir[Min]);
 302   do_force(log,cr,parm,nsb,my_vir[Min],pme_vir[Min],mdstep,nrnb,
 303            top,grps,x,v,force[Min],buf,md,ener,bVerbose && !PAR(cr),
 304            lambda,graph,bDoNS,FALSE,fr,mu_tot,FALSE);
 305   sum_lrforces(force[Min],fr,start,homenr);
 306
 307   df[Min]=rms_force(cr,force[Min],nshell,shells);
 308   df[Try]=0;
 309   if (debug) {
 310     fprintf(debug,"df = %g  %g\n",df[Min],df[Try]);
 311     sprintf(cbuf,"myvir step %d",0);
 312     pr_rvecs(debug,0,cbuf,my_vir[Min],DIM);
 313   }
 314
 315 #ifdef DEBUG
 316   if (debug) {
 317     pr_rvecs(debug,0,"force0",force[Min],md->nr);
 318   }
 319 #endif
 320   if (nshell > 0) {
 321     /* Copy x to pos[Min] & pos[Try]: during minimization only the
 322      * shell positions are updated, therefore the other particles must
 323      * be set here.
 324      */
 325     memcpy(pos[Min],x,nsb->natoms*sizeof(x[0]));
 326     memcpy(pos[Try],x,nsb->natoms*sizeof(x[0]));
 327   }
 328   /* Sum the potential energy terms from group contributions */
 329   sum_epot(&(parm->ir.opts),grps,ener);
 330   Epot[Min]=ener[F_EPOT];
 331
 332   if (PAR(cr))
 333     gmx_sum(NEPOT,Epot,cr);
 334
 335   step=step0;
 336
 337   if (bVerbose && MASTER(cr) && (nshell > 0))
 338     print_epot(stdout,mdstep,0,step,Epot[Min],df[Min],FALSE);
 339
 340   if (debug) {
 341     fprintf(debug,"%17s: %14.10e\n",
 342             interaction_function[F_EKIN].longname, ener[F_EKIN]);
 343     fprintf(debug,"%17s: %14.10e\n",
 344             interaction_function[F_EPOT].longname, ener[F_EPOT]);
 345     fprintf(debug,"%17s: %14.10e\n",
 346             interaction_function[F_ETOT].longname, ener[F_ETOT]);
 347     fprintf(debug,"SHELLSTEP %d\n",mdstep);
 348   }
 349
 350   /* First check whether we should do shells, or whether the force is
 351    * low enough even without minimization.
 352    */
 353   *bConverged = bDone = (df[Min] < ftol) || (nshell == 0);
 354
 355   for(count=1; (!bDone && (count < number_steps)); count++) {
 356
 357     /* New positions, Steepest descent */
 358     shell_pos_sd(log,step,pos[Min],pos[Try],force[Min],nshell,shells);
 359
 360     if (fr->k_dirmin) {
 361       if (min_f_dir == NULL)
 362         snew(min_f_dir,homenr);
 363       constrain(log,top,&(parm->ir),mdstep,md,start,end,
 364                 pos[Min],force[Min],min_f_dir,parm->box,
 365                 lambda,&dum,nrnb,FALSE);
 366       directional_sd(log,step,pos[Min],pos[Try],min_f_dir,start,homenr,
 367                      fr->k_dirmin);
 368     }
 369
 370     if (debug) {
 371       pr_rvecs(debug,0,"pos[Try] b4 do_force",pos[Try] + start,homenr);
 372       pr_rvecs(debug,0,"pos[Min] b4 do_force",pos[Min] + start,homenr);
 373     }
 374     /* Try the new positions */
 375     clear_mat(my_vir[Try]);
 376     clear_mat(pme_vir[Try]);
 377     do_force(log,cr,parm,nsb,my_vir[Try],pme_vir[Try],1,nrnb,
 378              top,grps,pos[Try],v,force[Try],buf,md,ener,bVerbose && !PAR(cr),
 379              lambda,graph,FALSE,FALSE,fr,mu_tot,FALSE);
 380     sum_lrforces(force[Try],fr,start,homenr);
 381     df[Try]=rms_force(cr,force[Try],nshell,shells);
 382     if (debug)
 383       fprintf(debug,"df = %g  %g\n",df[Min],df[Try]);
 384
 385     if (debug) {
 386       /*pr_rvecs(debug,0,"F na do_force",force[Try] + start,homenr);*/
 387       sprintf(cbuf,"myvir step %d",count);
 388       pr_rvecs(debug,0,cbuf,my_vir[Try],DIM);
 389       /*fprintf(debug,"SHELL ITER %d\n",count);
 390         dump_shells(debug,pos[Try],force[Try],ftol,nshell,shells);*/
 391     }
 392     /* Sum the potential energy terms from group contributions */
 393     sum_epot(&(parm->ir.opts),grps,ener);
 394     Epot[Try]=ener[F_EPOT];
 395
 396     if (PAR(cr))
 397       gmx_sum(1,&Epot[Try],cr);
 398
 399     if (bVerbose && MASTER(cr))
 400       print_epot(stdout,mdstep,count,step,Epot[Try],df[Try],FALSE);
 401
 402     *bConverged = (df[Try] < ftol);
 403     bDone       = *bConverged || (step < 0.01);
 404
 405     /* if ((Epot[Try] < Epot[Min])) { */
 406     if ((df[Try] < df[Min])) {
 407       if (debug)
 408         fprintf(debug,"Swapping Min and Try\n");
 409       Min  = Try;
 410       step = step0;
 411     }
 412     else
 413       step *= 0.8;
 414   }
 415   if (MASTER(cr) && !bDone)
 416     fprintf(stderr,"EM did not converge in %d steps\n",number_steps);
 417
 418   /*if (bDone)
 419     zero_shell_forces(log,force[Min],nshell,shells);
 420   */
 421   /* Parallelise this one! */
 422   if (EEL_LR(fr->eeltype)) {
 423     for(i=start; (i<end); i++)
 424       rvec_dec(force[Min][i],fr->f_pme[i]);
 425   }
 426   memcpy(f,force[Min],nsb->natoms*sizeof(f[0]));
 427
 428   /* CHECK VIRIAL */
 429   copy_mat(my_vir[Min],vir_part);
 430   copy_mat(pme_vir[Min],pme_vir_part);
 431
 432   if (debug) {
 433     sprintf(cbuf,"myvir step %d",count);
 434     pr_rvecs(debug,0,cbuf,vir_part,DIM);
 435   }
 436   if (nshell > 0)
 437     memcpy(x,pos[Min],nsb->natoms*sizeof(x[0]));
 438
 439   return count;
 440 }
 441