Commit 53dd4a57 authored by Francois Gygi's avatar Francois Gygi
Browse files

Add IonicStepper.cpp file

parent f8949045
////////////////////////////////////////////////////////////////////////////////
//
// 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
public:
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();
}
IonicStepper(Sample& s);
double r0(int is, int i) const { return r0_[is][i]; }
double v0(int is, int i) const { return v0_[is][i]; }
......@@ -82,80 +57,15 @@ class IonicStepper
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;
}
D3vector rcm(const std::vector<std::vector<double> >& r);
// 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;
}
}
}
std::vector<std::vector<double> >& r2);
// 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;
}
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)
{
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_);
}
void reset_vcm(std::vector<std::vector<double> >& v);
void setup_constraints(void) { constraints_.setup(atoms_); }
virtual void compute_r(double e0,
const std::vector<std::vector< double> >& f0) = 0;
virtual void compute_v(double e0,
......@@ -166,6 +76,5 @@ class IonicStepper
virtual double temp(void) const { return 0.0; }
virtual ~IonicStepper() {}
};
#endif
......@@ -50,8 +50,8 @@ OBJECTS=MPIdata.o qb.o AtomSet.o Atom.o Species.o \
SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
XMLGFPreprocessor.o Base64Transcoder.o \
CPSampleStepper.o BOSampleStepper.o \
MDWavefunctionStepper.o SDIonicStepper.o MDIonicStepper.o \
BMDIonicStepper.o \
MDWavefunctionStepper.o IonicStepper.o \
SDIonicStepper.o MDIonicStepper.o BMDIonicStepper.o \
CellStepper.o SDCellStepper.o \
CGCellStepper.o ConfinementPotential.o Preconditioner.o \
AndersonMixer.o SDAIonicStepper.o CGIonicStepper.o \
......
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