Merge branch 'colvars-2016-10-24' of https://github.com/giacomofiorin/lammps into colvars-update (ba6d1528) · Commits · 郑智淋 / lammps

lib/colvars/colvar.cpp

+13 −7

Original line number	Diff line number	Diff line
		@@ -1144,9 +1144,9 @@ cvm::real colvar::update_forces_energy()
		// For a periodic colvar, both walls may be applicable at the same time
		// in which case we pick the closer one
		if ( (!is_enabled(f_cv_upper_wall)) \|\|
		(this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {
		(this->dist2(x_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) {

		cvm::real const grad = this->dist2_lgrad(x, lower_wall);
		cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall);
		if (grad < 0.0) {
		fw = -0.5 * lower_wall_k * grad;
		f += fw;
		@@ -1157,7 +1157,7 @@ cvm::real colvar::update_forces_energy()

		} else {

		cvm::real const grad = this->dist2_lgrad(x, upper_wall);
		cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall);
		if (grad > 0.0) {
		fw = -0.5 * upper_wall_k * grad;
		f += fw;
		@@ -1177,15 +1177,21 @@ cvm::real colvar::update_forces_energy()
		// atoms only feel the harmonic force
		// fr: bias force on extended variable (without harmonic spring), for output in trajectory
		// f_ext: total force on extended variable (including harmonic spring)
		// f: - initially, external biasing force
		// f: - initially, external biasing force (including wall forces)
		// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
		fr = f;
		f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
		f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);

		if (is_enabled(f_cv_subtract_applied_force)) {
		// Report a "system" force without the biases on this colvar
		// that is, just the spring force
		ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
		} else {
		// The total force acting on the extended variable is f_ext
		// This will be used in the next timestep
		ft_reported = f_ext;
		}

		// leapfrog: starting from x_i, f_i, v_(i-1/2)
		vr += (0.5 * dt) * f_ext / ext_mass;

lib/colvars/colvaratoms.cpp

+39 −6

Original line number	Diff line number	Diff line
		@@ -367,6 +367,7 @@ int cvm::atom_group::parse(std::string const &conf)
		cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR);
		}

		// NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
		cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
		}
		}
		@@ -403,11 +404,21 @@ int cvm::atom_group::parse(std::string const &conf)
		}
		}

		// We need to know the fitting options to decide whether the group is scalable
		parse_error \|= parse_fitting_options(group_conf);

		if (is_available(f_ag_scalable_com) && !b_rotate) {
		enable(f_ag_scalable_com);
		enable(f_ag_scalable);
		}

		if (is_enabled(f_ag_scalable) && !b_dummy) {
		cvm::log("Enabling scalable calculation for group \""+this->key+"\".\n");
		index = (cvm::proxy)->init_atom_group(atoms_ids);
		}

		parse_error \|= parse_fitting_options(group_conf);
		bool b_print_atom_ids = false;
		get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false, colvarparse::parse_silent);

		// TODO move this to colvarparse object
		check_keywords(group_conf, key.c_str());
		@@ -427,6 +438,10 @@ int cvm::atom_group::parse(std::string const &conf)
		cvm::to_str(total_mass)+", total charge = "+
		cvm::to_str(total_charge)+".\n");

		if (b_print_atom_ids) {
		cvm::log("Internal definition of the atom group:\n");
		}

		cvm::decrease_depth();

		return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
		@@ -583,6 +598,21 @@ int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
		}


		std::string const cvm::atom_group::print_atom_ids() const
		{
		size_t line_count = 0;
		std::ostringstream os("");
		for (size_t i = 0; i < atoms_ids.size(); i++) {
		os << " " << std::setw(9) << atoms_ids[i];
		if (++line_count == 7) {
		os << "\n";
		line_count = 0;
		}
		}
		return os.str();
		}


		int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
		{
		bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
		@@ -1118,8 +1148,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
		log("Communicating a colvar force from atom group to the MD engine.\n");
		}

		if (b_dummy)
		return;
		if (b_dummy) return;

		if (noforce) {
		cvm::error("Error: sending a force to a group that has "
		@@ -1161,17 +1190,21 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)

		void cvm::atom_group::apply_force(cvm::rvector const &force)
		{
		if (b_dummy)
		return;
		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 "
		"\"disableForces\" defined.\n");
		"\"enableForces\" set to off.\n");
		return;
		}

		if (is_enabled(f_ag_scalable)) {
		(cvm::proxy)->apply_atom_group_force(index, force);
		return;
		}

		if (b_rotate) {

lib/colvars/colvaratoms.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -253,6 +253,8 @@ public:
		return atoms.size();
		}

		std::string const print_atom_ids() const;

		/// \brief If this option is on, this group merely acts as a wrapper
		/// for a fixed position; any calls to atoms within or to
		/// functions that return disaggregated data will fail

lib/colvars/colvarbias_abf.cpp

+28 −29

Original line number	Diff line number	Diff line
		@@ -80,6 +80,7 @@ int colvarbias_abf::init(std::string const &conf)

		if (update_bias) {
		// Request calculation of total force (which also checks for availability)
		// TODO - change this to a dependency - needs ABF-specific features
		if(enable(f_cvb_get_total_force)) return cvm::get_error();
		}

		@@ -133,6 +134,10 @@ int colvarbias_abf::init(std::string const &conf)

		// Data for eABF z-based estimator
		if (b_extended) {
		// CZAR output files for stratified eABF
		get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false,
		colvarparse::parse_silent);

		z_bin.assign(colvars.size(), 0);
		z_samples = new colvar_grid_count(colvars);
		z_samples->request_actual_value();
		@@ -241,7 +246,7 @@ int colvarbias_abf::update()

		for (size_t i = 0; i < colvars.size(); i++) {
		// get total forces (lagging by 1 timestep) from colvars
		// and subtract previous ABF force
		// and subtract previous ABF force if necessary
		update_system_force(i);
		}
		gradients->acc_force(force_bin, system_force);
		@@ -457,28 +462,30 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app

		if (z_gradients) {
		// Write eABF-related quantities

		std::string z_samples_out_name = prefix + ".zcount";
		std::string z_gradients_out_name = prefix + ".zgrad";
		std::string czar_gradients_out_name = prefix + ".czar";
		cvm::ofstream z_samples_os;
		cvm::ofstream z_gradients_os;
		cvm::ofstream czar_gradients_os;

		if (!append) cvm::backup_file(z_samples_out_name.c_str());
		z_samples_os.open(z_samples_out_name.c_str(), mode);
		if (!z_samples_os.is_open()) {
		cvm::error("Error opening ABF z sample file " + z_samples_out_name + " for writing");
		cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
		}
		z_samples->write_multicol(z_samples_os);
		z_samples_os.close();

		if (b_czar_window_file) {
		std::string z_gradients_out_name = prefix + ".zgrad";
		cvm::ofstream z_gradients_os;

		if (!append) cvm::backup_file(z_gradients_out_name.c_str());
		z_gradients_os.open(z_gradients_out_name.c_str(), mode);
		if (!z_gradients_os.is_open()) {
		cvm::error("Error opening ABF z gradient file " + z_gradients_out_name + " for writing");
		cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
		}
		z_gradients->write_multicol(z_gradients_os);
		z_gradients_os.close();
		}

		// Calculate CZAR estimator of gradients
		for (std::vector<int> ix = czar_gradients->new_index();
		@@ -490,6 +497,9 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
		}
		}

		std::string czar_gradients_out_name = prefix + ".czar.grad";
		cvm::ofstream czar_gradients_os;

		if (!append) cvm::backup_file(czar_gradients_out_name.c_str());
		czar_gradients_os.open(czar_gradients_out_name.c_str(), mode);
		if (!czar_gradients_os.is_open()) {
		@@ -499,17 +509,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
		czar_gradients_os.close();

		if (colvars.size() == 1) {
		std::string z_pmf_out_name = prefix + ".zpmf";
		if (!append) cvm::backup_file(z_pmf_out_name.c_str());
		cvm::ofstream z_pmf_os;
		// Do numerical integration and output a PMF
		z_pmf_os.open(z_pmf_out_name.c_str(), mode);
		if (!z_pmf_os.is_open()) cvm::error("Error opening z pmf file " + z_pmf_out_name + " for writing");
		z_gradients->write_1D_integral(z_pmf_os);
		z_pmf_os << std::endl;
		z_pmf_os.close();

		std::string czar_pmf_out_name = prefix + ".czarpmf";
		std::string czar_pmf_out_name = prefix + ".czar.pmf";
		if (!append) cvm::backup_file(czar_pmf_out_name.c_str());
		cvm::ofstream czar_pmf_os;
		// Do numerical integration and output a PMF
		@@ -520,8 +520,6 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
		czar_pmf_os.close();
		}
		}


		return;
		}

		@@ -559,7 +557,7 @@ void colvarbias_abf::read_gradients_samples()

		std::ifstream is;

		cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name);
		cvm::log("Reading sample count from " + samples_in_name + " and gradient from " + gradients_in_name);
		is.open(samples_in_name.c_str());
		if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");
		samples->read_multicol(is, true);
		@@ -572,17 +570,18 @@ void colvarbias_abf::read_gradients_samples()
		is.close();

		if (z_gradients) {
		cvm::log("Reading z sample count from " + z_samples_in_name + " and z gradients from " + z_gradients_in_name);
		// Read eABF z-averaged data for CZAR
		cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name);

		is.clear();
		is.open(z_samples_in_name.c_str());
		if (!is.is_open()) cvm::error("Error opening ABF z samples file " + z_samples_in_name + " for reading");
		if (!is.is_open()) cvm::error("Error opening eABF z-histogram file " + z_samples_in_name + " for reading");
		z_samples->read_multicol(is, true);
		is.close();
		is.clear();

		is.open(z_gradients_in_name.c_str());
		if (!is.is_open()) cvm::error("Error opening ABF z gradient file " + z_gradients_in_name + " for reading");
		if (!is.is_open()) cvm::error("Error opening eABF z-gradient file " + z_gradients_in_name + " for reading");
		z_gradients->read_multicol(is, true);
		is.close();
		}

lib/colvars/colvarbias_abf.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -40,6 +40,8 @@ private:
		int output_freq;
		/// Write combined files with a history of all output data?
		bool b_history_files;
		/// Write CZAR output file for stratified eABF (.zgrad)
		bool b_czar_window_file;
		size_t history_freq;

		/// Cap applied biasing force?

Admin message