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