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 "gromacs/utility/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 homenr,
 116               t_mdatoms *md)
 117 {
 118     gmx_bool              bLJPME;
 119     gmx_mtop_atomlookup_t alook;
 120     int                   i;
 121     t_grpopts            *opts;
 122     gmx_groups_t         *groups;
 123     gmx_molblock_t       *molblock;
 124
 125     bLJPME = EVDW_PME(ir->vdwtype);
 126
 127     opts = &ir->opts;
 128
 129     groups = &mtop->groups;
 130
 131     molblock = mtop->molblock;
 132
 133     /* Index==NULL indicates no DD (unless we have a DD node with no
 134      * atoms), so also check for homenr. This should be
 135      * signaled properly with an extra parameter or nindex==-1.
 136      */
 137     if (index == NULL && (homenr > 0))
 138     {
 139         md->nr = mtop->natoms;
 140     }
 141     else
 142     {
 143         md->nr = nindex;
 144     }
 145
 146     if (md->nr > md->nalloc)
 147     {
 148         md->nalloc = over_alloc_dd(md->nr);
 149
 150         if (md->nMassPerturbed)
 151         {
 152             srenew(md->massA, md->nalloc);
 153             srenew(md->massB, md->nalloc);
 154         }
 155         srenew(md->massT, md->nalloc);
 156         srenew(md->invmass, md->nalloc);
 157         srenew(md->chargeA, md->nalloc);
 158         if (bLJPME)
 159         {
 160             srenew(md->sqrt_c6A, md->nalloc);
 161             srenew(md->sigmaA, md->nalloc);
 162             srenew(md->sigma3A, md->nalloc);
 163         }
 164         if (md->nPerturbed)
 165         {
 166             srenew(md->chargeB, md->nalloc);
 167             if (bLJPME)
 168             {
 169                 srenew(md->sqrt_c6B, md->nalloc);
 170                 srenew(md->sigmaB, md->nalloc);
 171                 srenew(md->sigma3B, md->nalloc);
 172             }
 173         }
 174         srenew(md->typeA, md->nalloc);
 175         if (md->nPerturbed)
 176         {
 177             srenew(md->typeB, md->nalloc);
 178         }
 179         srenew(md->ptype, md->nalloc);
 180         if (opts->ngtc > 1)
 181         {
 182             srenew(md->cTC, md->nalloc);
 183             /* We always copy cTC with domain decomposition */
 184         }
 185         srenew(md->cENER, md->nalloc);
 186         if (opts->ngacc > 1)
 187         {
 188             srenew(md->cACC, md->nalloc);
 189         }
 190         if (opts->nFreeze &&
 191             (opts->ngfrz > 1 ||
 192              opts->nFreeze[0][XX] || opts->nFreeze[0][YY] || opts->nFreeze[0][ZZ]))
 193         {
 194             srenew(md->cFREEZE, md->nalloc);
 195         }
 196         if (md->bVCMgrps)
 197         {
 198             srenew(md->cVCM, md->nalloc);
 199         }
 200         if (md->bOrires)
 201         {
 202             srenew(md->cORF, md->nalloc);
 203         }
 204         if (md->nPerturbed)
 205         {
 206             srenew(md->bPerturbed, md->nalloc);
 207         }
 208
 209         /* Note that these user t_mdatoms array pointers are NULL
 210          * when there is only one group present.
 211          * Therefore, when adding code, the user should use something like:
 212          * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
 213          */
 214         if (mtop->groups.grpnr[egcUser1] != NULL)
 215         {
 216             srenew(md->cU1, md->nalloc);
 217         }
 218         if (mtop->groups.grpnr[egcUser2] != NULL)
 219         {
 220             srenew(md->cU2, md->nalloc);
 221         }
 222
 223         if (ir->bQMMM)
 224         {
 225             srenew(md->bQM, md->nalloc);
 226         }
 227         if (ir->bAdress)
 228         {
 229             srenew(md->wf, md->nalloc);
 230             srenew(md->tf_table_index, md->nalloc);
 231         }
 232     }
 233
 234     alook = gmx_mtop_atomlookup_init(mtop);
 235
 236 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntDefault)) schedule(static)
 237     for (i = 0; i < md->nr; i++)
 238     {
 239         int      g, ag, molb;
 240         real     mA, mB, fac;
 241         real     c6, c12;
 242         t_atom  *atom;
 243
 244         if (index == NULL)
 245         {
 246             ag = i;
 247         }
 248         else
 249         {
 250             ag   = index[i];
 251         }
 252         gmx_mtop_atomnr_to_atom(alook, ag, &atom);
 253
 254         if (md->cFREEZE)
 255         {
 256             md->cFREEZE[i] = ggrpnr(groups, egcFREEZE, ag);
 257         }
 258         if (EI_ENERGY_MINIMIZATION(ir->eI))
 259         {
 260             /* Displacement is proportional to F, masses used for constraints */
 261             mA = 1.0;
 262             mB = 1.0;
 263         }
 264         else if (ir->eI == eiBD)
 265         {
 266             /* With BD the physical masses are irrelevant.
 267              * To keep the code simple we use most of the normal MD code path
 268              * for BD. Thus for constraining the masses should be proportional
 269              * to the friction coefficient. We set the absolute value such that
 270              * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
 271              * Then if we set the (meaningless) velocity to v=dx/dt, we get the
 272              * correct kinetic energy and temperature using the usual code path.
 273              * Thus with BD v*dt will give the displacement and the reported
 274              * temperature can signal bad integration (too large time step).
 275              */
 276             if (ir->bd_fric > 0)
 277             {
 278                 mA = 0.5*ir->bd_fric*ir->delta_t;
 279                 mB = 0.5*ir->bd_fric*ir->delta_t;
 280             }
 281             else
 282             {
 283                 /* The friction coefficient is mass/tau_t */
 284                 fac = ir->delta_t/opts->tau_t[md->cTC ? groups->grpnr[egcTC][ag] : 0];
 285                 mA  = 0.5*atom->m*fac;
 286                 mB  = 0.5*atom->mB*fac;
 287             }
 288         }
 289         else
 290         {
 291             mA = atom->m;
 292             mB = atom->mB;
 293         }
 294         if (md->nMassPerturbed)
 295         {
 296             md->massA[i]  = mA;
 297             md->massB[i]  = mB;
 298         }
 299         md->massT[i]    = mA;
 300         if (mA == 0.0)
 301         {
 302             md->invmass[i]    = 0;
 303         }
 304         else if (md->cFREEZE)
 305         {
 306             g = md->cFREEZE[i];
 307             if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] && opts->nFreeze[g][ZZ])
 308             {
 309                 /* Set the mass of completely frozen particles to ALMOST_ZERO iso 0
 310                  * to avoid div by zero in lincs or shake.
 311                  * Note that constraints can still move a partially frozen particle.
 312                  */
 313                 md->invmass[i]  = ALMOST_ZERO;
 314             }
 315             else
 316             {
 317                 md->invmass[i]  = 1.0/mA;
 318             }
 319         }
 320         else
 321         {
 322             md->invmass[i]    = 1.0/mA;
 323         }
 324         md->chargeA[i]      = atom->q;
 325         md->typeA[i]        = atom->type;
 326         if (bLJPME)
 327         {
 328             c6                = mtop->ffparams.iparams[atom->type*(mtop->ffparams.atnr+1)].lj.c6;
 329             c12               = mtop->ffparams.iparams[atom->type*(mtop->ffparams.atnr+1)].lj.c12;
 330             md->sqrt_c6A[i]   = sqrt(c6);
 331             if (c6 == 0.0 || c12 == 0)
 332             {
 333                 md->sigmaA[i] = 1.0;
 334             }
 335             else
 336             {
 337                 md->sigmaA[i] = pow(c12/c6, 1.0/6.0);
 338             }
 339             md->sigma3A[i]    = 1/(md->sigmaA[i]*md->sigmaA[i]*md->sigmaA[i]);
 340         }
 341         if (md->nPerturbed)
 342         {
 343             md->bPerturbed[i] = PERTURBED(*atom);
 344             md->chargeB[i]    = atom->qB;
 345             md->typeB[i]      = atom->typeB;
 346             if (bLJPME)
 347             {
 348                 c6                = mtop->ffparams.iparams[atom->typeB*(mtop->ffparams.atnr+1)].lj.c6;
 349                 c12               = mtop->ffparams.iparams[atom->typeB*(mtop->ffparams.atnr+1)].lj.c12;
 350                 md->sqrt_c6B[i]   = sqrt(c6);
 351                 if (c6 == 0.0 || c12 == 0)
 352                 {
 353                     md->sigmaB[i] = 1.0;
 354                 }
 355                 else
 356                 {
 357                     md->sigmaB[i] = pow(c12/c6, 1.0/6.0);
 358                 }
 359                 md->sigma3B[i]    = 1/(md->sigmaB[i]*md->sigmaB[i]*md->sigmaB[i]);
 360             }
 361         }
 362         md->ptype[i]    = atom->ptype;
 363         if (md->cTC)
 364         {
 365             md->cTC[i]    = groups->grpnr[egcTC][ag];
 366         }
 367         md->cENER[i]    =
 368             (groups->grpnr[egcENER] ? groups->grpnr[egcENER][ag] : 0);
 369         if (md->cACC)
 370         {
 371             md->cACC[i]   = groups->grpnr[egcACC][ag];
 372         }
 373         if (md->cVCM)
 374         {
 375             md->cVCM[i]       = groups->grpnr[egcVCM][ag];
 376         }
 377         if (md->cORF)
 378         {
 379             md->cORF[i]       = groups->grpnr[egcORFIT][ag];
 380         }
 381
 382         if (md->cU1)
 383         {
 384             md->cU1[i]        = groups->grpnr[egcUser1][ag];
 385         }
 386         if (md->cU2)
 387         {
 388             md->cU2[i]        = groups->grpnr[egcUser2][ag];
 389         }
 390
 391         if (ir->bQMMM)
 392         {
 393             if (groups->grpnr[egcQMMM] == 0 ||
 394                 groups->grpnr[egcQMMM][ag] < groups->grps[egcQMMM].nr-1)
 395             {
 396                 md->bQM[i]      = TRUE;
 397             }
 398             else
 399             {
 400                 md->bQM[i]      = FALSE;
 401             }
 402         }
 403         /* Initialize AdResS weighting functions to adressw */
 404         if (ir->bAdress)
 405         {
 406             md->wf[i]           = 1.0;
 407             /* if no tf table groups specified, use default table */
 408             md->tf_table_index[i] = DEFAULT_TF_TABLE;
 409             if (ir->adress->n_tf_grps > 0)
 410             {
 411                 /* if tf table groups specified, tf is only applied to thoose energy groups*/
 412                 md->tf_table_index[i] = NO_TF_TABLE;
 413                 /* check wether atom is in one of the relevant energy groups and assign a table index */
 414                 for (g = 0; g < ir->adress->n_tf_grps; g++)
 415                 {
 416                     if (md->cENER[i] == ir->adress->tf_table_index[g])
 417                     {
 418                         md->tf_table_index[i] = g;
 419                     }
 420                 }
 421             }
 422         }
 423     }
 424
 425     gmx_mtop_atomlookup_destroy(alook);
 426
 427     md->homenr = homenr;
 428     md->lambda = 0;
 429 }
 430
 431 void update_mdatoms(t_mdatoms *md, real lambda)
 432 {
 433     int    al, end;
 434     real   L1 = 1.0-lambda;
 435
 436     end = md->nr;
 437
 438     if (md->nMassPerturbed)
 439     {
 440         for (al = 0; (al < end); al++)
 441         {
 442             if (md->bPerturbed[al])
 443             {
 444                 md->massT[al] = L1*md->massA[al]+ lambda*md->massB[al];
 445                 if (md->invmass[al] > 1.1*ALMOST_ZERO)
 446                 {
 447                     md->invmass[al] = 1.0/md->massT[al];
 448                 }
 449             }
 450         }
 451         md->tmass = L1*md->tmassA + lambda*md->tmassB;
 452     }
 453     else
 454     {
 455         md->tmass = md->tmassA;
 456     }
 457     md->lambda = lambda;
 458 }