MDIonicStepper.C 8.33 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 18 19
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "MDIonicStepper.h"
20
#include "sampling.h"
21
#include <cstdlib> // drand48
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
    atoms_.set_velocities(v0_);
124 125 126 127
  }
  else if ( thermostat_ == "LOWE" )
  {
    const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
Francois Gygi committed
128
    const int nat = atoms_.size();
129 130 131 132 133 134 135 136 137
    double collision_probability = 0.0;
    if ( th_time_ == 0 )
    {
      collision_probability = 1.0;
    }
    else
    {
      if ( nat > 1 )
      {
Francois Gygi committed
138
        collision_probability = min(1.0,fabs(dt_)/(0.5*(nat-1)*th_time_));
139 140
      }
    }
141

142 143
    // scan all atom pairs in the space (is,ia)
    //int npairs = 0;
144
    if ( s_.ctxt_.onpe0() )
145
    {
146 147 148
      // compute collision only on task 0 and synchronize later
      // since calculation involves drand48
      for ( int is1 = 0; is1 < nsp_; is1++ )
149
      {
150
        for ( int is2 = is1; is2 < nsp_; is2++ )
151
        {
152 153 154 155 156
          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++ )
157
          {
158 159 160
            int ia2min = 0;
            if ( is1 == is2 ) ia2min = ia1 + 1;
            for ( int ia2 = ia2min; ia2 < na_[is2]; ia2++ )
161
            {
162 163 164 165 166 167 168 169 170 171 172 173 174 175
              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
176 177 178
                double xi1, xi2; // xi2 not used
                normal_dev(&xi1,&xi2);
                double lambda = xi1 * width;
179 180 181
                D3vector dv12 = mu * ( lambda - v12*e12 ) * e12;
                D3vector v1p = v1 + (1.0/m1) * dv12;
                D3vector v2p = v2 - (1.0/m2) * dv12;
Francois Gygi committed
182

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

187 188 189 190 191
                v0_[is2][3*ia2+0] = v2p.x;
                v0_[is2][3*ia2+1] = v2p.y;
                v0_[is2][3*ia2+2] = v2p.z;
              }
              //npairs++;
192 193 194 195 196
            }
          }
        }
      }
    }
197
    atoms_.set_velocities(v0_);
198 199
    //cout << " npairs: " << npairs << endl;
  }
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 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247
  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_);
  }
Francois Gygi committed
248
  constraints_.enforce_v(r0_,v0_);
249 250
  // recompute ekin as velocities may be affected by constraints
  compute_ekin();
Francois Gygi committed
251
  atoms_.set_velocities(v0_);
252 253 254
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
255
void MDIonicStepper::compute_ekin(void)
256
{
257
  ekin_ = 0.0;
Francois Gygi committed
258
  for ( int is = 0; is < v0_.size(); is++ )
259
  {
Francois Gygi committed
260
    for ( int i = 0; i < v0_[is].size(); i++ )
261
    {
Francois Gygi committed
262
      const double v = v0_[is][i];
263
      ekin_ += 0.5 * pmass_[is] * v * v;
264 265 266
    }
  }
}