src/mdlib/mdatom.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Getting the Right Output Means no Artefacts in Calculating Stuff
  31  */
  32 static char *SRCID_mdatom_c = "$Id$";
  33 #include "typedefs.h"
  34 #include "mdatoms.h"
  35 #include "smalloc.h"
  36 #include "main.h"
  37
  38 #define ALMOST_ZERO 1e-30
  39
  40 t_mdatoms *atoms2md(FILE *fp,t_atoms *atoms,ivec nFreeze[],
  41                     bool bBD,real delta_t,real fric,real tau_t[],
  42                     bool bPert,bool bFree)
  43 {
  44   int       i,np,g;
  45   real      fac;
  46   double    tm;
  47   t_mdatoms *md;
  48
  49   snew(md,1);
  50   md->nr = atoms->nr;
  51   snew(md->massA,md->nr);
  52   snew(md->massB,md->nr);
  53   snew(md->massT,md->nr);
  54   snew(md->invmass,md->nr);
  55   snew(md->chargeA,md->nr);
  56   snew(md->chargeB,md->nr);
  57   snew(md->chargeT,md->nr);
  58   snew(md->resnr,md->nr);
  59   snew(md->typeA,md->nr);
  60   snew(md->typeB,md->nr);
  61   snew(md->ptype,md->nr);
  62   snew(md->cTC,md->nr);
  63   snew(md->cENER,md->nr);
  64   snew(md->cACC,md->nr);
  65   snew(md->cFREEZE,md->nr);
  66   snew(md->cXTC,md->nr);
  67   snew(md->cVCM,md->nr);
  68   snew(md->cORF,md->nr);
  69   snew(md->bPerturbed,md->nr);
  70
  71   snew(md->cU1,md->nr);
  72   snew(md->cU2,md->nr);
  73
  74   np=0;
  75   tm=0.0;
  76   for(i=0; (i<md->nr); i++) {
  77     if (bBD) {
  78       /* Make the mass proportional to the friction coefficient for BD.
  79        * This is necessary for the constraint algorithms.
  80        */
  81       if (fric) {
  82         md->massA[i]    = fric*delta_t;
  83         md->massB[i]    = fric*delta_t;
  84       } else {
  85         fac = delta_t/tau_t[atoms->atom[i].grpnr[egcTC]];
  86         md->massA[i]    = atoms->atom[i].m*fac;
  87         md->massB[i]    = atoms->atom[i].mB*fac;
  88       }
  89     } else {
  90       md->massA[i]      = atoms->atom[i].m;
  91       md->massB[i]      = atoms->atom[i].mB;
  92     }
  93     md->massT[i]        = md->massA[i];
  94     md->chargeA[i]      = atoms->atom[i].q;
  95     md->chargeB[i]      = atoms->atom[i].qB;
  96     md->resnr[i]        = atoms->atom[i].resnr;
  97     md->typeA[i]        = atoms->atom[i].type;
  98     md->typeB[i]        = atoms->atom[i].typeB;
  99     md->ptype[i]        = atoms->atom[i].ptype;
 100     md->cTC[i]          = atoms->atom[i].grpnr[egcTC];
 101     md->cENER[i]        = atoms->atom[i].grpnr[egcENER];
 102     md->cACC[i]         = atoms->atom[i].grpnr[egcACC];
 103     md->cFREEZE[i]      = atoms->atom[i].grpnr[egcFREEZE];
 104     md->cXTC[i]         = atoms->atom[i].grpnr[egcXTC];
 105     md->cVCM[i]         = atoms->atom[i].grpnr[egcVCM];
 106     md->cORF[i]         = atoms->atom[i].grpnr[egcORFIT];
 107     if (md->massA[i] != 0.0) {
 108       tm               += md->massT[i];
 109       g = md->cFREEZE[i];
 110       if (nFreeze[g][XX] && nFreeze[g][YY] && nFreeze[g][ZZ])
 111         /* Set the mass of completely frozen particles to ALMOST_ZERO iso 0
 112            to avoid div by zero in lincs or shake.
 113            Note that constraints can still move a partially frozen particle. */
 114         md->invmass[i]  = ALMOST_ZERO;
 115       else if (md->massT[i] == 0)
 116         md->invmass[i]  = 0;
 117       else
 118         md->invmass[i]  = 1.0/md->massT[i];
 119     }
 120     if (bPert) {
 121       md->bPerturbed[i]   = PERTURBED(atoms->atom[i]);
 122       if (md->bPerturbed[i])
 123         np++;
 124     }
 125
 126     md->cU1[i]          = atoms->atom[i].grpnr[egcUser1];
 127     md->cU2[i]          = atoms->atom[i].grpnr[egcUser2];
 128   }
 129   md->tmass  = tm;
 130
 131   if (bFree) {
 132     sfree(atoms->atom);
 133     atoms->atom=NULL;
 134   }
 135
 136   if (fp)
 137     fprintf(fp,"There are %d atoms for free energy perturbation\n",np);
 138
 139   return md;
 140 }
 141
 142 void md2atoms(t_mdatoms *md,t_atoms *atoms,bool bFree)
 143 {
 144   int i;
 145
 146   snew(atoms->atom,md->nr);
 147   for(i=0; (i<md->nr); i++) {
 148     atoms->atom[i].m                = md->massT[i];
 149     atoms->atom[i].q                = md->chargeT[i];
 150     atoms->atom[i].resnr            = md->resnr[i];
 151     atoms->atom[i].type             = md->typeA[i];
 152     atoms->atom[i].ptype            = md->ptype[i];
 153     atoms->atom[i].grpnr[egcTC]     = md->cTC[i];
 154     atoms->atom[i].grpnr[egcENER]   = md->cENER[i];
 155     atoms->atom[i].grpnr[egcACC]    = md->cACC[i];
 156     atoms->atom[i].grpnr[egcFREEZE] = md->cFREEZE[i];
 157     atoms->atom[i].grpnr[egcVCM]    = md->cVCM[i];
 158     atoms->atom[i].grpnr[egcXTC]    = md->cXTC[i];
 159     atoms->atom[i].grpnr[egcORFIT]  = md->cORF[i];
 160
 161     atoms->atom[i].grpnr[egcUser1]  = md->cU1[i];
 162     atoms->atom[i].grpnr[egcUser2]  = md->cU2[i];
 163
 164   }
 165   if (bFree) {
 166     sfree(md->massA);
 167     sfree(md->massB);
 168     sfree(md->massT);
 169     sfree(md->invmass);
 170     sfree(md->chargeA);
 171     sfree(md->chargeB);
 172     sfree(md->chargeT);
 173     sfree(md->resnr);
 174     sfree(md->typeA);
 175     sfree(md->typeB);
 176     sfree(md->ptype);
 177     sfree(md->cTC);
 178     sfree(md->cENER);
 179     sfree(md->cACC);
 180     sfree(md->cFREEZE);
 181     sfree(md->cVCM);
 182     sfree(md->cXTC);
 183     sfree(md->cORF);
 184
 185     sfree(md->cU1);
 186     sfree(md->cU2);
 187   }
 188 }
 189
 190 void init_mdatoms(t_mdatoms *md,real lambda,bool bFirst)
 191 {
 192   static real lambda0;
 193   int    i,end;
 194   real   L1=1.0-lambda;
 195
 196   if (bFirst)
 197     lambda0 = lambda;
 198   end=md->nr;
 199
 200   /* Only do this loop the first time, or when lambda has changed.
 201    * One could also check whether there is any perturbed atom at all,
 202    * but if you don't have perturbed atoms, it does not make sense to modify lambda.
 203    * In principle this has to be parallellized, although it would mean extra
 204    * communication. Basically only the charges are used on other nodes...
 205    */
 206   if (bFirst || (lambda0 != lambda)) {
 207     for(i=0; (i<end); i++) {
 208       if (md->bPerturbed[i] || bFirst) {
 209         md->massT[i]=L1*md->massA[i]+lambda*md->massB[i];
 210         if (md->invmass[i] > 1.1*ALMOST_ZERO)
 211           md->invmass[i]=1.0/md->massT[i];
 212         md->chargeT[i]=L1*md->chargeA[i]+lambda*md->chargeB[i];
 213       }
 214     }
 215   }
 216   lambda0 = lambda;
 217 }
 218
 219
 220
 221