Commit 95148d96 authored by Francois Gygi's avatar Francois Gygi
Browse files

Add lock_cm variable

parent c530e147
......@@ -83,6 +83,9 @@ void ANDIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
i++;
}
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
......
......@@ -31,6 +31,10 @@ void BMDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
for ( int is = 0; is < r0_.size(); is++ )
for ( int i = 0; i < r0_[is].size(); i++ )
rp_[is][i] = r0_[is][i] + v0_[is][i] + bmd_fac_ * f0[is][i];
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
fm_ = f0;
......
......@@ -510,6 +510,7 @@ void BOSampleStepper::step(int niter)
}
}
cout << "</atomset>" << endl;
cout << "<rcm> " << atoms.rcm() << "</rcm>" << endl;
cout << setprecision(6);
cout << "<unit_cell_a_norm> " << atoms.cell().a_norm(0)
<< " </unit_cell_a_norm>" << endl;
......
......@@ -88,6 +88,9 @@ void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
for ( int j = 0; j < r0_[is].size(); j++ )
rp_[is][j] = xp[i++];
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
......
......@@ -196,6 +196,7 @@ void CPSampleStepper::step(int niter)
}
}
cout << "</atomset>" << endl;
cout << "<rcm> " << atoms.rcm() << "</rcm>" << endl;
}
if ( s_.constraints.size() > 0 )
......
......@@ -47,6 +47,8 @@ struct Control
std::string thermostat;
double th_temp,th_time, th_width; // thermostat control
bool lock_cm; // lock center of mass
std::string stress;
std::string cell_dyn;
std::string cell_lock;
......
......@@ -81,6 +81,41 @@ class IonicStepper
void get_velocities(void) { atoms_.get_velocities(v0_); }
void set_positions(void) { atoms_.set_positions(r0_); }
void set_velocities(void) { atoms_.set_velocities(v0_); }
// center of mass position
D3vector rcm(const std::vector<std::vector<double> >& r)
{
D3vector mrsum;
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++ )
{
mrsum += mass * D3vector(r[is][3*ia],r[is][3*ia+1],r[is][3*ia+2]);
msum += mass;
}
}
if ( msum == 0.0 ) return D3vector(0,0,0);
return mrsum / msum;
}
// reset center of mass position of r2 to be equal to that of r1
void reset_rcm(const std::vector<std::vector<double> >& r1,
std::vector<std::vector<double> >& r2)
{
D3vector drcm = rcm(r2) - rcm(r1);
for ( int is = 0; is < nsp_; is++ )
{
int nais = na_[is];
for ( int ia = 0; ia < nais; ia++ )
{
r2[is][3*ia] -= drcm.x;
r2[is][3*ia+1] -= drcm.y;
r2[is][3*ia+2] -= drcm.z;
}
}
}
void setup_constraints(void)
{
......
////////////////////////////////////////////////////////////////////////////////
//
// 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 <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// LockCm.h
//
////////////////////////////////////////////////////////////////////////////////
#ifndef LOCKCM_H
#define LOCKCM_H
#include<iostream>
#include<iomanip>
#include<sstream>
#include<stdlib.h>
#include "Sample.h"
class LockCm : public Var
{
Sample *s;
public:
const char *name ( void ) const { return "lock_cm"; };
int set ( int argc, char **argv )
{
if ( argc != 2 )
{
if ( ui->onpe0() )
cout << " lock_cm takes only one value" << endl;
return 1;
}
string v = argv[1];
if ( !( v == "ON" || v == "OFF" ) )
{
if ( ui->onpe0() )
cout << " lock_cm must be ON or OFF" << endl;
return 1;
}
s->ctrl.lock_cm = ( v == "ON" );
return 0;
}
string print (void) const
{
ostringstream st;
st.setf(ios::left,ios::adjustfield);
st << setw(10) << name() << " = ";
st.setf(ios::right,ios::adjustfield);
if ( s->ctrl.lock_cm )
st << setw(10) << "ON";
else
st << setw(10) << "OFF";
return st.str();
}
LockCm(Sample *sample) : s(sample) { s->ctrl.lock_cm = false; }
};
#endif
......@@ -40,6 +40,9 @@ void MDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
}
}
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
......
......@@ -51,6 +51,10 @@ void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
rp_[is][i] = r0_[is][i] + alpha_sd * f0[is][i];
}
}
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
......
......@@ -31,6 +31,10 @@ void SDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
rp_[is][i] = r0_[is][i] + dt2bym * f0[is][i];
}
}
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
......
......@@ -105,6 +105,7 @@ using namespace std;
#include "FermiTemp.h"
#include "IterCmd.h"
#include "IterCmdPeriod.h"
#include "LockCm.h"
#include "Dt.h"
#include "MuRSH.h"
#include "Nempty.h"
......@@ -456,6 +457,7 @@ int main(int argc, char **argv, char **envp)
ui.addVar(new FermiTemp(s));
ui.addVar(new IterCmd(s));
ui.addVar(new IterCmdPeriod(s));
ui.addVar(new LockCm(s));
ui.addVar(new MuRSH(s));
ui.addVar(new Nempty(s));
ui.addVar(new NetCharge(s));
......
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