Skip to content
Snippets Groups Projects
Commit 7d94b05f authored by Davide Branduardi's avatar Davide Branduardi
Browse files

libmolfile plugin in driver

parent 8ff491c1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 153 additions and 5 deletions
src/cltools/Driver.cpp
+ 153
− 5
View file @ 7d94b05f
......@@ -28,9 +28,16 @@
#include <cstdio>
#include <string>
#include <vector>
#include <map>
#include "tools/Units.h"
#include "tools/PDB.h"
// when using molfile plugin
#ifdef __PLUMED_HAS_MOLFILE
#include "libmolfile_plugin.h"
#include "molfile_plugin.h"
#endif
using namespace std;
namespace PLMD {
......@@ -66,8 +73,47 @@ had we run the calculation we are running with driver when the MD simulation was
plumed driver --plumed plumed.dat --ixyz trajectory.xyz --trajectory-stride 100 --timestep 0.001
\endverbatim
Additionally, when you recompile with the -D__PLUMED_HAS_MOLFILE you have access to all the types of file read
by vmd. In order to do so, you should add
\verbatim
DYNAMIC_LIBS= [all the usual libs] \\
-lmolfile_plugin -ltcl8.5 -L/mypathtomolfilelibrary/ -L/mypathtotcl
CPPFLAGS= [all the usual flags] \\
-D__PLUMED_HAS_MOLFILE -I/mypathtolibmolfile_plugin.h/ -I/mypathtomolfile_plugin.h/
\endverbatim
The easy way to get this is to download the SOURCE of VMD and you find a plugins directory. Just adapt build.sh and
compile it. At the end, you should get the molfile plugin compiled. Just locate libmolfile_plugin.a ibmolfile_plugin.h
and molfile_plugin.h.
When you succeed you can use the following
\verbatim
plumed driver --plumed plumed.dat --pdb diala.pdb --mf_xtc traj.xtc --trajectory-stride 100 --timestep 0.001
\endverbatim
where the prefix --mf_ stands for the molfile plugin format accepted. There are many available. Check them in VMD.
Limitations therein applies.
*/
//+ENDPLUMEDOC
//
#ifdef __PLUMED_HAS_MOLFILE
static vector<molfile_plugin_t *> plugins;
static map <string, unsigned> pluginmap;
static int register_cb(void *v, vmdplugin_t *p){
//const char *key = p->name;
std::pair<std::map<string,unsigned>::iterator,bool> ret;
ret = pluginmap.insert ( std::pair<string,unsigned>(string(p->name),plugins.size()) );
if (ret.second==false) {
//cerr<<"MOLFILE: found duplicate plugin for "<<key<<" : not inserted "<<endl;
}else{
//cerr<<"MOLFILE: loading plugin "<<key<<" number "<<plugins.size()-1<<endl;
plugins.push_back((molfile_plugin_t *)p);
}
return VMDPLUGIN_SUCCESS;
}
#endif
template<typename real>
class Driver : public CLTool {
......@@ -98,15 +144,23 @@ void Driver<real>::registerKeywords( Keywords& keys ){
keys.add("hidden","--debug-pd","use a fake particle decomposition");
keys.add("hidden","--debug-grex","use a fake gromacs-like replica exchange");
keys.add("hidden","--debug-grex-log","log file for debug=grex");
#ifdef __PLUMED_HAS_MOLFILE
MOLFILE_INIT_ALL
MOLFILE_REGISTER_ALL(NULL, register_cb)
for(int i=0;i<plugins.size();i++){
string kk="--mf_"+string(plugins[i]->name);
string mm=" molfile: the trajectory in "+string(plugins[i]->name)+" format " ;
cerr<<"REGISTERING "<<kk<<mm<<endl;
keys.add("atoms",kk,mm);
}
#endif
}
template<typename real>
Driver<real>::Driver(const CLToolOptions& co ):
CLTool(co)
{
inputdata=commandline;
}
template<typename real>
string Driver<real>::description()const{ return "analyze trajectories with plumed"; }
......@@ -188,12 +242,33 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
string trajectory_fmt;
#ifdef __PLUMED_HAS_MOLFILE
molfile_plugin_t *api=NULL;
bool use_molfile=false;
void *h_in=NULL;
molfile_timestep_t ts_in; // this is the structure that has the timestep
#endif
// Read in an xyz file
string trajectoryFile(""), pdbfile("");
bool pbc_cli_given=false; vector<double> pbc_cli_box(9,0.0);
if(!noatoms){
std::string traj_xyz; parse("--ixyz",traj_xyz);
std::string traj_gro; parse("--igro",traj_gro);
#ifdef __PLUMED_HAS_MOLFILE
for(int i=0;i<plugins.size();i++){
string molfile_key="--mf_"+string(plugins[i]->name);
string traj_molfile;
parse(molfile_key,traj_molfile);
if(traj_molfile.length()>0){
fprintf(out,"\nDRIVER: Found molfile format trajectory %s with name %s\n",plugins[i]->name,traj_molfile.c_str());
trajectoryFile=traj_molfile;
trajectory_fmt=string(plugins[i]->name);
use_molfile=true;
api = plugins[i];
}
}
#endif
if(traj_xyz.length()>0 && traj_gro.length()>0){
fprintf(stderr,"ERROR: cannot provide more than one trajectory file\n");
if(multi_log)fclose(multi_log);
......@@ -270,19 +345,29 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
trajectoryFile+="."+n;
}
int natoms;
FILE* fp=NULL; FILE* fp_forces=NULL;
if(!noatoms){
if (trajectoryFile=="-")
fp=in;
else {
#ifdef __PLUMED_HAS_MOLFILE
if(use_molfile==true){
h_in = api->open_file_read(trajectoryFile.c_str(), trajectory_fmt.c_str(), &natoms);
}else{
#endif
fp=fopen(trajectoryFile.c_str(),"r");
if(!fp){
string msg="ERROR: Error opening XYZ file "+trajectoryFile;
string msg="ERROR: Error opening trajectory file "+trajectoryFile;
fprintf(stderr,"%s\n",msg.c_str());
return 1;
}
#ifdef __PLUMED_HAS_MOLFILE
}
#endif
}
if(dumpforces.length()>0){
if(Communicator::initialized() && pc.Get_size()>1){
string n;
......@@ -317,14 +402,36 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
while(true){
if(!noatoms){
#ifdef __PLUMED_HAS_MOLFILE
if(use_molfile==true){
int rc;
rc = api->read_next_timestep(h_in, natoms, &ts_in);
//if(rc==MOLFILE_SUCCESS){
// printf(" read this one :success \n");
//}
if(rc==MOLFILE_EOF){
printf(" read this one :eof or error \n");
break;
}
}else{
#endif
if(!Tools::getline(fp,line)) break;
#ifdef __PLUMED_HAS_MOLFILE
}
#endif
} else if ( plumedStopCondition ) break;
// add the condition to close the file and die
int natoms;
bool first_step=false;
if(!noatoms){
#ifdef __PLUMED_HAS_MOLFILE
if(use_molfile==false){
#endif
if(trajectory_fmt=="gro") if(!Tools::getline(fp,line)) error("premature end of trajectory file");
sscanf(line.c_str(),"%100d",&natoms);
#ifdef __PLUMED_HAS_MOLFILE
}
#endif
}
if(checknatoms<0 && !noatoms){
pd_nlocal=natoms;
......@@ -332,6 +439,7 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
first_step=true;
masses.assign(natoms,real(1.0));
charges.assign(natoms,real(0.0));
//case pdb: structure
if(pdbfile.length()>0){
for(unsigned i=0;i<pdb.size();++i){
AtomNumber an=pdb.getAtomNumbers()[i];
......@@ -414,6 +522,42 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
p.cmd("setStep",&step);
if(!noatoms){
#ifdef __PLUMED_HAS_MOLFILE
if(use_molfile){
if(pbc_cli_given==false) {
// info on the cell: convert using pbcset.tcl from pbctools in vmd distribution
double cosBC=cos(ts_in.alpha*pi/180.);
//double sinBC=sin(ts_in.alpha*pi/180.);
double cosAC=cos(ts_in.beta*pi/180.);
double cosAB=cos(ts_in.gamma*pi/180.);
double sinAB=sin(ts_in.gamma*pi/180.);
double Ax=ts_in.A;
double Bx=ts_in.B*cosAB;
double By=ts_in.B*sinAB;
// If sinAB is zero, then we can't determine C uniquely since it's defined
// in terms of the angle between A and B.
double Cx,Cy,Cz;
if (sinAB>0.){
Cx=ts_in.C*cosAC;
Cy=ts_in.C*cosBC-cosAC*cosAB/sinAB;
Cz=ts_in.C*pow(1.-Cx*Cx-Cy*Cy,0.5);
}else{Cx=0.;Cy=0.;Cz=0.;}
//convert to nm
cell[0]=Ax/10.;cell[1]=0.;cell[2]=0.;
cell[3]=Bx/10.;cell[4]=By/10.;cell[5]=0.;
cell[6]=Cx/10.;cell[7]=Cy/10.;cell[8]=Cz/10.;
//cerr<<"CELL "<<cell[0]<<" "<<cell[1]<<" "<<cell[2]<<" "<<cell[3]<<" "<<cell[4]<<" "<<cell[5]<<" "<<cell[6]<<" "<<cell[7]<<" "<<cell[8]<<endl;
}else{
for(unsigned i=0;i<9;i++)cell[i]=pbc_cli_box[i];
}
// info on coords
// the order is xyzxyz...
for(int i=0;i<3*natoms;i++){
coordinates[i]=real(ts_in.coords[i])/10.; //convert to nm
//cerr<<"COOR "<<coordinates[i]<<endl;
}
}else{
#endif
if(trajectory_fmt=="xyz"){
if(!Tools::getline(fp,line)) error("premature end of trajectory file");
......@@ -469,6 +613,10 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
if(words.size()>8) Tools::convert(words[8],cell[7]);
}
#ifdef __PLUMED_HAS_MOLFILE
}
#endif
if(debug_dd){
for(int i=0;i<dd_nlocal;++i){
int kk=dd_gatindex[i];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment