Commit f8949045 authored by Francois Gygi's avatar Francois Gygi
Browse files

Add reset vcm if lock_cm is ON

parent 95148d96
......@@ -74,6 +74,10 @@ void BMDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
for ( int i = 0; i < r0_[is].size(); i++ )
v0_[is][i] *= 1.05;
}
if ( s_.ctrl.lock_cm )
reset_vcm(v0_);
constraints_.enforce_v(r0_,v0_);
compute_ekin();
}
......
......@@ -117,6 +117,41 @@ class IonicStepper
}
}
// center of mass velocity
D3vector vcm(const std::vector<std::vector<double> >& v)
{
D3vector mvsum;
double msum = 0.0;
for ( int is = 0; is < nsp_; is++ )
{
int nais = na_[is];
double mass = pmass_[is];
for ( int ia = 0; ia < nais; ia++ )
{
mvsum += mass * D3vector(v[is][3*ia],v[is][3*ia+1],v[is][3*ia+2]);
msum += mass;
}
}
if ( msum == 0.0 ) return D3vector(0,0,0);
return mvsum / msum;
}
// reset center of mass velocity to zero
void reset_vcm(std::vector<std::vector<double> >& v)
{
D3vector dvcm = vcm(v);
for ( int is = 0; is < nsp_; is++ )
{
int nais = na_[is];
for ( int ia = 0; ia < nais; ia++ )
{
v[is][3*ia] -= dvcm.x;
v[is][3*ia+1] -= dvcm.y;
v[is][3*ia+2] -= dvcm.z;
}
}
}
void setup_constraints(void)
{
constraints_.setup(atoms_);
......
......@@ -248,6 +248,10 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
atoms_.set_velocities(v0_);
}
if ( s_.ctrl.lock_cm )
reset_vcm(v0_);
constraints_.enforce_v(r0_,v0_);
// recompute ekin as velocities may be affected by constraints
compute_ekin();
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment