MDIonicStepper.C 8.55 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
#include <cstdlib> // drand48
23
using namespace std;
24 25

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

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

Francois Gygi committed
43 44 45
  constraints_.enforce_r(r0_,rp_);
  rm_ = r0_;
  r0_ = rp_;
46
  atoms_.sync_positions(r0_);
Francois Gygi committed
47 48 49 50 51 52 53 54 55
  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
56
  // compute ekin
57

58
  assert(dt_ != 0.0);
Francois Gygi committed
59 60 61 62 63 64 65 66 67 68
  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];
    }
  }
69 70 71 72 73 74
  // if using a thermostat, compute ekin for use in the thermostat
  if ( thermostat_ != "OFF" )
  {
    constraints_.enforce_v(r0_,v0_);
    compute_ekin();
  }
75

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

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

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

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

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

190 191 192 193 194
                v0_[is2][3*ia2+0] = v2p.x;
                v0_[is2][3*ia2+1] = v2p.y;
                v0_[is2][3*ia2+2] = v2p.z;
              }
              //npairs++;
195 196 197 198 199
            }
          }
        }
      }
    }
200
    atoms_.sync_velocities(v0_);
201
    atoms_.set_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
  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);

250
    atoms_.sync_velocities(v0_);
251 252
    atoms_.set_velocities(v0_);
  }
Francois Gygi committed
253
  constraints_.enforce_v(r0_,v0_);
254 255
  // recompute ekin as velocities may be affected by constraints
  compute_ekin();
256
  atoms_.sync_velocities(v0_);
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
    }
  }
}