src/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  * 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 "typedefs.h"
  43 #include "mdatoms.h"
  44 #include "smalloc.h"
  45 #include "main.h"
  46 #include "qmmm.h"
  47 #include "mtop_util.h"
  48 #include "gmx_omp_nthreads.h"
  49
  50 #define ALMOST_ZERO 1e-30
  51
  52 t_mdatoms *init_mdatoms(FILE *fp,gmx_mtop_t *mtop,gmx_bool bFreeEnergy)
  53 {
  54   int    mb,a,g,nmol;
  55   double tmA,tmB;
  56   t_atom *atom;
  57   t_mdatoms *md;
  58   gmx_mtop_atomloop_all_t aloop;
  59   t_ilist *ilist;
  60
  61   snew(md,1);
  62
  63   md->nenergrp = mtop->groups.grps[egcENER].nr;
  64   md->bVCMgrps = FALSE;
  65   tmA = 0.0;
  66   tmB = 0.0;
  67
  68   aloop = gmx_mtop_atomloop_all_init(mtop);
  69   while(gmx_mtop_atomloop_all_next(aloop,&a,&atom)) {
  70     if (ggrpnr(&mtop->groups,egcVCM,a) > 0)
  71       md->bVCMgrps = TRUE;
  72
  73     if (bFreeEnergy && PERTURBED(*atom)) {
  74       md->nPerturbed++;
  75       if (atom->mB != atom->m)
  76         md->nMassPerturbed++;
  77       if (atom->qB != atom->q)
  78         md->nChargePerturbed++;
  79     }
  80
  81     tmA += atom->m;
  82     tmB += atom->mB;
  83   }
  84
  85   md->tmassA = tmA;
  86   md->tmassB = tmB;
  87
  88   if (bFreeEnergy && fp)
  89     fprintf(fp,
  90             "There are %d atoms and %d charges for free energy perturbation\n",
  91             md->nPerturbed,md->nChargePerturbed);
  92
  93   md->bOrires = gmx_mtop_ftype_count(mtop,F_ORIRES);
  94
  95   return md;
  96 }
  97
  98 void atoms2md(gmx_mtop_t *mtop,t_inputrec *ir,
  99               int nindex,int *index,
 100               int start,int homenr,
 101               t_mdatoms *md)
 102 {
 103   gmx_mtop_atomlookup_t alook;
 104   int       i;
 105   t_grpopts *opts;
 106   gmx_groups_t *groups;
 107   gmx_molblock_t *molblock;
 108
 109   opts = &ir->opts;
 110
 111   groups = &mtop->groups;
 112
 113   molblock = mtop->molblock;
 114
 115   /* Index==NULL indicates particle decomposition,
 116    * unless we have an empty DD node, so also check for homenr and start.
 117    * This should be signaled properly with an extra parameter or nindex==-1.
 118    */
 119   if (index == NULL && (homenr > 0 || start > 0)) {
 120     md->nr = mtop->natoms;
 121   } else {
 122     md->nr = nindex;
 123   }
 124
 125   if (md->nr > md->nalloc) {
 126     md->nalloc = over_alloc_dd(md->nr);
 127
 128     if (md->nMassPerturbed) {
 129       srenew(md->massA,md->nalloc);
 130       srenew(md->massB,md->nalloc);
 131     }
 132     srenew(md->massT,md->nalloc);
 133     srenew(md->invmass,md->nalloc);
 134     srenew(md->chargeA,md->nalloc);
 135     if (md->nPerturbed) {
 136       srenew(md->chargeB,md->nalloc);
 137     }
 138     srenew(md->typeA,md->nalloc);
 139     if (md->nPerturbed) {
 140       srenew(md->typeB,md->nalloc);
 141     }
 142     srenew(md->ptype,md->nalloc);
 143     if (opts->ngtc > 1) {
 144       srenew(md->cTC,md->nalloc);
 145       /* We always copy cTC with domain decomposition */
 146     }
 147     srenew(md->cENER,md->nalloc);
 148     if (opts->ngacc > 1)
 149       srenew(md->cACC,md->nalloc);
 150     if (opts->nFreeze &&
 151         (opts->ngfrz > 1 ||
 152          opts->nFreeze[0][XX] || opts->nFreeze[0][YY] || opts->nFreeze[0][ZZ]))
 153       srenew(md->cFREEZE,md->nalloc);
 154     if (md->bVCMgrps)
 155       srenew(md->cVCM,md->nalloc);
 156     if (md->bOrires)
 157       srenew(md->cORF,md->nalloc);
 158     if (md->nPerturbed)
 159       srenew(md->bPerturbed,md->nalloc);
 160
 161     /* Note that these user t_mdatoms array pointers are NULL
 162      * when there is only one group present.
 163      * Therefore, when adding code, the user should use something like:
 164      * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
 165      */
 166     if (mtop->groups.grpnr[egcUser1] != NULL)
 167       srenew(md->cU1,md->nalloc);
 168     if (mtop->groups.grpnr[egcUser2] != NULL)
 169       srenew(md->cU2,md->nalloc);
 170
 171     if (ir->bQMMM)
 172       srenew(md->bQM,md->nalloc);
 173     if (ir->bAdress)
 174       srenew(md->wf,md->nalloc);
 175       srenew(md->tf_table_index,md->nalloc);
 176
 177       md->purecg = FALSE;
 178       md->pureex = FALSE;
 179   }
 180
 181   alook = gmx_mtop_atomlookup_init(mtop);
 182
 183 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntDefault)) schedule(static)
 184   for(i=0; i<md->nr; i++) {
 185     int     g,ag,molb;
 186     real    mA,mB,fac;
 187     t_atom  *atom;
 188
 189     if (index == NULL) {
 190       ag = i;
 191     } else {
 192       ag   = index[i];
 193     }
 194     gmx_mtop_atomnr_to_atom(alook,ag,&atom);
 195
 196     if (md->cFREEZE) {
 197       md->cFREEZE[i] = ggrpnr(groups,egcFREEZE,ag);
 198     }
 199         if (EI_ENERGY_MINIMIZATION(ir->eI))
 200         {
 201             /* Displacement is proportional to F, masses used for constraints */
 202             mA = 1.0;
 203             mB = 1.0;
 204         }
 205         else if (ir->eI == eiBD)
 206         {
 207             /* With BD the physical masses are irrelevant.
 208              * To keep the code simple we use most of the normal MD code path
 209              * for BD. Thus for constraining the masses should be proportional
 210              * to the friction coefficient. We set the absolute value such that
 211              * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
 212              * Then if we set the (meaningless) velocity to v=dx/dt, we get the
 213              * correct kinetic energy and temperature using the usual code path.
 214              * Thus with BD v*dt will give the displacement and the reported
 215              * temperature can signal bad integration (too large time step).
 216              */
 217             if (ir->bd_fric > 0)
 218             {
 219                 mA = 0.5*ir->bd_fric*ir->delta_t;
 220                 mB = 0.5*ir->bd_fric*ir->delta_t;
 221             }
 222             else
 223             {
 224                 /* The friction coefficient is mass/tau_t */
 225                 fac = ir->delta_t/opts->tau_t[md->cTC ? groups->grpnr[egcTC][ag] : 0];
 226                 mA = 0.5*atom->m*fac;
 227                 mB = 0.5*atom->mB*fac;
 228             }
 229         }
 230         else
 231         {
 232             mA = atom->m;
 233             mB = atom->mB;
 234         }
 235     if (md->nMassPerturbed) {
 236       md->massA[i]      = mA;
 237       md->massB[i]      = mB;
 238     }
 239     md->massT[i]        = mA;
 240     if (mA == 0.0) {
 241       md->invmass[i]    = 0;
 242     } else if (md->cFREEZE) {
 243       g = md->cFREEZE[i];
 244       if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] && opts->nFreeze[g][ZZ])
 245         /* Set the mass of completely frozen particles to ALMOST_ZERO iso 0
 246          * to avoid div by zero in lincs or shake.
 247          * Note that constraints can still move a partially frozen particle.
 248          */
 249         md->invmass[i]  = ALMOST_ZERO;
 250       else
 251         md->invmass[i]  = 1.0/mA;
 252     } else {
 253       md->invmass[i]    = 1.0/mA;
 254     }
 255     md->chargeA[i]      = atom->q;
 256     md->typeA[i]        = atom->type;
 257     if (md->nPerturbed) {
 258       md->chargeB[i]    = atom->qB;
 259       md->typeB[i]      = atom->typeB;
 260       md->bPerturbed[i] = PERTURBED(*atom);
 261     }
 262     md->ptype[i]        = atom->ptype;
 263     if (md->cTC)
 264       md->cTC[i]        = groups->grpnr[egcTC][ag];
 265     md->cENER[i]        =
 266       (groups->grpnr[egcENER] ? groups->grpnr[egcENER][ag] : 0);
 267     if (md->cACC)
 268       md->cACC[i]       = groups->grpnr[egcACC][ag];
 269     if (md->cVCM)
 270       md->cVCM[i]       = groups->grpnr[egcVCM][ag];
 271     if (md->cORF)
 272       md->cORF[i]       = groups->grpnr[egcORFIT][ag];
 273
 274     if (md->cU1)
 275       md->cU1[i]        = groups->grpnr[egcUser1][ag];
 276     if (md->cU2)
 277       md->cU2[i]        = groups->grpnr[egcUser2][ag];
 278
 279     if (ir->bQMMM) {
 280       if (groups->grpnr[egcQMMM] == 0 ||
 281           groups->grpnr[egcQMMM][ag] < groups->grps[egcQMMM].nr-1) {
 282         md->bQM[i]      = TRUE;
 283       } else {
 284         md->bQM[i]      = FALSE;
 285       }
 286     }
 287     /* Initialize AdResS weighting functions to adressw */
 288     if (ir->bAdress){
 289        md->wf[i]           = 1.0;
 290         /* if no tf table groups specified, use default table */
 291        md->tf_table_index[i] = DEFAULT_TF_TABLE;
 292        if (ir->adress->n_tf_grps > 0){
 293             /* if tf table groups specified, tf is only applied to thoose energy groups*/
 294             md->tf_table_index[i] = NO_TF_TABLE;
 295             /* check wether atom is in one of the relevant energy groups and assign a table index */
 296             for (g=0; g<ir->adress->n_tf_grps; g++){
 297                 if (md->cENER[i] == ir->adress->tf_table_index[g]){
 298                    md->tf_table_index[i] = g;
 299                 }
 300             }
 301         }
 302     }
 303   }
 304
 305   gmx_mtop_atomlookup_destroy(alook);
 306
 307   md->start  = start;
 308   md->homenr = homenr;
 309   md->lambda = 0;
 310 }
 311
 312 void update_mdatoms(t_mdatoms *md,real lambda)
 313 {
 314   int    al,end;
 315   real   L1=1.0-lambda;
 316
 317   end=md->nr;
 318
 319   if (md->nMassPerturbed) {
 320     for(al=0; (al<end); al++) {
 321       if (md->bPerturbed[al]) {
 322         md->massT[al] = L1*md->massA[al]+ lambda*md->massB[al];
 323         if (md->invmass[al] > 1.1*ALMOST_ZERO)
 324           md->invmass[al] = 1.0/md->massT[al];
 325       }
 326     }
 327     md->tmass = L1*md->tmassA + lambda*md->tmassB;
 328   } else {
 329     md->tmass = md->tmassA;
 330   }
 331   md->lambda = lambda;
 332 }