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