update colvars to version 2016-08-10 (3b476d91) · Commits · 郑智淋 / lammps

lib/colvars/colvar.cpp

+30 −21

Original line number	Diff line number	Diff line
		@@ -374,10 +374,19 @@ colvar::colvar(std::string const &conf)
		}

		{
		bool b_output_system_force;
		get_keyval(conf, "outputSystemForce", b_output_system_force, false);
		if (b_output_system_force) {
		enable(f_cv_output_system_force);
		bool temp;
		if (get_keyval(conf, "outputSystemForce", temp, false)) {
		cvm::error("Colvar option outputSystemForce is deprecated.\n"
		"Please use outputTotalForce, or outputSystemForce within an ABF bias.");
		return;
		}
		}

		{
		bool b_output_total_force;
		get_keyval(conf, "outputTotalForce", b_output_total_force, false);
		if (b_output_total_force) {
		enable(f_cv_output_total_force);
		}
		}

		@@ -555,8 +564,8 @@ int colvar::init_components(std::string const &conf)
		int colvar::refresh_deps()
		{
		// If enabled features are changed upstream, the features below should be refreshed
		if (is_enabled(f_cv_system_force_calc)) {
		cvm::request_system_force();
		if (is_enabled(f_cv_total_force_calc)) {
		cvm::request_total_force();
		}
		if (is_enabled(f_cv_collect_gradient) && atom_ids.size() == 0) {
		build_atom_list();
		@@ -972,17 +981,17 @@ int colvar::calc_cvc_sys_forces(int first_cvc, size_t num_cvcs)
		size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
		size_t i, cvc_count;

		if (is_enabled(f_cv_system_force_calc)) {
		if (is_enabled(f_cv_total_force_calc)) {
		if (cvm::debug())
		cvm::log("Calculating system force of colvar \""+this->name+"\".\n");
		cvm::log("Calculating total force of colvar \""+this->name+"\".\n");

		// if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
		// Disabled check to allow for explicit system force calculation
		// Disabled check to allow for explicit total force calculation
		// even with extended Lagrangian

		if (cvm::step_relative() > 0) {
		cvm::increase_depth();
		// get from the cvcs the system forces from the PREVIOUS step
		// get from the cvcs the total forces from the PREVIOUS step
		for (i = first_cvc, cvc_count = 0;
		(i < cvcs.size()) && (cvc_count < cvc_max_count);
		i++) {
		@@ -994,7 +1003,7 @@ int colvar::calc_cvc_sys_forces(int first_cvc, size_t num_cvcs)
		}

		if (cvm::debug())
		cvm::log("Done calculating system force of colvar \""+this->name+"\".\n");
		cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
		}

		return COLVARS_OK;
		@@ -1003,20 +1012,20 @@ int colvar::calc_cvc_sys_forces(int first_cvc, size_t num_cvcs)

		int colvar::collect_cvc_sys_forces()
		{
		if (is_enabled(f_cv_system_force_calc)) {
		if (is_enabled(f_cv_total_force_calc)) {
		ft.reset();

		if (cvm::step_relative() > 0) {
		// get from the cvcs the system forces from the PREVIOUS step
		// get from the cvcs the total forces from the PREVIOUS step
		for (size_t i = 0; i < cvcs.size(); i++) {
		if (!cvcs[i]->is_enabled()) continue;
		// linear combination is assumed
		ft += (cvcs[i])->system_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
		ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
		}
		}

		if (!is_enabled(f_cv_hide_Jacobian)) {
		// add the Jacobian force to the system force, and don't apply any silent
		// add the Jacobian force to the total force, and don't apply any silent
		// correction internally: biases such as colvarbias_abf will handle it
		ft += fj;
		}
		@@ -1088,7 +1097,7 @@ int colvar::calc_colvar_properties()
		// report the restraint center as "value"
		x_reported = xr;
		v_reported = vr;
		// the "system force" with the extended Lagrangian is just the
		// the "total force" with the extended Lagrangian is just the
		// harmonic term acting on the extended coordinate
		// Note: this is the force for current timestep
		ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
		@@ -1117,7 +1126,7 @@ cvm::real colvar::update_forces_energy()
		f += fb;

		if (is_enabled(f_cv_Jacobian)) {
		// the instantaneous Jacobian force was not included in the reported system force;
		// the instantaneous Jacobian force was not included in the reported total force;
		// instead, it is subtracted from the applied force (silent Jacobian correction)
		if (is_enabled(f_cv_hide_Jacobian))
		f -= fj;
		@@ -1482,7 +1491,7 @@ std::istream & colvar::read_traj(std::istream &is)
		}
		}

		if (is_enabled(f_cv_output_system_force)) {
		if (is_enabled(f_cv_output_total_force)) {
		is >> ft;
		ft_reported = ft;
		}
		@@ -1567,8 +1576,8 @@ std::ostream & colvar::write_traj_label(std::ostream & os)
		<< cvm::wrap_string(this->name, this_cv_width-3);
		}

		if (is_enabled(f_cv_output_system_force)) {
		os << " fs_"
		if (is_enabled(f_cv_output_total_force)) {
		os << " ft_"
		<< cvm::wrap_string(this->name, this_cv_width-3);
		}

		@@ -1620,7 +1629,7 @@ std::ostream & colvar::write_traj(std::ostream &os)
		}


		if (is_enabled(f_cv_output_system_force)) {
		if (is_enabled(f_cv_output_total_force)) {
		os << " "
		<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
		<< ft_reported;

lib/colvars/colvar.h

+7 −7

Original line number	Diff line number	Diff line
		@@ -54,7 +54,7 @@ public:
		/// \brief Current velocity (previously set by calc() or by read_traj())
		colvarvalue const & velocity() const;

		/// \brief Current system force (previously obtained from calc() or
		/// \brief Current total force (previously obtained from calc() or
		/// read_traj()). Note: this is calculated using the atomic forces
		/// from the last simulation step.
		///
		@@ -62,7 +62,7 @@ public:
		/// acting on the collective variable is calculated summing those
		/// from all colvar components, the bias and walls forces are
		/// subtracted.
		colvarvalue const & system_force() const;
		colvarvalue const & total_force() const;

		/// \brief Typical fluctuation amplitude for this collective
		/// variable (e.g. local width of a free energy basin)
		@@ -160,7 +160,7 @@ protected:
		/// \brief Jacobian force, when Jacobian_force is enabled
		colvarvalue fj;

		/// Cached reported system force
		/// Cached reported total force
		colvarvalue ft_reported;

		public:
		@@ -176,7 +176,7 @@ public:
		colvarvalue f;

		/// \brief Total force, as derived from the atomic trajectory;
		/// should equal the total system force plus \link f \endlink
		/// should equal the system force plus \link f \endlink
		colvarvalue ft;


		@@ -253,7 +253,7 @@ public:
		int calc_cvc_values(int first, size_t num_cvcs);
		/// \brief Same as \link colvar::calc_cvc_values \endlink but for gradients
		int calc_cvc_gradients(int first, size_t num_cvcs);
		/// \brief Same as \link colvar::calc_cvc_values \endlink but for system forces
		/// \brief Same as \link colvar::calc_cvc_values \endlink but for total forces
		int calc_cvc_sys_forces(int first, size_t num_cvcs);
		/// \brief Same as \link colvar::calc_cvc_values \endlink but for Jacobian derivatives/forces
		int calc_cvc_Jacobians(int first, size_t num_cvcs);
		@@ -265,7 +265,7 @@ public:
		int collect_cvc_values();
		/// \brief Same as \link colvar::collect_cvc_values \endlink but for gradients
		int collect_cvc_gradients();
		/// \brief Same as \link colvar::collect_cvc_values \endlink but for system forces
		/// \brief Same as \link colvar::collect_cvc_values \endlink but for total forces
		int collect_cvc_sys_forces();
		/// \brief Same as \link colvar::collect_cvc_values \endlink but for Jacobian derivatives/forces
		int collect_cvc_Jacobians();
		@@ -555,7 +555,7 @@ inline colvarvalue const & colvar::velocity() const
		}


		inline colvarvalue const & colvar::system_force() const
		inline colvarvalue const & colvar::total_force() const
		{
		return ft_reported;
		}

lib/colvars/colvaratoms.cpp

+12 −12

Original line number	Diff line number	Diff line
		@@ -870,21 +870,21 @@ void cvm::atom_group::read_velocities()


		// TODO make this a calc function
		void cvm::atom_group::read_system_forces()
		void cvm::atom_group::read_total_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);
		ai->read_total_force();
		ai->total_force = rot.rotate(ai->total_force);
		}

		} else {

		for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
		ai->read_system_force();
		ai->read_total_force();
		}
		}
		}
		@@ -1070,15 +1070,15 @@ std::vector<cvm::rvector> cvm::atom_group::velocities() const
		return v;
		}

		std::vector<cvm::rvector> cvm::atom_group::system_forces() const
		std::vector<cvm::rvector> cvm::atom_group::total_forces() const
		{
		if (b_dummy) {
		cvm::error("Error: system forces are not available "
		cvm::error("Error: total 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 "
		cvm::error("Error: atomic total forces are not available "
		"from a scalable atom group.\n", INPUT_ERROR);
		}

		@@ -1086,27 +1086,27 @@ std::vector<cvm::rvector> cvm::atom_group::system_forces() const
		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;
		*fi = ai->total_force;
		}
		return f;
		}


		// TODO make this an accessor
		cvm::rvector cvm::atom_group::system_force() const
		cvm::rvector cvm::atom_group::total_force() const
		{
		if (b_dummy) {
		cvm::error("Error: total system forces are not available "
		cvm::error("Error: total total 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);
		return (cvm::proxy)->get_atom_group_total_force(index);
		}

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

lib/colvars/colvaratoms.h

+10 −10

Original line number	Diff line number	Diff line
		@@ -45,7 +45,7 @@ public:

		/// \brief System force at the previous step (copied from the
		/// program, can be modified if necessary)
		cvm::rvector system_force;
		cvm::rvector total_force;

		/// \brief Gradient of a scalar collective variable with respect
		/// to this atom
		@@ -86,7 +86,7 @@ public:
		inline void reset_data()
		{
		pos = cvm::atom_pos(0.0);
		vel = grad = system_force = cvm::rvector(0.0);
		vel = grad = total_force = cvm::rvector(0.0);
		}

		/// Get the latest value of the mass
		@@ -113,10 +113,10 @@ public:
		vel = (cvm::proxy)->get_atom_velocity(index);
		}

		/// Get the system force
		inline void read_system_force()
		/// Get the total force
		inline void read_total_force()
		{
		system_force = (cvm::proxy)->get_atom_system_force(index);
		total_force = (cvm::proxy)->get_atom_total_force(index);
		}

		/// \brief Apply a force to the atom
		@@ -336,10 +336,10 @@ public:
		/// rotation applied to the coordinates will be used
		void read_velocities();

		/// \brief Get the current system_forces; this must be called always
		/// \brief Get the current total_forces; this must be called always
		/// after read_positions(); if b_rotate is defined, the same
		/// rotation applied to the coordinates will be used
		void read_system_forces();
		void read_total_forces();

		/// Call reset_data() for each atom
		inline void reset_atoms_data()
		@@ -410,11 +410,11 @@ public:
		return dip;
		}

		/// \brief Return a copy of the system forces
		std::vector<cvm::rvector> system_forces() const;
		/// \brief Return a copy of the total forces
		std::vector<cvm::rvector> total_forces() const;

		/// \brief Return a copy of the aggregated total force on the group
		cvm::rvector system_force() const;
		cvm::rvector total_force() const;


		/// \brief Shorthand: save the specified gradient on each atom,

lib/colvars/colvarbias_abf.cpp

+24 −20

Original line number	Diff line number	Diff line
		@@ -7,7 +7,7 @@

		colvarbias_abf::colvarbias_abf(char const *key)
		: colvarbias(key),
		force(NULL),
		system_force(NULL),
		gradients(NULL),
		samples(NULL),
		z_gradients(NULL),
		@@ -81,8 +81,8 @@ int colvarbias_abf::init(std::string const &conf)
		}

		if (update_bias) {
		// Request calculation of system force (which also checks for availability)
		if(enable(f_cvb_get_system_force)) return cvm::get_error();
		// Request calculation of total force (which also checks for availability)
		if(enable(f_cvb_get_total_force)) return cvm::get_error();
		}
		if (apply_bias) {
		if(enable(f_cvb_apply_force)) return cvm::get_error();
		@@ -124,7 +124,7 @@ int colvarbias_abf::init(std::string const &conf)

		bin.assign(colvars.size(), 0);
		force_bin.assign(colvars.size(), 0);
		force = new cvm::real [colvars.size()];
		system_force = new cvm::real [colvars.size()];

		// Construct empty grids based on the colvars
		if (cvm::debug()) {
		@@ -208,9 +208,9 @@ colvarbias_abf::~colvarbias_abf()
		last_gradients = NULL;
		}

		if (force) {
		delete [] force;
		force = NULL;
		if (system_force) {
		delete [] system_force;
		system_force = NULL;
		}

		if (cvm::n_abf_biases > 0)
		@@ -244,10 +244,13 @@ int colvarbias_abf::update()
		if ( update_bias && samples->index_ok(force_bin) ) {
		// Only if requested and within bounds of the grid...

		for (size_t i=0; i<colvars.size(); i++) { // get forces(lagging by 1 timestep) from colvars
		force[i] = colvars[i]->system_force();
		for (size_t i = 0; i < colvars.size(); i++) {
		// get total forces (lagging by 1 timestep) from colvars
		// and subtract previous ABF force
		system_force[i] = colvars[i]->total_force().real_value
		- colvar_forces[i].real_value;
		}
		gradients->acc_force(force_bin, force);
		gradients->acc_force(force_bin, system_force);
		}
		if ( z_gradients && update_bias ) {
		for (size_t i = 0; i < colvars.size(); i++) {
		@@ -257,9 +260,10 @@ int colvarbias_abf::update()
		for (size_t i = 0; i < colvars.size(); i++) {
		// If we are outside the range of xi, the force has not been obtained above
		// the function is just an accessor, so cheap to call again anyway
		force[i] = colvars[i]->system_force();
		system_force[i] = colvars[i]->total_force().real_value
		- colvar_forces[i].real_value;
		}
		z_gradients->acc_force(z_bin, force);
		z_gradients->acc_force(z_bin, system_force);
		}
		}
		}

Admin message