colvaratoms.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Nov 18, 06:46

colvaratoms.cpp
View Options

	/// -- c++ --

	#include "colvarmodule.h"
	#include "colvarparse.h"
	#include "colvaratoms.h"


	// member functions of the "atom" class depend tightly on the MD interface, and are
	// thus defined in colvarproxy_xxx.cpp

	// in this file only atom_group functions are defined


	// Note: "conf" is the configuration of the cvc who is using this atom group;
	// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
	cvm::atom_group::atom_group(std::string const &conf,
	char const *key)
	: b_center(false), b_rotate(false), b_user_defined_fit(false),
	b_fit_gradients(false),
	ref_pos_group(NULL),
	noforce(false)
	{
	cvm::log("Defining atom group \""+
	std::string(key)+"\".\n");
	// real work is done by parse
	parse(conf, key);
	}


	cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms)
	: b_dummy(false), b_center(false), b_rotate(false),
	b_fit_gradients(false), ref_pos_group(NULL),
	noforce(false)
	{
	this->reserve(atoms.size());
	for (size_t i = 0; i < atoms.size(); i++) {
	this->push_back(atoms[i]);
	}
	total_mass = 0.0;
	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	total_mass += ai->mass;
	}
	}


	cvm::atom_group::atom_group()
	: b_dummy(false), b_center(false), b_rotate(false),
	b_user_defined_fit(false), b_fit_gradients(false),
	ref_pos_group(NULL), noforce(false)
	{
	total_mass = 0.0;
	}


	cvm::atom_group::~atom_group()
	{
	if (ref_pos_group) {
	delete ref_pos_group;
	ref_pos_group = NULL;
	}
	}


	void cvm::atom_group::add_atom(cvm::atom const &a)
	{
	if (b_dummy) {
	cvm::error("Error: cannot add atoms to a dummy group.\n", INPUT_ERROR);
	} else {
	this->push_back(a);
	total_mass += a.mass;
	}
	}


	void cvm::atom_group::reset_mass(std::string &name, int i, int j)
	{
	total_mass = 0.0;
	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	total_mass += ai->mass;
	}
	cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+
	cvm::to_str(j)+". "+ cvm::to_str(this->size())+
	" atoms: total mass = "+cvm::to_str(this->total_mass)+".\n");
	}

	int cvm::atom_group::parse(std::string const &conf,
	char const *key)
	{
	std::string group_conf;

	// save_delimiters is set to false for this call, because "conf" is
	// not the config string of this group, but of its parent object
	// (which has already taken care of the delimiters)
	save_delimiters = false;
	key_lookup(conf, key, group_conf, dummy_pos);
	// restoring the normal value, because we do want keywords checked
	// inside "group_conf"
	save_delimiters = true;

	if (group_conf.size() == 0) {
	cvm::error("Error: atom group \""+std::string(key)+
	"\" is set, but has no definition.\n",
	INPUT_ERROR);
	return COLVARS_ERROR;
	}

	cvm::increase_depth();

	cvm::log("Initializing atom group \""+std::string(key)+"\".\n");

	// whether or not to include messages in the log
	// colvarparse::Parse_Mode mode = parse_silent;
	// {
	// bool b_verbose;
	// get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
	// if (b_verbose) mode = parse_normal;
	// }
	colvarparse::Parse_Mode mode = parse_normal;

	{
	// std::vector<int> atom_indexes;
	std::string numbers_conf = "";
	size_t pos = 0;
	std::vector<int> atom_indexes;
	while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {
	if (numbers_conf.size()) {
	std::istringstream is(numbers_conf);
	int ia;
	while (is >> ia) {
	atom_indexes.push_back(ia);
	}
	}

	if (atom_indexes.size()) {
	this->reserve(this->size()+atom_indexes.size());
	for (size_t i = 0; i < atom_indexes.size(); i++) {
	this->push_back(cvm::atom(atom_indexes[i]));
	}
	if (cvm::get_error()) return COLVARS_ERROR;
	} else {
	cvm::error("Error: no numbers provided for \""
	"atomNumbers\".\n", INPUT_ERROR);
	return COLVARS_ERROR;
	}

	atom_indexes.clear();
	}

	std::string index_group_name;
	if (get_keyval(group_conf, "indexGroup", index_group_name)) {
	// use an index group from the index file read globally
	std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
	std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();
	for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) {
	if (*names_i == index_group_name)
	break;
	}
	if (names_i == cvm::index_group_names.end()) {
	cvm::error("Error: could not find index group "+
	index_group_name+" among those provided by the index file.\n",
	INPUT_ERROR);
	return COLVARS_ERROR;
	}
	this->reserve(index_groups_i->size());
	for (size_t i = 0; i < index_groups_i->size(); i++) {
	this->push_back(cvm::atom((*index_groups_i)[i]));
	}
	if (cvm::get_error()) return COLVARS_ERROR;
	}
	}

	{
	std::string range_conf = "";
	size_t pos = 0;
	while (key_lookup(group_conf, "atomNumbersRange",
	range_conf, pos)) {

	if (range_conf.size()) {
	std::istringstream is(range_conf);
	int initial, final;
	char dash;
	if ( (is >> initial) && (initial > 0) &&
	(is >> dash) && (dash == '-') &&
	(is >> final) && (final > 0) ) {
	for (int anum = initial; anum <= final; anum++) {
	this->push_back(cvm::atom(anum));
	}
	if (cvm::get_error()) return COLVARS_ERROR;
	range_conf = "";
	continue;
	}
	}

	cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+
	range_conf+"\".\n", INPUT_ERROR);
	}
	}

	std::vector<std::string> psf_segids;
	get_keyval(group_conf, "psfSegID", psf_segids, std::vector<std::string> (), mode);
	for (std::vector<std::string>::iterator psii = psf_segids.begin();
	psii < psf_segids.end(); ++psii) {

	if ( (psii->size() == 0) \|\| (psii->size() > 4) ) {
	cvm::error("Error: invalid segmend identifier provided, \""+
	(*psii)+"\".\n", INPUT_ERROR);
	}
	}

	{
	std::string range_conf = "";
	size_t pos = 0;
	size_t range_count = 0;
	std::vector<std::string>::iterator psii = psf_segids.begin();
	while (key_lookup(group_conf, "atomNameResidueRange",
	range_conf, pos)) {
	range_count++;

	if (range_count > psf_segids.size()) {
	cvm::error("Error: more instances of \"atomNameResidueRange\" than "
	"values of \"psfSegID\".\n", INPUT_ERROR);
	}

	std::string const &psf_segid = psf_segids.size() ? *psii : std::string("");

	if (range_conf.size()) {

	std::istringstream is(range_conf);
	std::string atom_name;
	int initial, final;
	char dash;
	if ( (is >> atom_name) && (atom_name.size()) &&
	(is >> initial) && (initial > 0) &&
	(is >> dash) && (dash == '-') &&
	(is >> final) && (final > 0) ) {
	for (int resid = initial; resid <= final; resid++) {
	this->push_back(cvm::atom(resid, atom_name, psf_segid));
	}
	if (cvm::get_error()) return COLVARS_ERROR;
	range_conf = "";
	} else {
	cvm::error("Error: cannot parse definition for \""
	"atomNameResidueRange\", \""+
	range_conf+"\".\n");
	}

	} else {
	cvm::error("Error: atomNameResidueRange with empty definition.\n");
	}

	if (psf_segid.size())
	++psii;
	}
	}

	{
	// read the atoms from a file
	std::string atoms_file_name;
	if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""), mode)) {

	std::string atoms_col;
	if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""), mode)) {
	cvm::error("Error: parameter atomsCol is required if atomsFile is set.\n",
	INPUT_ERROR);
	}

	double atoms_col_value;
	bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0, mode);
	if (atoms_col_value_defined && (!atoms_col_value)) {
	cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR);
	}

	cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
	}
	}

	// Catch any errors from all the initialization steps above
	if (cvm::get_error()) return COLVARS_ERROR;

	for (std::vector<cvm::atom>::iterator a1 = this->begin();
	a1 != this->end(); ++a1) {
	std::vector<cvm::atom>::iterator a2 = a1;
	++a2;
	for ( ; a2 != this->end(); ++a2) {
	if (a1->id == a2->id) {
	if (cvm::debug())
	cvm::log("Discarding doubly counted atom with number "+
	cvm::to_str(a1->id+1)+".\n");
	a2 = this->erase(a2);
	if (a2 == this->end())
	break;
	}
	}
	}

	if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) {
	b_dummy = true;
	this->total_mass = 1.0;
	} else
	b_dummy = false;

	if (b_dummy && (this->size())) {
	cvm::error("Error: cannot set up group \""+
	std::string(key)+"\" as a dummy atom "
	"and provide it with atom definitions.\n", INPUT_ERROR);
	}

	#if(! defined(COLVARS_STANDALONE))
	if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) {
	cvm::error("Error: no atoms defined for atom group \""+
	std::string(key)+"\".\n");
	}
	#endif

	if (!b_dummy) {

	// calculate total mass (TODO: this is the step that most needs deferred re-initialization)
	this->total_mass = 0.0;
	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	this->total_mass += ai->mass;
	}

	// whether these atoms will ever receive forces or not
	bool enable_forces = true;
	// disableForces is deprecated
	if (get_keyval(group_conf, "enableForces", enable_forces, true, mode)) {
	noforce = !enable_forces;
	} else {
	get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent);
	}
	}

	// FITTING OPTIONS

	bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false, mode);
	bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false, mode);
	// is the user setting explicit options?
	b_user_defined_fit = b_defined_center \|\| b_defined_rotate;

	get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true, mode);

	if (b_center \|\| b_rotate) {

	if (b_dummy)
	cvm::error("Error: centerReference or rotateReference "
	"cannot be defined for a dummy atom.\n");

	if (key_lookup(group_conf, "refPositionsGroup")) {
	// instead of this group, define another group to compute the fit
	if (ref_pos_group) {
	cvm::error("Error: the atom group \""+
	std::string(key)+"\" has already a reference group "
	"for the rototranslational fit, which was communicated by the "
	"colvar component. You should not use refPositionsGroup "
	"in this case.\n");
	}
	cvm::log("Within atom group \""+std::string(key)+"\":\n");
	ref_pos_group = new atom_group(group_conf, "refPositionsGroup");

	// regardless of the configuration, fit gradients must be calculated by refPositionsGroup
	ref_pos_group->b_fit_gradients = this->b_fit_gradients;
	this->b_fit_gradients = false;
	}

	atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;

	get_keyval(group_conf, "refPositions", ref_pos, ref_pos, mode);

	std::string ref_pos_file;
	if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""), mode)) {

	if (ref_pos.size()) {
	cvm::error("Error: cannot specify \"refPositionsFile\" and "
	"\"refPositions\" at the same time.\n");
	}

	std::string ref_pos_col;
	double ref_pos_col_value=0.0;

	if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""), mode)) {
	// if provided, use PDB column to select coordinates
	bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode);
	if (found && ref_pos_col_value == 0.0)
	cvm::error("Error: refPositionsColValue, "
	"if provided, must be non-zero.\n");
	} else {
	// if not, rely on existing atom indices for the group
	group_for_fit->create_sorted_ids();
	ref_pos.resize(group_for_fit->size());
	}

	cvm::load_coords(ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids,
	ref_pos_col, ref_pos_col_value);
	}

	if (ref_pos.size()) {

	if (b_rotate) {
	if (ref_pos.size() != group_for_fit->size())
	cvm::error("Error: the number of reference positions provided("+
	cvm::to_str(ref_pos.size())+
	") does not match the number of atoms within \""+
	std::string(key)+
	"\" ("+cvm::to_str(group_for_fit->size())+
	"): to perform a rotational fit, "+
	"these numbers should be equal.\n", INPUT_ERROR);
	}

	// save the center of geometry of ref_pos and subtract it
	center_ref_pos();

	} else {
	#if(! defined(COLVARS_STANDALONE))
	cvm::error("Error: no reference positions provided.\n");
	#endif
	}

	if (b_fit_gradients) {
	group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0));
	rot.request_group1_gradients(group_for_fit->size());
	}

	if (b_rotate && !noforce) {
	cvm::log("Warning: atom group \""+std::string(key)+
	"\" will be aligned to a fixed orientation given by the reference positions provided. "
	"If the internal structure of the group changes too much (i.e. its RMSD is comparable "
	"to its radius of gyration), the optimal rotation and its gradients may become discontinuous. "
	"If that happens, use refPositionsGroup (or a different definition for it if already defined) "
	"to align the coordinates.\n");
	// initialize rot member data
	rot.request_group1_gradients(this->size());
	}

	}

	if (cvm::debug())
	cvm::log("Done initializing atom group with name \""+
	std::string(key)+"\".\n");

	this->check_keywords(group_conf, key);
	if (cvm::get_error()) {
	cvm::error("Error setting up atom group \""+std::string(key)+"\".");
	return COLVARS_ERROR;
	}

	cvm::log("Atom group \""+std::string(key)+"\" defined, "+
	cvm::to_str(this->size())+" atoms initialized: total mass = "+
	cvm::to_str(this->total_mass)+".\n");

	cvm::decrease_depth();

	return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
	}


	int cvm::atom_group::create_sorted_ids(void)
	{
	// Only do the work if the vector is not yet populated
	if (sorted_ids.size())
	return COLVARS_OK;

	std::list<int> temp_id_list;
	cvm::atom_iter ai;
	for (ai = this->begin(); ai != this->end(); ai++) {
	temp_id_list.push_back(ai->id);
	}
	temp_id_list.sort();
	temp_id_list.unique();
	if (temp_id_list.size() != this->size()) {
	cvm::error("Error: duplicate atom IDs in atom group? (found " +
	cvm::to_str(temp_id_list.size()) +
	" unique atom IDs instead of" +
	cvm::to_str(this->size()) + ").\n");
	return COLVARS_ERROR;
	}
	sorted_ids = std::vector<int> (temp_id_list.size());
	unsigned int id_i = 0;
	std::list<int>::iterator li;
	for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
	sorted_ids[id_i] = *li;
	id_i++;
	}
	return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
	}

	void cvm::atom_group::center_ref_pos()
	{
	// save the center of geometry of ref_pos and then subtract it from
	// them; in this way it will be possible to use ref_pos also for
	// the rotational fit
	// This is called either by atom_group::parse or by CVCs that set
	// reference positions (eg. RMSD, eigenvector)
	ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0);
	std::vector<cvm::atom_pos>::iterator pi;
	for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
	ref_pos_cog += *pi;
	}
	ref_pos_cog /= (cvm::real) ref_pos.size();
	for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
	(*pi) -= ref_pos_cog;
	}
	}

	void cvm::atom_group::read_positions()
	{
	if (b_dummy) return;

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->read_position();
	}

	if (ref_pos_group)
	ref_pos_group->read_positions();
	}

	void cvm::atom_group::calc_apply_roto_translation()
	{
	atom_group *fit_group = ref_pos_group ? ref_pos_group : this;

	if (b_center) {
	// center on the origin first
	cvm::atom_pos const cog = fit_group->center_of_geometry();
	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->pos -= cog;
	}
	}

	if (b_rotate) {
	// rotate the group (around the center of geometry if b_center is
	// true, around the origin otherwise)
	rot.calc_optimal_rotation(fit_group->positions(), ref_pos);

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->pos = rot.rotate(ai->pos);
	}
	}

	if (b_center) {
	// align with the center of geometry of ref_pos
	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->pos += ref_pos_cog;
	}
	}
	}

	void cvm::atom_group::apply_translation(cvm::rvector const &t)
	{
	if (b_dummy) return;

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->pos += t;
	}
	}

	void cvm::atom_group::apply_rotation(cvm::rotation const &rot)
	{
	if (b_dummy) return;

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->pos = rot.rotate(ai->pos);
	}
	}


	void cvm::atom_group::read_velocities()
	{
	if (b_dummy) return;

	if (b_rotate) {

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->read_velocity();
	ai->vel = rot.rotate(ai->vel);
	}

	} else {

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->read_velocity();
	}
	}
	}

	void cvm::atom_group::read_system_forces()
	{
	if (b_dummy) return;

	if (b_rotate) {

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->read_system_force();
	ai->system_force = rot.rotate(ai->system_force);
	}

	} else {

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->read_system_force();
	}
	}
	}

	cvm::atom_pos cvm::atom_group::center_of_geometry() const
	{
	if (b_dummy)
	return dummy_atom_pos;

	cvm::atom_pos cog(0.0, 0.0, 0.0);
	for (cvm::atom_const_iter ai = this->begin();
	ai != this->end(); ai++) {
	cog += ai->pos;
	}
	cog /= this->size();
	return cog;
	}

	cvm::atom_pos cvm::atom_group::center_of_mass() const
	{
	if (b_dummy)
	return dummy_atom_pos;

	cvm::atom_pos com(0.0, 0.0, 0.0);
	for (cvm::atom_const_iter ai = this->begin();
	ai != this->end(); ai++) {
	com += ai->mass * ai->pos;
	}
	com /= this->total_mass;
	return com;
	}


	void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
	{
	if (b_dummy) return;

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->grad = (ai->mass/this->total_mass) * grad;
	}
	}


	void cvm::atom_group::calc_fit_gradients()
	{
	if (b_dummy) return;

	if ((!b_center) && (!b_rotate)) return; // no fit

	if (cvm::debug())
	cvm::log("Calculating fit gradients.\n");

	atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
	group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::rvector(0.0, 0.0, 0.0));

	if (b_center) {
	// add the center of geometry contribution to the gradients
	for (size_t i = 0; i < this->size(); i++) {
	// need to bring the gradients in original frame first
	cvm::rvector const atom_grad = b_rotate ?
	(rot.inverse()).rotate((*this)[i].grad) :
	(*this)[i].grad;
	for (size_t j = 0; j < group_for_fit->size(); j++) {
	group_for_fit->fit_gradients[j] +=
	(-1.0)/(cvm::real(group_for_fit->size())) *
	atom_grad;
	}
	}
	}

	if (b_rotate) {

	// add the rotation matrix contribution to the gradients
	cvm::rotation const rot_inv = rot.inverse();
	cvm::atom_pos const cog = this->center_of_geometry();

	for (size_t i = 0; i < this->size(); i++) {

	cvm::atom_pos const pos_orig = rot_inv.rotate((b_center ? ((this)[i].pos - cog) : ((this)[i].pos)));

	for (size_t j = 0; j < group_for_fit->size(); j++) {
	// calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
	cvm::quaternion const dxdq =
	rot.q.position_derivative_inner(pos_orig, (*this)[i].grad);
	// multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients
	for (size_t iq = 0; iq < 4; iq++) {
	group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
	}
	}
	}
	}
	if (cvm::debug())
	cvm::log("Done calculating fit gradients.\n");
	}


	std::vector<cvm::atom_pos> cvm::atom_group::positions() const
	{
	if (b_dummy) {
	cvm::error("Error: positions are not available "
	"from a dummy atom group.\n");
	}

	std::vector<cvm::atom_pos> x(this->size(), 0.0);
	cvm::atom_const_iter ai = this->begin();
	std::vector<cvm::atom_pos>::iterator xi = x.begin();
	for ( ; ai != this->end(); ++xi, ++ai) {
	*xi = ai->pos;
	}
	return x;
	}

	std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted(cvm::rvector const &shift) const
	{
	if (b_dummy) {
	cvm::error("Error: positions are not available "
	"from a dummy atom group.\n");
	}

	std::vector<cvm::atom_pos> x(this->size(), 0.0);
	cvm::atom_const_iter ai = this->begin();
	std::vector<cvm::atom_pos>::iterator xi = x.begin();
	for ( ; ai != this->end(); ++xi, ++ai) {
	*xi = (ai->pos + shift);
	}
	return x;
	}

	std::vector<cvm::rvector> cvm::atom_group::velocities() const
	{
	if (b_dummy) {
	cvm::error("Error: velocities are not available "
	"from a dummy atom group.\n");
	}

	std::vector<cvm::rvector> v(this->size(), 0.0);
	cvm::atom_const_iter ai = this->begin();
	std::vector<cvm::atom_pos>::iterator vi = v.begin();
	for ( ; ai != this->end(); vi++, ai++) {
	*vi = ai->vel;
	}
	return v;
	}

	std::vector<cvm::rvector> cvm::atom_group::system_forces() const
	{
	if (b_dummy) {
	cvm::error("Error: system forces are not available "
	"from a dummy atom group.\n");
	}

	std::vector<cvm::rvector> f(this->size(), 0.0);
	cvm::atom_const_iter ai = this->begin();
	std::vector<cvm::atom_pos>::iterator fi = f.begin();
	for ( ; ai != this->end(); ++fi, ++ai) {
	*fi = ai->system_force;
	}
	return f;
	}

	cvm::rvector cvm::atom_group::system_force() const
	{
	if (b_dummy) {
	cvm::error("Error: system forces are not available "
	"from a dummy atom group.\n");
	}

	cvm::rvector f(0.0);
	for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
	f += ai->system_force;
	}
	return f;
	}


	void cvm::atom_group::apply_colvar_force(cvm::real const &force)
	{
	if (b_dummy)
	return;

	if (noforce) {
	cvm::error("Error: sending a force to a group that has "
	"\"enableForces\" set to off.\n");
	return;
	}

	if (b_rotate) {

	// rotate forces back to the original frame
	cvm::rotation const rot_inv = rot.inverse();
	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->apply_force(rot_inv.rotate(force * ai->grad));
	}

	} else {

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->apply_force(force * ai->grad);
	}
	}

	if ((b_center \|\| b_rotate) && b_fit_gradients) {

	atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;

	// add the contribution from the roto-translational fit to the gradients
	if (b_rotate) {
	// rotate forces back to the original frame
	cvm::rotation const rot_inv = rot.inverse();
	for (size_t j = 0; j < group_for_fit->size(); j++) {
	(group_for_fit)[j].apply_force(rot_inv.rotate(force group_for_fit->fit_gradients[j]));
	}
	} else {
	for (size_t j = 0; j < group_for_fit->size(); j++) {
	(group_for_fit)[j].apply_force(force group_for_fit->fit_gradients[j]);
	}
	}
	}

	}


	void cvm::atom_group::apply_force(cvm::rvector const &force)
	{
	if (b_dummy)
	return;

	if (noforce) {
	cvm::error("Error: sending a force to a group that has "
	"\"disableForces\" defined.\n");
	return;
	}

	if (b_rotate) {

	cvm::rotation const rot_inv = rot.inverse();
	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->apply_force(rot_inv.rotate((ai->mass/this->total_mass) * force));
	}

	} else {

	for (cvm::atom_iter ai = this->begin();
	ai != this->end(); ai++) {
	ai->apply_force((ai->mass/this->total_mass) * force);
	}
	}
	}


	void cvm::atom_group::apply_forces(std::vector<cvm::rvector> const &forces)
	{
	if (b_dummy)
	return;

	if (noforce)
	cvm::error("Error: sending a force to a group that has "
	"\"disableForces\" defined.\n");

	if (forces.size() != this->size()) {
	cvm::error("Error: trying to apply an array of forces to an atom "
	"group which does not have the same length.\n");
	}

	if (b_rotate) {

	cvm::rotation const rot_inv = rot.inverse();
	cvm::atom_iter ai = this->begin();
	std::vector<cvm::rvector>::const_iterator fi = forces.begin();
	for ( ; ai != this->end(); ++fi, ++ai) {
	ai->apply_force(rot_inv.rotate(*fi));
	}

	} else {

	cvm::atom_iter ai = this->begin();
	std::vector<cvm::rvector>::const_iterator fi = forces.begin();
	for ( ; ai != this->end(); ++fi, ++ai) {
	ai->apply_force(*fi);
	}
	}
	}

colvaratoms.cppNo OneTemporaryActions

File Metadata

colvaratoms.cppView Options

Event Timeline

colvaratoms.cpp
No OneTemporary
Actions

colvaratoms.cpp
View Options