src/gmxlib/coupling.c

   1 #include <stdio.h>
   2 #include "typedefs.h"
   3 #include "smalloc.h"
   4 #include "update.h"
   5 #include "vec.h"
   6 #include "macros.h"
   7 #include "physics.h"
   8 #include "names.h"
   9 #include "fatal.h"
  10 #include "txtdump.h"
  11 #include "nrnb.h"
  12
  13 /*
  14  * This file implements temperature and pressure coupling algorithms:
  15  * For now only the Weak coupling and the modified weak coupling.
  16  *
  17  * Furthermore computation of pressure and temperature is done here
  18  *
  19  */
  20
  21 void calc_pres(matrix box,tensor ekin,tensor vir,tensor pres,real Elr)
  22 {
  23   int  n,m;
  24   real fac,Plr;
  25
  26   /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E. */
  27   /* Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in het */
  28   /* systeem...       */
  29
  30   /* Long range correctie for periodic systems, see
  31    * Neumann et al. JCP
  32    * divide by 6 because it is multiplied by fac later on.
  33    * If Elr = 0, no correction is made.
  34    */
  35   Plr = Elr/6.0;
  36
  37   fac=PRESFAC*2.0/det(box);
  38   for(n=0; (n<DIM); n++)
  39     for(m=0; (m<DIM); m++)
  40       pres[n][m]=(ekin[n][m]-vir[n][m]+Plr)*fac;
  41
  42   if (debug) {
  43     pr_rvecs(debug,0,"PC: pres",pres,DIM);
  44     pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
  45     pr_rvecs(debug,0,"PC: vir ",vir, DIM);
  46     pr_rvecs(debug,0,"PC: box ",box, DIM);
  47   }
  48 }
  49
  50 real calc_temp(real ekin,int nrdf)
  51 {
  52   return (2.0*ekin)/(nrdf*BOLTZ);
  53 }
  54
  55 real run_aver(real old,real cur,int step,int nmem)
  56 {
  57   nmem   = max(1,nmem);
  58
  59   return ((nmem-1)*old+cur)/nmem;
  60 }
  61
  62
  63 void do_pcoupl(t_inputrec *ir,int step,tensor pres,
  64                matrix box,int start,int nr_atoms,
  65                rvec x[],ushort cFREEZE[],
  66                t_nrnb *nrnb,rvec freezefac[])
  67 {
  68   static bool bFirst=TRUE;
  69   static rvec PPP;
  70   int    n,d,m,g,ncoupl=0;
  71   real   scalar_pressure, xy_pressure, p_corr_z;
  72   real   X,Y,Z,dx,dy,dz;
  73   rvec   factor;
  74   tensor mu;
  75   real   muxx,muxy,muxz,muyx,muyy,muyz,muzx,muzy,muzz;
  76   real   fgx,fgy,fgz;
  77
  78   /*
  79    *  PRESSURE SCALING
  80    *  Step (2P)
  81    */
  82   if (bFirst) {
  83     /* Initiate the pressure to the reference one */
  84     for(m=0; m<DIM; m++)
  85       PPP[m] = ir->ref_p[m];
  86     bFirst=FALSE;
  87   }
  88   scalar_pressure=0;
  89   xy_pressure=0;
  90   for(m=0; m<DIM; m++) {
  91     PPP[m]           = run_aver(PPP[m],pres[m][m],step,ir->npcmemory);
  92     scalar_pressure += PPP[m]/DIM;
  93     if (m != ZZ)
  94       xy_pressure += PPP[m]/(DIM-1);
  95   }
  96
  97   /* Pressure is now in bar, everywhere. */
  98   if ((ir->epc != epcNO) && (scalar_pressure != 0.0)) {
  99     for(m=0; m<DIM; m++)
 100       factor[m] = ir->compress[m]*ir->delta_t/ir->tau_p;
 101     clear_mat(mu);
 102     switch (ir->epc) {
 103     case epcISOTROPIC:
 104       for(m=0; m<DIM; m++)
 105         mu[m][m] = pow(1.0-factor[m]*(ir->ref_p[m]-scalar_pressure),1.0/DIM);
 106       break;
 107     case epcSEMIISOTROPIC:
 108       for(m=0; m<ZZ; m++)
 109         mu[m][m] = pow(1.0-factor[m]*(ir->ref_p[m]-xy_pressure),1.0/DIM);
 110       mu[ZZ][ZZ] = pow(1.0-factor[ZZ]*(ir->ref_p[ZZ] - PPP[ZZ]),1.0/DIM);
 111       break;
 112     case epcANISOTROPIC:
 113       for (m=0; m<DIM; m++)
 114         mu[m][m] = pow(1.0-factor[m]*(ir->ref_p[m] - PPP[m]),1.0/DIM);
 115       break;
 116     case epcSURFACETENSION:
 117       /* ir->ref_p[0/1] is the reference surface-tension times *
 118        * the number of surfaces                                */
 119       if (ir->compress[ZZ])
 120         p_corr_z = ir->delta_t/ir->tau_p*(ir->ref_p[ZZ] - PPP[ZZ]);
 121       else
 122         /* when the compressibity is zero, set the pressure correction   *
 123          * in the z-direction to zero to get the correct surface tension */
 124         p_corr_z = 0;
 125       mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ]*p_corr_z;
 126       for(m=0; m<ZZ; m++)
 127         mu[m][m] = sqrt(1.0+factor[m]*(ir->ref_p[m]/(mu[ZZ][ZZ]*box[ZZ][ZZ]) -
 128         (PPP[ZZ]+p_corr_z - xy_pressure)));
 129       break;
 130     case epcTRICLINIC:
 131     default:
 132       fatal_error(0,"Pressure coupling type %s not supported yet\n",
 133                   EPCOUPLTYPE(ir->epc));
 134     }
 135     if (debug) {
 136       pr_rvecs(debug,0,"PC: PPP ",(rvec *)&PPP,1);
 137       pr_rvecs(debug,0,"PC: fac ",(rvec *)&factor,1);
 138       pr_rvecs(debug,0,"PC: mu  ",mu,DIM);
 139     }
 140     /* Scale the positions using matrix operation */
 141     nr_atoms+=start;
 142
 143     muxx=mu[XX][XX],muxy=mu[XX][YY],muxz=mu[XX][ZZ];
 144     muyx=mu[YY][XX],muyy=mu[YY][YY],muyz=mu[YY][ZZ];
 145     muzx=mu[ZZ][XX],muzy=mu[ZZ][YY],muzz=mu[ZZ][ZZ];
 146     for (n=start; n<nr_atoms; n++) {
 147       g=cFREEZE[n];
 148       fgx=freezefac[g][XX];
 149       fgy=freezefac[g][YY];
 150       fgz=freezefac[g][ZZ];
 151
 152       X=x[n][XX];
 153       Y=x[n][YY];
 154       Z=x[n][ZZ];
 155       dx=muxx*X+muxy*Y+muxz*Z;
 156       dy=muyx*X+muyy*Y+muyz*Z;
 157       dz=muzx*X+muzy*Y+muzz*Z;
 158       x[n][XX]=X+fgx*(dx-X);
 159       x[n][YY]=Y+fgy*(dy-Y);
 160       x[n][ZZ]=Z+fgz*(dz-Z);
 161
 162       ncoupl++;
 163     }
 164     /* compute final boxlengths */
 165     for (d=0; d<DIM; d++)
 166       for (m=0; m<DIM; m++)
 167         box[d][m] *= mu[d][m];
 168   }
 169   inc_nrnb(nrnb,eNR_PCOUPL,ncoupl);
 170 }
 171
 172 void tcoupl(bool bTC,t_grpopts *opts,t_groups *grps,
 173             real dt,real SAfactor,int step,int nmem)
 174 {
 175   static real *Told=NULL;
 176   int    i;
 177   real   T,reft,lll;
 178
 179   if (!Told) {
 180     snew(Told,opts->ngtc);
 181     for(i=0; (i<opts->ngtc); i++)
 182       Told[i]=opts->ref_t[i]*SAfactor;
 183   }
 184
 185   for(i=0; (i<opts->ngtc); i++) {
 186     reft=opts->ref_t[i]*SAfactor;
 187     if (reft < 0)
 188       reft=0;
 189
 190     Told[i] = run_aver(Told[i],grps->tcstat[i].T,step,nmem);
 191     T       = Told[i];
 192
 193     if ((bTC) && (T != 0.0)) {
 194       lll=sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
 195       grps->tcstat[i].lambda=max(min(lll,1.25),0.8);
 196     }
 197     else
 198       grps->tcstat[i].lambda=1.0;
 199 #ifdef DEBUGTC
 200     fprintf(stdlog,"group %d: T: %g, Lambda: %g\n",
 201             i,T,grps->tcstat[i].lambda);
 202 #endif
 203   }
 204 }
 205
 206