MDIonicStepper.C 8.51 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
18
// $Id: MDIonicStepper.C,v 1.23 2010-04-16 22:41:55 fgygi Exp $
19 20

#include "MDIonicStepper.h"
21
#include "sampling.h"
22
using namespace std;
23 24

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
25
void MDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
26
{
27
  // f0 contains forces at r0
Francois Gygi committed
28 29 30
  // Compute new positions rp using the velocity Verlet algorithm
  // enforce constraints for rp
  // update rm <- r0, r0 <- rp, and update atomset
31

Francois Gygi committed
32 33
  // compute rp
  for ( int is = 0; is < r0_.size(); is++ )
34
  {
Francois Gygi committed
35 36
    const double dt2by2m = dt_ * dt_ / ( 2.0 * pmass_[is] );
    for ( int i = 0; i < r0_[is].size(); i++ )
37
    {
Francois Gygi committed
38
      rp_[is][i] = r0_[is][i] + v0_[is][i] * dt_ + dt2by2m * f0[is][i];
39
    }
40
  }
41

Francois Gygi committed
42 43 44 45 46 47 48 49 50 51 52 53
  constraints_.enforce_r(r0_,rp_);
  rm_ = r0_;
  r0_ = rp_;
  atoms_.set_positions(r0_);
}

////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
{
  // compute velocities v0_ using r0_, rm_ and f0(r0)
  // enforce constraints for vp
  // adjust velocities with the thermostat
54
  // compute ekin
55

56
  assert(dt_ != 0.0);
Francois Gygi committed
57 58 59 60 61 62 63 64 65 66
  for ( int is = 0; is < v0_.size(); is++ )
  {
    assert(pmass_[is] > 0.0);
    const double dtby2m = dt_ / ( 2.0 * pmass_[is] );
    for ( int i = 0; i < v0_[is].size(); i++ )
    {
      const double vhalf = ( r0_[is][i] - rm_[is][i] ) / dt_;
      v0_[is][i] = vhalf + dtby2m * f0[is][i];
    }
  }
67 68 69 70 71 72
  // if using a thermostat, compute ekin for use in the thermostat
  if ( thermostat_ != "OFF" )
  {
    constraints_.enforce_v(r0_,v0_);
    compute_ekin();
  }
73

74
  if ( thermostat_ == "SCALING" )
75 76
  {
    eta_ = tanh ( ( temp() - th_temp_ ) / th_width_ ) / th_time_;
77 78
    if ( s_.ctxt_.onpe0() )
    {
79 80 81
      cout << "  thermostat: temp=" << temp() << endl;
      cout << "  thermostat: tref=" << th_temp_ << endl;
      cout << "  thermostat: eta=" << eta_ << endl;
82
    }
83

84
    const double fac = (1.0 - eta_ * fabs(dt_));
85 86 87 88
    for ( int is = 0; is < r0_.size(); is++ )
    {
      for ( int i = 0; i < r0_[is].size(); i++ )
      {
Francois Gygi committed
89
        v0_[is][i] *= fac;
90 91
      }
    }
92
  }
93
  else if ( thermostat_ == "ANDERSEN" )
94 95
  {
    const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
96
    // if th_time_ is zero or less than dt, collision probability is one
Francois Gygi committed
97
    const double collision_probability = th_time_ == 0 ? 1.0 :
98
                 min(fabs(dt_) / th_time_, 1.0);
Francois Gygi committed
99

100
    if ( s_.ctxt_.onpe0() )
101
    {
102 103
      // compute collision on task 0 and synchronize later
      for ( int is = 0; is < nsp_; is++ )
104
      {
105 106 107
        const double m = pmass_[is];
        const double width = sqrt( boltz * th_temp_ / m );
        for ( int ia = 0; ia < na_[is]; ia++ )
108
        {
109 110 111
          if ( drand48() < collision_probability )
          {
            // cout << " collision: atom is=" << is << " ia=" << ia << endl;
112 113 114 115
            // draw pairs of unit variance gaussian random variables
            double xi0, xi1, xi2, xi3; // xi3 not used
            normal_dev(&xi0,&xi1);
            normal_dev(&xi2,&xi3);
116 117 118 119
            v0_[is][3*ia+0] = width * xi0;
            v0_[is][3*ia+1] = width * xi1;
            v0_[is][3*ia+2] = width * xi2;
          }
120 121 122
        }
      }
    }
123 124 125
    atoms_.set_velocities(v0_);
    atoms_.sync();
    atoms_.get_velocities(v0_);
126 127 128 129
  }
  else if ( thermostat_ == "LOWE" )
  {
    const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
Francois Gygi committed
130
    const int nat = atoms_.size();
131 132 133 134 135 136 137 138 139
    double collision_probability = 0.0;
    if ( th_time_ == 0 )
    {
      collision_probability = 1.0;
    }
    else
    {
      if ( nat > 1 )
      {
Francois Gygi committed
140
        collision_probability = min(1.0,fabs(dt_)/(0.5*(nat-1)*th_time_));
141 142
      }
    }
143

144 145
    // scan all atom pairs in the space (is,ia)
    //int npairs = 0;
146
    if ( s_.ctxt_.onpe0() )
147
    {
148 149 150
      // compute collision only on task 0 and synchronize later
      // since calculation involves drand48
      for ( int is1 = 0; is1 < nsp_; is1++ )
151
      {
152
        for ( int is2 = is1; is2 < nsp_; is2++ )
153
        {
154 155 156 157 158
          const double m1 = pmass_[is1];
          const double m2 = pmass_[is2];
          const double mu = m1*m2/(m1+m2);
          const double width = sqrt(boltz * th_temp_ / mu);
          for ( int ia1 = 0; ia1 < na_[is1]; ia1++ )
159
          {
160 161 162
            int ia2min = 0;
            if ( is1 == is2 ) ia2min = ia1 + 1;
            for ( int ia2 = ia2min; ia2 < na_[is2]; ia2++ )
163
            {
164 165 166 167 168 169 170 171 172 173 174 175 176 177
              if ( drand48() < collision_probability )
              {
                // cout << " collision: pair " << is1 << " " << ia1 << " "
                //      << is2 << " " << ia2 << endl;
                D3vector r1(&r0_[is1][3*ia1]);
                s_.wf.cell().fold_in_ws(r1);
                D3vector r2(&r0_[is2][3*ia2]);
                s_.wf.cell().fold_in_ws(r2);
                D3vector r12 = r1 - r2;
                D3vector e12 = normalized(r12);
                D3vector v1(&v0_[is1][3*ia1]);
                D3vector v2(&v0_[is2][3*ia2]);
                D3vector v12 = v1 - v2;
                // draw a gaussian random variable
178 179 180
                double xi1, xi2; // xi2 not used
                normal_dev(&xi1,&xi2);
                double lambda = xi1 * width;
181 182 183
                D3vector dv12 = mu * ( lambda - v12*e12 ) * e12;
                D3vector v1p = v1 + (1.0/m1) * dv12;
                D3vector v2p = v2 - (1.0/m2) * dv12;
Francois Gygi committed
184

185 186 187
                v0_[is1][3*ia1+0] = v1p.x;
                v0_[is1][3*ia1+1] = v1p.y;
                v0_[is1][3*ia1+2] = v1p.z;
Francois Gygi committed
188

189 190 191 192 193
                v0_[is2][3*ia2+0] = v2p.x;
                v0_[is2][3*ia2+1] = v2p.y;
                v0_[is2][3*ia2+2] = v2p.z;
              }
              //npairs++;
194 195 196 197 198
            }
          }
        }
      }
    }
199 200 201
    atoms_.set_velocities(v0_);
    atoms_.sync();
    atoms_.get_velocities(v0_);
202 203
    //cout << " npairs: " << npairs << endl;
  }
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 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
  else if ( thermostat_ == "BDP" )
  {
    // stochastic velocity rescaling following Bussi,Donadio,Parrinello,
    // J.Chem.Phys. 126, 014101 (2007)
    // Choice of sign of alpha following
    // Bussi, Parrinello, Comput. Phys. Comm. 179, 26 (2008).
    if ( s_.ctxt_.onpe0() )
    {
      // rescale velocities on task 0 and synchronize later
      // since calculation involves drand48
      const double c = (th_time_ > 0.0) ? exp(-dt_/th_time_) : 0.0;
      // generate a pair of normal deviates
      double r1, r2;
      normal_dev(&r1, &r2);
      // compute the chi square deviate for the sum of ndofs_-1 squared normal
      // deviates
      double chi2;
      if ( (ndofs_-1)%2 == 0 )
        chi2 = 2.0 * gamma_dev((ndofs_-1)/2);
      else
        chi2 = 2.0 * gamma_dev((ndofs_-2)/2) + r2*r2;
      double alpha = 1.0, alpha2 = 1.0;
      if ( temp() != 0.0 )
      {
        double z = th_temp_ / ( ndofs_ * temp() );
        alpha2 = c + (1.0-c) * z * (chi2 + r1*r1)
          + 2.0 * r1 * sqrt( c * (1.0-c) * z );
        alpha = sqrt(alpha2);
        if ( r1 + sqrt(c / ((1.0-c)*z)) < 0.0 )
          alpha = -alpha;
      }
      for ( int is = 0; is < r0_.size(); is++ )
      {
        for ( int i = 0; i < r0_[is].size(); i++ )
        {
          v0_[is][i] *= alpha;
        }
      }
      // accumulate the change in kinetic energy in ekin_stepper_
      // ekin(t+dt) = alpha2 * ekin(t)
      ekin_stepper_ += ( 1.0 - alpha2 ) * ekin_;
      s_.ctxt_.dbcast_send(1,1,&ekin_stepper_,1);
    }
    if ( !s_.ctxt_.onpe0() )
      s_.ctxt_.dbcast_recv(1,1,&ekin_stepper_,1,0,0);

    atoms_.set_velocities(v0_);
    atoms_.sync();
    atoms_.get_velocities(v0_);
  }
Francois Gygi committed
254
  constraints_.enforce_v(r0_,v0_);
255 256
  // recompute ekin as velocities may be affected by constraints
  compute_ekin();
Francois Gygi committed
257
  atoms_.set_velocities(v0_);
258 259 260
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
261
void MDIonicStepper::compute_ekin(void)
262
{
263
  ekin_ = 0.0;
Francois Gygi committed
264
  for ( int is = 0; is < v0_.size(); is++ )
265
  {
Francois Gygi committed
266
    for ( int i = 0; i < v0_[is].size(); i++ )
267
    {
Francois Gygi committed
268
      const double v = v0_[is][i];
269
      ekin_ += 0.5 * pmass_[is] * v * v;
270 271 272
    }
  }
}