forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8165 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
aa8aeb0294
commit
5ea3c4d23a
|
|
@ -0,0 +1,45 @@
|
|||
# Install/unInstall package files in LAMMPS
|
||||
# edit 2 Makefile.package files to include/exclude CUDA info
|
||||
# do not install child files if parent does not exist
|
||||
|
||||
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_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
|
||||
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(colvars_SYSPATH) |' ../Makefile.package
|
||||
fi
|
||||
|
||||
if (test -e ../Makefile.package.settings) then
|
||||
sed -i -e '/^include.*colvars.*$/d' ../Makefile.package.settings
|
||||
# multiline form needed for BSD sed on Macs
|
||||
sed -i -e '4 i \
|
||||
include ..\/..\/lib\/colvars\/Makefile.lammps\
|
||||
' ../Makefile.package.settings
|
||||
|
||||
fi
|
||||
|
||||
cp colvarproxy_lammps.h ..
|
||||
cp colvarproxy_lammps.cpp ..
|
||||
cp fix_colvars.h ..
|
||||
cp fix_colvars.cpp ..
|
||||
|
||||
elif (test $1 = 0) then
|
||||
|
||||
if (test -e ../Makefile.package) then
|
||||
sed -i -e 's/[^ \t]*colvars[^ \t]* //g' ../Makefile.package
|
||||
fi
|
||||
|
||||
if (test -e ../Makefile.package.settings) then
|
||||
sed -i -e '/^include.*colvars.*$/d' ../Makefile.package.settings
|
||||
fi
|
||||
|
||||
rm -f ../colvarproxy_lammps.h
|
||||
rm -f ../colvarproxy_lammps.cpp
|
||||
rm -f ../fix_colvars.h
|
||||
rm -f ../fix_colvars.cpp
|
||||
fi
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
This package implements the "fix colvars" command which can be used
|
||||
in a LAMMPS input script.
|
||||
|
||||
This fix allows to use "collective variables" to implement
|
||||
Adaptive Biasing Force, Metadynamics, Steered MD, Umbrella
|
||||
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).
|
||||
|
||||
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,
|
||||
J. Chem. Theory Comput., 6, 35-47 (2010).
|
||||
|
||||
- the colvars fix and a thin interface layer, which exchanges
|
||||
information between LAMMPS and the collective variable module.
|
||||
|
||||
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
|
||||
|
|
@ -0,0 +1,425 @@
|
|||
|
||||
#include "mpi.h"
|
||||
#include "lammps.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "output.h"
|
||||
#include "random_park.h"
|
||||
|
||||
#include "fix_colvars.h"
|
||||
|
||||
#include "colvarmodule.h"
|
||||
#include "colvaratoms.h"
|
||||
#include "colvarproxy.h"
|
||||
#include "colvarproxy_lammps.h"
|
||||
|
||||
#include <sys/types.h>
|
||||
#include <sys/stat.h>
|
||||
#include <unistd.h>
|
||||
|
||||
#include <cerrno>
|
||||
#include <cstdio>
|
||||
#include <cstring>
|
||||
|
||||
#include <iostream>
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
|
||||
#define HASH_FAIL -1
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// local helper functions
|
||||
|
||||
// safely move filename to filename.extension
|
||||
static void my_backup_file(const char *filename, const char *extension)
|
||||
{
|
||||
struct stat sbuf;
|
||||
if (stat(filename, &sbuf) == 0) {
|
||||
if ( ! extension ) extension = ".BAK";
|
||||
char *backup = new char[strlen(filename)+strlen(extension)+1];
|
||||
strcpy(backup, filename);
|
||||
strcat(backup, extension);
|
||||
#if defined(_WIN32) && !defined(__CYGWIN__)
|
||||
remove(backup);
|
||||
#endif
|
||||
if ( rename(filename,backup) ) {
|
||||
char *sys_err_msg = strerror(errno);
|
||||
if ( !sys_err_msg ) sys_err_msg = (char *) "(unknown error)";
|
||||
fprintf(stderr,"Error renaming file %s to %s: %s\n",
|
||||
filename, backup, sys_err_msg);
|
||||
}
|
||||
delete [] backup;
|
||||
}
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
|
||||
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)
|
||||
{
|
||||
if (cvm::debug())
|
||||
log("Initializing the colvars proxy object.\n");
|
||||
|
||||
_random = new LAMMPS_NS::RanPark(lmp,seed);
|
||||
|
||||
first_timestep=true;
|
||||
system_force_requested=false;
|
||||
previous_step=-1;
|
||||
t_target=temp;
|
||||
do_exit=false;
|
||||
|
||||
// set input restart name and strip the extension, if present
|
||||
input_prefix_str = std::string(inp_name ? inp_name : "");
|
||||
if (input_prefix_str.rfind(".colvars.state") != std::string::npos)
|
||||
input_prefix_str.erase(input_prefix_str.rfind(".colvars.state"),
|
||||
std::string(".colvars.state").size());
|
||||
|
||||
// output prefix is always given
|
||||
output_prefix_str = std::string(out_name);
|
||||
// not so for restarts
|
||||
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);
|
||||
}
|
||||
// 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);
|
||||
|
||||
// create the colvarmodule instance
|
||||
colvars = new colvarmodule(conf_file,this);
|
||||
|
||||
if (_lmp->update->ntimestep != 0) {
|
||||
cvm::log ("Initializing step number as firstTimestep.\n");
|
||||
colvars->it = colvars->it_restart = _lmp->update->ntimestep;
|
||||
}
|
||||
|
||||
if (cvm::debug()) {
|
||||
log("colvars_atoms = "+cvm::to_str(colvars_atoms)+"\n");
|
||||
log("colvars_atoms_ncopies = "+cvm::to_str(colvars_atoms_ncopies)+"\n");
|
||||
log("positions = "+cvm::to_str(positions)+"\n");
|
||||
log("applied_forces = "+cvm::to_str(applied_forces)+"\n");
|
||||
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;
|
||||
colvars = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
// trigger colvars computation
|
||||
double colvarproxy_lammps::compute()
|
||||
{
|
||||
if (first_timestep) {
|
||||
first_timestep = false;
|
||||
} else {
|
||||
// Use the time step number inherited from LAMMPS
|
||||
if ( _lmp->update->ntimestep - previous_step == 1 )
|
||||
colvars->it++;
|
||||
// Other cases could mean:
|
||||
// - run 0
|
||||
// - beginning of a new run statement
|
||||
// then the internal counter should not be incremented
|
||||
}
|
||||
previous_step = _lmp->update->ntimestep;
|
||||
|
||||
if (cvm::debug()) {
|
||||
cvm::log(cvm::line_marker+
|
||||
"colvarproxy_lammps, step no. "+cvm::to_str(colvars->it)+"\n"+
|
||||
"Updating internal data.\n");
|
||||
}
|
||||
|
||||
// call the collective variable module
|
||||
colvars->calc();
|
||||
|
||||
#if 0
|
||||
for (int i=0; i < colvars_atoms.size(); ++i) {
|
||||
fprintf(stderr,"CV: atom %d/%d/%d pos: %g %g %g for: %g %g %g\n",
|
||||
colvars_atoms[i], colvars_atoms_ncopies[i],
|
||||
positions[i].type, positions[i].x, positions[i].y, positions[i].z,
|
||||
applied_forces[i].x, applied_forces[i].y, applied_forces[i].z);
|
||||
}
|
||||
#endif
|
||||
|
||||
return bias_energy;
|
||||
}
|
||||
|
||||
cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1,
|
||||
cvm::atom_pos const &pos2)
|
||||
{
|
||||
double xtmp = pos2.x - pos1.x;
|
||||
double ytmp = pos2.y - pos1.y;
|
||||
double ztmp = pos2.z - pos1.z;
|
||||
_lmp->domain->minimum_image(xtmp,ytmp,ztmp);
|
||||
return cvm::rvector(xtmp, ytmp, ztmp);
|
||||
}
|
||||
|
||||
cvm::real colvarproxy_lammps::position_dist2(cvm::atom_pos const &pos1,
|
||||
cvm::atom_pos const &pos2)
|
||||
{
|
||||
double xtmp = pos2.x - pos1.x;
|
||||
double ytmp = pos2.y - pos1.y;
|
||||
double ztmp = pos2.z - pos1.z;
|
||||
_lmp->domain->minimum_image(xtmp,ytmp,ztmp);
|
||||
return cvm::real(xtmp*xtmp + ytmp*ytmp + ztmp*ztmp);
|
||||
}
|
||||
|
||||
|
||||
inline void colvarproxy_lammps::select_closest_image(cvm::atom_pos &pos,
|
||||
cvm::atom_pos const &ref)
|
||||
{
|
||||
double xtmp = pos.x - ref.x;
|
||||
double ytmp = pos.y - ref.y;
|
||||
double ztmp = pos.z - ref.z;
|
||||
_lmp->domain->minimum_image(xtmp,ytmp,ztmp);
|
||||
pos.x = ref.x + xtmp;
|
||||
pos.y = ref.y + ytmp;
|
||||
pos.z = ref.z + ztmp;
|
||||
}
|
||||
|
||||
void colvarproxy_lammps::log(std::string const &message)
|
||||
{
|
||||
std::istringstream is(message);
|
||||
std::string line;
|
||||
while (std::getline(is, line)) {
|
||||
if (_lmp->screen)
|
||||
fprintf(_lmp->screen,"colvars: %s\n",line.c_str());
|
||||
if (_lmp->logfile)
|
||||
fprintf(_lmp->logfile,"colvars: %s\n",line.c_str());
|
||||
}
|
||||
}
|
||||
|
||||
void colvarproxy_lammps::fatal_error(std::string const &message)
|
||||
{
|
||||
log(message);
|
||||
if (!cvm::debug())
|
||||
log("If this error message is unclear, try recompiling the "
|
||||
"colvars library and LAMMPS with -DCOLVARS_DEBUG.\n");
|
||||
|
||||
_lmp->error->one(FLERR,
|
||||
"Fatal error in the collective variables module.\n");
|
||||
}
|
||||
|
||||
void colvarproxy_lammps::exit(std::string const &message)
|
||||
{
|
||||
log(message);
|
||||
log("Request to exit the simulation made.\n");
|
||||
do_exit=true;
|
||||
}
|
||||
|
||||
enum e_pdb_field {
|
||||
e_pdb_none,
|
||||
e_pdb_occ,
|
||||
e_pdb_beta,
|
||||
e_pdb_x,
|
||||
e_pdb_y,
|
||||
e_pdb_z,
|
||||
e_pdb_ntot
|
||||
};
|
||||
|
||||
|
||||
e_pdb_field pdb_field_str2enum(std::string const &pdb_field_str)
|
||||
{
|
||||
e_pdb_field pdb_field = e_pdb_none;
|
||||
|
||||
if (colvarparse::to_lower_cppstr(pdb_field_str) ==
|
||||
colvarparse::to_lower_cppstr("O")) {
|
||||
pdb_field = e_pdb_occ;
|
||||
}
|
||||
|
||||
if (colvarparse::to_lower_cppstr(pdb_field_str) ==
|
||||
colvarparse::to_lower_cppstr("B")) {
|
||||
pdb_field = e_pdb_beta;
|
||||
}
|
||||
|
||||
if (colvarparse::to_lower_cppstr(pdb_field_str) ==
|
||||
colvarparse::to_lower_cppstr("X")) {
|
||||
pdb_field = e_pdb_x;
|
||||
}
|
||||
|
||||
if (colvarparse::to_lower_cppstr(pdb_field_str) ==
|
||||
colvarparse::to_lower_cppstr("Y")) {
|
||||
pdb_field = e_pdb_y;
|
||||
}
|
||||
|
||||
if (colvarparse::to_lower_cppstr(pdb_field_str) ==
|
||||
colvarparse::to_lower_cppstr("Z")) {
|
||||
pdb_field = e_pdb_z;
|
||||
}
|
||||
|
||||
if (pdb_field == e_pdb_none) {
|
||||
cvm::fatal_error("Error: unsupported PDB field, \""+
|
||||
pdb_field_str+"\".\n");
|
||||
}
|
||||
|
||||
return pdb_field;
|
||||
}
|
||||
|
||||
void colvarproxy_lammps::load_coords(char const *pdb_filename,
|
||||
std::vector<cvm::atom_pos> &pos,
|
||||
const std::vector<int> &indices,
|
||||
std::string const pdb_field_str,
|
||||
double const pdb_field_value)
|
||||
{
|
||||
|
||||
cvm::fatal_error("Reading collective variable coordinates "
|
||||
"from a PDB file is currently not supported.\n");
|
||||
}
|
||||
|
||||
void colvarproxy_lammps::load_atoms(char const *pdb_filename,
|
||||
std::vector<cvm::atom> &atoms,
|
||||
std::string const pdb_field_str,
|
||||
double const pdb_field_value)
|
||||
{
|
||||
cvm::fatal_error("Selecting collective variable atoms "
|
||||
"from a PDB file is currently not supported.\n");
|
||||
}
|
||||
|
||||
void colvarproxy_lammps::backup_file(char const *filename)
|
||||
{
|
||||
if (std::string(filename).rfind(std::string(".colvars.state"))
|
||||
!= std::string::npos) {
|
||||
my_backup_file(filename, ".old");
|
||||
} else {
|
||||
my_backup_file(filename, ".BAK");
|
||||
}
|
||||
}
|
||||
|
||||
int colvarproxy_lammps::init_lammps_atom(const int &aid, cvm::atom *atom)
|
||||
{
|
||||
atom->id = aid;
|
||||
atom->mass = _lmp->atom->mass[_typemap[aid]];
|
||||
|
||||
for (size_t i = 0; i < colvars_atoms.size(); i++) {
|
||||
if (colvars_atoms[i] == aid) {
|
||||
// this atom id was already recorded
|
||||
colvars_atoms_ncopies[i] += 1;
|
||||
return i;
|
||||
}
|
||||
}
|
||||
|
||||
// allocate a new slot for this atom
|
||||
colvars_atoms_ncopies.push_back(1);
|
||||
colvars_atoms.push_back(aid);
|
||||
struct commdata c;
|
||||
c.tag = aid;
|
||||
c.type = _typemap[aid];
|
||||
c.x = c.y = c.z = 0.0;
|
||||
positions.push_back(c);
|
||||
total_forces.push_back(c);
|
||||
applied_forces.push_back(c);
|
||||
|
||||
return colvars_atoms.size()-1;
|
||||
}
|
||||
|
||||
// atom member functions, LAMMPS specific implementations
|
||||
|
||||
cvm::atom::atom(const int &id)
|
||||
{
|
||||
|
||||
if (cvm::debug())
|
||||
cvm::log("Adding atom "+cvm::to_str(id)+
|
||||
" for collective variables calculation.\n");
|
||||
|
||||
if (id < 0)
|
||||
cvm::fatal_error("Error: invalid atom ID specified, "+
|
||||
cvm::to_str(id)+"\n");
|
||||
|
||||
int idx = ((colvarproxy_lammps *) cvm::proxy)->init_lammps_atom(id,this);
|
||||
if (idx < 0)
|
||||
cvm::fatal_error("Error: atom ID , "+cvm::to_str(id)+" does not exist.\n");
|
||||
|
||||
this->index = idx;
|
||||
if (cvm::debug())
|
||||
cvm::log("The index of this atom in the colvarproxy_lammps arrays is "+
|
||||
cvm::to_str(this->index)+".\n");
|
||||
|
||||
this->reset_data();
|
||||
}
|
||||
|
||||
|
||||
/// For AMBER topologies, the segment id is automatically set to
|
||||
/// "MAIN" (the segment id assigned by NAMD's AMBER topology parser),
|
||||
/// and is therefore optional when an AMBER topology is used
|
||||
cvm::atom::atom(cvm::residue_id const &residue,
|
||||
std::string const &atom_name,
|
||||
std::string const &segment_id)
|
||||
{
|
||||
cvm::fatal_error("Creating collective variable atoms "
|
||||
"from a PDB file is currently not supported.\n");
|
||||
}
|
||||
|
||||
|
||||
// copy constructor
|
||||
cvm::atom::atom(cvm::atom const &a)
|
||||
: index(a.index), id(a.id), mass(a.mass)
|
||||
{
|
||||
// init_lammps_atom() has already been called by a's constructor, no
|
||||
// need to call it again
|
||||
|
||||
// need to increment the counter anyway
|
||||
colvarproxy_lammps *cp = (colvarproxy_lammps *) cvm::proxy;
|
||||
cp->colvars_atoms_ncopies[this->index] += 1;
|
||||
}
|
||||
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
void cvm::atom::read_velocity()
|
||||
{
|
||||
cvm::fatal_error("Error: read_velocity is not yet implemented.\n");
|
||||
}
|
||||
|
||||
void cvm::atom::read_system_force()
|
||||
{
|
||||
colvarproxy_lammps const * const cp = (colvarproxy_lammps *) cvm::proxy;
|
||||
this->system_force.x = cp->total_forces[this->index].x
|
||||
- cp->applied_forces[this->index].x;
|
||||
this->system_force.y = cp->total_forces[this->index].y
|
||||
- cp->applied_forces[this->index].y;
|
||||
this->system_force.z = cp->total_forces[this->index].z
|
||||
- cp->applied_forces[this->index].z;
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,141 @@
|
|||
#ifndef COLVARPROXY_LAMMPS_H
|
||||
#define COLVARPROXY_LAMMPS_H
|
||||
|
||||
#include "colvarmodule.h"
|
||||
#include "colvarproxy.h"
|
||||
|
||||
#include "lammps.h"
|
||||
#include "domain.h"
|
||||
#include "force.h"
|
||||
#include "random_park.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <iostream>
|
||||
|
||||
/* struct for packed data communication of coordinates and forces. */
|
||||
struct commdata {
|
||||
int tag,type;
|
||||
double x,y,z;
|
||||
};
|
||||
|
||||
inline std::ostream & operator<< (std::ostream &out, const commdata &cd)
|
||||
{
|
||||
out << " (" << cd.tag << "/" << cd.type << ": "
|
||||
<< cd.x << ", " << cd.y << ", " << cd.z << ") ";
|
||||
return out;
|
||||
};
|
||||
|
||||
/// \brief Communication between colvars and LAMMPS
|
||||
/// (implementation of \link colvarproxy \endlink)
|
||||
class colvarproxy_lammps : public colvarproxy {
|
||||
|
||||
// LAMMPS specific data objects and flags
|
||||
protected:
|
||||
|
||||
// pointers to LAMMPS class instances
|
||||
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;
|
||||
int restart_every;
|
||||
int previous_step;
|
||||
|
||||
std::string input_prefix_str;
|
||||
std::string output_prefix_str;
|
||||
std::string restart_prefix_str;
|
||||
|
||||
bool first_timestep;
|
||||
bool system_force_requested;
|
||||
bool do_exit;
|
||||
|
||||
std::vector<int> colvars_atoms;
|
||||
std::vector<size_t> colvars_atoms_ncopies;
|
||||
std::vector<struct commdata> positions;
|
||||
std::vector<struct commdata> total_forces;
|
||||
std::vector<struct commdata> applied_forces;
|
||||
|
||||
public:
|
||||
friend class cvm::atom;
|
||||
colvarproxy_lammps (LAMMPS_NS::LAMMPS *lmp, const char *, const char *,
|
||||
const char *, const int, const double, const int *);
|
||||
virtual ~colvarproxy_lammps();
|
||||
|
||||
// disable default and copy constructor
|
||||
private:
|
||||
colvarproxy_lammps() {};
|
||||
colvarproxy_lammps(const colvarproxy_lammps &) {};
|
||||
|
||||
// methods for lammps to move data or trigger actions in the proxy
|
||||
public:
|
||||
void set_temperature(double t) { t_target = t; };
|
||||
bool need_system_forces() const { return system_force_requested; };
|
||||
bool want_exit() const { return do_exit; };
|
||||
std::vector<int> * get_tags() { return &colvars_atoms; };
|
||||
std::vector<struct commdata> * get_coords() { return &positions; };
|
||||
std::vector<struct commdata> * get_forces() { return &applied_forces; };
|
||||
std::vector<struct commdata> * get_oforce() { return &total_forces; };
|
||||
|
||||
// initialize atom structure
|
||||
int init_lammps_atom(const int &, cvm::atom *);
|
||||
|
||||
// perform colvars computation. returns biasing energy
|
||||
double compute();
|
||||
|
||||
// implementation of pure methods from base class
|
||||
public:
|
||||
|
||||
inline cvm::real unit_angstrom() { return _lmp->force->angstrom; };
|
||||
cvm::real boltzmann() { return _lmp->force->boltz; };
|
||||
cvm::real temperature() { return t_target; };
|
||||
cvm::real dt() { return _lmp->update->dt * _lmp->force->femtosecond; };
|
||||
|
||||
inline std::string input_prefix() { return input_prefix_str; };
|
||||
inline std::string output_prefix() { return output_prefix_str; };
|
||||
inline std::string restart_output_prefix() { return restart_prefix_str; };
|
||||
inline size_t restart_frequency() { return restart_every; };
|
||||
|
||||
void add_energy (cvm::real energy) { bias_energy = energy; };
|
||||
void request_system_force (bool yesno) { system_force_requested = yesno; };
|
||||
|
||||
void log (std::string const &message);
|
||||
void fatal_error (std::string const &message);
|
||||
void exit (std::string const &message);
|
||||
|
||||
cvm::rvector position_distance (cvm::atom_pos const &pos1,
|
||||
cvm::atom_pos const &pos2);
|
||||
cvm::real position_dist2 (cvm::atom_pos const &pos1,
|
||||
cvm::atom_pos const &pos2);
|
||||
void select_closest_image (cvm::atom_pos &pos,
|
||||
cvm::atom_pos const &ref_pos);
|
||||
|
||||
void load_atoms(char const *filename,
|
||||
std::vector<cvm::atom> &atoms,
|
||||
std::string const pdb_field,
|
||||
double const pdb_field_value = 0.0);
|
||||
|
||||
void load_coords(char const *filename,
|
||||
std::vector<cvm::atom_pos> &pos,
|
||||
const std::vector<int> &indices,
|
||||
std::string const pdb_field,
|
||||
double const pdb_field_value = 0.0);
|
||||
|
||||
void backup_file(char const *filename);
|
||||
|
||||
cvm::real rand_gaussian(void) { return _random->gaussian(); };
|
||||
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
// Emacs
|
||||
// Local Variables:
|
||||
// mode: C++
|
||||
// End:
|
||||
|
|
@ -0,0 +1,864 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (TempleU)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <errno.h>
|
||||
|
||||
#include "fix_colvars.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "group.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "respa.h"
|
||||
#include "update.h"
|
||||
|
||||
#include "colvarproxy_lammps.h"
|
||||
|
||||
/* re-usable integer hash table code with static linkage. */
|
||||
|
||||
/** hash table top level data structure */
|
||||
typedef struct inthash_t {
|
||||
struct inthash_node_t **bucket; /* array of hash nodes */
|
||||
int size; /* size of the array */
|
||||
int entries; /* number of entries in table */
|
||||
int downshift; /* shift cound, used in hash function */
|
||||
int mask; /* used to select bits for hashing */
|
||||
} inthash_t;
|
||||
|
||||
/** hash table node data structure */
|
||||
typedef struct inthash_node_t {
|
||||
int data; /* data in hash node */
|
||||
int key; /* key for hash lookup */
|
||||
struct inthash_node_t *next; /* next node in hash chain */
|
||||
} inthash_node_t;
|
||||
|
||||
#define HASH_FAIL -1
|
||||
#define HASH_LIMIT 0.5
|
||||
|
||||
/* initialize new hash table */
|
||||
static void inthash_init(inthash_t *tptr, int buckets);
|
||||
/* lookup entry in hash table */
|
||||
static int inthash_lookup(void *tptr, int key);
|
||||
/* insert an entry into hash table. */
|
||||
static int inthash_insert(inthash_t *tptr, int key, int data);
|
||||
/* delete the hash table */
|
||||
static void inthash_destroy(inthash_t *tptr);
|
||||
/* adapted sort for in-place sorting of map indices. */
|
||||
static void id_sort(int *idmap, int left, int right);
|
||||
|
||||
/************************************************************************
|
||||
* integer hash code:
|
||||
************************************************************************/
|
||||
|
||||
/* inthash() - Hash function returns a hash number for a given key.
|
||||
* tptr: Pointer to a hash table, key: The key to create a hash number for */
|
||||
static int inthash(const inthash_t *tptr, int key) {
|
||||
int hashvalue;
|
||||
|
||||
hashvalue = (((key*1103515249)>>tptr->downshift) & tptr->mask);
|
||||
if (hashvalue < 0) {
|
||||
hashvalue = 0;
|
||||
}
|
||||
|
||||
return hashvalue;
|
||||
}
|
||||
|
||||
/*
|
||||
* rebuild_table_int() - Create new hash table when old one fills up.
|
||||
*
|
||||
* tptr: Pointer to a hash table
|
||||
*/
|
||||
static void rebuild_table_int(inthash_t *tptr) {
|
||||
inthash_node_t **old_bucket, *old_hash, *tmp;
|
||||
int old_size, h, i;
|
||||
|
||||
old_bucket=tptr->bucket;
|
||||
old_size=tptr->size;
|
||||
|
||||
/* create a new table and rehash old buckets */
|
||||
inthash_init(tptr, old_size<<1);
|
||||
for (i=0; i<old_size; i++) {
|
||||
old_hash=old_bucket[i];
|
||||
while(old_hash) {
|
||||
tmp=old_hash;
|
||||
old_hash=old_hash->next;
|
||||
h=inthash(tptr, tmp->key);
|
||||
tmp->next=tptr->bucket[h];
|
||||
tptr->bucket[h]=tmp;
|
||||
tptr->entries++;
|
||||
} /* while */
|
||||
} /* for */
|
||||
|
||||
/* free memory used by old table */
|
||||
free(old_bucket);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
/*
|
||||
* inthash_init() - Initialize a new hash table.
|
||||
*
|
||||
* tptr: Pointer to the hash table to initialize
|
||||
* buckets: The number of initial buckets to create
|
||||
*/
|
||||
void inthash_init(inthash_t *tptr, int buckets) {
|
||||
|
||||
/* make sure we allocate something */
|
||||
if (buckets==0)
|
||||
buckets=16;
|
||||
|
||||
/* initialize the table */
|
||||
tptr->entries=0;
|
||||
tptr->size=2;
|
||||
tptr->mask=1;
|
||||
tptr->downshift=29;
|
||||
|
||||
/* ensure buckets is a power of 2 */
|
||||
while (tptr->size<buckets) {
|
||||
tptr->size<<=1;
|
||||
tptr->mask=(tptr->mask<<1)+1;
|
||||
tptr->downshift--;
|
||||
} /* while */
|
||||
|
||||
/* allocate memory for table */
|
||||
tptr->bucket=(inthash_node_t **) calloc(tptr->size, sizeof(inthash_node_t *));
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
/*
|
||||
* inthash_lookup() - Lookup an entry in the hash table and return a pointer to
|
||||
* it or HASH_FAIL if it wasn't found.
|
||||
*
|
||||
* tptr: Pointer to the hash table
|
||||
* key: The key to lookup
|
||||
*/
|
||||
int inthash_lookup(void *ptr, int key) {
|
||||
const inthash_t *tptr = (const inthash_t *) ptr;
|
||||
int h;
|
||||
inthash_node_t *node;
|
||||
|
||||
|
||||
/* find the entry in the hash table */
|
||||
h=inthash(tptr, key);
|
||||
for (node=tptr->bucket[h]; node!=NULL; node=node->next) {
|
||||
if (node->key == key)
|
||||
break;
|
||||
}
|
||||
|
||||
/* return the entry if it exists, or HASH_FAIL */
|
||||
return(node ? node->data : HASH_FAIL);
|
||||
}
|
||||
|
||||
/*
|
||||
* inthash_insert() - Insert an entry into the hash table. If the entry already
|
||||
* exists return a pointer to it, otherwise return HASH_FAIL.
|
||||
*
|
||||
* tptr: A pointer to the hash table
|
||||
* key: The key to insert into the hash table
|
||||
* data: A pointer to the data to insert into the hash table
|
||||
*/
|
||||
int inthash_insert(inthash_t *tptr, int key, int data) {
|
||||
int tmp;
|
||||
inthash_node_t *node;
|
||||
int h;
|
||||
|
||||
/* check to see if the entry exists */
|
||||
if ((tmp=inthash_lookup(tptr, key)) != HASH_FAIL)
|
||||
return(tmp);
|
||||
|
||||
/* expand the table if needed */
|
||||
while (tptr->entries>=HASH_LIMIT*tptr->size)
|
||||
rebuild_table_int(tptr);
|
||||
|
||||
/* insert the new entry */
|
||||
h=inthash(tptr, key);
|
||||
node=(struct inthash_node_t *) malloc(sizeof(inthash_node_t));
|
||||
node->data=data;
|
||||
node->key=key;
|
||||
node->next=tptr->bucket[h];
|
||||
tptr->bucket[h]=node;
|
||||
tptr->entries++;
|
||||
|
||||
return HASH_FAIL;
|
||||
}
|
||||
|
||||
/*
|
||||
* inthash_destroy() - Delete the entire table, and all remaining entries.
|
||||
*
|
||||
*/
|
||||
void inthash_destroy(inthash_t *tptr) {
|
||||
inthash_node_t *node, *last;
|
||||
int i;
|
||||
|
||||
for (i=0; i<tptr->size; i++) {
|
||||
node = tptr->bucket[i];
|
||||
while (node != NULL) {
|
||||
last = node;
|
||||
node = node->next;
|
||||
free(last);
|
||||
}
|
||||
}
|
||||
|
||||
/* free the entire array of buckets */
|
||||
if (tptr->bucket != NULL) {
|
||||
free(tptr->bucket);
|
||||
memset(tptr, 0, sizeof(inthash_t));
|
||||
}
|
||||
}
|
||||
|
||||
/************************************************************************
|
||||
* integer list sort code:
|
||||
************************************************************************/
|
||||
|
||||
/* sort for integer map. initial call id_sort(idmap, 0, natoms - 1); */
|
||||
static void id_sort(int *idmap, int left, int right)
|
||||
{
|
||||
int pivot, l_hold, r_hold;
|
||||
|
||||
l_hold = left;
|
||||
r_hold = right;
|
||||
pivot = idmap[left];
|
||||
|
||||
while (left < right) {
|
||||
while ((idmap[right] >= pivot) && (left < right))
|
||||
right--;
|
||||
if (left != right) {
|
||||
idmap[left] = idmap[right];
|
||||
left++;
|
||||
}
|
||||
while ((idmap[left] <= pivot) && (left < right))
|
||||
left++;
|
||||
if (left != right) {
|
||||
idmap[right] = idmap[left];
|
||||
right--;
|
||||
}
|
||||
}
|
||||
idmap[left] = pivot;
|
||||
pivot = left;
|
||||
left = l_hold;
|
||||
right = r_hold;
|
||||
|
||||
if (left < pivot)
|
||||
id_sort(idmap, left, pivot-1);
|
||||
if (right > pivot)
|
||||
id_sort(idmap, pivot+1, right);
|
||||
}
|
||||
|
||||
/***************************************************************/
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
// initialize static class members
|
||||
int FixColvars::instances=0;
|
||||
|
||||
/***************************************************************
|
||||
create class and parse arguments in LAMMPS script. Syntax:
|
||||
|
||||
fix ID group-ID colvars <config_file> [optional flags...]
|
||||
|
||||
optional keyword value pairs:
|
||||
|
||||
input <input prefix> (for restarting/continuing, defaults to
|
||||
NULL, but set to <output prefix> at end)
|
||||
output <output prefix> (defaults to 'out')
|
||||
seed <integer> (seed for RNG, defaults to '1966')
|
||||
tstat <fix label> (label of thermostatting fix)
|
||||
|
||||
TODO: add (optional) arguments for RNG seed, temperature compute
|
||||
***************************************************************/
|
||||
FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, 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)
|
||||
error->all(FLERR,"Only one fix colvars can be active at a time");
|
||||
++instances;
|
||||
|
||||
scalar_flag = 1;
|
||||
global_freq = 1;
|
||||
nevery = 1;
|
||||
extscalar = 1;
|
||||
|
||||
me = comm->me;
|
||||
|
||||
conf_file = strdup(arg[3]);
|
||||
rng_seed = 1966;
|
||||
|
||||
inp_name = NULL;
|
||||
out_name = NULL;
|
||||
tmp_name = NULL;
|
||||
|
||||
/* parse optional arguments */
|
||||
int argsdone = 4;
|
||||
while (argsdone+1 < narg) {
|
||||
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]);
|
||||
} else if (0 == strcmp(arg[argsdone], "tstat")) {
|
||||
tmp_name = strdup(arg[argsdone+1]);
|
||||
} else {
|
||||
error->all(FLERR,"Unknown fix colvars parameter");
|
||||
}
|
||||
++argsdone; ++argsdone;
|
||||
}
|
||||
|
||||
if (!out_name) out_name = strdup("out");
|
||||
|
||||
/* initialize various state variables. */
|
||||
tstat_id = -1;
|
||||
energy = 0.0;
|
||||
nlevels_respa = 0;
|
||||
num_coords = 0;
|
||||
coords = forces = oforce = NULL;
|
||||
comm_buf = NULL;
|
||||
force_buf = NULL;
|
||||
proxy = NULL;
|
||||
idmap = NULL;
|
||||
|
||||
/* storage required to communicate a single coordinate or force. */
|
||||
size_one = sizeof(struct commdata);
|
||||
}
|
||||
|
||||
/*********************************
|
||||
* Clean up on deleting the fix. *
|
||||
*********************************/
|
||||
FixColvars::~FixColvars()
|
||||
{
|
||||
memory->sfree(conf_file);
|
||||
memory->sfree(inp_name);
|
||||
memory->sfree(out_name);
|
||||
memory->sfree(tmp_name);
|
||||
deallocate();
|
||||
--instances;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixColvars::deallocate()
|
||||
{
|
||||
memory->destroy(comm_buf);
|
||||
|
||||
if (proxy) {
|
||||
delete proxy;
|
||||
inthash_t *hashtable = (inthash_t *)idmap;
|
||||
inthash_destroy(hashtable);
|
||||
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);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixColvars::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= THERMO_ENERGY;
|
||||
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.
|
||||
|
||||
void FixColvars::init()
|
||||
{
|
||||
if (strstr(update->integrate_style,"respa"))
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
|
||||
|
||||
int i,nme,tmp,ndata;
|
||||
|
||||
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. :-(
|
||||
|
||||
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;
|
||||
|
||||
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) {
|
||||
typemap = new int[tag_max+1];
|
||||
memset(typemap,0,sizeof(int)*tag_max);
|
||||
}
|
||||
type_buf = new int[2*nlocal_max];
|
||||
|
||||
if (me == 0) {
|
||||
for (i=0; i<nlocal; ++i)
|
||||
typemap[tag[i]] = type[i];
|
||||
|
||||
// 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];
|
||||
}
|
||||
} else { // me != 0
|
||||
|
||||
// copy tag/type data into communication buffer
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
// now create and initialize the colvar proxy
|
||||
|
||||
if (me == 0) {
|
||||
|
||||
if (inp_name) {
|
||||
if (strcmp(inp_name,"NULL") == 0) {
|
||||
memory->sfree(inp_name);
|
||||
inp_name = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
|
||||
// determine size of comm buffer
|
||||
nme=0;
|
||||
for (i=0; i < num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal))
|
||||
++nme;
|
||||
}
|
||||
|
||||
MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world);
|
||||
memory->create(comm_buf,nmax,"colvars:comm_buf");
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
|
||||
if (me == 0) {
|
||||
|
||||
std::vector<struct commdata> &cd = *coords;
|
||||
std::vector<struct commdata> &of = *oforce;
|
||||
|
||||
// store coordinate data in holding array, clear old forces
|
||||
|
||||
for (i=0; i<num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal)) {
|
||||
of[i].tag = cd[i].tag = tag[k];
|
||||
of[i].type = cd[i].type = type[k];
|
||||
cd[i].x = x[k][0];
|
||||
cd[i].y = x[k][1];
|
||||
cd[i].z = x[k][2];
|
||||
of[i].x = of[i].y = of[i].z = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// loop over procs to receive and apply remote data
|
||||
|
||||
for (i=1; i < comm->nprocs; ++i) {
|
||||
int maxbuf = nmax*size_one;
|
||||
MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request);
|
||||
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
|
||||
MPI_Wait(&request, &status);
|
||||
MPI_Get_count(&status, MPI_BYTE, &ndata);
|
||||
ndata /= size_one;
|
||||
|
||||
for (int k=0; k<ndata; ++k) {
|
||||
const int j = inthash_lookup(idmap, comm_buf[k].tag);
|
||||
if (j != HASH_FAIL) {
|
||||
of[j].tag = cd[j].tag = comm_buf[k].tag;
|
||||
of[j].type = cd[j].type = comm_buf[k].type;
|
||||
cd[j].x = comm_buf[k].x;
|
||||
cd[j].y = comm_buf[k].y;
|
||||
cd[j].z = comm_buf[k].z;
|
||||
of[j].x = of[j].y = of[j].z = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else { // me != 0
|
||||
|
||||
// copy coordinate data into communication buffer
|
||||
|
||||
nme = 0;
|
||||
for (i=0; i<num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal)) {
|
||||
comm_buf[nme].tag = tag[k];
|
||||
comm_buf[nme].type = type[k];
|
||||
comm_buf[nme].x = x[k][0];
|
||||
comm_buf[nme].y = x[k][1];
|
||||
comm_buf[nme].z = x[k][2];
|
||||
++nme;
|
||||
}
|
||||
}
|
||||
/* blocking receive to wait until it is our turn to send data. */
|
||||
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
|
||||
MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
|
||||
}
|
||||
|
||||
// clear temporary storage
|
||||
if (me == 0) delete typemap;
|
||||
delete type_buf;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixColvars::setup(int vflag)
|
||||
{
|
||||
if (strstr(update->integrate_style,"verlet"))
|
||||
post_force(vflag);
|
||||
else
|
||||
post_force_respa(vflag,0,0);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* Main colvars handler:
|
||||
* Send coodinates and add colvar forces to atoms. */
|
||||
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");
|
||||
|
||||
if (tstat_id < 0) {
|
||||
proxy->set_temperature(0.0);
|
||||
} else {
|
||||
int tmp;
|
||||
// get thermostat target temperature from corresponding fix,
|
||||
// if the fix supports extraction.
|
||||
double *tt = (double *) modify->fix[tstat_id]->extract("t_target",tmp);
|
||||
if (tt)
|
||||
proxy->set_temperature(*tt);
|
||||
else
|
||||
proxy->set_temperature(0.0);
|
||||
}
|
||||
}
|
||||
|
||||
const int * const tag = atom->tag;
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = atom->f;
|
||||
const int nlocal = atom->nlocal;
|
||||
|
||||
/* check and potentially grow local communication buffers. */
|
||||
int i,nmax_new,nme=0;
|
||||
for (i=0; i < num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal))
|
||||
++nme;
|
||||
}
|
||||
|
||||
MPI_Allreduce(&nme,&nmax_new,1,MPI_INT,MPI_MAX,world);
|
||||
if (nmax_new > nmax) {
|
||||
nmax = nmax_new;
|
||||
memory->grow(comm_buf,nmax,"colvars:comm_buf");
|
||||
}
|
||||
|
||||
MPI_Status status;
|
||||
MPI_Request request;
|
||||
int tmp, ndata;
|
||||
|
||||
if (me == 0) {
|
||||
std::vector<struct commdata> &cd = *coords;
|
||||
|
||||
// store coordinate data
|
||||
|
||||
for (i=0; i<num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal)) {
|
||||
cd[i].x = x[k][0];
|
||||
cd[i].y = x[k][1];
|
||||
cd[i].z = x[k][2];
|
||||
}
|
||||
}
|
||||
|
||||
/* loop over procs to receive remote data */
|
||||
for (i=1; i < comm->nprocs; ++i) {
|
||||
int maxbuf = nmax*size_one;
|
||||
MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request);
|
||||
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
|
||||
MPI_Wait(&request, &status);
|
||||
MPI_Get_count(&status, MPI_BYTE, &ndata);
|
||||
ndata /= size_one;
|
||||
|
||||
for (int k=0; k<ndata; ++k) {
|
||||
const int j = inthash_lookup(idmap, comm_buf[k].tag);
|
||||
if (j != HASH_FAIL) {
|
||||
cd[j].x = comm_buf[k].x;
|
||||
cd[j].y = comm_buf[k].y;
|
||||
cd[j].z = comm_buf[k].z;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} else { // me != 0
|
||||
/* copy coordinate data into communication buffer */
|
||||
nme = 0;
|
||||
for (i=0; i<num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal)) {
|
||||
comm_buf[nme].tag = tag[k];
|
||||
comm_buf[nme].x = x[k][0];
|
||||
comm_buf[nme].y = x[k][1];
|
||||
comm_buf[nme].z = x[k][2];
|
||||
++nme;
|
||||
}
|
||||
}
|
||||
/* blocking receive to wait until it is our turn to send data. */
|
||||
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
|
||||
MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// call our workhorse and retrieve additional information.
|
||||
if (me == 0) {
|
||||
energy = proxy->compute();
|
||||
store_forces = proxy->need_system_forces();
|
||||
}
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
|
||||
// broadcast store_forces flag and energy data to all processors
|
||||
MPI_Bcast(&energy, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&store_forces, 1, MPI_INT, 0, world);
|
||||
|
||||
// broadcast and apply biasing forces
|
||||
|
||||
if (me == 0) {
|
||||
std::vector<struct commdata> &fo = *forces;
|
||||
double *fbuf = force_buf;
|
||||
for (int j=0; j < num_coords; ++j) {
|
||||
*fbuf++ = fo[j].x;
|
||||
*fbuf++ = fo[j].y;
|
||||
*fbuf++ = fo[j].z;
|
||||
}
|
||||
}
|
||||
MPI_Bcast(force_buf, 3*num_coords, MPI_DOUBLE, 0, world);
|
||||
|
||||
for (int i=0; i < num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal)) {
|
||||
f[k][0] += force_buf[3*i+0];
|
||||
f[k][1] += force_buf[3*i+1];
|
||||
f[k][2] += force_buf[3*i+2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixColvars::min_post_force(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixColvars::post_force_respa(int vflag, int ilevel, int iloop)
|
||||
{
|
||||
/* only process colvar forces on the outmost RESPA level. */
|
||||
if (ilevel == nlevels_respa-1) post_force(vflag);
|
||||
return;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixColvars::end_of_step()
|
||||
{
|
||||
if (store_forces) {
|
||||
|
||||
const int * const tag = atom->tag;
|
||||
double * const * const f = atom->f;
|
||||
const int nlocal = atom->nlocal;
|
||||
|
||||
/* check and potentially grow local communication buffers. */
|
||||
int i,nmax_new,nme=0;
|
||||
for (i=0; i < num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal))
|
||||
++nme;
|
||||
}
|
||||
|
||||
MPI_Allreduce(&nme,&nmax_new,1,MPI_INT,MPI_MAX,world);
|
||||
if (nmax_new > nmax) {
|
||||
nmax = nmax_new;
|
||||
memory->grow(comm_buf,nmax,"colvars:comm_buf");
|
||||
}
|
||||
|
||||
MPI_Status status;
|
||||
MPI_Request request;
|
||||
int tmp, ndata;
|
||||
|
||||
if (me == 0) {
|
||||
|
||||
// store old force data
|
||||
std::vector<struct commdata> &of = *oforce;
|
||||
|
||||
for (i=0; i<num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal)) {
|
||||
|
||||
const int j = inthash_lookup(idmap, tag[k]);
|
||||
if (j != HASH_FAIL) {
|
||||
of[j].x = f[k][0];
|
||||
of[j].y = f[k][1];
|
||||
of[j].z = f[k][2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* loop over procs to receive remote data */
|
||||
for (i=1; i < comm->nprocs; ++i) {
|
||||
int maxbuf = nmax*size_one;
|
||||
MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request);
|
||||
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
|
||||
MPI_Wait(&request, &status);
|
||||
MPI_Get_count(&status, MPI_BYTE, &ndata);
|
||||
ndata /= size_one;
|
||||
|
||||
for (int k=0; k<ndata; ++k) {
|
||||
const int j = inthash_lookup(idmap, comm_buf[k].tag);
|
||||
if (j != HASH_FAIL) {
|
||||
of[j].x = comm_buf[k].x;
|
||||
of[j].y = comm_buf[k].y;
|
||||
of[j].z = comm_buf[k].z;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} else { // me != 0
|
||||
/* copy total force data into communication buffer */
|
||||
nme = 0;
|
||||
for (i=0; i<num_coords; ++i) {
|
||||
const int k = atom->map(taglist[i]);
|
||||
if ((k >= 0) && (k < nlocal)) {
|
||||
comm_buf[nme].tag = tag[k];
|
||||
comm_buf[nme].x = f[k][0];
|
||||
comm_buf[nme].y = f[k][1];
|
||||
comm_buf[nme].z = f[k][2];
|
||||
++nme;
|
||||
}
|
||||
}
|
||||
/* blocking receive to wait until it is our turn to send data. */
|
||||
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
|
||||
MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixColvars::compute_scalar()
|
||||
{
|
||||
return energy;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* local memory usage. approximately. */
|
||||
double FixColvars::memory_usage(void)
|
||||
{
|
||||
double bytes = (double) (num_coords * (2*sizeof(int)+3*sizeof(double)));
|
||||
bytes += (double) (nmax*size_one) + sizeof(this);
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,98 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(colvars,FixColvars)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_COLVARS_H
|
||||
#define LMP_FIX_COLVARS_H
|
||||
|
||||
#include "fix.h"
|
||||
#include <vector>
|
||||
|
||||
// forward declaration
|
||||
class colvarproxy_lammps;
|
||||
struct commdata;
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixColvars : public Fix {
|
||||
|
||||
public:
|
||||
FixColvars(class LAMMPS *, int, char **);
|
||||
virtual ~FixColvars();
|
||||
|
||||
virtual int setmask();
|
||||
virtual void init();
|
||||
virtual void setup(int);
|
||||
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
|
||||
|
||||
protected:
|
||||
class colvarproxy_lammps *proxy; // pointer to the colvars proxy class
|
||||
char *conf_file; // name of colvars config file
|
||||
char *inp_name; // name/prefix of colvars restart file
|
||||
char *out_name; // prefix string for all output files
|
||||
char *tmp_name; // name of thermostat fix.
|
||||
int rng_seed; // seed to initialize random number generator
|
||||
int tstat_id; // id of the thermostat fix
|
||||
double energy; // biasing energy of the fix
|
||||
|
||||
int me; // my MPI rank in this "world".
|
||||
int num_coords; // total number of atoms controlled by this fix
|
||||
int *taglist; // list of all atom IDs referenced by colvars.
|
||||
std::vector<struct commdata> *coords; // coordinates of colvar atoms
|
||||
std::vector<struct commdata> *forces; // received forces of colvar atoms
|
||||
std::vector<struct commdata> *oforce; // old total forces of colvar atoms
|
||||
|
||||
int nmax; // size of atom communication buffer.
|
||||
int size_one; // bytes per atom in communication buffer.
|
||||
struct commdata *comm_buf; // communication buffer
|
||||
double *force_buf; // communication buffer
|
||||
|
||||
void *idmap; // hash for mapping atom indices to consistent order.
|
||||
int *rev_idmap; // list of the hash keys for reverse mapping.
|
||||
|
||||
int nlevels_respa; // flag to determine respa levels.
|
||||
int store_forces; // flag to determine whether to store total forces
|
||||
static int instances; // count fix instances, since colvars currently
|
||||
// only supports one instance at a time
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
// Local Variables:
|
||||
// mode: c++
|
||||
// compile-command: "make -j4 openmpi"
|
||||
// c-basic-offset: 2
|
||||
// fill-column: 76
|
||||
// indent-tabs-mode: nil
|
||||
// End:
|
||||
Loading…
Reference in New Issue