src/gromacs/mdlib/mdatom.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-2004, The GROMACS development team.
   6  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42
  43 #include "typedefs.h"
  44 #include "mdatoms.h"
  45 #include "smalloc.h"
  46 #include "main.h"
  47 #include "qmmm.h"
  48 #include "mtop_util.h"
  49 #include "gmx_omp_nthreads.h"
  50
  51 #define ALMOST_ZERO 1e-30
  52
  53 t_mdatoms *init_mdatoms(FILE *fp, gmx_mtop_t *mtop, gmx_bool bFreeEnergy)
  54 {
  55     int                     mb, a, g, nmol;
  56     double                  tmA, tmB;
  57     t_atom                 *atom;
  58     t_mdatoms              *md;
  59     gmx_mtop_atomloop_all_t aloop;
  60     t_ilist                *ilist;
  61
  62     snew(md, 1);
  63
  64     md->nenergrp = mtop->groups.grps[egcENER].nr;
  65     md->bVCMgrps = FALSE;
  66     tmA          = 0.0;
  67     tmB          = 0.0;
  68
  69     aloop = gmx_mtop_atomloop_all_init(mtop);
  70     while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
  71     {
  72         if (ggrpnr(&mtop->groups, egcVCM, a) > 0)
  73         {
  74             md->bVCMgrps = TRUE;
  75         }
  76
  77         if (bFreeEnergy && PERTURBED(*atom))
  78         {
  79             md->nPerturbed++;
  80             if (atom->mB != atom->m)
  81             {
  82                 md->nMassPerturbed++;
  83             }
  84             if (atom->qB != atom->q)
  85             {
  86                 md->nChargePerturbed++;
  87             }
  88             if (atom->typeB != atom->type)
  89             {
  90                 md->nTypePerturbed++;
  91             }
  92         }
  93
  94         tmA += atom->m;
  95         tmB += atom->mB;
  96     }
  97
  98     md->tmassA = tmA;
  99     md->tmassB = tmB;
 100
 101     if (bFreeEnergy && fp)
 102     {
 103         fprintf(fp,
 104                 "There are %d atoms and %d charges for free energy perturbation\n",
 105                 md->nPerturbed, md->nChargePerturbed);
 106     }
 107
 108     md->bOrires = gmx_mtop_ftype_count(mtop, F_ORIRES);
 109
 110     return md;
 111 }
 112
 113 void atoms2md(gmx_mtop_t *mtop, t_inputrec *ir,
 114               int nindex, int *index,
 115               int start, int homenr,
 116               t_mdatoms *md)
 117 {
 118     gmx_mtop_atomlookup_t alook;
 119     int                   i;
 120     t_grpopts            *opts;
 121     real                  c6, c12;
 122     gmx_groups_t         *groups;
 123     gmx_molblock_t       *molblock;
 124
 125     opts = &ir->opts;
 126
 127     groups = &mtop->groups;
 128
 129     molblock = mtop->molblock;
 130
 131     /* Index==NULL indicates particle decomposition,
 132      * unless we have an empty DD node, so also check for homenr and start.
 133      * This should be signaled properly with an extra parameter or nindex==-1.
 134      */
 135     if (index == NULL && (homenr > 0 || start > 0))
 136     {
 137         md->nr = mtop->natoms;
 138     }
 139     else
 140     {
 141         md->nr = nindex;
 142     }
 143
 144     if (md->nr > md->nalloc)
 145     {
 146         md->nalloc = over_alloc_dd(md->nr);
 147
 148         if (md->nMassPerturbed)
 149         {
 150             srenew(md->massA, md->nalloc);
 151             srenew(md->massB, md->nalloc);
 152         }
 153         srenew(md->massT, md->nalloc);
 154         srenew(md->invmass, md->nalloc);
 155         srenew(md->chargeA, md->nalloc);
 156         srenew(md->c6A, md->nalloc);
 157         srenew(md->sigmaA, md->nalloc);
 158         srenew(md->sigma3A, md->nalloc);
 159         if (md->nPerturbed)
 160         {
 161             srenew(md->chargeB, md->nalloc);
 162             srenew(md->c6B, md->nalloc);
 163             srenew(md->sigmaB, md->nalloc);
 164             srenew(md->sigma3B, md->nalloc);
 165         }
 166         srenew(md->typeA, md->nalloc);
 167         if (md->nPerturbed)
 168         {
 169             srenew(md->typeB, md->nalloc);
 170         }
 171         srenew(md->ptype, md->nalloc);
 172         if (opts->ngtc > 1)
 173         {
 174             srenew(md->cTC, md->nalloc);
 175             /* We always copy cTC with domain decomposition */
 176         }
 177         srenew(md->cENER, md->nalloc);
 178         if (opts->ngacc > 1)
 179         {
 180             srenew(md->cACC, md->nalloc);
 181         }
 182         if (opts->nFreeze &&
 183             (opts->ngfrz > 1 ||
 184              opts->nFreeze[0][XX] || opts->nFreeze[0][YY] || opts->nFreeze[0][ZZ]))
 185         {
 186             srenew(md->cFREEZE, md->nalloc);
 187         }
 188         if (md->bVCMgrps)
 189         {
 190             srenew(md->cVCM, md->nalloc);
 191         }
 192         if (md->bOrires)
 193         {
 194             srenew(md->cORF, md->nalloc);
 195         }
 196         if (md->nPerturbed)
 197         {
 198             srenew(md->bPerturbed, md->nalloc);
 199         }
 200
 201         /* Note that these user t_mdatoms array pointers are NULL
 202          * when there is only one group present.
 203          * Therefore, when adding code, the user should use something like:
 204          * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
 205          */
 206         if (mtop->groups.grpnr[egcUser1] != NULL)
 207         {
 208             srenew(md->cU1, md->nalloc);
 209         }
 210         if (mtop->groups.grpnr[egcUser2] != NULL)
 211         {
 212             srenew(md->cU2, md->nalloc);
 213         }
 214
 215         if (ir->bQMMM)
 216         {
 217             srenew(md->bQM, md->nalloc);
 218         }
 219         if (ir->bAdress)
 220         {
 221             srenew(md->wf, md->nalloc);
 222             srenew(md->tf_table_index, md->nalloc);
 223         }
 224     }
 225
 226     alook = gmx_mtop_atomlookup_init(mtop);
 227
 228 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntDefault)) schedule(static)
 229     for (i = 0; i < md->nr; i++)
 230     {
 231         int      g, ag, molb;
 232         real     mA, mB, fac;
 233         t_atom  *atom;
 234
 235         if (index == NULL)
 236         {
 237             ag = i;
 238         }
 239         else
 240         {
 241             ag   = index[i];
 242         }
 243         gmx_mtop_atomnr_to_atom(alook, ag, &atom);
 244
 245         if (md->cFREEZE)
 246         {
 247             md->cFREEZE[i] = ggrpnr(groups, egcFREEZE, ag);
 248         }
 249         if (EI_ENERGY_MINIMIZATION(ir->eI))
 250         {
 251             /* Displacement is proportional to F, masses used for constraints */
 252             mA = 1.0;
 253             mB = 1.0;
 254         }
 255         else if (ir->eI == eiBD)
 256         {
 257             /* With BD the physical masses are irrelevant.
 258              * To keep the code simple we use most of the normal MD code path
 259              * for BD. Thus for constraining the masses should be proportional
 260              * to the friction coefficient. We set the absolute value such that
 261              * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
 262              * Then if we set the (meaningless) velocity to v=dx/dt, we get the
 263              * correct kinetic energy and temperature using the usual code path.
 264              * Thus with BD v*dt will give the displacement and the reported
 265              * temperature can signal bad integration (too large time step).
 266              */
 267             if (ir->bd_fric > 0)
 268             {
 269                 mA = 0.5*ir->bd_fric*ir->delta_t;
 270                 mB = 0.5*ir->bd_fric*ir->delta_t;
 271             }
 272             else
 273             {
 274                 /* The friction coefficient is mass/tau_t */
 275                 fac = ir->delta_t/opts->tau_t[md->cTC ? groups->grpnr[egcTC][ag] : 0];
 276                 mA  = 0.5*atom->m*fac;
 277                 mB  = 0.5*atom->mB*fac;
 278             }
 279         }
 280         else
 281         {
 282             mA = atom->m;
 283             mB = atom->mB;
 284         }
 285         if (md->nMassPerturbed)
 286         {
 287             md->massA[i]  = mA;
 288             md->massB[i]  = mB;
 289         }
 290         md->massT[i]    = mA;
 291         if (mA == 0.0)
 292         {
 293             md->invmass[i]    = 0;
 294         }
 295         else if (md->cFREEZE)
 296         {
 297             g = md->cFREEZE[i];
 298             if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] && opts->nFreeze[g][ZZ])
 299             {
 300                 /* Set the mass of completely frozen particles to ALMOST_ZERO iso 0
 301                  * to avoid div by zero in lincs or shake.
 302                  * Note that constraints can still move a partially frozen particle.
 303                  */
 304                 md->invmass[i]  = ALMOST_ZERO;
 305             }
 306             else
 307             {
 308                 md->invmass[i]  = 1.0/mA;
 309             }
 310         }
 311         else
 312         {
 313             md->invmass[i]    = 1.0/mA;
 314         }
 315         md->chargeA[i]      = atom->q;
 316         md->typeA[i]        = atom->type;
 317         c6                  = mtop->ffparams.iparams[atom->type*(mtop->ffparams.atnr+1)].lj.c6;
 318         c12                 = mtop->ffparams.iparams[atom->type*(mtop->ffparams.atnr+1)].lj.c12;
 319         md->c6A[i]          = sqrt(c6);
 320         if (c6 == 0.0 || c12 == 0.0)
 321         {
 322             md->sigmaA[i]     = 1.0;
 323         }
 324         else
 325         {
 326             md->sigmaA[i]     = pow(c12/c6, 1.0/6.0);
 327         }
 328         md->sigma3A[i] = 1/(md->sigmaA[i]*md->sigmaA[i]*md->sigmaA[i]);
 329
 330         if (md->nPerturbed)
 331         {
 332             md->chargeB[i]    = atom->qB;
 333             md->typeB[i]      = atom->typeB;
 334             c6                = mtop->ffparams.iparams[atom->typeB*(mtop->ffparams.atnr+1)].lj.c6;
 335             c12               = mtop->ffparams.iparams[atom->typeB*(mtop->ffparams.atnr+1)].lj.c12;
 336             md->c6B[i]        = sqrt(c6);
 337             if (c6 == 0.0 || c12 == 0.0)
 338             {
 339                 md->sigmaB[i]   = 1.0;
 340             }
 341             else
 342             {
 343                 md->sigmaB[i]   = pow(c12/c6, 1.0/6.0);
 344             }
 345             md->bPerturbed[i] = PERTURBED(*atom);
 346             md->sigma3B[i]    = 1/(md->sigmaB[i]*md->sigmaB[i]*md->sigmaB[i]);
 347         }
 348         md->ptype[i]    = atom->ptype;
 349         if (md->cTC)
 350         {
 351             md->cTC[i]    = groups->grpnr[egcTC][ag];
 352         }
 353         md->cENER[i]    =
 354             (groups->grpnr[egcENER] ? groups->grpnr[egcENER][ag] : 0);
 355         if (md->cACC)
 356         {
 357             md->cACC[i]   = groups->grpnr[egcACC][ag];
 358         }
 359         if (md->cVCM)
 360         {
 361             md->cVCM[i]       = groups->grpnr[egcVCM][ag];
 362         }
 363         if (md->cORF)
 364         {
 365             md->cORF[i]       = groups->grpnr[egcORFIT][ag];
 366         }
 367
 368         if (md->cU1)
 369         {
 370             md->cU1[i]        = groups->grpnr[egcUser1][ag];
 371         }
 372         if (md->cU2)
 373         {
 374             md->cU2[i]        = groups->grpnr[egcUser2][ag];
 375         }
 376
 377         if (ir->bQMMM)
 378         {
 379             if (groups->grpnr[egcQMMM] == 0 ||
 380                 groups->grpnr[egcQMMM][ag] < groups->grps[egcQMMM].nr-1)
 381             {
 382                 md->bQM[i]      = TRUE;
 383             }
 384             else
 385             {
 386                 md->bQM[i]      = FALSE;
 387             }
 388         }
 389         /* Initialize AdResS weighting functions to adressw */
 390         if (ir->bAdress)
 391         {
 392             md->wf[i]           = 1.0;
 393             /* if no tf table groups specified, use default table */
 394             md->tf_table_index[i] = DEFAULT_TF_TABLE;
 395             if (ir->adress->n_tf_grps > 0)
 396             {
 397                 /* if tf table groups specified, tf is only applied to thoose energy groups*/
 398                 md->tf_table_index[i] = NO_TF_TABLE;
 399                 /* check wether atom is in one of the relevant energy groups and assign a table index */
 400                 for (g = 0; g < ir->adress->n_tf_grps; g++)
 401                 {
 402                     if (md->cENER[i] == ir->adress->tf_table_index[g])
 403                     {
 404                         md->tf_table_index[i] = g;
 405                     }
 406                 }
 407             }
 408         }
 409     }
 410
 411     gmx_mtop_atomlookup_destroy(alook);
 412
 413     md->start  = start;
 414     md->homenr = homenr;
 415     md->lambda = 0;
 416 }
 417
 418 void update_mdatoms(t_mdatoms *md, real lambda)
 419 {
 420     int    al, end;
 421     real   L1 = 1.0-lambda;
 422
 423     end = md->nr;
 424
 425     if (md->nMassPerturbed)
 426     {
 427         for (al = 0; (al < end); al++)
 428         {
 429             if (md->bPerturbed[al])
 430             {
 431                 md->massT[al] = L1*md->massA[al]+ lambda*md->massB[al];
 432                 if (md->invmass[al] > 1.1*ALMOST_ZERO)
 433                 {
 434                     md->invmass[al] = 1.0/md->massT[al];
 435                 }
 436             }
 437         }
 438         md->tmass = L1*md->tmassA + lambda*md->tmassB;
 439     }
 440     else
 441     {
 442         md->tmass = md->tmassA;
 443     }
 444     md->lambda = lambda;
 445 }