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

This commit is contained in:
sjplimp 2013-06-27 22:50:07 +00:00
parent 6dc2447695
commit 738d900762
6 changed files with 259 additions and 181 deletions

View File

@ -35,7 +35,7 @@ if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*colvars[^ \t]* //g' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-I..\/..\/lib\/colvars |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L..\/..\/lib\/colvars |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L..\/..\/lib\/colvars$(LIBOBJDIR) |' ../Makefile.package
sed -i -e 's|^PKG_LIB =[ \t]*|&-lcolvars |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(colvars_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(colvars_SYSLIB) |' ../Makefile.package

View File

@ -7,12 +7,20 @@ Sampling and Restraints. This code consists of two parts:
- a portable collective variable module library written
and maintained by Giacomo Fiorin (ICMS, Temple University,
Philadelphia, PA, USA) and Jerome Henin (LISM, CNRS,
Marseille, France).
Philadelphia, PA, USA) and Jerome Henin (IBPC, CNRS,
Paris, France).
A copy of this code is located in the directory lib/colvars
and needs to be compiled first. More info about this code
can be found at: http://colvars.github.io/
and in the publications:
Using collective variables to drive molecular dynamics
simulations,
Giacomo Fiorin , Michael L. Klein & Jérôme Hénin (2013):
Molecular Physics DOI:10.1080/00268976.2013.813594
This code is located in the directory lib/colvars and
needs to be compiled first. More info about this code
can be found in the publication:
Exploring Multidimensional Free Energy Landscapes Using
Time-Dependent Biases on Collective Variables,
J. Hénin, G. Fiorin, C. Chipot, and M. L. Klein,
@ -20,23 +28,17 @@ Sampling and Restraints. This code consists of two parts:
- the colvars fix and a thin interface layer, which exchanges
information between LAMMPS and the collective variable module.
This interface was written and is maintained by
Axel Kohlmeyer (akohlmey@gmail.com)
See the doc page of fix colvars for more details
There are example scripts for using this package in
examples/USER/colvars
This is a very new interface that does not yet support all
features in the module and will see future optimizations
and improvements. The colvars module library is also available
in NAMD has been thoroughly used and tested there. Bugs and
problems are likely due to the interface layers code.
Thus the current version of this package should be considered
beta quality.
The person who created this package is Axel Kohlmeyer at Temple U
(akohlmey at gmail.com). Contact him directly if you have questions.
---------------------------------
version: 0.1 / 2012-03-03
Version: 2013-06-27

View File

@ -56,13 +56,11 @@ static void my_backup_file(const char *filename, const char *extension)
////////////////////////////////////////////////////////////////////////
colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
const char *conf_file,
const char *inp_name,
const char *out_name,
const int seed,
const double temp,
const int *typemap)
: _lmp(lmp), _typemap(typemap)
const double temp)
: _lmp(lmp)
{
if (cvm::debug())
log("Initializing the colvars proxy object.\n");
@ -87,15 +85,20 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
restart_prefix_str = std::string("rest");
// try to extract a restart prefix from a potential restart command.
if (_lmp->output->restart_every_single > 0) {
restart_prefix_str = std::string(_lmp->output->restart1);
} else if (_lmp->output->restart_every_double > 0) {
restart_prefix_str = std::string(_lmp->output->restart2a);
LAMMPS_NS::Output *outp = _lmp->output;
if ((outp->restart_every_single > 0) && (outp->restart1 != 0)) {
restart_prefix_str = std::string(outp->restart1);
} else if ((outp->restart_every_double > 0) && (outp->restart2a != 0)) {
restart_prefix_str = std::string(outp->restart2a);
}
// trim off unwanted stuff from the restart prefix
if (restart_prefix_str.rfind(".*") != std::string::npos)
restart_prefix_str.erase(restart_prefix_str.rfind(".*"),2);
}
void colvarproxy_lammps::init(const char *conf_file)
{
// create the colvarmodule instance
colvars = new colvarmodule(conf_file,this);
@ -112,14 +115,11 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
log(cvm::line_marker);
log("Info: done initializing the colvars proxy object.\n");
}
// this is only valid through the constructor.
_typemap = NULL;
}
colvarproxy_lammps::~colvarproxy_lammps()
{
delete _random;
if (colvars != NULL) {
colvars->write_output_files();
delete colvars;
@ -127,6 +127,12 @@ colvarproxy_lammps::~colvarproxy_lammps()
}
}
// re-initialize data where needed
void colvarproxy_lammps::setup()
{
colvars->setup();
}
// trigger colvars computation
double colvarproxy_lammps::compute()
{
@ -149,6 +155,11 @@ double colvarproxy_lammps::compute()
"Updating internal data.\n");
}
// zero the forces on the atoms, so that they can be accumulated by the colvars
for (size_t i = 0; i < applied_forces.size(); i++) {
applied_forces[i].x = applied_forces[i].y = applied_forces[i].z = 0.0;
}
// call the collective variable module
colvars->calc();
@ -164,6 +175,26 @@ double colvarproxy_lammps::compute()
return bias_energy;
}
void colvarproxy_lammps::serialize_status(std::string &rst)
{
std::ostringstream os;
colvars->write_restart(os);
rst = os.str();
}
// set status from string
bool colvarproxy_lammps::deserialize_status(std::string &rst)
{
std::istringstream is;
is.str(rst);
if (!colvars->read_restart(is)) {
return false;
} else {
return true;
}
}
cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
@ -308,7 +339,7 @@ void colvarproxy_lammps::backup_file(char const *filename)
int colvarproxy_lammps::init_lammps_atom(const int &aid, cvm::atom *atom)
{
atom->id = aid;
atom->mass = _lmp->atom->mass[_typemap[aid]];
atom->mass = 0.0;
for (size_t i = 0; i < colvars_atoms.size(); i++) {
if (colvars_atoms[i] == aid) {
@ -323,7 +354,7 @@ int colvarproxy_lammps::init_lammps_atom(const int &aid, cvm::atom *atom)
colvars_atoms.push_back(aid);
struct commdata c;
c.tag = aid;
c.type = _typemap[aid];
c.type = 0;
c.x = c.y = c.z = 0.0;
positions.push_back(c);
total_forces.push_back(c);
@ -385,18 +416,20 @@ cvm::atom::atom(cvm::atom const &a)
cvm::atom::~atom()
{
colvarproxy_lammps *cp = (colvarproxy_lammps *) cvm::proxy;
if (cp->colvars_atoms_ncopies[this->index] > 0)
cp->colvars_atoms_ncopies[this->index] -= 1;
if (this->index >= 0) {
colvarproxy_lammps *cp = (colvarproxy_lammps *) cvm::proxy;
if (cp->colvars_atoms_ncopies[this->index] > 0)
cp->colvars_atoms_ncopies[this->index] -= 1;
}
}
void cvm::atom::read_position()
{
colvarproxy_lammps const * const cp = (colvarproxy_lammps *) cvm::proxy;
this->pos.x = cp->positions[this->index].x;
this->pos.y = cp->positions[this->index].y;
this->pos.z = cp->positions[this->index].z;
this->mass = cp->positions[this->index].m;
}
void cvm::atom::read_velocity()
@ -418,7 +451,7 @@ void cvm::atom::read_system_force()
void cvm::atom::apply_force(cvm::rvector const &new_force)
{
colvarproxy_lammps *cp = (colvarproxy_lammps *) cvm::proxy;
cp->applied_forces[this->index].x = new_force.x;
cp->applied_forces[this->index].y = new_force.y;
cp->applied_forces[this->index].z = new_force.z;
cp->applied_forces[this->index].x += new_force.x;
cp->applied_forces[this->index].y += new_force.y;
cp->applied_forces[this->index].z += new_force.z;
}

View File

@ -17,7 +17,7 @@
/* struct for packed data communication of coordinates and forces. */
struct commdata {
int tag,type;
double x,y,z;
double x,y,z,m;
};
inline std::ostream & operator<< (std::ostream &out, const commdata &cd)
@ -38,9 +38,6 @@ class colvarproxy_lammps : public colvarproxy {
class LAMMPS_NS::LAMMPS *_lmp;
class LAMMPS_NS::RanPark *_random;
// pointers to LAMMPS provided storage
const int *_typemap;
// state of LAMMPS properties
double t_target;
double bias_energy;
@ -63,9 +60,11 @@ class colvarproxy_lammps : public colvarproxy {
public:
friend class cvm::atom;
colvarproxy_lammps (LAMMPS_NS::LAMMPS *lmp, const char *, const char *,
const char *, const int, const double, const int *);
colvarproxy_lammps (LAMMPS_NS::LAMMPS *lmp, const char *,
const char *, const int, const double);
virtual ~colvarproxy_lammps();
void init(const char*);
void setup();
// disable default and copy constructor
private:
@ -88,6 +87,12 @@ class colvarproxy_lammps : public colvarproxy {
// perform colvars computation. returns biasing energy
double compute();
// dump status to string
void serialize_status(std::string &);
// set status from string
bool deserialize_status(std::string &);
// implementation of pure methods from base class
public:

View File

@ -293,20 +293,15 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
if (narg < 4)
error->all(FLERR,"Illegal fix colvars command: too few arguments");
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use fix colvars without atom IDs defined");
if (atom->rmass_flag)
error->all(FLERR,"Cannot use fix colvars for atoms with rmass attribute");
if (instances > 0)
error->all(FLERR,"Only one fix colvars can be active at a time");
error->all(FLERR,"Only one colvars fix can be active at a time");
++instances;
scalar_flag = 1;
global_freq = 1;
nevery = 1;
extscalar = 1;
restart_global = 1;
me = comm->me;
@ -320,13 +315,17 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
/* parse optional arguments */
int argsdone = 4;
while (argsdone+1 < narg) {
while (argsdone < narg) {
// we have keyword/value pairs. check if value is missing
if (argsdone+1 == narg)
error->all(FLERR,"Missing argument to keyword");
if (0 == strcmp(arg[argsdone], "input")) {
inp_name = strdup(arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "output")) {
out_name = strdup(arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "seed")) {
rng_seed = atoi(arg[argsdone+1]);
rng_seed = force->inumeric(FLERR,arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "unwrap")) {
if (0 == strcmp(arg[argsdone+1], "yes")) {
unwrap_flag = 1;
@ -349,6 +348,7 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
tstat_id = -1;
energy = 0.0;
nlevels_respa = 0;
init_flag = 0;
num_coords = 0;
coords = forces = oforce = NULL;
comm_buf = NULL;
@ -369,15 +369,7 @@ FixColvars::~FixColvars()
memory->sfree(inp_name);
memory->sfree(out_name);
memory->sfree(tmp_name);
deallocate();
--instances;
}
/* ---------------------------------------------------------------------- */
void FixColvars::deallocate()
{
memory->destroy(comm_buf);
memory->sfree(comm_buf);
if (proxy) {
delete proxy;
@ -386,21 +378,7 @@ void FixColvars::deallocate()
delete hashtable;
}
proxy = NULL;
idmap = NULL;
coords = NULL;
forces = NULL;
oforce = NULL;
comm_buf = NULL;
}
/* ---------------------------------------------------------------------- */
void FixColvars::post_run()
{
deallocate();
memory->sfree(inp_name);
inp_name = strdup(out_name);
--instances;
}
/* ---------------------------------------------------------------------- */
@ -412,127 +390,99 @@ int FixColvars::setmask()
mask |= MIN_POST_FORCE;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= POST_RUN;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
// initial setup of colvars run.
// initial checks for colvars run.
void FixColvars::init()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use fix colvars without atom IDs");
if (atom->map_style == 0)
error->all(FLERR,"Fix colvars requires an atom map");
if ((me == 0) && (update->whichflag == 2))
error->warning(FLERR,"Using fix colvars with minimization");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
int i,nme,tmp,ndata;
void FixColvars::setup(int vflag)
{
const int * const tag = atom->tag;
const int * const type = atom->type;
int i,nme,tmp,ndata,nlocal_max,tag_max,max;
int nlocal = atom->nlocal;
MPI_Status status;
MPI_Request request;
// collect a list of atom type by atom id for the entire system.
// the colvar module requires this information to set masses. :-(
// one time initialization
if (init_flag == 0) {
init_flag = 1;
int *typemap,*type_buf;
int nlocal_max,tag_max,max;
const int * const tag = atom->tag;
const int * const type = atom->type;
int nlocal = atom->nlocal;
// now create and initialize the colvars proxy
max=0;
for (i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
MPI_Allreduce(&max,&tag_max,1,MPI_INT,MPI_MAX,world);
MPI_Allreduce(&nlocal,&nlocal_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
if (me == 0) {
typemap = new int[tag_max+1];
memset(typemap,0,sizeof(int)*tag_max);
}
type_buf = new int[2*nlocal_max];
// input (= restart) name == "NULL" means, no restart.
if (inp_name) {
if (strcmp(inp_name,"NULL") == 0) {
memory->sfree(inp_name);
inp_name = NULL;
}
}
if (me == 0) {
for (i=0; i<nlocal; ++i)
typemap[tag[i]] = type[i];
// try to determine thermostat target temperature
double t_target = 0.0;
if (tmp_name) {
if (strcmp(tmp_name,"NULL") == 0)
tstat_id = -1;
else {
tstat_id = modify->find_fix(tmp_name);
if (tstat_id < 0) error->one(FLERR,"Could not find tstat fix ID");
double *tt = (double*)modify->fix[tstat_id]->extract("t_target",tmp);
if (tt) t_target = *tt;
}
}
// loop over procs to receive and apply remote data
for (i=1; i < comm->nprocs; ++i) {
MPI_Irecv(type_buf, 2*nlocal_max, MPI_INT, i, 0, world, &request);
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
MPI_Wait(&request, &status);
MPI_Get_count(&status, MPI_INT, &ndata);
for (int k=0; k<ndata; k+=2)
typemap[type_buf[k]] = type_buf[k+1];
proxy = new colvarproxy_lammps(lmp,inp_name,out_name,rng_seed,t_target);
proxy->init(conf_file);
coords = proxy->get_coords();
forces = proxy->get_forces();
oforce = proxy->get_oforce();
num_coords = coords->size();
}
} else { // me != 0
// copy tag/type data into communication buffer
// send the list of all colvar atom IDs to all nodes.
// also initialize and build hashtable on master.
nme = 0;
for (i=0; i<nlocal; ++i) {
type_buf[nme] = tag[i];
type_buf[nme+1] = type[i];
nme +=2;
}
/* blocking receive to wait until it is our turn to send data. */
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
MPI_Rsend(type_buf, nme, MPI_INT, 0, 0, world);
}
MPI_Bcast(&num_coords, 1, MPI_INT, 0, world);
memory->create(taglist,num_coords,"colvars:taglist");
memory->create(force_buf,3*num_coords,"colvars:force_buf");
// now create and initialize the colvar proxy
if (me == 0) {
std::vector<int> *tags_list = proxy->get_tags();
std::vector<int> &tl = *tags_list;
inthash_t *hashtable=new inthash_t;
inthash_init(hashtable, num_coords);
idmap = (void *)hashtable;
if (me == 0) {
if (inp_name) {
if (strcmp(inp_name,"NULL") == 0) {
memory->sfree(inp_name);
inp_name = NULL;
for (i=0; i < num_coords; ++i) {
taglist[i] = tl[i];
inthash_insert(hashtable, tl[i], i);
}
}
double t_target = 0.0;
if (tmp_name) {
if (strcmp(tmp_name,"NULL") == 0)
tstat_id = -1;
else {
tstat_id = modify->find_fix(tmp_name);
if (tstat_id < 0) error->one(FLERR,"Could not find tstat fix ID");
double *tt = (double*)modify->fix[tstat_id]->extract("t_target",tmp);
if (tt) t_target = *tt;
}
}
proxy = new colvarproxy_lammps(lmp,conf_file,inp_name,out_name,
rng_seed,t_target,typemap);
coords = proxy->get_coords();
forces = proxy->get_forces();
oforce = proxy->get_oforce();
num_coords = coords->size();
}
// send the list of all colvar atom IDs to all nodes.
// also initialize and build hashtable on master.
MPI_Bcast(&num_coords, 1, MPI_INT, 0, world);
memory->create(taglist,num_coords,"colvars:taglist");
memory->create(force_buf,3*num_coords,"colvars:force_buf");
if (me == 0) {
std::vector<int> *tags_list = proxy->get_tags();
std::vector<int> &tl = *tags_list;
inthash_t *hashtable=new inthash_t;
inthash_init(hashtable, num_coords);
idmap = (void *)hashtable;
for (i=0; i < num_coords; ++i) {
taglist[i] = tl[i];
inthash_insert(hashtable, tl[i], i);
}
}
MPI_Bcast(taglist, num_coords, MPI_INT, 0, world);
MPI_Bcast(taglist, num_coords, MPI_INT, 0, world);
} // end of one time initialization
// determine size of comm buffer
nme=0;
@ -583,6 +533,11 @@ void FixColvars::init()
cd[i].y = x[k][1];
cd[i].z = x[k][2];
}
if (atom->rmass_flag) {
cd[i].m = atom->rmass[k];
} else {
cd[i].m = atom->mass[type[k]];
}
}
}
@ -604,6 +559,7 @@ void FixColvars::init()
cd[j].x = comm_buf[k].x;
cd[j].y = comm_buf[k].y;
cd[j].z = comm_buf[k].z;
cd[j].m = comm_buf[k].m;
of[j].x = of[j].y = of[j].z = 0.0;
}
}
@ -634,6 +590,12 @@ void FixColvars::init()
comm_buf[nme].z = x[k][2];
}
if (atom->rmass_flag) {
comm_buf[nme].m = atom->rmass[k];
} else {
comm_buf[nme].m = atom->mass[type[k]];
}
++nme;
}
}
@ -642,19 +604,17 @@ void FixColvars::init()
MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
}
// clear temporary storage
if (me == 0) delete typemap;
delete type_buf;
}
// run pre-run setup in colvarproxy
proxy->setup();
/* ---------------------------------------------------------------------- */
void FixColvars::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
// initialize forces
if (strstr(update->integrate_style,"verlet") || (update->whichflag == 2))
post_force(vflag);
else
post_force_respa(vflag,0,0);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
@ -665,7 +625,7 @@ void FixColvars::post_force(int vflag)
// some housekeeping: update status of the proxy as needed.
if (me == 0) {
if (proxy->want_exit())
error->one(FLERR,"Run halted on request from colvars module.\n");
error->one(FLERR,"Run aborted on request from colvars module.\n");
if (tstat_id < 0) {
proxy->set_temperature(0.0);
@ -921,6 +881,31 @@ void FixColvars::end_of_step()
/* ---------------------------------------------------------------------- */
void FixColvars::write_restart(FILE *fp)
{
if (me == 0) {
std::string rest_text("");
proxy->serialize_status(rest_text);
const char *cvm_state = rest_text.c_str();
int len = strlen(cvm_state) + 1; // need to include terminating NULL byte.
fwrite(&len,sizeof(int),1,fp);
fwrite(cvm_state,1,len,fp);
}
}
/* ---------------------------------------------------------------------- */
void FixColvars::restart(char *buf)
{
init();
if (me == 0) {
std::string rest_text(buf);
proxy->deserialize_status(rest_text);
}
}
/* ---------------------------------------------------------------------- */
double FixColvars::compute_scalar()
{
return energy;

View File

@ -42,16 +42,16 @@ class FixColvars : public Fix {
virtual int setmask();
virtual void init();
virtual void setup(int);
virtual void min_setup(int vflag) {setup(vflag);};
virtual void min_post_force(int);
virtual void post_force(int);
virtual void post_force_respa(int, int, int);
virtual void post_run();
virtual void end_of_step();
virtual double compute_scalar();
virtual double memory_usage();
protected:
void deallocate(); // free internal buffers
virtual void write_restart(FILE *);
virtual void restart(char *);
protected:
class colvarproxy_lammps *proxy; // pointer to the colvars proxy class
@ -81,6 +81,7 @@ class FixColvars : public Fix {
int nlevels_respa; // flag to determine respa levels.
int store_forces; // flag to determine whether to store total forces
int unwrap_flag; // 1 if atom coords are unwrapped, 0 if not
int init_flag; // 1 if initialized, 0 if not
static int instances; // count fix instances, since colvars currently
// only supports one instance at a time
};
@ -90,3 +91,55 @@ class FixColvars : public Fix {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot use fix colvars for atoms with rmass attribute
The colvars library assigns atom masses per atom type, thus atom styles
which allow setting individual per atom masses are not supported.
E: Missing argument to keyword
Self-explanatory. A keyword was recognized, but no corresponding value
found. Check the input script syntax and compare to the documentation
for the command.
E: Incorrect fix colvars unwrap flag
Self-explanatory. Check the input script syntax.
E: Unknown fix colvars parameter
Self-explanatory. Check your input script syntax.
E: Cannot use fix colvars without atom IDs
Atom IDs are not defined, but fix colvars needs them to identify an atom.
E: Fix colvars requires an atom map
Use the atom_modify command to create an atom map.
W: Using fix colvars with minimization
Some of the functionality supported with the colvars library is not
meaningful with minimization calculations.
E: Could not find tstat fix ID
Self-explanatory. The thermostat fix ID provided with the tstat keyword
is not defined (anymore). Check your input file.
E: Run aborted on request from colvars module
Some error condition happened inside the colvars library that prohibits
it from continuing. Please examine the output for additional information.
*/