src/mdlib/coupling.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_coupling_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "smalloc.h"
  33 #include "update.h"
  34 #include "vec.h"
  35 #include "macros.h"
  36 #include "physics.h"
  37 #include "names.h"
  38 #include "fatal.h"
  39 #include "txtdump.h"
  40 #include "nrnb.h"
  41
  42 /*
  43  * This file implements temperature and pressure coupling algorithms:
  44  * For now only the Weak coupling and the modified weak coupling.
  45  *
  46  * Furthermore computation of pressure and temperature is done here
  47  *
  48  */
  49
  50 void calc_pres(int eBox,matrix box,tensor ekin,tensor vir,tensor pres,real Elr)
  51 {
  52   int  n,m;
  53   real fac,Plr;
  54
  55   if (eBox == ebtNONE)
  56     clear_mat(pres);
  57   else {
  58     /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
  59      * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
  60      * het systeem...
  61      */
  62
  63     /* Long range correction for periodic systems, see
  64      * Neumann et al. JCP
  65      * divide by 6 because it is multiplied by fac later on.
  66      * If Elr = 0, no correction is made.
  67      */
  68
  69     /* This formula should not be used with Ewald or PME,
  70      * where the full long-range virial is calculated. EL 990823
  71      */
  72     Plr = Elr/6.0;
  73
  74     fac=PRESFAC*2.0/det(box);
  75     for(n=0; (n<DIM); n++)
  76       for(m=0; (m<DIM); m++)
  77         pres[n][m]=(ekin[n][m]-vir[n][m]+Plr)*fac;
  78
  79     if (debug) {
  80       pr_rvecs(debug,0,"PC: pres",pres,DIM);
  81       pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
  82       pr_rvecs(debug,0,"PC: vir ",vir, DIM);
  83       pr_rvecs(debug,0,"PC: box ",box, DIM);
  84     }
  85   }
  86 }
  87
  88 real calc_temp(real ekin,int nrdf)
  89 {
  90   return (2.0*ekin)/(nrdf*BOLTZ);
  91 }
  92
  93 real run_aver(real old,real cur,int step,int nmem)
  94 {
  95   nmem   = max(1,nmem);
  96
  97   return ((nmem-1)*old+cur)/nmem;
  98 }
  99
 100
 101 void do_pcoupl(t_inputrec *ir,int step,tensor pres,
 102                matrix box,int start,int nr_atoms,
 103                rvec x[],unsigned short cFREEZE[],
 104                t_nrnb *nrnb,ivec nFreeze[])
 105 {
 106   static bool bFirst=TRUE;
 107   static rvec PPP;
 108   int    n,d,m,g,ncoupl=0;
 109   real   scalar_pressure, xy_pressure, p_corr_z;
 110   real   X,Y,Z;
 111   rvec   factor;
 112   tensor mu;
 113   real   muxx,muxy,muxz,muyx,muyy,muyz,muzx,muzy,muzz;
 114   real   fgx,fgy,fgz;
 115
 116   /*
 117    *  PRESSURE SCALING
 118    *  Step (2P)
 119    */
 120   if (bFirst) {
 121     /* Initiate the pressure to the reference one */
 122     for(m=0; m<DIM; m++)
 123       PPP[m] = ir->ref_p[m];
 124     bFirst=FALSE;
 125   }
 126   scalar_pressure=0;
 127   xy_pressure=0;
 128   for(m=0; m<DIM; m++) {
 129     PPP[m]           = run_aver(PPP[m],pres[m][m],step,ir->npcmemory);
 130     scalar_pressure += PPP[m]/DIM;
 131     if (m != ZZ)
 132       xy_pressure += PPP[m]/(DIM-1);
 133   }
 134
 135   /* Pressure is now in bar, everywhere. */
 136   if ((ir->epc != epcNO) && (scalar_pressure != 0.0)) {
 137     for(m=0; m<DIM; m++)
 138       factor[m] = ir->compress[m]*ir->delta_t/ir->tau_p;
 139     clear_mat(mu);
 140     switch (ir->epc) {
 141     case epcISOTROPIC:
 142       for(m=0; m<DIM; m++)
 143         mu[m][m] = pow(1.0-factor[m]*(ir->ref_p[m]-scalar_pressure),1.0/DIM);
 144       break;
 145     case epcSEMIISOTROPIC:
 146       for(m=0; m<ZZ; m++)
 147         mu[m][m] = pow(1.0-factor[m]*(ir->ref_p[m]-xy_pressure),1.0/DIM);
 148       mu[ZZ][ZZ] = pow(1.0-factor[ZZ]*(ir->ref_p[ZZ] - PPP[ZZ]),1.0/DIM);
 149       break;
 150     case epcANISOTROPIC:
 151       for (m=0; m<DIM; m++)
 152         mu[m][m] = pow(1.0-factor[m]*(ir->ref_p[m] - PPP[m]),1.0/DIM);
 153       break;
 154     case epcSURFACETENSION:
 155       /* ir->ref_p[0/1] is the reference surface-tension times *
 156        * the number of surfaces                                */
 157       if (ir->compress[ZZ])
 158         p_corr_z = ir->delta_t/ir->tau_p*(ir->ref_p[ZZ] - PPP[ZZ]);
 159       else
 160         /* when the compressibity is zero, set the pressure correction   *
 161          * in the z-direction to zero to get the correct surface tension */
 162         p_corr_z = 0;
 163       mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ]*p_corr_z;
 164       for(m=0; m<ZZ; m++)
 165         mu[m][m] = sqrt(1.0+factor[m]*(ir->ref_p[m]/(mu[ZZ][ZZ]*box[ZZ][ZZ]) -
 166         (PPP[ZZ]+p_corr_z - xy_pressure)));
 167       break;
 168     case epcTRICLINIC:
 169     default:
 170       fatal_error(0,"Pressure coupling type %s not supported yet\n",
 171                   EPCOUPLTYPE(ir->epc));
 172     }
 173     if (debug) {
 174       pr_rvecs(debug,0,"PC: PPP ",&PPP,1);
 175       pr_rvecs(debug,0,"PC: fac ",&factor,1);
 176       pr_rvecs(debug,0,"PC: mu  ",mu,DIM);
 177     }
 178     /* Scale the positions using matrix operation */
 179     nr_atoms+=start;
 180     muxx=mu[XX][XX],muxy=mu[XX][YY],muxz=mu[XX][ZZ];
 181     muyx=mu[YY][XX],muyy=mu[YY][YY],muyz=mu[YY][ZZ];
 182     muzx=mu[ZZ][XX],muzy=mu[ZZ][YY],muzz=mu[ZZ][ZZ];
 183     for (n=start; n<nr_atoms; n++) {
 184       g=cFREEZE[n];
 185
 186       X=x[n][XX];
 187       Y=x[n][YY];
 188       Z=x[n][ZZ];
 189
 190       if (!nFreeze[g][XX])
 191         x[n][XX] = muxx*X+muxy*Y+muxz*Z;
 192       if (!nFreeze[g][YY])
 193         x[n][YY] = muyx*X+muyy*Y+muyz*Z;
 194       if (!nFreeze[g][ZZ])
 195         x[n][ZZ] = muzx*X+muzy*Y+muzz*Z;
 196
 197       ncoupl++;
 198     }
 199     /* compute final boxlengths */
 200     for (d=0; d<DIM; d++)
 201       for (m=0; m<DIM; m++)
 202         box[d][m] *= mu[d][m];
 203   }
 204   inc_nrnb(nrnb,eNR_PCOUPL,ncoupl);
 205 }
 206
 207 void tcoupl(bool bTC,t_grpopts *opts,t_groups *grps,
 208             real dt,real SAfactor,int step,int nmem)
 209 {
 210   static real *Told=NULL;
 211   int    i;
 212   real   T,reft,lll;
 213
 214   if (!Told) {
 215     snew(Told,opts->ngtc);
 216     for(i=0; (i<opts->ngtc); i++)
 217       Told[i]=opts->ref_t[i]*SAfactor;
 218   }
 219
 220   for(i=0; (i<opts->ngtc); i++) {
 221     reft=opts->ref_t[i]*SAfactor;
 222     if (reft < 0)
 223       reft=0;
 224
 225     Told[i] = run_aver(Told[i],grps->tcstat[i].T,step,nmem);
 226     T       = Told[i];
 227
 228     if ((bTC) && (T != 0.0)) {
 229       lll=sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
 230       grps->tcstat[i].lambda=max(min(lll,1.25),0.8);
 231     }
 232     else
 233       grps->tcstat[i].lambda=1.0;
 234 #ifdef DEBUGTC
 235     fprintf(stdlog,"group %d: T: %g, Lambda: %g\n",
 236             i,T,grps->tcstat[i].lambda);
 237 #endif
 238   }
 239 }
 240
 241