Merge pull request #309 from giacomofiorin/colvars-2016-12-22 (41a6a307) · Commits · 郑智淋 / lammps

doc/src/PDF/colvars-refman-lammps.pdf

+1.89 KiB (553 KiB)

File changed.

No diff preview for this file type.

doc/src/fix_colvars.txt

+5 −12

Original line number	Diff line number	Diff line
		@@ -31,21 +31,19 @@ fix abf all colvars colvars.inp tstat 1 :pre

		[Description:]

		This fix interfaces LAMMPS to a "collective variables" or "colvars"
		module library which allows to calculate potentials of mean force
		This fix interfaces LAMMPS to the collective variables "Colvars"
		library, which allows to calculate potentials of mean force
		(PMFs) for any set of colvars, using different sampling methods:
		currently implemented are the Adaptive Biasing Force (ABF) method,
		metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
		(US) via a flexible harmonic restraint bias. The colvars library is
		hosted at "http://colvars.github.io/"_http://colvars.github.io/
		(US) via a flexible harmonic restraint bias.

		This documentation describes only the fix colvars command itself and
		LAMMPS specific parts of the code. The full documentation of the
		colvars library is available as "this supplementary PDF document"_PDF/colvars-refman-lammps.pdf

		A detailed discussion of the implementation of the portable collective
		variable library is in "(Fiorin)"_#Fiorin. Additional information can
		be found in "(Henin)"_#Henin.
		The Colvars library is developed at "https://github.com/colvars/colvars"_https://github.com/colvars/colvars
		A detailed discussion of its implementation is in "(Fiorin)"_#Fiorin.

		There are some example scripts for using this package with LAMMPS in the
		examples/USER/colvars directory.
		@@ -129,8 +127,3 @@ and tstat = NULL.

		:link(Fiorin)
		[(Fiorin)] Fiorin , Klein, Henin, Mol. Phys., DOI:10.1080/00268976.2013.813594

		:link(Henin)
		[(Henin)] Henin, Fiorin, Chipot, Klein, J. Chem. Theory Comput., 6,
		35-47 (2010)

lib/colvars/Makefile.colvars

0 → 100644

+119 −0

Original line number	Diff line number	Diff line
		# library build -- makefile -- for colvars module

		# which file will be copied to Makefile.lammps

		EXTRAMAKE = Makefile.lammps.empty

		# ------ SETTINGS ------

		CXX = g++
		CXXFLAGS = -O2 -g -Wall -fPIC -funroll-loops # -DCOLVARS_DEBUG
		ARCHIVE = ar
		ARCHFLAG = -rscv
		SHELL = /bin/sh

		# ------ DEFINITIONS ------

		SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
		colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
		colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
		colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
		colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
		colvarscript.cpp colvartypes.cpp colvarvalue.cpp

		LIB = libcolvars.a
		OBJ = $(SRC:.cpp=.o)
		EXE = #colvars_standalone

		# ------ MAKE PROCEDURE ------

		default: $(LIB) $(EXE) Makefile.lammps

		Makefile.lammps:
		@cp $(EXTRAMAKE) Makefile.lammps

		$(LIB): $(OBJ)
		$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)

		colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
		$(CXX) -o $@ $(CXXFLAGS) $^

		# ------ MAKE FLAGS ------

		.SUFFIXES:
		.SUFFIXES: .cpp .o

		.PHONY: default clean

		# ------ COMPILE RULES ------

		.cpp.o:
		$(CXX) $(CXXFLAGS) -c $<

		# ------ DEPENDENCIES ------
		#
		colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h colvarparse.h colvardeps.h colvaratoms.h
		colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
		colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
		colvarbias_abf.h colvarbias.h colvargrid.h
		colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \
		colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.h colvarparse.h \
		colvardeps.h colvarbias_restraint.h colvarbias.h
		colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
		colvarbias_histogram.o: colvarbias_histogram.cpp colvarmodule.h \
		colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \
		colvardeps.h colvarbias_histogram.h colvarbias.h colvargrid.h
		colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
		colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
		colvarbias_meta.h colvarbias.h colvargrid.h
		colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \
		colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \
		colvarbias.h colvar.h colvarparse.h colvardeps.h
		colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
		colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
		colvarcomp.h colvaratoms.h
		colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
		colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
		colvaratoms.h colvar.h colvarcomp.h
		colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
		colvaratoms.h
		colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
		colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
		colvar.h colvarcomp.h colvaratoms.h
		colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
		colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
		colvarcomp.h colvaratoms.h
		colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
		colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
		colvar.h colvarcomp.h colvaratoms.h
		colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
		colvaratoms.h colvarscript.h colvarbias.h
		colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
		colvarproxy.h colvarvalue.h colvarparse.h
		colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
		colvaratoms.h colvargrid.h
		colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
		colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
		colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
		colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
		colvarscript.h
		colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h colvarparse.h
		colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
		colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
		colvarparse.h colvardeps.h
		colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h colvarparse.h
		colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
		colvarvalue.h

		# ------ CLEAN ------

		clean:
		-rm .o ~ $(LIB)

lib/colvars/colvar.cpp

+114 −56

Original line number	Diff line number	Diff line
		// -- c++ --

		// This file is part of the Collective Variables module (Colvars).
		// The original version of Colvars and its updates are located at:
		// https://github.com/colvars/colvars
		// Please update all Colvars source files before making any changes.
		// If you wish to distribute your changes, please submit them to the
		// Colvars repository at GitHub.

		#include "colvarmodule.h"
		#include "colvarvalue.h"
		#include "colvarparse.h"
		@@ -163,30 +170,30 @@ colvar::colvar(std::string const &conf)
		}
		feature_states[f_cv_homogeneous]->enabled = homogeneous;
		}
		// Colvar is deemed periodic iff:

		// Colvar is deemed periodic if:
		// - it is homogeneous
		// - all cvcs are periodic
		// - all cvcs have the same period

		b_periodic = cvcs[0]->b_periodic && is_enabled(f_cv_homogeneous);
		if (cvcs[0]->b_periodic) { // TODO make this a CVC feature
		bool b_periodic = true;
		period = cvcs[0]->period;
		for (i = 1; i < cvcs.size(); i++) {
		if (!cvcs[i]->b_periodic \|\| cvcs[i]->period != period) {
		b_periodic = false;
		period = 0.0;
		cvm::log("Warning: although one component is periodic, this colvar will "
		"not be treated as periodic, either because the exponent is not "
		"1, or because components of different periodicity are defined. "
		"Make sure that you know what you are doing!");
		}
		}
		feature_states[f_cv_periodic]->enabled = b_periodic;
		}

		// check that cvcs are compatible

		for (i = 0; i < cvcs.size(); i++) {
		if ((cvcs[i])->b_periodic && !b_periodic) {
		cvm::log("Warning: although this component is periodic, the colvar will "
		"not be treated as periodic, either because the exponent is not "
		"1, or because multiple components are present. Make sure that "
		"you know what you are doing!");
		}

		// components may have different types only for scripted functions
		if (!is_enabled(f_cv_scripted) && (colvarvalue::check_types(cvcs[i]->value(),
		@@ -207,16 +214,15 @@ colvar::colvar(std::string const &conf)
		// at this point, the colvar's type is defined
		f.type(value());
		f_accumulated.type(value());
		fb.type(value());

		reset_bias_force();

		get_keyval(conf, "width", width, 1.0);
		if (width <= 0.0) {
		cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
		return;
		}

		// NOTE: not porting wall stuff to new deps, as this will change to a separate bias
		// the grid functions will wait a little as well

		lower_boundary.type(value());
		lower_wall.type(value());

		@@ -308,6 +314,9 @@ colvar::colvar(std::string const &conf)
		enable(f_cv_extended_Lagrangian);
		provide(f_cv_Langevin);

		// The extended mass will apply forces
		enable(f_cv_gradient);

		xr.type(value());
		vr.type(value());
		fr.type(value());
		@@ -400,6 +409,9 @@ colvar::colvar(std::string const &conf)
		f_old.type(value());
		f_old.reset();

		x_restart.type(value());
		after_restart = false;

		if (cvm::b_analysis)
		parse_analysis(conf);

		@@ -761,9 +773,13 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
		return error_code;
		}

		if (cvm::step_relative() > 0) {
		// Total force depends on Jacobian derivative from previous timestep
		error_code \|= calc_cvc_total_force(first_cvc, num_cvcs);
		}
		// atom coordinates are updated by the next line
		error_code \|= calc_cvc_values(first_cvc, num_cvcs);
		error_code \|= calc_cvc_gradients(first_cvc, num_cvcs);
		error_code \|= calc_cvc_total_force(first_cvc, num_cvcs);
		error_code \|= calc_cvc_Jacobians(first_cvc, num_cvcs);

		if (cvm::debug())
		@@ -780,9 +796,12 @@ int colvar::collect_cvc_data()

		int error_code = COLVARS_OK;

		if (cvm::step_relative() > 0) {
		// Total force depends on Jacobian derivative from previous timestep
		error_code \|= collect_cvc_total_forces();
		}
		error_code \|= collect_cvc_values();
		error_code \|= collect_cvc_gradients();
		error_code \|= collect_cvc_total_forces();
		error_code \|= collect_cvc_Jacobians();
		error_code \|= calc_colvar_properties();

		@@ -872,6 +891,22 @@ int colvar::collect_cvc_values()
		cvm::log("Colvar \""+this->name+"\" has value "+
		cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");

		if (after_restart) {
		after_restart = false;
		if (cvm::proxy->simulation_running()) {
		cvm::real const jump2 = dist2(x, x_restart) / (width*width);
		if (jump2 > 0.25) {
		cvm::error("Error: the calculated value of colvar \""+name+
		"\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
		"last read from the state file:\n"+cvm::to_str(x_restart)+
		"\nPossible causes are changes in configuration, "
		"wrong state file, or how PBC wrapping is handled.\n",
		INPUT_ERROR);
		return INPUT_ERROR;
		}
		}
		}

		return COLVARS_OK;
		}

		@@ -979,13 +1014,8 @@ int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
		if (cvm::debug())
		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 total force calculation
		// even with extended Lagrangian

		if (cvm::step_relative() > 0) {
		cvm::increase_depth();
		// 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 +1024,7 @@ int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
		(cvcs[i])->calc_force_invgrads();
		}
		cvm::decrease_depth();
		}


		if (cvm::debug())
		cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
		@@ -1013,6 +1043,11 @@ int colvar::collect_cvc_total_forces()
		// 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;
		if (cvm::debug())
		cvm::log("Colvar component no. "+cvm::to_str(i+1)+
		" within colvar \""+this->name+"\" has total force "+
		cvm::to_str((cvcs[i])->total_force(),
		cvm::cv_width, cvm::cv_prec)+".\n");
		// linear combination is assumed
		ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
		}
		@@ -1056,6 +1091,11 @@ int colvar::collect_cvc_Jacobians()
		fj.reset();
		for (size_t i = 0; i < cvcs.size(); i++) {
		if (!cvcs[i]->is_enabled()) continue;
		if (cvm::debug())
		cvm::log("Colvar component no. "+cvm::to_str(i+1)+
		" within colvar \""+this->name+"\" has Jacobian derivative"+
		cvm::to_str((cvcs[i])->Jacobian_derivative(),
		cvm::cv_width, cvm::cv_prec)+".\n");
		// linear combination is assumed
		fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
		}
		@@ -1085,7 +1125,7 @@ int colvar::calc_colvar_properties()
		// TODO: put it in the restart information
		if (cvm::step_relative() == 0) {
		xr = x;
		vr = 0.0; // (already 0; added for clarity)
		vr.reset(); // (already 0; added for clarity)
		}

		// report the restraint center as "value"
		@@ -1128,28 +1168,28 @@ cvm::real colvar::update_forces_energy()
		if (is_enabled(f_cv_Jacobian)) {
		// the instantaneous Jacobian force was not included in the reported total force;
		// instead, it is subtracted from the applied force (silent Jacobian correction)
		// This requires the Jacobian term for the current timestep
		if (is_enabled(f_cv_hide_Jacobian))
		f -= fj;
		}

		if (is_enabled(f_cv_lower_wall) \|\| is_enabled(f_cv_upper_wall)) {

		// Wall force
		colvarvalue fw(x);
		fw.reset();

		if (is_enabled(f_cv_lower_wall) \|\| is_enabled(f_cv_upper_wall)) {

		if (cvm::debug())
		cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n");

		// 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_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) {
		(this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {

		cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall);
		cvm::real const grad = this->dist2_lgrad(x, lower_wall);
		if (grad < 0.0) {
		fw = -0.5 * lower_wall_k * grad;
		f += fw;
		if (cvm::debug())
		cvm::log("Applying a lower wall force("+
		cvm::to_str(fw)+") to \""+this->name+"\".\n");
		@@ -1157,10 +1197,9 @@ cvm::real colvar::update_forces_energy()

		} else {

		cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall);
		cvm::real const grad = this->dist2_lgrad(x, upper_wall);
		if (grad > 0.0) {
		fw = -0.5 * upper_wall_k * grad;
		f += fw;
		if (cvm::debug())
		cvm::log("Applying an upper wall force("+
		cvm::to_str(fw)+") to \""+this->name+"\".\n");
		@@ -1168,17 +1207,26 @@ cvm::real colvar::update_forces_energy()
		}
		}

		// At this point f is the force f from external biases that will be applied to the
		// extended variable if there is one

		if (is_enabled(f_cv_extended_Lagrangian)) {

		if (cvm::debug()) {
		cvm::log("Updating extended-Lagrangian degrees of freedom.\n");
		}

		cvm::real dt = cvm::dt();
		cvm::real f_ext;
		colvarvalue f_ext(fr.type());
		f_ext.reset();

		// the total force is applied to the fictitious mass, while the
		// atoms only feel the harmonic force
		// atoms only feel the harmonic force + wall 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 (including wall forces)
		// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
		// f: - initially, external biasing force
		// - after this code block, colvar force to be applied to atomic coordinates
		// ie. spring force + wall 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);
		@@ -1200,15 +1248,24 @@ cvm::real colvar::update_forces_energy()
		potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
		// leap to v_(i+1/2)
		if (is_enabled(f_cv_Langevin)) {
		vr -= dt * ext_gamma * vr.real_value;
		vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
		vr -= dt * ext_gamma * vr;
		colvarvalue rnd(x);
		rnd.set_random();
		vr += dt * ext_sigma * rnd / ext_mass;
		}
		vr += (0.5 * dt) * f_ext / ext_mass;
		xr += dt * vr;
		xr.apply_constraints();
		if (this->b_periodic) this->wrap(xr);
		if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
		}

		// TODO remove the wall force
		f += fw;
		// Now adding the force on the actual colvar (for those biases who
		// bypass the extended Lagrangian mass)
		f += fb_actual;

		// Store force to be applied, possibly summed over several timesteps
		f_accumulated += f;

		if (is_enabled(f_cv_fdiff_velocity)) {
		@@ -1425,14 +1482,15 @@ std::istream & colvar::read_restart(std::istream &is)
		}
		}

		if ( !(get_keyval(conf, "x", x,
		colvarvalue(x.type()), colvarparse::parse_silent)) ) {
		if ( !(get_keyval(conf, "x", x, x, colvarparse::parse_silent)) ) {
		cvm::log("Error: restart file does not contain "
		"the value of the colvar \""+
		name+"\" .\n");
		} else {
		cvm::log("Restarting collective variable \""+name+"\" from value: "+
		cvm::to_str(x)+"\n");
		x_restart = x;
		after_restart = true;
		}

		if (is_enabled(f_cv_extended_Lagrangian)) {

lib/colvars/colvar.h

+31 −5

Original line number	Diff line number	Diff line
		// -- c++ --

		// This file is part of the Collective Variables module (Colvars).
		// The original version of Colvars and its updates are located at:
		// https://github.com/colvars/colvars
		// Please update all Colvars source files before making any changes.
		// If you wish to distribute your changes, please submit them to the
		// Colvars repository at GitHub.

		#ifndef COLVAR_H
		#define COLVAR_H

		@@ -170,6 +177,9 @@ public:
		/// the biases are updated
		colvarvalue fb;

		/// \brief Bias force to the actual value (only useful with extended Lagrangian)
		colvarvalue fb_actual;

		/// \brief Total \em applied force; fr (if extended_lagrangian
		/// is defined), fb (if biases are applied) and the walls' forces
		/// (if defined) contribute to it
		@@ -183,13 +193,9 @@ public:
		colvarvalue ft;


		/// Period, if it is a constant
		/// Period, if this variable is periodic
		cvm::real period;

		/// \brief Same as above, but also takes into account components
		/// with a variable period, such as distanceZ
		bool b_periodic;


		/// \brief Expand the boundaries of multiples of width, to keep the
		/// value always within range
		@@ -290,6 +296,9 @@ public:
		/// Add to the total force from biases
		void add_bias_force(colvarvalue const &force);

		/// Apply a force to the actual value (only meaningful with extended Lagrangian)
		void add_bias_force_actual_value(colvarvalue const &force);

		/// \brief Collect all forces on this colvar, integrate internal
		/// equations of motion of internal degrees of freedom; see also
		/// colvar::communicate_forces()
		@@ -386,6 +395,12 @@ protected:
		/// Previous value (to calculate velocities during analysis)
		colvarvalue x_old;

		/// Value read from the most recent state file (if any)
		colvarvalue x_restart;

		/// True if a state file was just read
		bool after_restart;

		/// Time series of values and velocities used in correlation
		/// functions
		std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
		@@ -577,9 +592,20 @@ inline void colvar::add_bias_force(colvarvalue const &force)
		}


		inline void colvar::add_bias_force_actual_value(colvarvalue const &force)
		{
		if (cvm::debug()) {
		cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
		}
		fb_actual += force;
		}


		inline void colvar::reset_bias_force() {
		fb.type(value());
		fb.reset();
		fb_actual.type(value());
		fb_actual.reset();
		}

		#endif

Admin message