Commit 8d3b1ab7 by Francois Gygi

position constraint to lock atoms


git-svn-id: http://qboxcode.org/svn/qb/trunk@692 cba15fb0-1239-40c8-b417-11db7ca47a34
parent bc5c618e
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// PositionConstraint.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PositionConstraint.C,v 1.1 2009-04-30 23:13:39 fgygi Exp $
#include "PositionConstraint.h"
#include "AtomSet.h"
#include "Atom.h"
#include "Species.h"
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void PositionConstraint::setup(const AtomSet& atoms)
{
// find position in tau array corresponding to atom name1
is1_ = atoms.is(name1_);
ia1_ = atoms.ia(name1_);
assert(is1_>=0);
assert(ia1_>=0);
}
////////////////////////////////////////////////////////////////////////////////
void PositionConstraint::update(double dt)
{
// nothing to update
}
////////////////////////////////////////////////////////////////////////////////
bool PositionConstraint::enforce_r(const vector<vector<double> > &r0,
vector<vector<double> > &rp) const
{
const double* pr1 = &r0[is1_][3*ia1_];
D3vector r1(pr1);
double* pr1p = &rp[is1_][3*ia1_];
D3vector r1p(pr1p);
double sigma = length(r1p-r1);
if ( sigma < tol_ ) return true;
pr1p[0] = pr1[0];
pr1p[1] = pr1[1];
pr1p[2] = pr1[2];
return false;
}
////////////////////////////////////////////////////////////////////////////////
bool PositionConstraint::enforce_v(const vector<vector<double> > &r0,
vector<vector<double> > &v0) const
{
const double* pr1 = &r0[is1_][3*ia1_];
D3vector r1(pr1);
double* pv1 = &v0[is1_][3*ia1_];
D3vector v1(pv1);
const double err = length(v1);
if ( err < tol_ ) return true;
pv1[0] = 0.0;
pv1[1] = 0.0;
pv1[2] = 0.0;
return false;
}
////////////////////////////////////////////////////////////////////////////////
void PositionConstraint::compute_force(const vector<vector<double> > &r0,
const vector<vector<double> > &f)
{
const double* pr1 = &r0[is1_][3*ia1_];
D3vector r1(pr1);
const double* pf1 = &f[is1_][3*ia1_];
D3vector f1(pf1);
force_ = length(f1);
}
////////////////////////////////////////////////////////////////////////////////
ostream& PositionConstraint::print( ostream &os )
{
os.setf(ios::left,ios::adjustfield);
os << " <constraint name=\"" << name();
os << "\" type=\"" << type();
os << "\" atoms=\"" << name1_ << "\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << " value=\"" << setprecision(6) << 0;
os << "\" velocity=\"" << setprecision(6) << 0 << "\"\n";
os << " force=\"" << setprecision(6) << force_;
os << "\" weight=\"" << setprecision(6) << weight_ << "\"/>";
return os;
}
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// PositionConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PositionConstraint.h,v 1.1 2009-04-30 23:13:39 fgygi Exp $
#ifndef POSITIONCONSTRAINT_H
#define POSITIONCONSTRAINT_H
#include "Constraint.h"
#include <cassert>
#include <cmath> // fabs
class AtomSet;
class PositionConstraint : public Constraint
{
std::string name1_;
int ia1_, is1_;
double force_, weight_, tol_;
public:
PositionConstraint(std::string name, std::string name1, double tolerance):
name1_(name1), tol_(tolerance)
{
name_ = name;
names_.resize(1);
names_[0] = name1_;
force_ = 0.0;
weight_ = 1.0;
}
~PositionConstraint(void) {}
std::string type(void) const { return "position"; }
double value(void) const { return 0.0; }
double velocity(void) const { return 0.0; }
double force(void) const { return force_; }
double weight(void) const { return weight_; }
double tolerance(void) const { return tol_; }
void set_value(double value)
{
// no value to set
}
void set_velocity(double velocity)
{
// no value to set
}
void setup(const AtomSet& atoms);
void update(double dt);
bool enforce_r(const std::vector<std::vector<double> > &r0,
std::vector<std::vector<double> > &rp) const;
bool enforce_v(const std::vector<std::vector<double> > &r0,
std::vector<std::vector<double> > &v0) const;
void compute_force(const std::vector<std::vector<double> > &r0,
const std::vector<std::vector<double> > &f);
std::ostream& print( std::ostream& os );
};
#endif
Markdown is supported
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