Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79180431
colvaratoms.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Aug 25, 03:55
Size
31 KB
Mime Type
text/x-c
Expires
Tue, Aug 27, 03:55 (2 d)
Engine
blob
Format
Raw Data
Handle
20135373
Attached To
rLAMMPS lammps
colvaratoms.cpp
View Options
// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvaratoms.h"
cvm::atom::atom()
{
index = -1;
id = -1;
reset_data();
}
cvm::atom::atom(int atom_number)
{
colvarproxy *p = cvm::proxy;
index = p->init_atom(atom_number);
if (cvm::debug()) {
cvm::log("The index of this atom in the colvarproxy arrays is "+
cvm::to_str(index)+".\n");
}
id = p->get_atom_id(index);
update_mass();
reset_data();
}
cvm::atom::atom(cvm::residue_id const &residue,
std::string const &atom_name,
std::string const &segment_id)
{
colvarproxy *p = cvm::proxy;
index = p->init_atom(residue, atom_name, segment_id);
if (cvm::debug()) {
cvm::log("The index of this atom in the colvarproxy_namd arrays is "+
cvm::to_str(index)+".\n");
}
id = p->get_atom_id(index);
update_mass();
reset_data();
}
cvm::atom::atom(atom const &a)
: index(a.index)
{
id = (cvm::proxy)->get_atom_id(index);
update_mass();
reset_data();
}
cvm::atom::~atom()
{
if (index >= 0) {
(cvm::proxy)->clear_atom(index);
}
}
// TODO change this arrangement
// 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_in)
{
key = key_in;
cvm::log("Defining atom group \"" + key + "\".\n");
init();
// real work is done by parse
parse(conf);
setup();
}
cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms_in)
{
init();
atoms = atoms_in;
setup();
}
cvm::atom_group::atom_group()
{
init();
}
cvm::atom_group::~atom_group()
{
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->clear_atom_group(index);
index = -1;
}
if (fitting_group) {
delete fitting_group;
fitting_group = NULL;
}
}
int cvm::atom_group::add_atom(cvm::atom const &a)
{
if (a.id < 0) {
return COLVARS_ERROR;
}
for (size_t i = 0; i < atoms_ids.size(); i++) {
if (atoms_ids[i] == a.id) {
if (cvm::debug())
cvm::log("Discarding doubly counted atom with number "+
cvm::to_str(a.id+1)+".\n");
return COLVARS_OK;
}
}
// for consistency with add_atom_id(), we update the list as well
atoms_ids.push_back(a.id);
atoms.push_back(a);
total_mass += a.mass;
total_charge += a.charge;
return COLVARS_OK;
}
int cvm::atom_group::add_atom_id(int aid)
{
if (aid < 0) {
return COLVARS_ERROR;
}
for (size_t i = 0; i < atoms_ids.size(); i++) {
if (atoms_ids[i] == aid) {
if (cvm::debug())
cvm::log("Discarding doubly counted atom with number "+
cvm::to_str(aid+1)+".\n");
return COLVARS_OK;
}
}
atoms_ids.push_back(aid);
return COLVARS_OK;
}
int cvm::atom_group::remove_atom(cvm::atom_iter ai)
{
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: cannot remove atoms from a scalable group.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
if (!this->size()) {
cvm::error("Error: trying to remove an atom from an empty group.\n", INPUT_ERROR);
return COLVARS_ERROR;
} else {
total_mass -= ai->mass;
total_charge -= ai->charge;
atoms_ids.erase(atoms_ids.begin() + (ai - atoms.begin()));
atoms.erase(ai);
}
return COLVARS_OK;
}
int cvm::atom_group::init()
{
if (!key.size()) key = "atoms";
atoms.clear();
// TODO: check with proxy whether atom forces etc are available
init_ag_requires();
index = -1;
b_center = false;
b_rotate = false;
b_user_defined_fit = false;
b_fit_gradients = false;
fitting_group = NULL;
noforce = false;
total_mass = 0.0;
total_charge = 0.0;
cog.reset();
com.reset();
return COLVARS_OK;
}
int cvm::atom_group::setup()
{
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
ai->update_mass();
ai->update_charge();
}
update_total_mass();
update_total_charge();
return COLVARS_OK;
}
void cvm::atom_group::update_total_mass()
{
if (b_dummy) {
total_mass = 1.0;
return;
}
if (is_enabled(f_ag_scalable)) {
total_mass = (cvm::proxy)->get_atom_group_mass(index);
} else {
total_mass = 0.0;
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
total_mass += ai->mass;
}
}
}
void cvm::atom_group::reset_mass(std::string &name, int i, int j)
{
update_total_mass();
cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+
cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+
" atoms: total mass = "+cvm::to_str(total_mass)+".\n");
}
void cvm::atom_group::update_total_charge()
{
if (b_dummy) {
total_charge = 0.0;
return;
}
if (is_enabled(f_ag_scalable)) {
total_charge = (cvm::proxy)->get_atom_group_charge(index);
} else {
total_charge = 0.0;
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
total_charge += ai->charge;
}
}
}
int cvm::atom_group::parse(std::string const &conf)
{
std::string group_conf;
// TODO move this to the cvc class constructor/init
// 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.c_str(), 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 \""+key+
"\" is set, but has no definition.\n",
INPUT_ERROR);
return COLVARS_ERROR;
}
cvm::increase_depth();
cvm::log("Initializing atom group \""+key+"\".\n");
description = "atom group " + key;
// 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;
int parse_error = COLVARS_OK;
{
std::string numbers_conf = "";
size_t pos = 0;
while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {
cvm::combine_errors(parse_error, add_atom_numbers(numbers_conf));
numbers_conf = "";
}
}
{
std::string index_group_name;
if (get_keyval(group_conf, "indexGroup", index_group_name)) {
// use an index group from the index file read globally
cvm::combine_errors(parse_error, add_index_group(index_group_name));
}
}
{
std::string range_conf = "";
size_t pos = 0;
while (key_lookup(group_conf, "atomNumbersRange",
range_conf, pos)) {
cvm::combine_errors(parse_error, add_atom_numbers_range(range_conf));
range_conf = "";
}
}
{
std::vector<std::string> psf_segids;
get_keyval(group_conf, "psfSegID", psf_segids, std::vector<std::string>());
std::vector<std::string>::iterator psii;
for (psii = psf_segids.begin(); psii < psf_segids.end(); ++psii) {
if ( (psii->size() == 0) || (psii->size() > 4) ) {
cvm::error("Error: invalid PSF segment identifier provided, \""+
(*psii)+"\".\n", INPUT_ERROR);
}
}
std::string range_conf = "";
size_t pos = 0;
size_t range_count = 0;
psii = psf_segids.begin();
while (key_lookup(group_conf, "atomNameResidueRange",
range_conf, pos)) {
range_count++;
if (psf_segids.size() && (range_count > psf_segids.size())) {
cvm::error("Error: more instances of \"atomNameResidueRange\" than "
"values of \"psfSegID\".\n", INPUT_ERROR);
} else {
cvm::combine_errors(parse_error, add_atom_name_residue_range(psf_segids.size() ?
*psii : std::string(""), range_conf));
if (psf_segids.size()) psii++;
}
range_conf = "";
}
}
{
// read the atoms from a file
std::string atoms_file_name;
if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""))) {
std::string atoms_col;
if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""))) {
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);
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 (parse_error || cvm::get_error()) return (parse_error || cvm::get_error());
// checks of doubly-counted atoms have been handled by add_atom() already
if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos())) {
b_dummy = true;
// note: atoms_ids.size() is used here in lieu of atoms.size(),
// which can be empty for scalable groups
if (atoms_ids.size()) {
cvm::error("Error: cannot set up group \""+
key+"\" as a dummy atom "
"and provide it with atom definitions.\n", INPUT_ERROR);
}
} else {
b_dummy = false;
if (!(atoms_ids.size())) {
cvm::error("Error: no atoms defined for atom group \""+
key+"\".\n", INPUT_ERROR);
}
// 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)) {
noforce = !enable_forces;
} else {
get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent);
}
}
if (is_enabled(f_ag_scalable) && !b_dummy) {
index = (cvm::proxy)->init_atom_group(atoms_ids);
}
cvm::combine_errors(parse_error, parse_fitting_options(group_conf));
// TODO move this to colvarparse object
check_keywords(group_conf, key.c_str());
if (cvm::get_error()) {
cvm::error("Error setting up atom group \""+key+"\".");
return COLVARS_ERROR;
}
// Calculate all required properties (such as total mass)
setup();
if (cvm::debug())
cvm::log("Done initializing atom group \""+key+"\".\n");
cvm::log("Atom group \""+key+"\" defined, "+
cvm::to_str(atoms_ids.size())+" atoms initialized: total mass = "+
cvm::to_str(total_mass)+", total charge = "+
cvm::to_str(total_charge)+".\n");
cvm::decrease_depth();
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
{
std::vector<int> atom_indexes;
if (numbers_conf.size()) {
std::istringstream is(numbers_conf);
int ia;
while (is >> ia) {
atom_indexes.push_back(ia);
}
}
if (atom_indexes.size()) {
atoms_ids.reserve(atoms_ids.size()+atom_indexes.size());
if (is_enabled(f_ag_scalable)) {
for (size_t i = 0; i < atom_indexes.size(); i++) {
add_atom_id((cvm::proxy)->check_atom_id(atom_indexes[i]));
}
} else {
// if we are handling the group on rank 0, better allocate the vector in one shot
atoms.reserve(atoms.size()+atom_indexes.size());
for (size_t i = 0; i < atom_indexes.size(); i++) {
add_atom(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;
}
return COLVARS_OK;
}
int cvm::atom_group::add_index_group(std::string const &index_group_name)
{
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;
}
atoms_ids.reserve(atoms_ids.size()+index_groups_i->size());
if (is_enabled(f_ag_scalable)) {
for (size_t i = 0; i < index_groups_i->size(); i++) {
add_atom_id((cvm::proxy)->check_atom_id((*index_groups_i)[i]));
}
} else {
atoms.reserve(atoms.size()+index_groups_i->size());
for (size_t i = 0; i < index_groups_i->size(); i++) {
add_atom(cvm::atom((*index_groups_i)[i]));
}
}
if (cvm::get_error())
return COLVARS_ERROR;
return COLVARS_OK;
}
int cvm::atom_group::add_atom_numbers_range(std::string const &range_conf)
{
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) ) {
atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
if (is_enabled(f_ag_scalable)) {
for (int anum = initial; anum <= final; anum++) {
add_atom_id((cvm::proxy)->check_atom_id(anum));
}
} else {
atoms.reserve(atoms.size() + (final - initial + 1));
for (int anum = initial; anum <= final; anum++) {
add_atom(cvm::atom(anum));
}
}
}
if (cvm::get_error()) return COLVARS_ERROR;
} else {
cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+
range_conf+"\".\n", INPUT_ERROR);
return COLVARS_ERROR;
}
return COLVARS_OK;
}
int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
std::string const &range_conf)
{
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) ) {
atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
if (is_enabled(f_ag_scalable)) {
for (int resid = initial; resid <= final; resid++) {
add_atom_id((cvm::proxy)->check_atom_id(resid, atom_name, psf_segid));
}
} else {
atoms.reserve(atoms.size() + (final - initial + 1));
for (int resid = initial; resid <= final; resid++) {
add_atom(cvm::atom(resid, atom_name, psf_segid));
}
}
if (cvm::get_error()) return COLVARS_ERROR;
} else {
cvm::error("Error: cannot parse definition for \""
"atomNameResidueRange\", \""+
range_conf+"\".\n");
return COLVARS_ERROR;
}
} else {
cvm::error("Error: atomNameResidueRange with empty definition.\n");
return COLVARS_ERROR;
}
return COLVARS_OK;
}
int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
{
bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false);
// 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);
if (b_center || b_rotate) {
if (b_dummy)
cvm::error("Error: centerReference or rotateReference "
"cannot be defined for a dummy atom.\n");
bool b_ref_pos_group = false;
if (key_lookup(group_conf, "refPositionsGroup")) {
b_ref_pos_group = true;
cvm::log("Warning: keyword \"refPositionsGroup\" is deprecated: please use \"fittingGroup\" instead.\n");
}
if (b_ref_pos_group || key_lookup(group_conf, "fittingGroup")) {
// instead of this group, define another group to compute the fit
if (fitting_group) {
cvm::error("Error: the atom group \""+
key+"\" has already a reference group "
"for the rototranslational fit, which was communicated by the "
"colvar component. You should not use fittingGroup "
"in this case.\n");
}
cvm::log("Within atom group \""+key+"\":\n");
fitting_group = b_ref_pos_group ?
new atom_group(group_conf, "refPositionsGroup") :
new atom_group(group_conf, "fittingGroup");
// regardless of the configuration, fit gradients must be calculated by fittingGroup
fitting_group->b_fit_gradients = this->b_fit_gradients;
}
atom_group *group_for_fit = fitting_group ? fitting_group : this;
get_keyval(group_conf, "refPositions", ref_pos, ref_pos);
std::string ref_pos_file;
if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""))) {
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(""))) {
// if provided, use PDB column to select coordinates
bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0);
if (found && ref_pos_col_value == 0.0) {
cvm::error("Error: refPositionsColValue, "
"if provided, must be non-zero.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
} 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 \""+
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 {
cvm::error("Error: no reference positions provided.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
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 \""+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 fittingGroup (or a different definition for it if already defined) "
"to align the coordinates.\n");
// initialize rot member data
rot.request_group1_gradients(group_for_fit->size());
}
}
return 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;
for (cvm::atom_iter 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()
{
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 (fitting_group)
fitting_group->read_positions();
}
int cvm::atom_group::calc_required_properties()
{
// TODO check if the com is needed?
calc_center_of_mass();
calc_center_of_geometry();
if (!is_enabled(f_ag_scalable)) {
if (b_center || b_rotate) {
if (fitting_group) {
fitting_group->calc_center_of_geometry();
}
calc_apply_roto_translation();
// update COM and COG after fitting
calc_center_of_geometry();
calc_center_of_mass();
if (fitting_group) {
fitting_group->calc_center_of_geometry();
}
}
}
// TODO calculate elements of scalable cvc's here before reduction
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
void cvm::atom_group::calc_apply_roto_translation()
{
// store the laborarory-frame COGs for when they are needed later
cog_orig = this->center_of_geometry();
if (fitting_group) {
fitting_group->cog_orig = fitting_group->center_of_geometry();
}
if (b_center) {
// center on the origin first
cvm::atom_pos const rpg_cog = fitting_group ?
fitting_group->center_of_geometry() : this->center_of_geometry();
apply_translation(-1.0 * rpg_cog);
if (fitting_group) {
fitting_group->apply_translation(-1.0 * rpg_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(fitting_group ?
fitting_group->positions() :
this->positions(),
ref_pos);
cvm::atom_iter ai;
for (ai = this->begin(); ai != this->end(); ai++) {
ai->pos = rot.rotate(ai->pos);
}
if (fitting_group) {
for (ai = fitting_group->begin(); ai != fitting_group->end(); ai++) {
ai->pos = rot.rotate(ai->pos);
}
}
}
if (b_center) {
// align with the center of geometry of ref_pos
apply_translation(ref_pos_cog);
if (fitting_group) {
fitting_group->apply_translation(ref_pos_cog);
}
}
// update of COM and COG is done from the calling routine
}
void cvm::atom_group::apply_translation(cvm::rvector const &t)
{
if (b_dummy) {
cvm::error("Error: cannot translate the coordinates of a dummy atom group.\n", INPUT_ERROR);
return;
}
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: cannot translate the coordinates of a scalable atom group.\n", INPUT_ERROR);
return;
}
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->pos += t;
}
}
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();
}
}
}
// TODO make this a calc function
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();
}
}
}
int cvm::atom_group::calc_center_of_geometry()
{
if (b_dummy) {
cog = dummy_atom_pos;
} else {
cog.reset();
for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
cog += ai->pos;
}
cog /= this->size();
}
return COLVARS_OK;
}
int cvm::atom_group::calc_center_of_mass()
{
if (b_dummy) {
com = dummy_atom_pos;
if (cvm::debug()) {
cvm::log("Dummy atom center of mass = "+cvm::to_str(com)+"\n");
}
} else if (is_enabled(f_ag_scalable)) {
com = (cvm::proxy)->get_atom_group_com(index);
} else {
com.reset();
for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
com += ai->mass * ai->pos;
}
com /= total_mass;
}
return COLVARS_OK;
}
int cvm::atom_group::calc_dipole(cvm::atom_pos const &com)
{
if (b_dummy) {
cvm::error("Error: trying to compute the dipole of an empty group.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
dip.reset();
for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
dip += ai->charge * (ai->pos - com);
}
return COLVARS_OK;
}
void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
{
if (b_dummy) return;
if (is_enabled(f_ag_scalable)) {
scalar_com_gradient = grad;
return;
}
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->grad = (ai->mass/total_mass) * grad;
}
}
void cvm::atom_group::calc_fit_gradients()
{
if (b_dummy) return;
if (cvm::debug())
cvm::log("Calculating fit gradients.\n");
atom_group *group_for_fit = fitting_group ? fitting_group : this;
if (b_center) {
// add the center of geometry contribution to the gradients
cvm::rvector atom_grad;
for (size_t i = 0; i < this->size(); i++) {
atom_grad += atoms[i].grad;
}
if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad);
atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
for (size_t j = 0; j < group_for_fit->size(); j++) {
group_for_fit->fit_gradients[j] = atom_grad;
}
}
if (b_rotate) {
// add the rotation matrix contribution to the gradients
cvm::rotation const rot_inv = rot.inverse();
for (size_t i = 0; i < this->size(); i++) {
// compute centered, unrotated position
cvm::atom_pos const pos_orig =
rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
// 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, atoms[i].grad);
for (size_t j = 0; j < group_for_fit->size(); j++) {
// multiply by {\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", INPUT_ERROR);
}
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic positions are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
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", INPUT_ERROR);
}
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic positions are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
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", INPUT_ERROR);
}
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic velocities are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
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", INPUT_ERROR);
}
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic system forces are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
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;
}
// TODO make this an accessor
cvm::rvector cvm::atom_group::system_force() const
{
if (b_dummy) {
cvm::error("Error: total system forces are not available "
"from a dummy atom group.\n", INPUT_ERROR);
}
if (is_enabled(f_ag_scalable)) {
return (cvm::proxy)->get_atom_group_system_force(index);
}
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 (cvm::debug()) {
log("Communicating a colvar force from atom group to the MD engine.\n");
}
if (b_dummy)
return;
if (noforce) {
cvm::error("Error: sending a force to a group that has "
"\"enableForces\" set to off.\n");
return;
}
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient);
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 = fitting_group ? fitting_group : this;
// Fit gradients are already calculated in "laboratory" frame
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 (is_enabled(f_ag_scalable)) {
(cvm::proxy)->apply_atom_group_force(index, force);
}
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/total_mass) * force));
}
} else {
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->apply_force((ai->mass/total_mass) * force);
}
}
}
// Static members
std::vector<cvm::deps::feature *> cvm::atom_group::ag_features;
Event Timeline
Log In to Comment