From 318ffa1ce5323e63ed2638a334e9ec63413dd3fe Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sat, 30 Aug 2014 17:18:24 +0200
Subject: [PATCH] Added SMAC collective variable

---
 src/crystallization/SMAC.cpp  | 82 +++++++++++++++++++++++++++++++++++
 src/tools/KernelFunctions.cpp | 24 ++++++++++
 src/tools/KernelFunctions.h   |  1 +
 3 files changed, 107 insertions(+)
 create mode 100644 src/crystallization/SMAC.cpp

diff --git a/src/crystallization/SMAC.cpp b/src/crystallization/SMAC.cpp
new file mode 100644
index 000000000..524737115
--- /dev/null
+++ b/src/crystallization/SMAC.cpp
@@ -0,0 +1,82 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2014 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/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_crystallization_SMAC_h
+#define __PLUMED_crystallization_SMAC_h
+#include "OrientationSphere.h"
+#include "core/ActionRegister.h"
+#include "tools/KernelFunctions.h"
+
+namespace PLMD {
+namespace crystallization {
+
+class SMAC : public OrientationSphere {
+private:
+  std::vector<KernelFunctions> kernels;
+  std::vector<double> deriv;
+  std::vector<Value*> pos;
+public:
+  static void registerKeywords( Keywords& keys ); 
+  SMAC(const ActionOptions& ao); 
+  ~SMAC();
+  double transformDotProduct( const double& dot, double& df ); 
+};
+
+PLUMED_REGISTER_ACTION(SMAC,"SMAC")
+
+void SMAC::registerKeywords( Keywords& keys ){
+  OrientationSphere::registerKeywords(keys);
+  keys.add("numbered","KERNEL","The kernels used in the function of the angle");
+  keys.reset_style("KERNEL","compulsory"); 
+}
+
+SMAC::SMAC(const ActionOptions& ao):
+Action(ao),
+OrientationSphere(ao),
+deriv(1)
+{
+   std::string kernelinpt;
+   for(int i=1;;i++){
+      if( !parseNumbered("KERNEL",i,kernelinpt) ) break;
+      KernelFunctions mykernel( kernelinpt, false );
+      kernels.push_back( mykernel ); 
+   }
+   if( kernels.size()==0 ) error("no kernels defined");
+
+   pos.push_back( new Value() ); 
+   pos[0]->setNotPeriodic();
+}
+
+SMAC::~SMAC(){
+   delete pos[0];
+}
+
+double SMAC::transformDotProduct( const double& dot, double& df ){
+  double ans=0; df=0; pos[0]->set( acos( dot ) ); double dcos=-1./sqrt( 1. - dot*dot );
+  for(unsigned i=0;i<kernels.size();++i){
+      ans += kernels[i].evaluate( pos, deriv );
+      df += deriv[0]*dcos;
+  }
+}
+
+}
+}
+#endif
diff --git a/src/tools/KernelFunctions.cpp b/src/tools/KernelFunctions.cpp
index 6c577858f..68b652827 100644
--- a/src/tools/KernelFunctions.cpp
+++ b/src/tools/KernelFunctions.cpp
@@ -86,6 +86,30 @@ the kernels we can use in this method.
 */
 //+ENDPLUMEDOC
 
+KernelFunctions::KernelFunctions( const std::string& input, const bool& normed ){
+  std::vector<std::string> data=Tools::getWords(input);
+  std::string name=data[0];
+  data.erase(data.begin());
+
+  std::vector<double> at; 
+  bool foundc = Tools::parseVector(data,"CENTER",at);
+  if(!foundc) plumed_merror("failed to find center keyword in definition of kernel");
+  std::vector<double> sig; 
+  bool founds = Tools::parseVector(data,"SIGMA",sig);
+  if(!foundc) plumed_merror("failed to find sigma keyword in definition of kernel");
+
+  bool multi; Tools::parseFlag(data,"MULTIVARIATE",multi);
+  if( center.size()==1 && multi ) plumed_merror("one dimensional kernel cannot be multivariate");
+  if( multi && sig.size()!=0.5*center.size()*(center.size()-1) ) plumed_merror("size mismatch between center size and sigma size");
+  if( !multi && sig.size()!=center.size() ) plumed_merror("size mismatch between center size and sigma size");
+
+  double h;
+  bool foundh = Tools::parse(data,"HEIGHT",h); 
+  if( !foundh) h=1.0;
+
+  KernelFunctions( at, sig, name, multi, h, normed );
+}
+
 KernelFunctions::KernelFunctions( const std::vector<double>& at, const std::vector<double>& sig, const std::string& type, const bool multivariate, const double& w, const bool norm ):
 center(at),
 width(sig)
diff --git a/src/tools/KernelFunctions.h b/src/tools/KernelFunctions.h
index 7e0817598..719601feb 100644
--- a/src/tools/KernelFunctions.h
+++ b/src/tools/KernelFunctions.h
@@ -46,6 +46,7 @@ private:
 /// Get the cutoff for a kernel
   double getCutoff( const double& width ) const ;
 public:
+  KernelFunctions( const std::string& input, const bool& normed );
   KernelFunctions( const std::vector<double>& at, const std::vector<double>& sig, const std::string& type, const bool multivariate ,const double& w, const bool norm );
 /// Get the dimensionality of the kernel
   unsigned ndim() const;
-- 
GitLab