src/mdlib/shellfc.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2008, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  *
  31  * And Hey:
  32  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  33  */
  34 #ifdef HAVE_CONFIG_H
  35 #include <config.h>
  36 #endif
  37
  38 #include <string.h>
  39 #include "typedefs.h"
  40 #include "smalloc.h"
  41 #include "gmx_fatal.h"
  42 #include "vec.h"
  43 #include "txtdump.h"
  44 #include "mdrun.h"
  45 #include "partdec.h"
  46 #include "mdatoms.h"
  47 #include "vsite.h"
  48 #include "network.h"
  49 #include "names.h"
  50 #include "constr.h"
  51 #include "domdec.h"
  52 #include "partdec.h"
  53 #include "physics.h"
  54 #include "copyrite.h"
  55 #include "shellfc.h"
  56 #include "mtop_util.h"
  57 #include "chargegroup.h"
  58
  59
  60 typedef struct {
  61   int     nnucl;
  62   atom_id shell;                /* The shell id                         */
  63   atom_id nucl1,nucl2,nucl3;    /* The nuclei connected to the shell    */
  64   /* gmx_bool    bInterCG; */       /* Coupled to nuclei outside cg?        */
  65   real    k;                    /* force constant                       */
  66   real    k_1;                  /* 1 over force constant                */
  67   rvec    xold;
  68   rvec    fold;
  69   rvec    step;
  70 } t_shell;
  71
  72 typedef struct gmx_shellfc {
  73   int     nshell_gl;       /* The number of shells in the system       */
  74   t_shell *shell_gl;       /* All the shells (for DD only)             */
  75   int     *shell_index_gl; /* Global shell index (for DD only)         */
  76   gmx_bool    bInterCG;        /* Are there inter charge-group shells?     */
  77   int     nshell;          /* The number of local shells               */
  78   t_shell *shell;          /* The local shells                         */
  79   int     shell_nalloc;    /* The allocation size of shell             */
  80   gmx_bool    bPredict;        /* Predict shell positions                  */
  81   gmx_bool    bForceInit;      /* Force initialization of shell positions  */
  82   int     nflexcon;        /* The number of flexible constraints       */
  83   rvec    *x[2];           /* Array for iterative minimization         */
  84   rvec    *f[2];           /* Array for iterative minimization         */
  85   int     x_nalloc;        /* The allocation size of x and f           */
  86   rvec    *acc_dir;        /* Acceleration direction for flexcon       */
  87   rvec    *x_old;          /* Old coordinates for flexcon              */
  88   int     flex_nalloc;     /* The allocation size of acc_dir and x_old */
  89   rvec    *adir_xnold;     /* Work space for init_adir                 */
  90   rvec    *adir_xnew;      /* Work space for init_adir                 */
  91   int     adir_nalloc;     /* Work space for init_adir                 */
  92 } t_gmx_shellfc;
  93
  94
  95 static void pr_shell(FILE *fplog,int ns,t_shell s[])
  96 {
  97   int i;
  98
  99   fprintf(fplog,"SHELL DATA\n");
 100   fprintf(fplog,"%5s  %8s  %5s  %5s  %5s\n",
 101           "Shell","Force k","Nucl1","Nucl2","Nucl3");
 102   for(i=0; (i<ns); i++) {
 103     fprintf(fplog,"%5d  %8.3f  %5d",s[i].shell,1.0/s[i].k_1,s[i].nucl1);
 104     if (s[i].nnucl == 2)
 105       fprintf(fplog,"  %5d\n",s[i].nucl2);
 106     else if (s[i].nnucl == 3)
 107       fprintf(fplog,"  %5d  %5d\n",s[i].nucl2,s[i].nucl3);
 108     else
 109       fprintf(fplog,"\n");
 110   }
 111 }
 112
 113 static void predict_shells(FILE *fplog,rvec x[],rvec v[],real dt,
 114                            int ns,t_shell s[],
 115                            real mass[],gmx_mtop_t *mtop,gmx_bool bInit)
 116 {
 117   int  i,m,s1,n1,n2,n3;
 118   real dt_1,dt_2,dt_3,fudge,tm,m1,m2,m3;
 119   rvec *ptr;
 120   t_atom *atom;
 121
 122   /* We introduce a fudge factor for performance reasons: with this choice
 123    * the initial force on the shells is about a factor of two lower than
 124    * without
 125    */
 126   fudge = 1.0;
 127
 128   if (bInit) {
 129     if (fplog)
 130       fprintf(fplog,"RELAX: Using prediction for initial shell placement\n");
 131     ptr  = x;
 132     dt_1 = 1;
 133   }
 134   else {
 135     ptr  = v;
 136     dt_1 = fudge*dt;
 137   }
 138
 139   for(i=0; (i<ns); i++) {
 140     s1 = s[i].shell;
 141     if (bInit)
 142       clear_rvec(x[s1]);
 143     switch (s[i].nnucl) {
 144     case 1:
 145       n1 = s[i].nucl1;
 146       for(m=0; (m<DIM); m++)
 147         x[s1][m]+=ptr[n1][m]*dt_1;
 148       break;
 149     case 2:
 150       n1 = s[i].nucl1;
 151       n2 = s[i].nucl2;
 152       if (mass) {
 153         m1 = mass[n1];
 154         m2 = mass[n2];
 155       } else {
 156         /* Not the correct masses with FE, but it is just a prediction... */
 157         m1 = atom[n1].m;
 158         m2 = atom[n2].m;
 159       }
 160       tm = dt_1/(m1+m2);
 161       for(m=0; (m<DIM); m++)
 162         x[s1][m]+=(m1*ptr[n1][m]+m2*ptr[n2][m])*tm;
 163       break;
 164     case 3:
 165       n1 = s[i].nucl1;
 166       n2 = s[i].nucl2;
 167       n3 = s[i].nucl3;
 168       if (mass) {
 169         m1 = mass[n1];
 170         m2 = mass[n2];
 171         m3 = mass[n3];
 172       } else {
 173         /* Not the correct masses with FE, but it is just a prediction... */
 174         gmx_mtop_atomnr_to_atom(mtop,n1,&atom);
 175         m1 = atom->m;
 176         gmx_mtop_atomnr_to_atom(mtop,n2,&atom);
 177         m2 = atom->m;
 178         gmx_mtop_atomnr_to_atom(mtop,n3,&atom);
 179         m3 = atom->m;
 180       }
 181       tm = dt_1/(m1+m2+m3);
 182       for(m=0; (m<DIM); m++)
 183         x[s1][m]+=(m1*ptr[n1][m]+m2*ptr[n2][m]+m3*ptr[n3][m])*tm;
 184       break;
 185     default:
 186       gmx_fatal(FARGS,"Shell %d has %d nuclei!",i,s[i].nnucl);
 187     }
 188   }
 189 }
 190
 191 gmx_shellfc_t init_shell_flexcon(FILE *fplog,
 192                                  gmx_mtop_t *mtop,int nflexcon,
 193                                  rvec *x)
 194 {
 195   struct gmx_shellfc *shfc;
 196   t_shell     *shell;
 197   int         *shell_index,*at2cg;
 198   t_atom      *atom;
 199   int         n[eptNR],ns,nshell,nsi;
 200   int         i,j,nmol,type,mb,mt,a_offset,cg,mol,ftype,nra;
 201   real        qS,alpha;
 202   int         aS,aN=0; /* Shell and nucleus */
 203   int         bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_WATER_POL };
 204 #define NBT asize(bondtypes)
 205   t_iatom     *ia;
 206   gmx_mtop_atomloop_all_t aloop;
 207   gmx_ffparams_t *ffparams;
 208   gmx_molblock_t *molb;
 209   gmx_moltype_t *molt;
 210   t_block     *cgs;
 211
 212   /* Count number of shells, and find their indices */
 213   for(i=0; (i<eptNR); i++) {
 214     n[i] = 0;
 215   }
 216   snew(shell_index,mtop->natoms);
 217
 218   aloop = gmx_mtop_atomloop_all_init(mtop);
 219   while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
 220     if (atom->ptype == eptShell) {
 221       shell_index[i] = n[atom->ptype];
 222     }
 223     n[atom->ptype]++;
 224   }
 225
 226   /*
 227   snew(shell_index,end-start);
 228   snew(at2cg,end-start);
 229   nsi = 0;
 230   for(cg=cg0; (cg<cg1); cg++) {
 231     for(i=cgindex[cg]; i<cgindex[cg+1]; i++) {
 232       n[atom[i].ptype]++;
 233       if (atom[i].ptype == eptShell)
 234         shell_index[i-start] = nsi++;
 235       at2cg[i-start] = cg;
 236     }
 237   }
 238   */
 239
 240   if (fplog) {
 241     /* Print the number of each particle type */
 242     for(i=0; (i<eptNR); i++) {
 243       if (n[i] != 0) {
 244         fprintf(fplog,"There are: %d %ss\n",n[i],ptype_str[i]);
 245       }
 246     }
 247   }
 248
 249   nshell = n[eptShell];
 250
 251   if (nshell == 0 && nflexcon == 0) {
 252     sfree(shell_index);
 253
 254     return NULL;
 255   }
 256
 257   snew(shfc,1);
 258   shfc->nflexcon = nflexcon;
 259
 260   if (nshell == 0) {
 261     sfree(shell_index);
 262
 263     return shfc;
 264   }
 265
 266   snew(shell,nshell);
 267   shfc->nflexcon = nflexcon;
 268
 269   /* Initiate the shell structures */
 270   for(i=0; (i<nshell); i++) {
 271     shell[i].shell = NO_ATID;
 272     shell[i].nnucl = 0;
 273     shell[i].nucl1 = NO_ATID;
 274     shell[i].nucl2 = NO_ATID;
 275     shell[i].nucl3 = NO_ATID;
 276     /* shell[i].bInterCG=FALSE; */
 277     shell[i].k_1   = 0;
 278     shell[i].k     = 0;
 279   }
 280
 281   ffparams = &mtop->ffparams;
 282
 283   /* Now fill the structures */
 284   shfc->bInterCG = FALSE;
 285   ns = 0;
 286   a_offset = 0;
 287   for(mb=0; mb<mtop->nmolblock; mb++) {
 288     molb = &mtop->molblock[mb];
 289     molt = &mtop->moltype[molb->type];
 290
 291     cgs = &molt->cgs;
 292     snew(at2cg,molt->atoms.nr);
 293     for(cg=0; cg<cgs->nr; cg++) {
 294       for(i=cgs->index[cg]; i<cgs->index[cg+1]; i++) {
 295         at2cg[i] = cg;
 296       }
 297     }
 298
 299     atom = molt->atoms.atom;
 300     for(mol=0; mol<molb->nmol; mol++) {
 301       for(j=0; (j<NBT); j++) {
 302         ia = molt->ilist[bondtypes[j]].iatoms;
 303         for(i=0; (i<molt->ilist[bondtypes[j]].nr); ) {
 304           type  = ia[0];
 305           ftype = ffparams->functype[type];
 306           nra   = interaction_function[ftype].nratoms;
 307
 308           /* Check whether we have a bond with a shell */
 309           aS = NO_ATID;
 310
 311           switch (bondtypes[j]) {
 312           case F_BONDS:
 313           case F_HARMONIC:
 314           case F_CUBICBONDS:
 315           case F_POLARIZATION:
 316             if (atom[ia[1]].ptype == eptShell) {
 317               aS = ia[1];
 318               aN = ia[2];
 319             }
 320             else if (atom[ia[2]].ptype == eptShell) {
 321               aS = ia[2];
 322               aN = ia[1];
 323             }
 324             break;
 325           case F_WATER_POL:
 326             aN    = ia[4];  /* Dummy */
 327             aS    = ia[5];  /* Shell */
 328             break;
 329           default:
 330             gmx_fatal(FARGS,"Death Horror: %s, %d",__FILE__,__LINE__);
 331           }
 332
 333           if (aS != NO_ATID) {
 334             qS = atom[aS].q;
 335
 336             /* Check whether one of the particles is a shell... */
 337             nsi = shell_index[a_offset+aS];
 338             if ((nsi < 0) || (nsi >= nshell))
 339               gmx_fatal(FARGS,"nsi is %d should be within 0 - %d. aS = %d",
 340                         nsi,nshell,aS);
 341             if (shell[nsi].shell == NO_ATID) {
 342               shell[nsi].shell = a_offset + aS;
 343               ns ++;
 344             }
 345             else if (shell[nsi].shell != a_offset+aS)
 346               gmx_fatal(FARGS,"Weird stuff in %s, %d",__FILE__,__LINE__);
 347
 348             if      (shell[nsi].nucl1 == NO_ATID) {
 349               shell[nsi].nucl1 = a_offset + aN;
 350             } else if (shell[nsi].nucl2 == NO_ATID) {
 351               shell[nsi].nucl2 = a_offset + aN;
 352             } else if (shell[nsi].nucl3 == NO_ATID) {
 353               shell[nsi].nucl3 = a_offset + aN;
 354             } else {
 355               if (fplog)
 356                 pr_shell(fplog,ns,shell);
 357               gmx_fatal(FARGS,"Can not handle more than three bonds per shell\n");
 358             }
 359             if (at2cg[aS] != at2cg[aN]) {
 360               /* shell[nsi].bInterCG = TRUE; */
 361               shfc->bInterCG = TRUE;
 362             }
 363
 364             switch (bondtypes[j]) {
 365             case F_BONDS:
 366             case F_HARMONIC:
 367               shell[nsi].k    += ffparams->iparams[type].harmonic.krA;
 368               break;
 369             case F_CUBICBONDS:
 370               shell[nsi].k    += ffparams->iparams[type].cubic.kb;
 371               break;
 372             case F_POLARIZATION:
 373               if (qS != atom[aS].qB)
 374                 gmx_fatal(FARGS,"polarize can not be used with qA != qB");
 375               shell[nsi].k    += sqr(qS)*ONE_4PI_EPS0/
 376                 ffparams->iparams[type].polarize.alpha;
 377               break;
 378             case F_WATER_POL:
 379               if (qS != atom[aS].qB)
 380                 gmx_fatal(FARGS,"water_pol can not be used with qA != qB");
 381               alpha          = (ffparams->iparams[type].wpol.al_x+
 382                                 ffparams->iparams[type].wpol.al_y+
 383                                 ffparams->iparams[type].wpol.al_z)/3.0;
 384               shell[nsi].k  += sqr(qS)*ONE_4PI_EPS0/alpha;
 385               break;
 386             default:
 387               gmx_fatal(FARGS,"Death Horror: %s, %d",__FILE__,__LINE__);
 388             }
 389             shell[nsi].nnucl++;
 390           }
 391           ia += nra+1;
 392           i  += nra+1;
 393         }
 394       }
 395       a_offset += molt->atoms.nr;
 396     }
 397     /* Done with this molecule type */
 398     sfree(at2cg);
 399   }
 400
 401   /* Verify whether it's all correct */
 402   if (ns != nshell)
 403     gmx_fatal(FARGS,"Something weird with shells. They may not be bonded to something");
 404
 405   for(i=0; (i<ns); i++)
 406     shell[i].k_1 = 1.0/shell[i].k;
 407
 408   if (debug)
 409     pr_shell(debug,ns,shell);
 410
 411
 412   shfc->nshell_gl      = ns;
 413   shfc->shell_gl       = shell;
 414   shfc->shell_index_gl = shell_index;
 415
 416   shfc->bPredict   = (getenv("GMX_NOPREDICT") == NULL);
 417   shfc->bForceInit = FALSE;
 418   if (!shfc->bPredict) {
 419     if (fplog)
 420       fprintf(fplog,"\nWill never predict shell positions\n");
 421   } else {
 422     shfc->bForceInit = (getenv("GMX_FORCEINIT") != NULL);
 423     if (shfc->bForceInit && fplog)
 424       fprintf(fplog,"\nWill always initiate shell positions\n");
 425   }
 426
 427   if (shfc->bPredict) {
 428     if (x) {
 429       predict_shells(fplog,x,NULL,0,shfc->nshell_gl,shfc->shell_gl,
 430                      NULL,mtop,TRUE);
 431     }
 432
 433     if (shfc->bInterCG) {
 434       if (fplog)
 435         fprintf(fplog,"\nNOTE: there all shells that are connected to particles outside thier own charge group, will not predict shells positions during the run\n\n");
 436       shfc->bPredict = FALSE;
 437     }
 438   }
 439
 440   return shfc;
 441 }
 442
 443 void make_local_shells(t_commrec *cr,t_mdatoms *md,
 444                        struct gmx_shellfc *shfc)
 445 {
 446   t_shell *shell;
 447   int a0,a1,*ind,nshell,i;
 448   gmx_domdec_t *dd=NULL;
 449
 450   if (PAR(cr)) {
 451     if (DOMAINDECOMP(cr)) {
 452       dd = cr->dd;
 453       a0 = 0;
 454       a1 = dd->nat_home;
 455     } else {
 456       pd_at_range(cr,&a0,&a1);
 457     }
 458   } else {
 459     /* Single node: we need all shells, just copy the pointer */
 460     shfc->nshell = shfc->nshell_gl;
 461     shfc->shell  = shfc->shell_gl;
 462
 463     return;
 464   }
 465
 466   ind = shfc->shell_index_gl;
 467
 468   nshell = 0;
 469   shell  = shfc->shell;
 470   for(i=a0; i<a1; i++) {
 471     if (md->ptype[i] == eptShell) {
 472       if (nshell+1 > shfc->shell_nalloc) {
 473         shfc->shell_nalloc = over_alloc_dd(nshell+1);
 474         srenew(shell,shfc->shell_nalloc);
 475       }
 476       if (dd) {
 477         shell[nshell] = shfc->shell_gl[ind[dd->gatindex[i]]];
 478       } else {
 479         shell[nshell] = shfc->shell_gl[ind[i]];
 480       }
 481       /* With inter-cg shells we can no do shell prediction,
 482        * so we do not need the nuclei numbers.
 483        */
 484       if (!shfc->bInterCG) {
 485         shell[nshell].nucl1   = i + shell[nshell].nucl1 - shell[nshell].shell;
 486         if (shell[nshell].nnucl > 1)
 487           shell[nshell].nucl2 = i + shell[nshell].nucl2 - shell[nshell].shell;
 488         if (shell[nshell].nnucl > 2)
 489           shell[nshell].nucl3 = i + shell[nshell].nucl3 - shell[nshell].shell;
 490       }
 491       shell[nshell].shell = i;
 492       nshell++;
 493     }
 494   }
 495
 496   shfc->nshell = nshell;
 497   shfc->shell  = shell;
 498 }
 499
 500 static void do_1pos(rvec xnew,rvec xold,rvec f,real step)
 501 {
 502   real xo,yo,zo;
 503   real dx,dy,dz;
 504
 505   xo=xold[XX];
 506   yo=xold[YY];
 507   zo=xold[ZZ];
 508
 509   dx=f[XX]*step;
 510   dy=f[YY]*step;
 511   dz=f[ZZ]*step;
 512
 513   xnew[XX]=xo+dx;
 514   xnew[YY]=yo+dy;
 515   xnew[ZZ]=zo+dz;
 516 }
 517
 518 static void do_1pos3(rvec xnew,rvec xold,rvec f,rvec step)
 519 {
 520   real xo,yo,zo;
 521   real dx,dy,dz;
 522
 523   xo=xold[XX];
 524   yo=xold[YY];
 525   zo=xold[ZZ];
 526
 527   dx=f[XX]*step[XX];
 528   dy=f[YY]*step[YY];
 529   dz=f[ZZ]*step[ZZ];
 530
 531   xnew[XX]=xo+dx;
 532   xnew[YY]=yo+dy;
 533   xnew[ZZ]=zo+dz;
 534 }
 535
 536 static void directional_sd(FILE *log,rvec xold[],rvec xnew[],rvec acc_dir[],
 537                            int start,int homenr,real step)
 538 {
 539   int  i;
 540
 541   for(i=start; i<homenr; i++)
 542     do_1pos(xnew[i],xold[i],acc_dir[i],step);
 543 }
 544
 545 static void shell_pos_sd(FILE *log,rvec xcur[],rvec xnew[],rvec f[],
 546                          int ns,t_shell s[],int count)
 547 {
 548   int  i,shell,d;
 549   real dx,df,k_est;
 550 #ifdef PRINT_STEP
 551   real step_min,step_max;
 552
 553   step_min = 1e30;
 554   step_max = 0;
 555 #endif
 556   for(i=0; (i<ns); i++) {
 557     shell = s[i].shell;
 558     if (count == 1) {
 559       for(d=0; d<DIM; d++) {
 560         s[i].step[d] = s[i].k_1;
 561 #ifdef PRINT_STEP
 562         step_min = min(step_min,s[i].step[d]);
 563         step_max = max(step_max,s[i].step[d]);
 564 #endif
 565       }
 566     } else {
 567       for(d=0; d<DIM; d++) {
 568         dx = xcur[shell][d] - s[i].xold[d];
 569         df =    f[shell][d] - s[i].fold[d];
 570         if (dx != 0 && df != 0) {
 571           k_est = -dx/df;
 572           if (k_est >= 2*s[i].step[d]) {
 573             s[i].step[d] *= 1.2;
 574           } else if (k_est <= 0) {
 575             s[i].step[d] *= 0.8;
 576           } else {
 577             s[i].step[d] = 0.8*s[i].step[d] + 0.2*k_est;
 578           }
 579         } else if (dx != 0) {
 580           s[i].step[d] *= 1.2;
 581         }
 582 #ifdef PRINT_STEP
 583         step_min = min(step_min,s[i].step[d]);
 584         step_max = max(step_max,s[i].step[d]);
 585 #endif
 586       }
 587     }
 588     copy_rvec(xcur[shell],s[i].xold);
 589     copy_rvec(f[shell],   s[i].fold);
 590
 591     do_1pos3(xnew[shell],xcur[shell],f[shell],s[i].step);
 592
 593     if (gmx_debug_at) {
 594       fprintf(debug,"shell[%d] = %d\n",i,shell);
 595       pr_rvec(debug,0,"fshell",f[shell],DIM,TRUE);
 596       pr_rvec(debug,0,"xold",xcur[shell],DIM,TRUE);
 597       pr_rvec(debug,0,"step",s[i].step,DIM,TRUE);
 598       pr_rvec(debug,0,"xnew",xnew[shell],DIM,TRUE);
 599     }
 600   }
 601 #ifdef PRINT_STEP
 602   printf("step %.3e %.3e\n",step_min,step_max);
 603 #endif
 604 }
 605
 606 static void decrease_step_size(int nshell,t_shell s[])
 607 {
 608   int i;
 609
 610   for(i=0; i<nshell; i++)
 611     svmul(0.8,s[i].step,s[i].step);
 612 }
 613
 614 static void print_epot(FILE *fp,gmx_large_int_t mdstep,int count,real epot,real df,
 615                        int ndir,real sf_dir)
 616 {
 617   char buf[22];
 618
 619   fprintf(fp,"MDStep=%5s/%2d EPot: %12.8e, rmsF: %6.2e",
 620           gmx_step_str(mdstep,buf),count,epot,df);
 621   if (ndir)
 622     fprintf(fp,", dir. rmsF: %6.2e\n",sqrt(sf_dir/ndir));
 623   else
 624     fprintf(fp,"\n");
 625 }
 626
 627
 628 static real rms_force(t_commrec *cr,rvec f[],int ns,t_shell s[],
 629                       int ndir,real *sf_dir,real *Epot)
 630 {
 631   int  i,shell,ntot;
 632   double buf[4];
 633
 634   buf[0] = *sf_dir;
 635   for(i=0; i<ns; i++) {
 636     shell = s[i].shell;
 637     buf[0]  += norm2(f[shell]);
 638   }
 639   ntot = ns;
 640
 641   if (PAR(cr)) {
 642     buf[1] = ntot;
 643     buf[2] = *sf_dir;
 644     buf[3] = *Epot;
 645     gmx_sumd(4,buf,cr);
 646     ntot = (int)(buf[1] + 0.5);
 647     *sf_dir = buf[2];
 648     *Epot   = buf[3];
 649   }
 650   ntot += ndir;
 651
 652   return (ntot ? sqrt(buf[0]/ntot) : 0);
 653 }
 654
 655 static void check_pbc(FILE *fp,rvec x[],int shell)
 656 {
 657   int m,now;
 658
 659   now = shell-4;
 660   for(m=0; (m<DIM); m++)
 661     if (fabs(x[shell][m]-x[now][m]) > 0.3) {
 662       pr_rvecs(fp,0,"SHELL-X",x+now,5);
 663       break;
 664     }
 665 }
 666
 667 static void dump_shells(FILE *fp,rvec x[],rvec f[],real ftol,int ns,t_shell s[])
 668 {
 669   int  i,shell;
 670   real ft2,ff2;
 671
 672   ft2 = sqr(ftol);
 673
 674   for(i=0; (i<ns); i++) {
 675     shell = s[i].shell;
 676     ff2   = iprod(f[shell],f[shell]);
 677     if (ff2 > ft2)
 678       fprintf(fp,"SHELL %5d, force %10.5f  %10.5f  %10.5f, |f| %10.5f\n",
 679               shell,f[shell][XX],f[shell][YY],f[shell][ZZ],sqrt(ff2));
 680     check_pbc(fp,x,shell);
 681   }
 682 }
 683
 684 static void init_adir(FILE *log,gmx_shellfc_t shfc,
 685                       gmx_constr_t constr,t_idef *idef,t_inputrec *ir,
 686                       t_commrec *cr,int dd_ac1,
 687                       gmx_large_int_t step,t_mdatoms *md,int start,int end,
 688                       rvec *x_old,rvec *x_init,rvec *x,
 689                       rvec *f,rvec *acc_dir,matrix box,
 690                       real lambda,real *dvdlambda,t_nrnb *nrnb)
 691 {
 692   rvec   *xnold,*xnew;
 693   double w_dt;
 694   int    gf,ga,gt;
 695   real   dt,scale;
 696   int    n,d;
 697   unsigned short *ptype;
 698   rvec   p,dx;
 699
 700   if (DOMAINDECOMP(cr))
 701     n = dd_ac1;
 702   else
 703     n = end - start;
 704   if (n > shfc->adir_nalloc) {
 705     shfc->adir_nalloc = over_alloc_dd(n);
 706     srenew(shfc->adir_xnold,shfc->adir_nalloc);
 707     srenew(shfc->adir_xnew ,shfc->adir_nalloc);
 708   }
 709   xnold = shfc->adir_xnold;
 710   xnew  = shfc->adir_xnew;
 711
 712   ptype = md->ptype;
 713
 714   dt = ir->delta_t;
 715
 716   /* Does NOT work with freeze or acceleration groups (yet) */
 717   for (n=start; n<end; n++) {
 718     w_dt = md->invmass[n]*dt;
 719
 720     for (d=0; d<DIM; d++) {
 721       if ((ptype[n] != eptVSite) && (ptype[n] != eptShell)) {
 722         xnold[n-start][d] = x[n][d] - (x_init[n][d] - x_old[n][d]);
 723         xnew[n-start][d] = 2*x[n][d] - x_old[n][d] + f[n][d]*w_dt*dt;
 724       } else {
 725         xnold[n-start][d] = x[n][d];
 726         xnew[n-start][d] = x[n][d];
 727       }
 728     }
 729   }
 730   constrain(log,FALSE,FALSE,constr,idef,ir,NULL,cr,step,0,md,
 731             x,xnold-start,NULL,box,
 732             lambda,dvdlambda,NULL,NULL,nrnb,econqCoord,FALSE,0,0);
 733   constrain(log,FALSE,FALSE,constr,idef,ir,NULL,cr,step,0,md,
 734             x,xnew-start,NULL,box,
 735             lambda,dvdlambda,NULL,NULL,nrnb,econqCoord,FALSE,0,0);
 736
 737   /* Set xnew to minus the acceleration */
 738   for (n=start; n<end; n++) {
 739     for(d=0; d<DIM; d++)
 740       xnew[n-start][d] =
 741         -(2*x[n][d]-xnold[n-start][d]-xnew[n-start][d])/sqr(dt)
 742         - f[n][d]*md->invmass[n];
 743     clear_rvec(acc_dir[n]);
 744   }
 745
 746   /* Project the acceleration on the old bond directions */
 747   constrain(log,FALSE,FALSE,constr,idef,ir,NULL,cr,step,0,md,
 748             x_old,xnew-start,acc_dir,box,
 749             lambda,dvdlambda,NULL,NULL,nrnb,econqDeriv_FlexCon,FALSE,0,0);
 750 }
 751
 752 int relax_shell_flexcon(FILE *fplog,t_commrec *cr,gmx_bool bVerbose,
 753                         gmx_large_int_t mdstep,t_inputrec *inputrec,
 754                         gmx_bool bDoNS,int force_flags,
 755                         gmx_bool bStopCM,
 756                         gmx_localtop_t *top,
 757                         gmx_mtop_t* mtop,
 758                         gmx_constr_t constr,
 759                         gmx_enerdata_t *enerd,t_fcdata *fcd,
 760                         t_state *state,rvec f[],
 761                         tensor force_vir,
 762                         t_mdatoms *md,
 763                         t_nrnb *nrnb,gmx_wallcycle_t wcycle,
 764                         t_graph *graph,
 765                         gmx_groups_t *groups,
 766                         struct gmx_shellfc *shfc,
 767                         t_forcerec *fr,
 768                         gmx_bool bBornRadii,
 769                         double t,rvec mu_tot,
 770                         int natoms,gmx_bool *bConverged,
 771                         gmx_vsite_t *vsite,
 772                         FILE *fp_field)
 773 {
 774   int    nshell;
 775   t_shell *shell;
 776   t_idef *idef;
 777   rvec   *pos[2],*force[2],*acc_dir=NULL,*x_old=NULL;
 778   real   Epot[2],df[2];
 779   rvec   dx;
 780   real   sf_dir,invdt;
 781   real   ftol,xiH,xiS,dum=0;
 782   char   sbuf[22];
 783   gmx_bool   bCont,bInit;
 784   int    nat,dd_ac0,dd_ac1=0,i;
 785   int    start=md->start,homenr=md->homenr,end=start+homenr,cg0,cg1;
 786   int    nflexcon,g,number_steps,d,Min=0,count=0;
 787 #define  Try (1-Min)             /* At start Try = 1 */
 788
 789   bCont        = (mdstep == inputrec->init_step) && inputrec->bContinuation;
 790   bInit        = (mdstep == inputrec->init_step) || shfc->bForceInit;
 791   ftol         = inputrec->em_tol;
 792   number_steps = inputrec->niter;
 793   nshell       = shfc->nshell;
 794   shell        = shfc->shell;
 795   nflexcon     = shfc->nflexcon;
 796
 797   idef = &top->idef;
 798
 799   if (DOMAINDECOMP(cr)) {
 800     nat = dd_natoms_vsite(cr->dd);
 801     if (nflexcon > 0) {
 802       dd_get_constraint_range(cr->dd,&dd_ac0,&dd_ac1);
 803       nat = max(nat,dd_ac1);
 804     }
 805   } else {
 806     nat = state->natoms;
 807   }
 808
 809   if (nat > shfc->x_nalloc) {
 810     /* Allocate local arrays */
 811     shfc->x_nalloc = over_alloc_dd(nat);
 812     for(i=0; (i<2); i++) {
 813       srenew(shfc->x[i],shfc->x_nalloc);
 814       srenew(shfc->f[i],shfc->x_nalloc);
 815     }
 816   }
 817   for(i=0; (i<2); i++) {
 818     pos[i]   = shfc->x[i];
 819     force[i] = shfc->f[i];
 820   }
 821
 822   /* With particle decomposition this code only works
 823    * when all particles involved with each shell are in the same cg.
 824    */
 825
 826   if (bDoNS && inputrec->ePBC != epbcNONE && !DOMAINDECOMP(cr)) {
 827     /* This is the only time where the coordinates are used
 828      * before do_force is called, which normally puts all
 829      * charge groups in the box.
 830      */
 831     if (PARTDECOMP(cr)) {
 832       pd_cg_range(cr,&cg0,&cg1);
 833     } else {
 834       cg0 = 0;
 835       cg1 = top->cgs.nr;
 836     }
 837     put_charge_groups_in_box(fplog,cg0,cg1,fr->ePBC,state->box,
 838                              &(top->cgs),state->x,fr->cg_cm);
 839     if (graph)
 840       mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
 841   }
 842
 843   /* After this all coordinate arrays will contain whole molecules */
 844   if (graph)
 845     shift_self(graph,state->box,state->x);
 846
 847   if (nflexcon) {
 848     if (nat > shfc->flex_nalloc) {
 849       shfc->flex_nalloc = over_alloc_dd(nat);
 850       srenew(shfc->acc_dir,shfc->flex_nalloc);
 851       srenew(shfc->x_old,shfc->flex_nalloc);
 852     }
 853     acc_dir = shfc->acc_dir;
 854     x_old   = shfc->x_old;
 855     for(i=0; i<homenr; i++) {
 856       for(d=0; d<DIM; d++)
 857         shfc->x_old[i][d] =
 858           state->x[start+i][d] - state->v[start+i][d]*inputrec->delta_t;
 859     }
 860   }
 861
 862   /* Do a prediction of the shell positions */
 863   if (shfc->bPredict && !bCont) {
 864     predict_shells(fplog,state->x,state->v,inputrec->delta_t,nshell,shell,
 865                    md->massT,NULL,bInit);
 866   }
 867
 868   /* do_force expected the charge groups to be in the box */
 869   if (graph)
 870     unshift_self(graph,state->box,state->x);
 871
 872   /* Calculate the forces first time around */
 873   if (gmx_debug_at) {
 874     pr_rvecs(debug,0,"x b4 do_force",state->x + start,homenr);
 875   }
 876   do_force(fplog,cr,inputrec,mdstep,nrnb,wcycle,top,mtop,groups,
 877            state->box,state->x,&state->hist,
 878            force[Min],force_vir,md,enerd,fcd,
 879            state->lambda,graph,
 880            fr,vsite,mu_tot,t,fp_field,NULL,bBornRadii,
 881            (bDoNS ? GMX_FORCE_NS : 0) | force_flags);
 882
 883   sf_dir = 0;
 884   if (nflexcon) {
 885     init_adir(fplog,shfc,
 886               constr,idef,inputrec,cr,dd_ac1,mdstep,md,start,end,
 887               shfc->x_old-start,state->x,state->x,force[Min],
 888               shfc->acc_dir-start,state->box,state->lambda,&dum,nrnb);
 889
 890     for(i=start; i<end; i++)
 891       sf_dir += md->massT[i]*norm2(shfc->acc_dir[i-start]);
 892   }
 893
 894   Epot[Min] = enerd->term[F_EPOT];
 895
 896   df[Min]=rms_force(cr,shfc->f[Min],nshell,shell,nflexcon,&sf_dir,&Epot[Min]);
 897   df[Try]=0;
 898   if (debug) {
 899     fprintf(debug,"df = %g  %g\n",df[Min],df[Try]);
 900   }
 901
 902   if (gmx_debug_at) {
 903     pr_rvecs(debug,0,"force0",force[Min],md->nr);
 904   }
 905
 906   if (nshell+nflexcon > 0) {
 907     /* Copy x to pos[Min] & pos[Try]: during minimization only the
 908      * shell positions are updated, therefore the other particles must
 909      * be set here.
 910      */
 911     memcpy(pos[Min],state->x,nat*sizeof(state->x[0]));
 912     memcpy(pos[Try],state->x,nat*sizeof(state->x[0]));
 913   }
 914
 915   if (bVerbose && MASTER(cr))
 916     print_epot(stdout,mdstep,0,Epot[Min],df[Min],nflexcon,sf_dir);
 917
 918   if (debug) {
 919     fprintf(debug,"%17s: %14.10e\n",
 920             interaction_function[F_EKIN].longname,enerd->term[F_EKIN]);
 921     fprintf(debug,"%17s: %14.10e\n",
 922             interaction_function[F_EPOT].longname,enerd->term[F_EPOT]);
 923     fprintf(debug,"%17s: %14.10e\n",
 924             interaction_function[F_ETOT].longname,enerd->term[F_ETOT]);
 925     fprintf(debug,"SHELLSTEP %s\n",gmx_step_str(mdstep,sbuf));
 926   }
 927
 928   /* First check whether we should do shells, or whether the force is
 929    * low enough even without minimization.
 930    */
 931   *bConverged = (df[Min] < ftol);
 932
 933   for(count=1; (!(*bConverged) && (count < number_steps)); count++) {
 934     if (vsite)
 935       construct_vsites(fplog,vsite,pos[Min],nrnb,inputrec->delta_t,state->v,
 936                        idef->iparams,idef->il,
 937                        fr->ePBC,fr->bMolPBC,graph,cr,state->box);
 938
 939     if (nflexcon) {
 940       init_adir(fplog,shfc,
 941                 constr,idef,inputrec,cr,dd_ac1,mdstep,md,start,end,
 942                 x_old-start,state->x,pos[Min],force[Min],acc_dir-start,
 943                 state->box,state->lambda,&dum,nrnb);
 944
 945       directional_sd(fplog,pos[Min],pos[Try],acc_dir-start,start,end,
 946                      fr->fc_stepsize);
 947     }
 948
 949     /* New positions, Steepest descent */
 950     shell_pos_sd(fplog,pos[Min],pos[Try],force[Min],nshell,shell,count);
 951
 952     /* do_force expected the charge groups to be in the box */
 953     if (graph)
 954       unshift_self(graph,state->box,pos[Try]);
 955
 956     if (gmx_debug_at) {
 957       pr_rvecs(debug,0,"RELAX: pos[Min]  ",pos[Min] + start,homenr);
 958       pr_rvecs(debug,0,"RELAX: pos[Try]  ",pos[Try] + start,homenr);
 959     }
 960     /* Try the new positions */
 961     do_force(fplog,cr,inputrec,1,nrnb,wcycle,
 962              top,mtop,groups,state->box,pos[Try],&state->hist,
 963              force[Try],force_vir,
 964              md,enerd,fcd,state->lambda,graph,
 965              fr,vsite,mu_tot,t,fp_field,NULL,bBornRadii,
 966              force_flags);
 967
 968     if (gmx_debug_at) {
 969       pr_rvecs(debug,0,"RELAX: force[Min]",force[Min] + start,homenr);
 970       pr_rvecs(debug,0,"RELAX: force[Try]",force[Try] + start,homenr);
 971     }
 972     sf_dir = 0;
 973     if (nflexcon) {
 974       init_adir(fplog,shfc,
 975                 constr,idef,inputrec,cr,dd_ac1,mdstep,md,start,end,
 976                 x_old-start,state->x,pos[Try],force[Try],acc_dir-start,
 977                 state->box,state->lambda,&dum,nrnb);
 978
 979       for(i=start; i<end; i++)
 980         sf_dir += md->massT[i]*norm2(acc_dir[i-start]);
 981     }
 982
 983     Epot[Try] = enerd->term[F_EPOT];
 984
 985     df[Try]=rms_force(cr,force[Try],nshell,shell,nflexcon,&sf_dir,&Epot[Try]);
 986
 987     if (debug)
 988       fprintf(debug,"df = %g  %g\n",df[Min],df[Try]);
 989
 990     if (debug) {
 991       if (gmx_debug_at)
 992         pr_rvecs(debug,0,"F na do_force",force[Try] + start,homenr);
 993       if (gmx_debug_at) {
 994         fprintf(debug,"SHELL ITER %d\n",count);
 995         dump_shells(debug,pos[Try],force[Try],ftol,nshell,shell);
 996       }
 997     }
 998
 999     if (bVerbose && MASTER(cr))
1000       print_epot(stdout,mdstep,count,Epot[Try],df[Try],nflexcon,sf_dir);
1001
1002     *bConverged = (df[Try] < ftol);
1003
1004     if ((df[Try] < df[Min])) {
1005       if (debug)
1006         fprintf(debug,"Swapping Min and Try\n");
1007       if (nflexcon) {
1008         /* Correct the velocities for the flexible constraints */
1009         invdt = 1/inputrec->delta_t;
1010         for(i=start; i<end; i++) {
1011           for(d=0; d<DIM; d++)
1012             state->v[i][d] += (pos[Try][i][d] - pos[Min][i][d])*invdt;
1013         }
1014       }
1015       Min  = Try;
1016     } else {
1017       decrease_step_size(nshell,shell);
1018     }
1019   }
1020   if (MASTER(cr) && !(*bConverged)) {
1021     /* Note that the energies and virial are incorrect when not converged */
1022     if (fplog)
1023       fprintf(fplog,
1024               "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1025               gmx_step_str(mdstep,sbuf),number_steps,df[Min]);
1026     fprintf(stderr,
1027             "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1028             gmx_step_str(mdstep,sbuf),number_steps,df[Min]);
1029   }
1030
1031   /* Copy back the coordinates and the forces */
1032   memcpy(state->x,pos[Min],nat*sizeof(state->x[0]));
1033   memcpy(f,force[Min],nat*sizeof(f[0]));
1034
1035   return count;
1036 }
1037