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

Merge branch 'develop'

parents 36bc0b41 3035809e
...@@ -83,6 +83,9 @@ void ANDIonicStepper::compute_r(double e0, const vector<vector<double> >& f0) ...@@ -83,6 +83,9 @@ void ANDIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
i++; i++;
} }
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
r0_ = rp_; r0_ = rp_;
......
...@@ -31,6 +31,10 @@ void BMDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) ...@@ -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 is = 0; is < r0_.size(); is++ )
for ( int i = 0; i < r0_[is].size(); i++ ) for ( int i = 0; i < r0_[is].size(); i++ )
rp_[is][i] = r0_[is][i] + v0_[is][i] + bmd_fac_ * f0[is][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_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
fm_ = f0; fm_ = f0;
...@@ -70,6 +74,10 @@ void BMDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0) ...@@ -70,6 +74,10 @@ void BMDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
for ( int i = 0; i < r0_[is].size(); i++ ) for ( int i = 0; i < r0_[is].size(); i++ )
v0_[is][i] *= 1.05; v0_[is][i] *= 1.05;
} }
if ( s_.ctrl.lock_cm )
reset_vcm(v0_);
constraints_.enforce_v(r0_,v0_); constraints_.enforce_v(r0_,v0_);
compute_ekin(); compute_ekin();
} }
......
...@@ -510,6 +510,7 @@ void BOSampleStepper::step(int niter) ...@@ -510,6 +510,7 @@ void BOSampleStepper::step(int niter)
} }
} }
cout << "</atomset>" << endl; cout << "</atomset>" << endl;
cout << "<rcm> " << atoms.rcm() << "</rcm>" << endl;
cout << setprecision(6); cout << setprecision(6);
cout << "<unit_cell_a_norm> " << atoms.cell().a_norm(0) cout << "<unit_cell_a_norm> " << atoms.cell().a_norm(0)
<< " </unit_cell_a_norm>" << endl; << " </unit_cell_a_norm>" << endl;
......
...@@ -88,6 +88,9 @@ void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0) ...@@ -88,6 +88,9 @@ void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
for ( int j = 0; j < r0_[is].size(); j++ ) for ( int j = 0; j < r0_[is].size(); j++ )
rp_[is][j] = xp[i++]; rp_[is][j] = xp[i++];
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
r0_ = rp_; r0_ = rp_;
......
...@@ -196,6 +196,7 @@ void CPSampleStepper::step(int niter) ...@@ -196,6 +196,7 @@ void CPSampleStepper::step(int niter)
} }
} }
cout << "</atomset>" << endl; cout << "</atomset>" << endl;
cout << "<rcm> " << atoms.rcm() << "</rcm>" << endl;
} }
if ( s_.constraints.size() > 0 ) if ( s_.constraints.size() > 0 )
......
...@@ -47,6 +47,8 @@ struct Control ...@@ -47,6 +47,8 @@ struct Control
std::string thermostat; std::string thermostat;
double th_temp,th_time, th_width; // thermostat control double th_temp,th_time, th_width; // thermostat control
bool lock_cm; // lock center of mass
std::string stress; std::string stress;
std::string cell_dyn; std::string cell_dyn;
std::string cell_lock; std::string cell_lock;
......
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// IonicStepper.cpp:
//
////////////////////////////////////////////////////////////////////////////////
#include "IonicStepper.h"
IonicStepper::IonicStepper (Sample& s) : s_(s), atoms_(s.atoms),
constraints_(s.constraints), dt_(s.ctrl.dt)
{
ndofs_ = 3 * atoms_.size() - constraints_.ndofs();
// if there are more constraints than dofs, set ndofs_ to zero
if ( ndofs_ < 0 ) ndofs_ = 0;
nsp_ = atoms_.nsp();
na_.resize(nsp_);
r0_.resize(nsp_);
rp_.resize(nsp_);
v0_.resize(nsp_);
pmass_.resize(nsp_);
for ( int is = 0; is < nsp_; is++ )
{
const int nais = atoms_.na(is);
na_[is] = nais;
r0_[is].resize(3*nais);
rp_[is].resize(3*nais);
v0_[is].resize(3*nais);
pmass_[is] = atoms_.species_list[is]->mass() * 1822.89;
}
natoms_ = atoms_.size();
// get positions and velocities from atoms_
get_positions();
get_velocities();
}
// center of mass position
D3vector IonicStepper::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 IonicStepper::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;
}
}
}
// center of mass velocity
D3vector IonicStepper::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 IonicStepper::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;
}
}
}
...@@ -45,32 +45,7 @@ class IonicStepper ...@@ -45,32 +45,7 @@ class IonicStepper
public: public:
IonicStepper (Sample& s) : s_(s), atoms_(s.atoms), IonicStepper(Sample& s);
constraints_(s.constraints), dt_(s.ctrl.dt)
{
ndofs_ = 3 * atoms_.size() - constraints_.ndofs();
// if there are more constraints than dofs, set ndofs_ to zero
if ( ndofs_ < 0 ) ndofs_ = 0;
nsp_ = atoms_.nsp();
na_.resize(nsp_);
r0_.resize(nsp_);
rp_.resize(nsp_);
v0_.resize(nsp_);
pmass_.resize(nsp_);
for ( int is = 0; is < nsp_; is++ )
{
const int nais = atoms_.na(is);
na_[is] = nais;
r0_[is].resize(3*nais);
rp_[is].resize(3*nais);
v0_[is].resize(3*nais);
pmass_[is] = atoms_.species_list[is]->mass() * 1822.89;
}
natoms_ = atoms_.size();
// get positions and velocities from atoms_
get_positions();
get_velocities();
}
double r0(int is, int i) const { return r0_[is][i]; } double r0(int is, int i) const { return r0_[is][i]; }
double v0(int is, int i) const { return v0_[is][i]; } double v0(int is, int i) const { return v0_[is][i]; }
...@@ -81,11 +56,16 @@ class IonicStepper ...@@ -81,11 +56,16 @@ class IonicStepper
void get_velocities(void) { atoms_.get_velocities(v0_); } void get_velocities(void) { atoms_.get_velocities(v0_); }
void set_positions(void) { atoms_.set_positions(r0_); } void set_positions(void) { atoms_.set_positions(r0_); }
void set_velocities(void) { atoms_.set_velocities(v0_); } void set_velocities(void) { atoms_.set_velocities(v0_); }
// center of mass position
void setup_constraints(void) D3vector rcm(const std::vector<std::vector<double> >& r);
{ // reset center of mass position of r2 to be equal to that of r1
constraints_.setup(atoms_); void reset_rcm(const std::vector<std::vector<double> >& r1,
} std::vector<std::vector<double> >& r2);
// center of mass velocity
D3vector vcm(const std::vector<std::vector<double> >& v);
// reset center of mass velocity to zero
void reset_vcm(std::vector<std::vector<double> >& v);
void setup_constraints(void) { constraints_.setup(atoms_); }
virtual void compute_r(double e0, virtual void compute_r(double e0,
const std::vector<std::vector< double> >& f0) = 0; const std::vector<std::vector< double> >& f0) = 0;
virtual void compute_v(double e0, virtual void compute_v(double e0,
...@@ -96,6 +76,5 @@ class IonicStepper ...@@ -96,6 +76,5 @@ class IonicStepper
virtual double temp(void) const { return 0.0; } virtual double temp(void) const { return 0.0; }
virtual ~IonicStepper() {} virtual ~IonicStepper() {}
}; };
#endif #endif
////////////////////////////////////////////////////////////////////////////////
//
// 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) ...@@ -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_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
r0_ = rp_; r0_ = rp_;
...@@ -245,6 +248,10 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0) ...@@ -245,6 +248,10 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
atoms_.set_velocities(v0_); atoms_.set_velocities(v0_);
} }
if ( s_.ctrl.lock_cm )
reset_vcm(v0_);
constraints_.enforce_v(r0_,v0_); constraints_.enforce_v(r0_,v0_);
// recompute ekin as velocities may be affected by constraints // recompute ekin as velocities may be affected by constraints
compute_ekin(); compute_ekin();
......
...@@ -50,8 +50,8 @@ OBJECTS=MPIdata.o qb.o AtomSet.o Atom.o Species.o \ ...@@ -50,8 +50,8 @@ OBJECTS=MPIdata.o qb.o AtomSet.o Atom.o Species.o \
SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \ SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
XMLGFPreprocessor.o Base64Transcoder.o \ XMLGFPreprocessor.o Base64Transcoder.o \
CPSampleStepper.o BOSampleStepper.o \ CPSampleStepper.o BOSampleStepper.o \
MDWavefunctionStepper.o SDIonicStepper.o MDIonicStepper.o \ MDWavefunctionStepper.o IonicStepper.o \
BMDIonicStepper.o \ SDIonicStepper.o MDIonicStepper.o BMDIonicStepper.o \
CellStepper.o SDCellStepper.o \ CellStepper.o SDCellStepper.o \
CGCellStepper.o ConfinementPotential.o Preconditioner.o \ CGCellStepper.o ConfinementPotential.o Preconditioner.o \
AndersonMixer.o SDAIonicStepper.o CGIonicStepper.o \ AndersonMixer.o SDAIonicStepper.o CGIonicStepper.o \
......
...@@ -51,6 +51,10 @@ void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) ...@@ -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]; 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_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
r0_ = rp_; r0_ = rp_;
......
...@@ -31,6 +31,10 @@ void SDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) ...@@ -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]; rp_[is][i] = r0_[is][i] + dt2bym * f0[is][i];
} }
} }
if ( s_.ctrl.lock_cm )
reset_rcm(r0_,rp_);
constraints_.enforce_r(r0_,rp_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
r0_ = rp_; r0_ = rp_;
......
...@@ -70,7 +70,10 @@ class StatusCmd : public Cmd ...@@ -70,7 +70,10 @@ class StatusCmd : public Cmd
if ( s->wfv != 0 ) if ( s->wfv != 0 )
s->wfv->info(cout,"wfv"); s->wfv->info(cout,"wfv");
if ( ui->onpe0() ) if ( ui->onpe0() )
{
cout << "<rcm> " << s->atoms.rcm() << " </rcm>" << endl;
cout << "<vcm> " << s->atoms.vcm() << " </vcm>" << endl; cout << "<vcm> " << s->atoms.vcm() << " </vcm>" << endl;
}
return 0; return 0;
} }
}; };
......
...@@ -105,6 +105,7 @@ using namespace std; ...@@ -105,6 +105,7 @@ using namespace std;
#include "FermiTemp.h" #include "FermiTemp.h"
#include "IterCmd.h" #include "IterCmd.h"
#include "IterCmdPeriod.h" #include "IterCmdPeriod.h"
#include "LockCm.h"
#include "Dt.h" #include "Dt.h"
#include "MuRSH.h" #include "MuRSH.h"
#include "Nempty.h" #include "Nempty.h"
...@@ -456,6 +457,7 @@ int main(int argc, char **argv, char **envp) ...@@ -456,6 +457,7 @@ int main(int argc, char **argv, char **envp)
ui.addVar(new FermiTemp(s)); ui.addVar(new FermiTemp(s));
ui.addVar(new IterCmd(s)); ui.addVar(new IterCmd(s));
ui.addVar(new IterCmdPeriod(s)); ui.addVar(new IterCmdPeriod(s));
ui.addVar(new LockCm(s));
ui.addVar(new MuRSH(s)); ui.addVar(new MuRSH(s));
ui.addVar(new Nempty(s)); ui.addVar(new Nempty(s));
ui.addVar(new NetCharge(s)); ui.addVar(new NetCharge(s));
......
...@@ -19,5 +19,5 @@ ...@@ -19,5 +19,5 @@
#include "release.h" #include "release.h"
std::string release(void) std::string release(void)
{ {
return std::string("rel1_73_5"); return std::string("rel1_73_5dev");
} }
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