ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8996 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2012-10-23 16:57:29 +00:00
parent c524604ec7
commit fdce177203
8 changed files with 263 additions and 310 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
lib/colvars/colvar.cpp
View File

@ -661,23 +661,6 @@ void colvar::disable (colvar::task const &t)
colvar::~colvar()
{
if (cvm::b_analysis) {
if (acf.size()) {
cvm::log ("Writing acf to file \""+acf_outfile+"\".\n");
std::ofstream acf_os (acf_outfile.c_str());
if (! acf_os.good())
cvm::fatal_error ("Cannot open file \""+acf_outfile+"\".\n");
write_acf (acf_os);
acf_os.close();
}
if (runave_os.good()) {
runave_os.close();
}
}
for (size_t i = 0; i < cvcs.size(); i++) {
delete cvcs[i];
}
@ -1330,6 +1313,26 @@ std::ostream & colvar::write_traj (std::ostream &os)
return os;
}
void colvar::write_output_files()
{
if (cvm::b_analysis) {
if (acf.size()) {
cvm::log ("Writing acf to file \""+acf_outfile+"\".\n");
std::ofstream acf_os (acf_outfile.c_str());
if (! acf_os.good())
cvm::fatal_error ("Cannot open file \""+acf_outfile+"\".\n");
write_acf (acf_os);
acf_os.close();
}
if (runave_os.good()) {
runave_os.close();
}
}
}
// ******************** ANALYSIS FUNCTIONS ********************

4
lib/colvars/colvar.h
View File

@ -370,6 +370,9 @@ public:
/// Write the collective variable to a restart file
std::ostream & write_restart (std::ostream &os);
/// Write output files (if defined, e.g. in analysis mode)
void write_output_files();
protected:
@ -478,7 +481,6 @@ public:
class selfcoordnum;
class h_bond;
class rmsd;
class logmsd;
class orientation_angle;
class tilt;
class spin_angle;

167
lib/colvars/colvaratoms.cpp
View File

@ -12,7 +12,7 @@
cvm::atom_group::atom_group (std::string const &conf,
char const *key,
atom_group *ref_pos_group_in)
: b_center (false), b_rotate (false),
: b_center (false), b_rotate (false), b_prevent_fitting (false),
ref_pos_group (NULL), // this is always set within parse(),
// regardless of ref_pos_group_in
noforce (false)
@ -49,12 +49,13 @@ void cvm::atom_group::parse (std::string const &conf,
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_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;
{
// get the atoms by numbers
@ -216,44 +217,49 @@ void cvm::atom_group::parse (std::string const &conf,
}
}
get_keyval (group_conf, "disableForces", noforce, false, mode);
if (!b_dummy)
get_keyval (group_conf, "disableForces", noforce, false, mode);
get_keyval (group_conf, "centerReference", b_center, false, mode);
get_keyval (group_conf, "rotateReference", b_rotate, false, mode);
if (b_center || b_rotate) {
// FITTING OPTIONS
if (b_dummy)
cvm::fatal_error ("Error: cannot set \"centerReference\" or "
"\"rotateReference\" with \"dummyAtom\".\n");
bool fit_defined_by_user =
( get_keyval (group_conf, "centerReference", b_center, false, mode) ||
get_keyval (group_conf, "rotateReference", b_rotate, false, mode) );
if ((!b_rotate) && (!b_center) && fit_defined_by_user)
b_prevent_fitting = true;
// use refPositionsGroup instead of this group as the one which is
// used to fit the coordinates
if (key_lookup (group_conf, "refPositionsGroup")) {
if (ref_pos_group) {
cvm::fatal_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");
// if ((b_center || b_rotate) && b_dummy)
// cvm::fatal_error ("Error: cannot set \"centerReference\" or "
// "\"rotateReference\" when \"dummyAtom\" is defined.\n");
// instead of this group, define another group (refPositionsGroup) to be the one
// used to fit the coordinates
if (key_lookup (group_conf, "refPositionsGroup")) {
if (ref_pos_group) {
cvm::fatal_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");
}
atom_group *ag = ref_pos_group ? ref_pos_group : this;
atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
if (get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode)) {
cvm::log ("Using reference positions from input file.\n");
if (ref_pos.size() != ag->size()) {
cvm::fatal_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 (ag->size())+").\n");
}
if (get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode)) {
if (ref_pos.size() != group_for_fit->size()) {
cvm::fatal_error ("Error: the number of reference positions provided ("+
cvm::to_str (ref_pos.size())+
") does not match the number of atoms of group \""+
std::string (key)+
"\" ("+cvm::to_str (group_for_fit->size())+").\n");
}
}
{
std::string ref_pos_file;
if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) {
@ -273,55 +279,56 @@ void cvm::atom_group::parse (std::string const &conf,
"if provided, must be non-zero.\n");
} else {
// if not, rely on existing atom indices for the group
ag->create_sorted_ids();
group_for_fit->create_sorted_ids();
}
cvm::load_coords (ref_pos_file.c_str(), ref_pos, ag->sorted_ids,
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() != ag->size())
cvm::fatal_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 (ag->size())+").\n");
}
// save the center of mass of ref_pos and then subtract it from
// them; in this way it is possible to use the coordinates for
// the rotational fit, if needed
ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0);
std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
for ( ; pi != ref_pos.end(); pi++) {
ref_pos_cog += *pi;
}
ref_pos_cog /= (cvm::real) ref_pos.size();
for (std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
pi != ref_pos.end(); pi++) {
(*pi) -= ref_pos_cog;
}
} else {
#if (! defined (COLVARS_STANDALONE))
if (!cvm::b_analysis)
cvm::fatal_error ("Error: no reference positions provided.\n");
#endif
}
if (b_rotate && !noforce) {
cvm::log ("Warning: atom group \""+std::string (key)+
"\" is set to be rotated to a reference orientation: "
"a torque different than zero on this group "
"could make the simulation unstable. "
"If this happens, set \"disableForces\" to yes "
"for this group.\n");
}
}
if (ref_pos.size()) {
if (b_rotate) {
if (ref_pos.size() != group_for_fit->size())
cvm::fatal_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");
}
// 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
ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0);
std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
for ( ; pi != ref_pos.end(); pi++) {
ref_pos_cog += *pi;
}
ref_pos_cog /= (cvm::real) ref_pos.size();
for (std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
pi != ref_pos.end(); pi++) {
(*pi) -= ref_pos_cog;
}
} else {
#if (! defined (COLVARS_STANDALONE))
cvm::fatal_error ("Error: no reference positions provided.\n");
#endif
}
if (b_rotate && !noforce) {
cvm::log ("Warning: atom group \""+std::string (key)+
"\" will be fitted automatically onto a fixed orientation: "
"in few cases, torques applied on this group "
"may make the simulation unstable. "
"If this happens, set \"disableForces\" to yes "
"for this group.\n");
// initialize rot member data
rot.request_group1_gradients (this->size());
}
if (cvm::debug())
cvm::log ("Done initializing atom group with name \""+

9
lib/colvars/colvaratoms.h
View File

@ -156,8 +156,13 @@ public:
/// Rotation between the group and its reference coordinates
cvm::rotation rot;
/// \brief In case b_center or b_rotate is true, use these reference
/// coordinates
/// \brief Indicates that b_center and b_rotate are not false by default,
/// but following the user's choice instead: cvc's will not override
/// the value of b_center or b_rotate
bool b_prevent_fitting;
/// \brief use these reference coordinates, for b_center, b_rotate, or to compute
/// certain colvar components (orientation, rmsd, etc)
std::vector<cvm::atom_pos> ref_pos;
/// \brief Center of geometry of the reference coordinates; regardless
/// of whether b_center is true, ref_pos is centered to zero at

44
lib/colvars/colvarcomp.h
View File

@ -1107,11 +1107,15 @@ public:
/// colvar::orientation \endlink to calculate the rotation matrix
/// (colvarvalue::type_scalar type, range [0:*))
class colvar::rmsd
: public colvar::orientation
: public colvar::cvc
{
protected:
/// Sum of the squares of ref_coords
cvm::real ref_pos_sum2;
/// Atom group
cvm::atom_group atoms;
/// Reference coordinates (for RMSD calculation only)
std::vector<cvm::atom_pos> ref_pos;
public:
@ -1134,39 +1138,6 @@ public:
};
/// \brief Colvar component: mean square deviation (RMSD) of a
/// group with respect to a set of reference coordinates; uses \link
/// colvar::orientation \endlink to calculate the rotation matrix
/// (colvarvalue::type_scalar type, range [0:*))
class colvar::logmsd
: public colvar::orientation
{
protected:
/// Sum of the squares of ref_coords
cvm::real ref_pos_sum2;
cvm::real MSD;
public:
/// Constructor
logmsd (std::string const &conf);
virtual inline ~logmsd() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
virtual void calc_Jacobian_derivative();
virtual void apply_force (colvarvalue const &force);
virtual cvm::real dist2 (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
// metrics functions for cvc implementations with a periodicity
@ -1305,7 +1276,6 @@ inline void colvar::spin_angle::wrap (colvarvalue &x) const
simple_scalar_dist_functions (inertia)
simple_scalar_dist_functions (inertia_z)
simple_scalar_dist_functions (rmsd)
simple_scalar_dist_functions (logmsd)
simple_scalar_dist_functions (orientation_angle)
simple_scalar_dist_functions (tilt)
simple_scalar_dist_functions (eigenvector)

300
lib/colvars/colvarcomp_distances.cpp
View File

@ -466,6 +466,13 @@ colvar::distance_inv::distance_inv (std::string const &conf)
cvm::fatal_error ("Error: negative or zero exponent provided.\n");
}
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
if (ai1->id == ai2->id)
cvm::fatal_error ("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n");
}
}
b_inverse_gradients = false;
b_Jacobian_derivative = false;
x.type (colvarvalue::type_scalar);
@ -716,17 +723,82 @@ void colvar::inertia_z::apply_force (colvarvalue const &force)
colvar::rmsd::rmsd (std::string const &conf)
: orientation (conf)
: cvc (conf)
{
b_inverse_gradients = true;
b_Jacobian_derivative = true;
function_type = "rmsd";
x.type (colvarvalue::type_scalar);
ref_pos_sum2 = 0.0;
for (size_t i = 0; i < ref_pos.size(); i++) {
ref_pos_sum2 += ref_pos[i].norm2();
parse_group (conf, "atoms", atoms);
atom_groups.push_back (&atoms);
if (atoms.b_dummy)
cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom.");
if (!atoms.b_prevent_fitting) {
// fit everything, unless the user made an explicit choice
cvm::log ("No value was specified within \"atoms\" for \"centerReference\" or "
"\"rotateReference\": enabling both.\n");
// NOTE: this won't work for class V cvc's
atoms.b_center = true;
atoms.b_rotate = true;
}
if (atoms.b_center || atoms.b_rotate) {
// request the calculation of the derivatives of the rotation defined by the atom group
atoms.rot.request_group1_gradients (atoms.size());
// request derivatives of optimal rotation wrt 2nd group for Jacobian: this is only
// required for ABF, but we do both groups here for better caching
atoms.rot.request_group2_gradients (atoms.size());
if (atoms.ref_pos_group != NULL) {
cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
"Jacobian derivatives of the RMSD will not be calculated.\n");
b_Jacobian_derivative = false;
}
}
// the following is a simplified version of the corresponding atom group options:
// to define the reference coordinates to compute this variable
if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
cvm::log ("Using reference positions from configuration file to calculate the variable.\n");
if (ref_pos.size() != atoms.size()) {
cvm::fatal_error ("Error: the number of reference positions provided ("+
cvm::to_str (ref_pos.size())+
") does not match the number of atoms of group \"atoms\" ("+
cvm::to_str (atoms.size())+").\n");
}
}
{
std::string ref_pos_file;
if (get_keyval (conf, "refPositionsFile", ref_pos_file, std::string (""))) {
if (ref_pos.size()) {
cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and "
"\"refPositions\" at the same time.\n");
}
std::string ref_pos_col;
double ref_pos_col_value;
if (get_keyval (conf, "refPositionsCol", ref_pos_col, std::string (""))) {
// if provided, use PDB column to select coordinates
bool found = get_keyval (conf, "refPositionsColValue", ref_pos_col_value, 0.0);
if (found && !ref_pos_col_value)
cvm::fatal_error ("Error: refPositionsColValue, "
"if provided, must be non-zero.\n");
} else {
// if not, rely on existing atom indices for the group
atoms.create_sorted_ids();
}
cvm::load_coords (ref_pos_file.c_str(), ref_pos, atoms.sorted_ids,
ref_pos_col, ref_pos_col_value);
}
}
}
@ -734,20 +806,28 @@ void colvar::rmsd::calc_value()
{
atoms.reset_atoms_data();
atoms.read_positions();
// rotational fit is done internally
atoms_cog = atoms.center_of_geometry();
rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
// atoms_cog = atoms.center_of_geometry();
// rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
cvm::real group_pos_sum2 = 0.0;
for (size_t i = 0; i < atoms.size(); i++) {
group_pos_sum2 += (atoms[i].pos - atoms_cog).norm2();
// cvm::real group_pos_sum2 = 0.0;
// for (size_t i = 0; i < atoms.size(); i++) {
// group_pos_sum2 += (atoms[i].pos - atoms_cog).norm2();
// }
// // value of the RMSD (Coutsias et al)
// cvm::real const MSD = 1.0/(cvm::real (atoms.size())) *
// ( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda );
// x.real_value = (MSD > 0.0) ? std::sqrt (MSD) : 0.0;
x.real_value = 0.0;
for (size_t ia = 0; ia < atoms.size(); ia++) {
x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2();
}
// value of the RMSD (Coutsias et al)
cvm::real const MSD = 1.0/(cvm::real (atoms.size())) *
( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda );
x.real_value = (MSD > 0.0) ? std::sqrt (MSD) : 0.0;
x.real_value /= cvm::real (atoms.size()); // MSD
x.real_value = (x.real_value > 0.0) ? std::sqrt (x.real_value) : 0.0;
}
@ -758,8 +838,7 @@ void colvar::rmsd::calc_gradients()
0.0;
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - atoms_cog -
rot.q.rotate (ref_pos[ia])));
atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - atoms.ref_pos[ia]));
}
}
@ -787,46 +866,47 @@ void colvar::rmsd::calc_force_invgrads()
void colvar::rmsd::calc_Jacobian_derivative()
{
// divergence of the back-rotated target coordinates
// divergence of the rotated coordinates (including only derivatives of the rotation matrix)
cvm::real divergence = 0.0;
// gradient of the rotation matrix
cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
if (atoms.b_rotate) {
// gradients of products of 2 quaternion components
cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
// gradient of the rotation matrix
cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
// gradients of products of 2 quaternion components
cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
for (size_t ia = 0; ia < atoms.size(); ia++) {
for (size_t ia = 0; ia < atoms.size(); ia++) {
// Gradient of optimal quaternion wrt current Cartesian position
cvm::vector1d< cvm::rvector, 4 > &dq = atoms.rot.dQ0_1[ia];
// Gradient of optimal quaternion wrt current Cartesian position
cvm::vector1d< cvm::rvector, 4 > &dq = rot.dQ0_2[ia];
g11 = 2.0 * (atoms.rot.q)[1]*dq[1];
g22 = 2.0 * (atoms.rot.q)[2]*dq[2];
g33 = 2.0 * (atoms.rot.q)[3]*dq[3];
g01 = (atoms.rot.q)[0]*dq[1] + (atoms.rot.q)[1]*dq[0];
g02 = (atoms.rot.q)[0]*dq[2] + (atoms.rot.q)[2]*dq[0];
g03 = (atoms.rot.q)[0]*dq[3] + (atoms.rot.q)[3]*dq[0];
g12 = (atoms.rot.q)[1]*dq[2] + (atoms.rot.q)[2]*dq[1];
g13 = (atoms.rot.q)[1]*dq[3] + (atoms.rot.q)[3]*dq[1];
g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2];
g11 = 2.0 * (rot.q)[1]*dq[1];
g22 = 2.0 * (rot.q)[2]*dq[2];
g33 = 2.0 * (rot.q)[3]*dq[3];
g01 = (rot.q)[0]*dq[1] + (rot.q)[1]*dq[0];
g02 = (rot.q)[0]*dq[2] + (rot.q)[2]*dq[0];
g03 = (rot.q)[0]*dq[3] + (rot.q)[3]*dq[0];
g12 = (rot.q)[1]*dq[2] + (rot.q)[2]*dq[1];
g13 = (rot.q)[1]*dq[3] + (rot.q)[3]*dq[1];
g23 = (rot.q)[2]*dq[3] + (rot.q)[3]*dq[2];
// Gradient of the rotation matrix wrt current Cartesian position
grad_rot_mat[0][0] = -2.0 * (g22 + g33);
grad_rot_mat[1][0] = 2.0 * (g12 + g03);
grad_rot_mat[2][0] = 2.0 * (g13 - g02);
grad_rot_mat[0][1] = 2.0 * (g12 - g03);
grad_rot_mat[1][1] = -2.0 * (g11 + g33);
grad_rot_mat[2][1] = 2.0 * (g01 + g23);
grad_rot_mat[0][2] = 2.0 * (g02 + g13);
grad_rot_mat[1][2] = 2.0 * (g23 - g01);
grad_rot_mat[2][2] = -2.0 * (g11 + g22);
// Gradient of the rotation matrix wrt current Cartesian position
grad_rot_mat[0][0] = -2.0 * (g22 + g33);
grad_rot_mat[1][0] = 2.0 * (g12 + g03);
grad_rot_mat[2][0] = 2.0 * (g13 - g02);
grad_rot_mat[0][1] = 2.0 * (g12 - g03);
grad_rot_mat[1][1] = -2.0 * (g11 + g33);
grad_rot_mat[2][1] = 2.0 * (g01 + g23);
grad_rot_mat[0][2] = 2.0 * (g02 + g13);
grad_rot_mat[1][2] = 2.0 * (g23 - g01);
grad_rot_mat[2][2] = -2.0 * (g11 + g22);
cvm::atom_pos &x = atoms[ia].pos;
cvm::atom_pos &y = ref_pos[ia];
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
divergence += grad_rot_mat[i][j][i] * y[j];
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
divergence += grad_rot_mat[i][j][i] * x[j];
}
}
}
}
@ -836,128 +916,6 @@ void colvar::rmsd::calc_Jacobian_derivative()
colvar::logmsd::logmsd (std::string const &conf)
: orientation (conf)
{
b_inverse_gradients = true;
b_Jacobian_derivative = true;
function_type = "logmsd";
x.type (colvarvalue::type_scalar);
ref_pos_sum2 = 0.0;
for (size_t i = 0; i < ref_pos.size(); i++) {
ref_pos_sum2 += ref_pos[i].norm2();
}
}
void colvar::logmsd::calc_value()
{
atoms.reset_atoms_data();
atoms.read_positions();
if (cvm::debug())
cvm::log ("colvar::logmsd: current com: "+
cvm::to_str (atoms.center_of_mass())+"\n");
atoms_cog = atoms.center_of_geometry();
rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
cvm::real group_pos_sum2 = 0.0;
for (size_t i = 0; i < atoms.size(); i++) {
group_pos_sum2 += (atoms[i].pos-atoms_cog).norm2();
}
// value of the MSD (Coutsias et al)
MSD = 1.0/(cvm::real (atoms.size())) *
( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda );
x.real_value = (MSD > 0.0) ? std::log(MSD) : 0.0;
}
void colvar::logmsd::calc_gradients()
{
cvm::real fact = (MSD > 0.0) ? 2.0/(cvm::real (atoms.size()) * MSD) : 0.0;
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = fact * (atoms[ia].pos - atoms_cog - rot.dL0_2[ia]);
}
}
void colvar::logmsd::apply_force (colvarvalue const &force)
{
if (!atoms.noforce)
atoms.apply_colvar_force (force.real_value);
}
void colvar::logmsd::calc_force_invgrads()
{
atoms.read_system_forces();
ft.real_value = 0.0;
// Note: gradient square norm is 4.0 / (N_atoms * E)
for (size_t ia = 0; ia < atoms.size(); ia++) {
ft.real_value += atoms[ia].grad * atoms[ia].system_force;
}
ft.real_value *= atoms.size() * MSD / 4.0;
}
void colvar::logmsd::calc_Jacobian_derivative()
{
// divergence of the back-rotated target coordinates
cvm::real divergence = 0.0;
// gradient of the rotation matrix
cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
// gradients of products of 2 quaternion components
cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
for (size_t ia = 0; ia < atoms.size(); ia++) {
// Gradient of optimal quaternion wrt current Cartesian position
cvm::vector1d< cvm::rvector, 4 > &dq = rot.dQ0_2[ia];
g11 = 2.0 * (rot.q)[1]*dq[1];
g22 = 2.0 * (rot.q)[2]*dq[2];
g33 = 2.0 * (rot.q)[3]*dq[3];
g01 = (rot.q)[0]*dq[1] + (rot.q)[1]*dq[0];
g02 = (rot.q)[0]*dq[2] + (rot.q)[2]*dq[0];
g03 = (rot.q)[0]*dq[3] + (rot.q)[3]*dq[0];
g12 = (rot.q)[1]*dq[2] + (rot.q)[2]*dq[1];
g13 = (rot.q)[1]*dq[3] + (rot.q)[3]*dq[1];
g23 = (rot.q)[2]*dq[3] + (rot.q)[3]*dq[2];
// Gradient of the rotation matrix wrt current Cartesian position
// Note: we are only going to use "diagonal" terms: grad_rot_mat[i][j][i]
grad_rot_mat[0][0] = -2.0 * (g22 + g33);
grad_rot_mat[1][0] = 2.0 * (g12 + g03);
grad_rot_mat[2][0] = 2.0 * (g13 - g02);
grad_rot_mat[0][1] = 2.0 * (g12 - g03);
grad_rot_mat[1][1] = -2.0 * (g11 + g33);
grad_rot_mat[2][1] = 2.0 * (g01 + g23);
grad_rot_mat[0][2] = 2.0 * (g02 + g13);
grad_rot_mat[1][2] = 2.0 * (g23 - g01);
grad_rot_mat[2][2] = -2.0 * (g11 + g22);
cvm::atom_pos &y = ref_pos[ia];
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
divergence += grad_rot_mat[i][j][i] * y[j];
}
}
}
jd.real_value = (3.0 * atoms.size() - 3.0 - divergence) / 2.0;
}
colvar::eigenvector::eigenvector (std::string const &conf)
: cvc (conf)

8
lib/colvars/colvarmodule.cpp
View File

@ -583,6 +583,14 @@ void colvarmodule::write_output_files()
this->write_restart (out);
out.close();
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_output_files();
}
cvm::decrease_depth();
// do not close to avoid problems with multiple NAMD runs
cv_traj_os.flush();
}

2
lib/colvars/colvarmodule.h
View File

@ -2,7 +2,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2012-08-08"
#define COLVARS_VERSION "2012-10-03"
#endif
#ifndef COLVARS_DEBUG