src/mdlib/shellfc.c

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