Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

file 232 lines (196 sloc) 7.115 kb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
/*
Copyright (C) 2010,2011,2012 The ESPResSo project
Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
Max-Planck-Institute for Polymer Research, Theory Group
This file is part of ESPResSo.
ESPResSo is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ESPResSo is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _THERMOSTAT_H
#define _THERMOSTAT_H
/** \file thermostat.h

*/

#include <math.h>
#include "utils.h"
#include "particle_data.h"
#include "random.h"
#include "global.h"
#include "integrate.h"
#include "cells.h"
#include "lb.h"
#include "dpd.h"
#include "virtual_sites.h"

/** \name Thermostat switches*/
/************************************************************/
/*@{*/

#define THERMO_OFF 0
#define THERMO_LANGEVIN 1
#define THERMO_DPD 2
#define THERMO_NPT_ISO 4
#define THERMO_LB 8
#define THERMO_INTER_DPD 16

/*@}*/

/************************************************
* exported variables
************************************************/

/** Switch determining which thermostat to use. This is a or'd value
of the different possible thermostats (defines: \ref THERMO_OFF,
\ref THERMO_LANGEVIN, \ref THERMO_DPD \ref THERMO_NPT_ISO). If it
is zero all thermostats are switched off and the temperature is
set to zero. */
extern int thermo_switch;

/** temperature. */
extern double temperature;

/** Langevin friction coefficient gamma. */
extern double langevin_gamma;

/** Friction coefficient for nptiso-thermostat's inline-function friction_therm0_nptiso */
extern double nptiso_gamma0;
/** Friction coefficient for nptiso-thermostat's inline-function friction_thermV_nptiso */
extern double nptiso_gammav;

/************************************************
* functions
************************************************/


/** initialize constants of the thermostat on
start of integration */
void thermo_init();

/** very nasty: if we recalculate force when leaving/reentering the integrator,
a(t) and a((t-dt)+dt) are NOT equal in the vv algorithm. The random
numbers are drawn twice, resulting in a different variance of the random force.
This is corrected by additional heat when restarting the integrator here.
Currently only works for the Langevin thermostat, although probably also others
are affected.
*/
void thermo_heat_up();

/** pendant to \ref thermo_heat_up */
void thermo_cool_down();

#ifdef NPT
/** add velocity-dependend noise and friction for NpT-sims to the particle's velocity
@param dt_vj j-component of the velocity scaled by time_step dt
@return j-component of the noise added to the velocity, also scaled by dt (contained in prefactors) */
MDINLINE double friction_therm0_nptiso(double dt_vj) {
  extern double nptiso_pref1, nptiso_pref2;
  if(thermo_switch & THERMO_NPT_ISO)
    return ( nptiso_pref1*dt_vj + nptiso_pref2*(d_random()-0.5) );
  return 0.0;
}

/** add p_diff-dependend noise and friction for NpT-sims to \ref nptiso_struct::p_diff */
MDINLINE double friction_thermV_nptiso(double p_diff) {
  extern double nptiso_pref3, nptiso_pref4;
  if(thermo_switch & THERMO_NPT_ISO)
    return ( nptiso_pref3*p_diff + nptiso_pref4*(d_random()-0.5) );
  return 0.0;
}
#endif

/** overwrite the forces of a particle with
the friction term, i.e. \f$ F_i= -\gamma v_i + \xi_i\f$.
*/
MDINLINE void friction_thermo_langevin(Particle *p)
{
  extern double langevin_pref1, langevin_pref2;
#ifdef LANGEVIN_PER_PARTICLE
  double langevin_pref1_temp, langevin_pref2_temp;
#endif
  
  int j;
#ifdef MASS
  double massf = sqrt(PMASS(*p));
#else
  double massf = 1;
#endif


#ifdef VIRTUAL_SITES
#ifndef VIRTUAL_SITES_THERMOSTAT
    if (ifParticleIsVirtual(p))
    {
     for (j=0;j<3;j++)
      p->f.f[j]=0;
    return;
   }
#endif

#ifdef THERMOSTAT_IGNORE_NON_VIRTUAL
    if (!ifParticleIsVirtual(p))
    {
     for (j=0;j<3;j++)
      p->f.f[j]=0;
    return;
   }
#endif
#endif

  for ( j = 0 ; j < 3 ; j++) {
#ifdef EXTERNAL_FORCES
    if (!(p->l.ext_flag & COORD_FIXED(j)))
#endif
    {
#ifdef LANGEVIN_PER_PARTICLE
      if(p->p.gamma >= 0.) {
        langevin_pref1_temp = -p->p.gamma/time_step;
        
        if(p->p.T >= 0.)
          langevin_pref2_temp = sqrt(24.0*p->p.T*p->p.gamma/time_step);
        else
          langevin_pref2_temp = sqrt(24.0*temperature*p->p.gamma/time_step);
        
        p->f.f[j] = langevin_pref1_temp*p->m.v[j]*PMASS(*p) + langevin_pref2_temp*(d_random()-0.5)*massf;
      }
      else {
        if(p->p.T >= 0.)
          langevin_pref2_temp = sqrt(24.0*p->p.T*langevin_gamma/time_step);
        else
          langevin_pref2_temp = langevin_pref2;
        
        p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) + langevin_pref2_temp*(d_random()-0.5)*massf;
      }
#else
      p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) + langevin_pref2*(d_random()-0.5)*massf;
#endif
    }
#ifdef EXTERNAL_FORCES
    else p->f.f[j] = 0;
#endif
  }
// printf("%d: %e %e %e %e %e %e\n",p->p.identity, p->f.f[0],p->f.f[1],p->f.f[2], p->m.v[0],p->m.v[1],p->m.v[2]);
  

  ONEPART_TRACE(if(p->p.identity==check_id) fprintf(stderr,"%d: OPT: LANG f = (%.3e,%.3e,%.3e)\n",this_node,p->f.f[0],p->f.f[1],p->f.f[2]));
  THERMO_TRACE(fprintf(stderr,"%d: Thermo: P %d: force=(%.3e,%.3e,%.3e)\n",this_node,p->p.identity,p->f.f[0],p->f.f[1],p->f.f[2]));
}

#ifdef ROTATION
/** set the particle torques to the friction term, i.e. \f$\tau_i=-\gamma w_i + \xi_i\f$.
The same friction coefficient \f$\gamma\f$ is used as that for translation.
*/
MDINLINE void friction_thermo_langevin_rotation(Particle *p)
{
  extern double langevin_pref2;

  int j;
#ifdef VIRTUAL_SITES
#ifndef VIRTUAL_SITES_THERMOSTAT
    if (ifParticleIsVirtual(p))
    {
     for (j=0;j<3;j++)
      p->f.torque[j]=0;
    return;
   }
#endif

#ifdef THERMOSTAT_IGNORE_NON_VIRTUAL
    if (!ifParticleIsVirtual(p))
    {
     for (j=0;j<3;j++)
      p->f.torque[j]=0;
    return;
   }
#endif
#endif
      for ( j = 0 ; j < 3 ; j++)
      {
#ifdef ROTATIONAL_INERTIA
p->f.torque[j] = -langevin_gamma*p->m.omega[j] *p->p.rinertia[j] + langevin_pref2*sqrt(p->p.rinertia[j]) * (d_random()-0.5);
#else
p->f.torque[j] = -langevin_gamma*p->m.omega[j] + langevin_pref2*(d_random()-0.5);
#endif
      }
      ONEPART_TRACE(if(p->p.identity==check_id) fprintf(stderr,"%d: OPT: LANG f = (%.3e,%.3e,%.3e)\n",this_node,p->f.f[0],p->f.f[1],p->f.f[2]));
      THERMO_TRACE(fprintf(stderr,"%d: Thermo: P %d: force=(%.3e,%.3e,%.3e)\n",this_node,p->p.identity,p->f.f[0],p->f.f[1],p->f.f[2]));
}
#endif

#endif
Something went wrong with that request. Please try again.