////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or .
//
////////////////////////////////////////////////////////////////////////////////
//
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
#include "MDIonicStepper.h"
#include "sampling.h"
#include // drand48
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::compute_r(double e0, const vector >& f0)
{
// f0 contains forces at r0
// Compute new positions rp using the velocity Verlet algorithm
// enforce constraints for rp
// update rm <- r0, r0 <- rp, and update atomset
// compute rp
for ( int is = 0; is < r0_.size(); is++ )
{
const double dt2by2m = dt_ * dt_ / ( 2.0 * pmass_[is] );
for ( int i = 0; i < r0_[is].size(); i++ )
{
rp_[is][i] = r0_[is][i] + v0_[is][i] * dt_ + dt2by2m * f0[is][i];
}
}
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::compute_v(double e0, const vector >& f0)
{
// compute velocities v0_ using r0_, rm_ and f0(r0)
// enforce constraints for vp
// adjust velocities with the thermostat
// compute ekin
assert(dt_ != 0.0);
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];
}
}
// if using a thermostat, compute ekin for use in the thermostat
if ( thermostat_ != "OFF" )
{
constraints_.enforce_v(r0_,v0_);
compute_ekin();
}
if ( thermostat_ == "SCALING" )
{
eta_ = tanh ( ( temp() - th_temp_ ) / th_width_ ) / th_time_;
if ( s_.ctxt_.onpe0() )
{
cout << " thermostat: temp=" << temp() << endl;
cout << " thermostat: tref=" << th_temp_ << endl;
cout << " thermostat: eta=" << eta_ << endl;
}
const double fac = (1.0 - eta_ * fabs(dt_));
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
v0_[is][i] *= fac;
}
}
}
else if ( thermostat_ == "ANDERSEN" )
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
// if th_time_ is zero or less than dt, collision probability is one
const double collision_probability = th_time_ == 0 ? 1.0 :
min(fabs(dt_) / th_time_, 1.0);
if ( s_.ctxt_.onpe0() )
{
// compute collision on task 0 and synchronize later
for ( int is = 0; is < nsp_; is++ )
{
const double m = pmass_[is];
const double width = sqrt( boltz * th_temp_ / m );
for ( int ia = 0; ia < na_[is]; ia++ )
{
if ( drand48() < collision_probability )
{
// cout << " collision: atom is=" << is << " ia=" << ia << endl;
// draw pairs of unit variance gaussian random variables
double xi0, xi1, xi2, xi3; // xi3 not used
normal_dev(&xi0,&xi1);
normal_dev(&xi2,&xi3);
v0_[is][3*ia+0] = width * xi0;
v0_[is][3*ia+1] = width * xi1;
v0_[is][3*ia+2] = width * xi2;
}
}
}
}
atoms_.set_velocities(v0_);
}
else if ( thermostat_ == "LOWE" )
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
const int nat = atoms_.size();
double collision_probability = 0.0;
if ( th_time_ == 0 )
{
collision_probability = 1.0;
}
else
{
if ( nat > 1 )
{
collision_probability = min(1.0,fabs(dt_)/(0.5*(nat-1)*th_time_));
}
}
// scan all atom pairs in the space (is,ia)
//int npairs = 0;
if ( s_.ctxt_.onpe0() )
{
// compute collision only on task 0 and synchronize later
// since calculation involves drand48
for ( int is1 = 0; is1 < nsp_; is1++ )
{
for ( int is2 = is1; is2 < nsp_; is2++ )
{
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++ )
{
int ia2min = 0;
if ( is1 == is2 ) ia2min = ia1 + 1;
for ( int ia2 = ia2min; ia2 < na_[is2]; ia2++ )
{
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
double xi1, xi2; // xi2 not used
normal_dev(&xi1,&xi2);
double lambda = xi1 * width;
D3vector dv12 = mu * ( lambda - v12*e12 ) * e12;
D3vector v1p = v1 + (1.0/m1) * dv12;
D3vector v2p = v2 - (1.0/m2) * dv12;
v0_[is1][3*ia1+0] = v1p.x;
v0_[is1][3*ia1+1] = v1p.y;
v0_[is1][3*ia1+2] = v1p.z;
v0_[is2][3*ia2+0] = v2p.x;
v0_[is2][3*ia2+1] = v2p.y;
v0_[is2][3*ia2+2] = v2p.z;
}
//npairs++;
}
}
}
}
}
atoms_.set_velocities(v0_);
//cout << " npairs: " << npairs << endl;
}
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_);
}
constraints_.enforce_v(r0_,v0_);
// recompute ekin as velocities may be affected by constraints
compute_ekin();
atoms_.set_velocities(v0_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::compute_ekin(void)
{
ekin_ = 0.0;
for ( int is = 0; is < v0_.size(); is++ )
{
for ( int i = 0; i < v0_[is].size(); i++ )
{
const double v = v0_[is][i];
ekin_ += 0.5 * pmass_[is] * v * v;
}
}
}