From eb7eed2a1c656d69bbd0de405ff56d428072b648 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Thu, 31 Dec 2015 12:04:54 +0000
Subject: [PATCH] Added find contour action to replace WillardChandler surface
action
---
regtest/adjmat/rt-dfg-wcsurf/plumed.dat | 4 +-
.../crystallization/rt-wcsurface/plumed.dat | 3 +-
src/gridtools/ActionWithInputGrid.cpp | 1 +
src/gridtools/FindContour.cpp | 169 +++++++++++
src/gridtools/GridVessel.cpp | 111 ++++---
src/gridtools/GridVessel.h | 15 +-
src/multicolvar/WillardChandlerSurface.cpp | 272 ------------------
src/tools/Minimise1DBrent.h | 6 +-
8 files changed, 251 insertions(+), 330 deletions(-)
create mode 100644 src/gridtools/FindContour.cpp
delete mode 100644 src/multicolvar/WillardChandlerSurface.cpp
diff --git a/regtest/adjmat/rt-dfg-wcsurf/plumed.dat b/regtest/adjmat/rt-dfg-wcsurf/plumed.dat
index 8a89ab0b3..08f843273 100644
--- a/regtest/adjmat/rt-dfg-wcsurf/plumed.dat
+++ b/regtest/adjmat/rt-dfg-wcsurf/plumed.dat
@@ -6,4 +6,6 @@ clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1
# DUMPMULTICOLVAR DATA=clust1 ORIGIN=1 UNITS=A FILE=atoms.xyz
# MULTICOLVARDENS DATA=clust1 ORIGIN=1 DIR=xyz NBINS=30,30,30 BANDWIDTH=1.0,1.0,1.0 RUN=1 OFILE=dens.cube
-WCSURFACE DATA=clust1 ORIGIN=1 FILE=surface.xyz NBINS=30,30,30 BANDWIDTH=1.0,1.0,1.0 UNITS=A CONTOUR=250.0 PRECISION=4
+dens: MULTICOLVARDENS DATA=clust1 ORIGIN=1 DIR=xyz NBINS=30,30,30 BANDWIDTH=1.0,1.0,1.0
+FIND_CONTOUR GRID=dens CONTOUR=250.0 FILE=surface.xyz PRECISION=4 UNITS=A
+# WCSURFACE DATA=clust1 ORIGIN=1 FILE=surface.xyz NBINS=30,30,30 BANDWIDTH=1.0,1.0,1.0 UNITS=A CONTOUR=250.0 PRECISION=4
diff --git a/regtest/crystallization/rt-wcsurface/plumed.dat b/regtest/crystallization/rt-wcsurface/plumed.dat
index 3f4be0507..452ac0e42 100644
--- a/regtest/crystallization/rt-wcsurface/plumed.dat
+++ b/regtest/crystallization/rt-wcsurface/plumed.dat
@@ -6,5 +6,4 @@ DUMPMULTICOLVAR DATA=fcc ORIGIN=1 FILE=atoms.xyz PRECISION=4
dens: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,28 BANDWIDTH=1.0,1.0,1.0
PRINT_CUBE GRID=dens STRIDE=1 FILE=dens.cube
-
-WCSURFACE DATA=fcc ORIGIN=1 NBINS=14,14,28 BANDWIDTH=1.0,1.0,1.0 CONTOUR=0.42 FILE=surface.xyz PRECISION=4
+FIND_CONTOUR GRID=dens STRIDE=1 CONTOUR=0.42 FILE=surface.xyz PRECISION=4
diff --git a/src/gridtools/ActionWithInputGrid.cpp b/src/gridtools/ActionWithInputGrid.cpp
index 7ee8323ba..7596d020b 100644
--- a/src/gridtools/ActionWithInputGrid.cpp
+++ b/src/gridtools/ActionWithInputGrid.cpp
@@ -82,6 +82,7 @@ void ActionWithInputGrid::update(){
}
void ActionWithInputGrid::runFinalJobs(){
+ if( !single_run ) return ;
if( unormalised ) norm = 1.0;
else norm=1.0/mygrid->getNorm();
performOperationsWithGrid( false );
diff --git a/src/gridtools/FindContour.cpp b/src/gridtools/FindContour.cpp
new file mode 100644
index 000000000..001ec6ce9
--- /dev/null
+++ b/src/gridtools/FindContour.cpp
@@ -0,0 +1,169 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2015 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "core/PlumedMain.h"
+#include "core/Atoms.h"
+#include "ActionWithInputGrid.h"
+#include "tools/OFile.h"
+#include "tools/RootFindingBase.h"
+
+namespace PLMD {
+namespace gridtools {
+
+class FindContour : public ActionWithInputGrid {
+private:
+ OFile of;
+ double lenunit;
+ std::string fmt_xyz;
+ unsigned mycomp;
+ double contour;
+public:
+ static void registerKeywords( Keywords& keys );
+ explicit FindContour(const ActionOptions&ao);
+ void performOperationsWithGrid( const bool& from_update );
+ double getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der );
+ unsigned getNumberOfDerivatives(){ return 0; }
+ void performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {}
+ bool isPeriodic(){ return false; }
+};
+
+PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
+
+void FindContour::registerKeywords( Keywords& keys ){
+ ActionWithInputGrid::registerKeywords( keys );
+ keys.add("compulsory","CONTOUR","the value we would like to draw the contour at in the space");
+// 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("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)
+{
+ if( mygrid->getNumberOfComponents()==1 ){
+ mycomp=0;
+ } else {
+ int tcomp=-1; parse("COMPONENT",tcomp);
+ if( tcomp<0 ) error("component of vector field was not specified - use COMPONENT keyword");
+ mycomp=tcomp;
+ }
+ if( mygrid->noDerivatives() ) error("cannot find contours if input grid has no derivatives");
+
+ parse("CONTOUR",contour);
+ log.printf(" calculating dividing surface along which function equals %f \n", contour);
+
+ // START OF BIT TO IMPROVE
+ std::string file; parse("FILE",file);
+ if(file.length()==0) error("name out output file was not specified");
+ std::string type=Tools::extension(file);
+ log<<" file name "<<file<<"\n";
+ if(type!="xyz") error("can only print xyz file type with DUMPMULTICOLVAR");
+
+ fmt_xyz="%f";
+ std::string precision; parse("PRECISION",precision);
+ if(precision.length()>0){
+ int p; Tools::convert(precision,p);
+ log<<" with precision "<<p<<"\n";
+ std::string a,b;
+ Tools::convert(p+5,a);
+ Tools::convert(p,b);
+ fmt_xyz="%"+a+"."+b+"f";
+ }
+
+ std::string unitname; parse("UNITS",unitname);
+ if(unitname!="PLUMED"){
+ Units myunit; myunit.setLength(unitname);
+ lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
+ }
+ else lenunit=1.0;
+ // END OF BIT TO IMPROVE
+
+ of.link(*this);
+ of.open(file);
+ checkRead();
+}
+
+double FindContour::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ){
+ return mygrid->getValueAndDerivatives( x, mycomp, der ) - contour;
+}
+
+void FindContour::performOperationsWithGrid( const bool& from_update ){
+ std::vector< std::vector<double> > contour_points( mygrid->getDimension()*mygrid->getNumberOfPoints() );
+
+ // Two points for search
+ std::vector<unsigned> ind( mygrid->getDimension() );
+ std::vector<double> dx( mygrid->getGridSpacing() );
+ std::vector<double> direction( mygrid->getDimension(), 0 );
+ // Retrieve nbin from grid (remember this returns nbin for restart)
+ std::vector<unsigned> nbin( mygrid->getNbin() );
+ for(unsigned i=0;i<nbin.size();++i){
+ if( !mygrid->isPeriodic(i) ) nbin[i]+=1;
+ }
+
+ // Run over whole grid
+ unsigned npoints=0; RootFindingBase<FindContour> mymin( this );
+ for(unsigned i=0;i<mygrid->getNumberOfPoints();++i){
+ // Get the index of the current grid point
+ mygrid->getIndices( i, ind );
+
+ // Get the value of a point on the grid
+ double val1=getGridElement( i, mycomp*(mygrid->getDimension()+1) ) - contour;
+
+ bool edge=false;
+ for(unsigned j=0;j<mygrid->getDimension();++j){
+ // 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; }
+ else ind[j]+=1;
+ double val2=getGridElement(ind,mycomp*(mygrid->getDimension()+1)) - contour;
+ if( val1*val2<0 ){
+ // Use initial point location as first guess for search
+ contour_points[npoints].resize( mygrid->getDimension() );
+ mygrid->getGridPointCoordinates( i, contour_points[npoints] );
+ // Setup direction vector
+ direction[j]=0.9999*dx[j];
+ // And do proper search for contour point
+ mymin.linesearch( direction, contour_points[npoints], &FindContour::getDifferenceFromContour );
+ direction[j]=0.0; npoints++;
+ }
+ if( mygrid->isPeriodic(j) && edge ){ edge=false; ind[j]=nbin[j]-1; }
+ else ind[j]-=1;
+ }
+
+ }
+
+ of.printf("%u\n",npoints);
+ of.printf("Points found on isocontour\n");
+ const char* nn="X"; Vector cpos;
+ for(unsigned i=0;i<npoints;++i){
+ for(unsigned j=0;j<3;++j) cpos[j]=contour_points[i][j];
+ // cpos=mycolv->getPbc().scaledToReal(fpos);
+ of.printf(("%s "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz).c_str(),nn,lenunit*cpos[0], lenunit*cpos[1], lenunit*cpos[2] );
+ of.printf("\n");
+ }
+}
+
+}
+}
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index e0496faac..9039b94d0 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -38,7 +38,6 @@ GridVessel::GridVessel( const vesselbase::VesselOptions& da ):
Vessel(da),
noderiv(false),
bounds_set(false),
-bold(0),
cube_units(1.0)
{
std::vector<std::string> compnames; parseVector("COMPONENTS",compnames);
@@ -48,9 +47,6 @@ cube_units(1.0)
str_min.resize( dimension); str_max.resize( dimension );
parseFlag("NOMEMORY",nomemory);
- tmp_indices.resize( str_min.size() );
- current_neigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
-
unsigned n=0; nper=compnames.size()*( 1 + coordnames.size() );
arg_names.resize( coordnames.size() + compnames.size()*( 1 + coordnames.size() ) );
for(unsigned i=0;i<coordnames.size();++i){ arg_names[n] = coordnames[i]; n++; }
@@ -136,13 +132,17 @@ unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
return index;
}
-unsigned GridVessel::getIndex( const std::vector<double>& point ) const {
- plumed_dbg_assert( bounds_set && point.size()==dimension );
- std::vector<unsigned> indices(dimension);
+void GridVessel::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
+ plumed_dbg_assert( bounds_set && point.size()==dimension && indices.size()==dimension );
for(unsigned i=0;i<dimension;++i){
- indices[i]=std::floor( (point[i] - min[i])/dx[i] );
+ indices[i]=std::floor( (point[i] - min[i])/dx[i] );
if( pbc[i] ) indices[i]=indices[i]%nbin[i];
}
+}
+
+unsigned GridVessel::getIndex( const std::vector<double>& point ) const {
+ plumed_dbg_assert( bounds_set && point.size()==dimension );
+ std::vector<unsigned> indices(dimension); getIndices( point, indices );
return getIndex( indices );
}
@@ -169,42 +169,22 @@ void GridVessel::getGridPointCoordinates( const unsigned& ipoint , std::vector<d
for(unsigned i=0;i<dimension;++i) x[i] = min[i] + dx[i]*tindices[i];
}
-unsigned GridVessel::getLocationOnGrid( const std::vector<double>& x, std::vector<double>& dd ){
- plumed_dbg_assert( bounds_set && x.size()==dimension && dd.size()==dimension );
- getIndices( bold, tmp_indices ); bool changebox=false;
- for(unsigned i=0;i<dimension;++i){
- double bb = x[i] - min[i];
- if ( bb<0.0 || bb>dx[i]*nbin[i] ){
- getAction()->error("Extrapolation of function is not allowed");
- } else if( bb<tmp_indices[i]*dx[i] || bb>(tmp_indices[i]+1)*dx[i] ){
- tmp_indices[i]=static_cast<unsigned>( std::floor(bb/dx[i]) );
- changebox=true;
+void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const {
+ mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
+
+ std::vector<unsigned> tmp_indices( dimension );
+ std::vector<unsigned> my_indices( dimension );
+ getIndices( mybox, my_indices );
+ for(unsigned i=0;i<mysneigh.size();++i){
+ unsigned tmp=i;
+ for(unsigned j=0;j<dimension;++j){
+ unsigned i0=tmp%2+my_indices[j]; tmp/=2;
+ if(!pbc[j] && i0==nbin[j]) getAction()->error("Extrapolating function on grid");
+ if( pbc[j] && i0==nbin[j]) i0=0;
+ tmp_indices[j]=i0;
}
- dd[i] = bb/dx[i] - static_cast<double>(tmp_indices[i]);
+ mysneigh[i]=getIndex( tmp_indices );
}
- if(changebox) bold = getIndex( tmp_indices );
- return bold;
-}
-
-void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ){
- if( mysneigh.size()!=current_neigh.size() ) mysneigh.resize( current_neigh.size() );
-
- if( bold!=mybox ){
- std::vector<unsigned> my_indices( dimension );
- getIndices( mybox, my_indices );
- for(unsigned i=0;i<current_neigh.size();++i){
- unsigned tmp=i;
- for(unsigned j=0;j<dimension;++j){
- unsigned i0=tmp%2+my_indices[j]; tmp/=2;
- if(!pbc[j] && i0==nbin[j]) getAction()->error("Extrapolating function on grid");
- if( pbc[j] && i0==nbin[j]) i0=0;
- tmp_indices[j]=i0;
- }
- current_neigh[i]=getIndex( tmp_indices );
- }
- bold=mybox;
- }
- for(unsigned i=0;i<current_neigh.size();++i) mysneigh[i]=current_neigh[i];
}
double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
@@ -311,6 +291,53 @@ std::string GridVessel::getInputString() const {
return mstring;
}
+double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const {
+ plumed_dbg_assert( der.size()==dimension && !noderiv && ind<getNumberOfComponents() );
+
+ double X,X2,X3,value=0; der.assign(der.size(),0.0);
+ std::vector<double> fd(dimension);
+ std::vector<double> C(dimension);
+ std::vector<double> D(dimension);
+ std::vector<double> dder(dimension);
+
+ std::vector<unsigned> nindices(dimension);
+ std::vector<unsigned> indices(dimension); getIndices( x, indices );
+ std::vector<unsigned> neigh; getSplineNeighbors( getIndex(indices), neigh );
+ std::vector<double> xfloor(dimension); getGridPointCoordinates( getIndex(x), xfloor );
+
+// loop over neighbors
+ for(unsigned int ipoint=0;ipoint<neigh.size();++ipoint){
+ double grid=getGridElement(neigh[ipoint], ind*(1+dimension) );
+ for(unsigned j=0;j<dimension;++j) dder[j] = getGridElement( neigh[ipoint], ind*(1+dimension) + 1 + j );
+
+ getIndices( neigh[ipoint], nindices );
+ double ff=1.0;
+
+ for(unsigned j=0;j<dimension;++j){
+ int x0=1;
+ if(nindices[j]==indices[j]) x0=0;
+ double ddx=dx[j];
+ X=fabs((x[j]-xfloor[j])/ddx-(double)x0);
+ X2=X*X;
+ X3=X2*X;
+ double yy;
+ if(fabs(grid)<0.0000001) yy=0.0;
+ else yy=-dder[j]/grid;
+ C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*ddx;
+ D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*ddx;
+ D[j]*=(x0?-1.0:1.0)/ddx;
+ ff*=C[j];
+ }
+ for(unsigned j=0;j<dimension;++j){
+ fd[j]=D[j];
+ for(unsigned i=0;i<dimension;++i) if(i!=j) fd[j]*=C[i];
+ }
+ value+=grid*ff;
+ for(unsigned j=0;j<dimension;++j) der[j]+=grid*fd[j];
+ }
+ return value;
+}
+
}
}
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index 84684c3f2..cf303df71 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -36,15 +36,10 @@ private:
bool noderiv;
/// Have the minimum and maximum for the grid been set
bool bounds_set;
-/// These two variables are used to
-/// remember the box we were in on the previous call
- unsigned bold;
/// The number of points in the grid
unsigned npoints;
/// Units for Gaussian Cube file
double cube_units;
-/// Remember the neighbors that were used last time
- std::vector<unsigned> current_neigh;
/// The names of the various columns in the grid file
std::vector<std::string> arg_names;
/// The minimum and maximum of the grid stored as doubles
@@ -53,8 +48,6 @@ private:
std::vector<unsigned> stride;
/// The number of bins in each grid direction
std::vector<unsigned> nbin;
-/// A tempory array that can be used to store indices for a point
- std::vector<unsigned> tmp_indices;
/// Flatten the grid and get the grid index for a point
unsigned getIndex( const std::vector<unsigned>& indices ) const ;
/// The grid point that was requested last by getGridPointCoordinates
@@ -78,6 +71,8 @@ protected:
/// Get the set of points neighouring a particular location in space
void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
+/// Get the indices of a particular point
+ void getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const ;
/// Convert a point in space the the correspoinding grid point
unsigned getIndex( const std::vector<double>& p ) const ;
public:
@@ -130,10 +125,8 @@ public:
double getCellVolume() const ;
/// Get the value of the ith grid element
double getGridElement( const unsigned&, const unsigned& ) const ;
-/// Get the numerical index for the box that contains a particular point
- unsigned getLocationOnGrid( const std::vector<double>& x, std::vector<double>& dd );
/// Get the points neighboring a particular spline point
- void getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh );
+ void getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const ;
/// Get the spacing between grid points
const std::vector<double>& getGridSpacing() const ;
/// Get the extent of the grid in one of the axis
@@ -151,6 +144,8 @@ public:
std::string getInputString() const ;
/// Does this have derivatives
bool noDerivatives() const ;
+/// Get the value and derivatives at a particular location using spline interpolation
+ double getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const ;
};
inline
diff --git a/src/multicolvar/WillardChandlerSurface.cpp b/src/multicolvar/WillardChandlerSurface.cpp
deleted file mode 100644
index f9dc0fd25..000000000
--- a/src/multicolvar/WillardChandlerSurface.cpp
+++ /dev/null
@@ -1,272 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- Copyright (c) 2015,2016 The plumed team
- (see the PEOPLE file at the root of the distribution for a list of names)
-
- See http://www.plumed.org for more information.
-
- This file is part of plumed, version 2.
-
- plumed is free software: you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- plumed is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
-
- You should have received a copy of the GNU Lesser General Public License
- along with plumed. If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "core/ActionAtomistic.h"
-#include "core/ActionPilot.h"
-#include "core/ActionRegister.h"
-#include "tools/Pbc.h"
-#include "tools/File.h"
-#include "core/PlumedMain.h"
-#include "core/Atoms.h"
-#include "tools/Units.h"
-#include <cstdio>
-#include "core/SetupMolInfo.h"
-#include "core/ActionSet.h"
-#include "MultiColvarBase.h"
-#include "tools/Grid.h"
-#include "tools/OFile.h"
-#include "tools/KernelFunctions.h"
-#include "vesselbase/ActionWithInputVessel.h"
-#include "vesselbase/StoreDataVessel.h"
-
-using namespace std;
-
-namespace PLMD {
-namespace multicolvar {
-
-//+PLUMEDOC ANALYSIS WCSURFACE
-/*
-Calculate the Willard Chanlder dividing surface
-
-\par Examples
-
-
-*/
-//+ENDPLUMEDOC
-
-class WillardChandlerSurface :
- public ActionPilot,
- public ActionAtomistic,
- public vesselbase::ActionWithInputVessel
-{
- std::string kerneltype;
- unsigned rstride;
- std::string filename;
- unsigned dir;
- double contour;
- OFile of;
- double lenunit;
- std::string fmt_xyz;
- MultiColvarBase* mycolv;
- std::vector<unsigned> nbins;
- std::vector<double> bw;
- std::vector<bool> nosearch_dirs, confined;
- std::vector<double> cmin, cmax;
-public:
- explicit WillardChandlerSurface(const ActionOptions&);
- static void registerKeywords( Keywords& keys );
- void calculate(){}
- void calculateNumericalDerivatives( ActionWithValue* a=NULL ){ plumed_error(); }
- void apply(){}
- void update();
-};
-
-PLUMED_REGISTER_ACTION(WillardChandlerSurface,"WCSURFACE")
-
-void WillardChandlerSurface::registerKeywords( Keywords& keys ){
- Action::registerKeywords( keys );
- ActionPilot::registerKeywords( keys );
- ActionAtomistic::registerKeywords( keys );
- ActionWithInputVessel::registerKeywords( keys );
- keys.add("compulsory","STRIDE","1","the frequency with which the data should be collected and added to the grid");
- keys.add("atoms","ORIGIN","we will use the position of this atom as the origin");
- keys.add("compulsory","NBINS","the number of bins to use to represent the density profile");
- keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
- keys.add("compulsory","KERNEL","gaussian","the kernel function you are using. More details on the kernels available "
- "in plumed plumed can be found in \\ref kernelfunctions.");
- keys.add("compulsory","CONTOUR","the value we would like to draw the contour at in the space");
-// 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("optional", "PRECISION","The number of digits in trajectory file");
- keys.addFlag("NO_X_SEARCH",false,"");
- keys.addFlag("XREDUCED",false,"");
- keys.add("optional","XLOWER","");
- keys.add("optional","XUPPER","");
- keys.addFlag("NO_Y_SEARCH",false,"");
- keys.addFlag("YREDUCED",false,"");
- keys.add("optional","YLOWER","");
- keys.add("optional","YUPPER","");
- keys.addFlag("NO_Z_SEARCH",false,"");
- keys.addFlag("ZREDUCED",false,"");
- keys.add("optional","ZLOWER","");
- keys.add("optional","ZUPPER","");
-// Better way for this bit
-}
-
-WillardChandlerSurface::WillardChandlerSurface(const ActionOptions&ao):
- Action(ao),
- ActionPilot(ao),
- ActionAtomistic(ao),
- ActionWithInputVessel(ao),
- nosearch_dirs(3,false),
- confined(3,false),
- cmin(3),
- cmax(3)
-{
- std::vector<AtomNumber> atom; parseAtomList("ORIGIN",atom);
- if( atom.size()!=1 ) error("should only be one atom specified");
- log.printf(" origin is at position of atom : %d\n",atom[0].serial() );
-
- readArgument("store");
- mycolv = dynamic_cast<MultiColvarBase*>( getDependencies()[0] );
- plumed_assert( getDependencies().size()==1 );
- if(!mycolv) error("action labeled " + mycolv->getLabel() + " is not a multicolvar");
-
- parse("KERNEL",kerneltype); parseVector("BANDWIDTH",bw); parseVector("NBINS",nbins); parse("CONTOUR",contour);
- if( bw.size()!=3 || nbins.size()!=3 ) error("wrong size found for bandwidth vector or number of bins");
- log.printf(" calculating dividing surface where average equals %f along ", contour);
- log.printf(" for colvars calculated by action %s \n",mycolv->getLabel().c_str() );
-
- // Read in confinement stuff
- bool tflag;
- parseFlag("XREDUCED",tflag); confined[0]=tflag;
- parseFlag("NO_X_SEARCH",tflag); nosearch_dirs[0]=tflag;
- if( confined[0] ){
- cmin[0]=cmax[0]=0.0;
- parse("XLOWER",cmin[0]); parse("XUPPER",cmax[0]);
- if( fabs(cmin[0]-cmax[0])<epsilon ) error("range set for x axis search makes no sense");
- log.printf(" confining search in x direction to between %f and %f \n",cmin[0],cmax[0]);
- }
- parseFlag("YREDUCED",tflag); confined[1]=tflag;
- parseFlag("NO_Y_SEARCH",tflag); nosearch_dirs[1]=tflag;
- if( confined[1] ){
- cmin[1]=cmax[1]=0.0;
- parse("YLOWER",cmin[1]); parse("YUPPER",cmax[1]);
- if( fabs(cmin[1]-cmax[1])<epsilon ) error("range set for y axis search makes no sense");
- log.printf(" confining search in y direction to between %f and %f \n",cmin[1],cmax[1]);
- }
- parseFlag("ZREDUCED",tflag); confined[2]=tflag;
- parseFlag("NO_Z_SEARCH",tflag); nosearch_dirs[2]=tflag;
- if( confined[2] ){
- cmin[2]=cmax[2]=0.0;
- parse("ZLOWER",cmin[2]); parse("ZUPPER",cmax[2]);
- if( fabs(cmin[2]-cmax[2])<epsilon ) error("range set for z axis search makes no sense");
- log.printf(" confining search in z direction to between %f and %f \n",cmin[2],cmax[2]);
- }
-
- // START OF BIT TO IMPROVE
- std::string file; parse("FILE",file);
- if(file.length()==0) error("name out output file was not specified");
- std::string type=Tools::extension(file);
- log<<" file name "<<file<<"\n";
- if(type!="xyz") error("can only print xyz file type with DUMPMULTICOLVAR");
-
- fmt_xyz="%f";
- std::string precision; parse("PRECISION",precision);
- if(precision.length()>0){
- int p; Tools::convert(precision,p);
- log<<" with precision "<<p<<"\n";
- string a,b;
- Tools::convert(p+5,a);
- Tools::convert(p,b);
- fmt_xyz="%"+a+"."+b+"f";
- }
-
- std::string unitname; parse("UNITS",unitname);
- if(unitname!="PLUMED"){
- Units myunit; myunit.setLength(unitname);
- lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
- }
- else lenunit=1.0;
- // END OF BIT TO IMPROVE
-
- of.link(*this);
- of.open(file);
-
- checkRead(); requestAtoms(atom);
- // Stupid dependencies cleared by requestAtoms - why GBussi why? That's got me so many times
- addDependency( mycolv );
-}
-
-void WillardChandlerSurface::update(){
- std::vector<bool> pbc(nbins.size()); std::vector<double> min(nbins.size()), max(nbins.size());
- std::vector<std::string> args(nbins.size()), gmin(nbins.size()), gmax(nbins.size());;
- args[0]="x"; args[1]="y"; args[2]="z";
- for(unsigned i=0;i<3;++i){ min[i]=-0.5; max[i]=0.5; pbc[i]=!confined[i]; }
-
- if( !mycolv->getPbc().isOrthorombic() ){
- error("I think that Willard-Chandler surfaces with non-orthorhombic cells don't work. If you want it have a look and see if you can work it out");
- }
-
- // Get the box (this is going to give us the extent of the grid)
- for(unsigned i=0;i<3;++i){
- if( !confined[i] ){
- min[i]*=mycolv->getBox()(i,i); max[i]*=mycolv->getBox()(i,i);
- } else {
- min[i]=cmin[i]; max[i]=cmax[i];
- }
- }
-
- // This will be the only thing we keep eventually
- for(unsigned i=0;i<3;++i){ Tools::convert(min[i],gmin[i]); Tools::convert(max[i],gmax[i]); }
-
- // Setup the grid
- std::string funcl=mycolv->getLabel() + ".dens";
- Grid gg(funcl,args,gmin,gmax,nbins,true,true,true,pbc,gmin,gmax);
-
- vesselbase::StoreDataVessel* stash=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToArgument() );
- plumed_dbg_assert( stash ); std::vector<double> cvals( mycolv->getNumberOfQuantities() );
- Vector origin = getPosition(0); std::vector<double> pp( 3 ); Vector fpos;
-
- unsigned rank=comm.Get_rank(), size=comm.Get_size();
-
- for(unsigned i=rank;i<mycolv->getFullNumberOfTasks();i+=size){
- // Skip if task was not active on last run through
- if( !mycolv->taskIsCurrentlyActive(i) ) continue ;
-
- stash->retrieveValue( i, false, cvals );
- Vector apos = pbcDistance( mycolv->getCentralAtomPos( mycolv->getTaskCode(i) ), origin );
-
- // fpos = getPbc().realToScaled( apos ); Ideally want to do with scaled coordinates eventually GAT
- for(unsigned j=0;j<3;++j) pp[j]=apos[j];
-
- KernelFunctions kernel( pp, bw, kerneltype, false, cvals[0]*cvals[1], true );
- gg.addKernel( kernel );
- }
- gg.mpiSumValuesAndDerivatives( comm );
-
- unsigned npoints=0; std::vector<std::vector<double> > contour_points;
- gg.findSetOfPointsOnContour( contour, nosearch_dirs, npoints, contour_points );
- if(npoints==0 ) warning("found no points on Willard-Chandler surface try changing the CONTOUR parameter");
-
- of.printf("%u\n",npoints);
- const Tensor & t(mycolv->getPbc().getBox());
- if(mycolv->getPbc().isOrthorombic()){
- of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2));
- }else{
- of.printf((" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz+"\n").c_str(),
- lenunit*t(0,0),lenunit*t(0,1),lenunit*t(0,2),
- lenunit*t(1,0),lenunit*t(1,1),lenunit*t(1,2),
- lenunit*t(2,0),lenunit*t(2,1),lenunit*t(2,2)
- );
- }
- const char* nn="X"; Vector cpos;
- for(unsigned i=0;i<npoints;++i){
- for(unsigned j=0;j<3;++j) cpos[j]=contour_points[i][j];
- // cpos=mycolv->getPbc().scaledToReal(fpos);
- of.printf(("%s "+fmt_xyz+" "+fmt_xyz+" "+fmt_xyz).c_str(),nn,lenunit*cpos[0], lenunit*cpos[1], lenunit*cpos[2] );
- of.printf("\n");
- }
-}
-
-}
-}
diff --git a/src/tools/Minimise1DBrent.h b/src/tools/Minimise1DBrent.h
index ee6761e32..b0f67c84c 100644
--- a/src/tools/Minimise1DBrent.h
+++ b/src/tools/Minimise1DBrent.h
@@ -45,7 +45,7 @@ private:
/// Use to prevent any possible division by zero
const double TINY;
/// Maximum number of interactions in line minimiser
- const int ITMAX;
+ const unsigned ITMAX;
/// The value of the golden ratio
const double CGOLD;
/// A small number that protects against trying to achieve fractional
@@ -78,10 +78,10 @@ TINY(1.0E-20),
ITMAX(100),
CGOLD(0.3819660),
ZEPS(epsilon*1.0E-3),
-myclass_func(pf),
ax(0),bx(0),cx(0),
fa(0),fb(0),fc(0),
-fmin(0)
+fmin(0),
+myclass_func(pf)
{
}
--
GitLab