From e89d66f623368fbaa79be3a840629b2df11d5263 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sat, 30 Jan 2016 17:09:59 +0000
Subject: [PATCH] Added functionality to calculate multicolvar density in
specific part of box.
---
src/gridtools/FindContour.cpp | 29 ++++++++++-
src/multicolvar/MultiColvarDensity.cpp | 66 ++++++++++++++++++++++++--
2 files changed, 89 insertions(+), 6 deletions(-)
diff --git a/src/gridtools/FindContour.cpp b/src/gridtools/FindContour.cpp
index 889f757e5..0a1c135ed 100644
--- a/src/gridtools/FindContour.cpp
+++ b/src/gridtools/FindContour.cpp
@@ -36,6 +36,7 @@ private:
std::string fmt_xyz;
unsigned mycomp;
double contour;
+ std::vector<bool> nosearch_dirs;
public:
static void registerKeywords( Keywords& keys );
explicit FindContour(const ActionOptions&ao);
@@ -54,13 +55,20 @@ void FindContour::registerKeywords( Keywords& keys ){
// We want a better way of doing this bit
keys.add("compulsory", "FILE", "file on which to output coordinates");
keys.add("compulsory", "UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
+ keys.add("compulsory","SEARCHDIR","all","In which directions do you wish to search for the contour. By default the code searches "
+ "for points on the dividing surface in all directions. You may whish, however to search along one direction "
+ "for points on the dividing surface. You can specify the a subset of directions as numbers so if your grid "
+ "is thre dimensional specifying SEARCHDIR=1,2 means that only points found along grid lines parallel to the "
+ "first and second axis of the grid will be looked for. The code will not search for the dividing surface "
+ "along the third axis of the grid");
keys.add("optional","COMPONENT","if your input is a vector field use this to specifiy the component of the input vector field for which you wish to find the contour");
keys.add("optional", "PRECISION","The number of digits in trajectory file");
}
FindContour::FindContour(const ActionOptions&ao):
Action(ao),
-ActionWithInputGrid(ao)
+ActionWithInputGrid(ao),
+nosearch_dirs( mygrid->getDimension() )
{
if( mygrid->getNumberOfComponents()==1 ){
mycomp=0;
@@ -74,6 +82,24 @@ ActionWithInputGrid(ao)
parse("CONTOUR",contour);
log.printf(" calculating dividing surface along which function equals %f \n", contour);
+ std::string searchdir_str; parse("SEARCHDIR",searchdir_str);
+ if( searchdir_str=="all" ){
+ nosearch_dirs.assign( nosearch_dirs.size(), false );
+ } else {
+ std::vector<std::string> searchdirs = Tools::getWords(searchdir_str,"\t\n ,");
+ if( searchdirs.size()>mygrid->getDimension() ) error("number of search directions is larger than number of input dimensions");
+ nosearch_dirs.assign( nosearch_dirs.size(), true );
+ unsigned nn; Tools::convert(searchdirs[0],nn); nosearch_dirs[nn-1]=false;
+ log.printf(" searching for contour along %d", nn);
+ for(unsigned i=1;i<searchdirs.size()-1;++i){
+ Tools::convert(searchdirs[i],nn);
+ nosearch_dirs[nn-1]=false; log.printf(", %d",nn);
+ }
+ if( searchdirs.size()>1 ){
+ Tools::convert(searchdirs[searchdirs.size()-1],nn);
+ nosearch_dirs[nn-1]=false; log.printf(" and %d directions \n",nn);
+ } else log.printf(" direction \n");
+ }
// START OF BIT TO IMPROVE
std::string file; parse("FILE",file);
if(file.length()==0) error("name out output file was not specified");
@@ -133,6 +159,7 @@ void FindContour::performOperationsWithGrid( const bool& from_update ){
bool edge=false;
for(unsigned j=0;j<mygrid->getDimension();++j){
+ if( nosearch_dirs[j] ) continue ;
// Make sure we don't search at the edge of the grid
if( !mygrid->isPeriodic(j) && (ind[j]+1)==nbin[j] ) continue;
else if( (ind[j]+1)==nbin[j] ){ edge=true; ind[j]=0; }
diff --git a/src/multicolvar/MultiColvarDensity.cpp b/src/multicolvar/MultiColvarDensity.cpp
index d76188970..6e9a6153c 100644
--- a/src/multicolvar/MultiColvarDensity.cpp
+++ b/src/multicolvar/MultiColvarDensity.cpp
@@ -64,6 +64,8 @@ class MultiColvarDensity :
MultiColvarBase* mycolv;
std::vector<unsigned> nbins;
std::vector<double> gspacing;
+ std::vector<bool> confined;
+ std::vector<double> cmin, cmax;
vesselbase::StoreDataVessel* stash;
gridtools::HistogramOnGrid* mygrid;
Vector origin;
@@ -99,6 +101,21 @@ void MultiColvarDensity::registerKeywords( Keywords& keys ){
"in plumed plumed can be found in \\ref kernelfunctions.");
keys.addFlag("FRACTIONAL",false,"use fractional coordinates on the x-axis");
keys.addFlag("NOMEMORY",false,"do a block averaging rather than a cumulative average");
+ keys.addFlag("XREDUCED",false,"limit the calculation of the density/average to a portion of the z-axis only");
+ keys.add("optional","XLOWER","this is required if you are using XREDUCED. It specifes the lower bound for the region of the x-axis that for "
+ "which you are calculating the density/average");
+ keys.add("optional","XUPPER","this is required if you are using XREDUCED. It specifes the upper bound for the region of the x-axis that for "
+ "which you are calculating the density/average");
+ keys.addFlag("YREDUCED",false,"limit the calculation of the density/average to a portion of the y-axis only");
+ keys.add("optional","YLOWER","this is required if you are using YREDUCED. It specifes the lower bound for the region of the y-axis that for "
+ "which you are calculating the density/average");
+ keys.add("optional","YUPPER","this is required if you are using YREDUCED. It specifes the upper bound for the region of the y-axis that for "
+ "which you are calculating the density/average");
+ keys.addFlag("ZREDUCED",false,"limit the calculation of the density/average to a portion of the z-axis only");
+ keys.add("optional","ZLOWER","this is required if you are using ZREDUCED. It specifes the lower bound for the region of the z-axis that for "
+ "which you are calculating the density/average");
+ keys.add("optional","ZUPPER","this is required if you are using ZREDUCED. It specifes the upper bound for the region of the z-axis that for "
+ "which you are calculating the density/average");
}
MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
@@ -152,8 +169,42 @@ MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
parseVector("NBINS",nbins); parseVector("SPACING",gspacing);
if( nbins.size()!=directions.size() && gspacing.size()!=directions.size() ) error("NBINS or SPACING must be set");
+ confined.resize( directions.size() ); cmin.resize( directions.size() ); cmax.resize( directions.size() );
+ for(unsigned i=0;i<directions.size();++i){
+ if( directions[i]==0 ){
+ bool tflag; parseFlag("XREDUCED",tflag); confined[i]=tflag;
+ if( confined[i] ){
+ cmin[i]=cmax[i]=0.0; parse("XLOWER",cmin[i]); parse("XUPPER",cmax[i]);
+ if( fractional ) error("XREDUCED is incompatible with FRACTIONAL");
+ if( fabs(cmin[i]-cmax[i])<epsilon ) error("range set for x axis makes no sense");
+ log.printf(" confining calculation in x direction to between %f and %f \n",cmin[i],cmax[i]);
+ }
+ } else if( directions[i]==1 ){
+ bool tflag; parseFlag("YREDUCED",tflag); confined[i]=tflag;
+ if( confined[i] ){
+ cmin[i]=cmax[i]=0.0; parse("YLOWER",cmin[i]); parse("YUPPER",cmax[i]);
+ if( fractional ) error("YREDUCED is incompatible with FRACTIONAL");
+ if( fabs(cmin[i]-cmax[i])<epsilon ) error("range set for y axis makes no sense");
+ log.printf(" confining calculation in y direction to between %f and %f \n",cmin[i],cmax[i]);
+ }
+ } else if( directions[i]==2 ){
+ bool tflag; parseFlag("ZREDUCED",tflag); confined[i]=tflag;
+ if( confined[i] ){
+ cmin[i]=cmax[i]=0.0; parse("ZLOWER",cmin[i]); parse("ZUPPER",cmax[i]);
+ if( fractional ) error("ZREDUCED is incompatible with FRACTIONAL");
+ if( fabs(cmin[i]-cmax[i])<epsilon ) error("range set for z axis search makes no sense");
+ log.printf(" confining calculation in z direction to between %f and %f \n",cmin[i],cmax[i]);
+ }
+ }
+ }
+
std::string vstring = getKeyword("KERNEL") + " " + getKeyword("BANDWIDTH");
- vstring += " PBC=T"; for(unsigned i=1;i<directions.size();++i) vstring+=",T";
+ if( confined[0] ) vstring +=" PBC=F";
+ else vstring += " PBC=T";
+ for(unsigned i=1;i<directions.size();++i){
+ if( confined[i] ) vstring += ",F";
+ else vstring += ",T";
+ }
vstring +=" COMPONENTS=" + getPntrToArgument()->getLabel() + ".dens";
vstring +=" COORDINATES=";
if( directions[0]==0 ) vstring+="x";
@@ -172,7 +223,6 @@ MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
Keywords keys; gridtools::HistogramOnGrid::registerKeywords( keys );
vesselbase::VesselOptions dar( da, keys );
mygrid = new gridtools::HistogramOnGrid(dar); addVessel( mygrid );
-
// Enusre units for cube files are set correctly
if( !fractional ){
if( plumed.getAtoms().usingNaturalUnits() ) mygrid->setCubeUnits( 1.0/0.5292 );
@@ -194,11 +244,17 @@ void MultiColvarDensity::update(){
std::vector<std::string> gmin(directions.size()), gmax(directions.size());;
for(unsigned i=0;i<directions.size();++i){ min[i]=-0.5; max[i]=0.5; }
if( !fractional ){
- if( !mycolv->getPbc().isOrthorombic() ) error("I think that density profiles with non-orthorhombic cells don't work. If you want it have a look and see if you can work it out");
+ if( !mycolv->getPbc().isOrthorombic() ){
+ error("I think that density profiles with non-orthorhombic cells don't work. If you want it have a look and see if you can work it out");
+ }
for(unsigned i=0;i<directions.size();++i){
- min[i]*=mycolv->getBox()(directions[i],directions[i]);
- max[i]*=mycolv->getBox()(directions[i],directions[i]);
+ if( !confined[i] ){
+ min[i]*=mycolv->getBox()(directions[i],directions[i]);
+ max[i]*=mycolv->getBox()(directions[i],directions[i]);
+ } else {
+ min[i]=cmin[i]; max[i]=cmax[i];
+ }
}
}
for(unsigned i=0;i<directions.size();++i){ Tools::convert(min[i],gmin[i]); Tools::convert(max[i],gmax[i]); }
--
GitLab