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