From 903b98b0af1b3eb64d74f95650b7901c1bb60a4a Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Wed, 25 Mar 2015 11:19:41 +0000
Subject: [PATCH] Added new metric and a new multicolvar
The new multicolvar measures the inplanedistance i.e. the
distane from a plane defined by a vector connecting two atoms.
The new refernce distance measures the dot product distance
between two points
---
src/multicolvar/InPlaneDistances.cpp | 133 +++++++++++++++++++++++++++
src/reference/DotProductDistance.cpp | 54 +++++++++++
2 files changed, 187 insertions(+)
create mode 100644 src/multicolvar/InPlaneDistances.cpp
create mode 100644 src/reference/DotProductDistance.cpp
diff --git a/src/multicolvar/InPlaneDistances.cpp b/src/multicolvar/InPlaneDistances.cpp
new file mode 100644
index 000000000..bf5ef3625
--- /dev/null
+++ b/src/multicolvar/InPlaneDistances.cpp
@@ -0,0 +1,133 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013,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/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "MultiColvar.h"
+#include "tools/SwitchingFunction.h"
+#include "core/ActionRegister.h"
+#include "vesselbase/LessThan.h"
+#include "vesselbase/Between.h"
+#include "tools/Angle.h"
+
+#include <string>
+#include <cmath>
+
+using namespace std;
+
+namespace PLMD{
+namespace multicolvar{
+
+//+PLUMEDOC MCOLVAR INPLANEDISTANCES
+/*
+Calculate distances in the plane perpendicular to an axis
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+class InPlaneDistances : public MultiColvar {
+public:
+ static void registerKeywords( Keywords& keys );
+ InPlaneDistances(const ActionOptions&);
+// active methods:
+ virtual double compute();
+ bool isPeriodic(){ return false; }
+ Vector getCentralAtom();
+};
+
+PLUMED_REGISTER_ACTION(InPlaneDistances,"INPLANEDISTANCES")
+
+void InPlaneDistances::registerKeywords( Keywords& keys ){
+ MultiColvar::registerKeywords( keys );
+ keys.use("MEAN"); keys.use("MIN"); keys.use("MAX"); keys.use("LESS_THAN");
+ keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
+ keys.add("atoms","VECTORSTART","The first atom position that is used to define the normal to the plane of interest");
+ keys.add("atoms","VECTOREND","The second atom position that is used to defin the normal to the plane of interest");
+ keys.add("atoms-2","GROUP","The set of atoms for which you wish to calculate the in plane distance ");
+}
+
+InPlaneDistances::InPlaneDistances(const ActionOptions&ao):
+PLUMED_MULTICOLVAR_INIT(ao)
+{
+ // Read in the atoms
+ std::vector<AtomNumber> all_atoms;
+ readThreeGroups("GROUP","VECTORSTART","VECTOREND",false, all_atoms);
+ if( all_atoms.size()>0 ) ActionAtomistic::requestAtoms( all_atoms );
+
+ // Check atoms are OK
+ if( getFullNumberOfTasks()!=getNumberOfAtoms()-2 ) error("you should specify one atom for VECTORSTART and one atom for VECTOREND only");
+
+ // Setup the multicolvar base
+ setupMultiColvarBase(); readVesselKeywords();
+ // And check everything has been read in correctly
+ checkRead();
+
+ // Now check if we can use link cells
+ bool use_link=false; double rcut;
+ if( getNumberOfVessels()>0 ){
+ vesselbase::LessThan* lt=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(0) );
+ if( lt ){
+ use_link=true; rcut=lt->getCutoff();
+ } else {
+ vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(0) );
+ if( bt ) use_link=true; rcut=bt->getCutoff();
+ }
+ if( use_link ){
+ for(unsigned i=1;i<getNumberOfVessels();++i){
+ vesselbase::LessThan* lt2=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(i) );
+ vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(i) );
+ if( lt2 ){
+ double tcut=lt2->getCutoff();
+ if( tcut>rcut ) rcut=tcut;
+ } else if( bt ){
+ double tcut=bt->getCutoff();
+ if( tcut>rcut ) rcut=tcut;
+ } else {
+ use_link=false;
+ }
+ }
+ }
+ if( use_link ) setLinkCellCutoff( rcut );
+ }
+}
+
+double InPlaneDistances::compute(){
+ Vector normal=getSeparation( getPosition(1), getPosition(2) );
+ Vector dir=getSeparation( getPosition(1), getPosition(0) );
+ PLMD::Angle a; Vector ddij, ddik; double angle=a.compute(normal,dir,ddij,ddik);
+ double sangle=sin(angle), cangle=cos(angle);
+ double dd=dir.modulo(), invdd=1.0/dd, val=dd*sangle;
+
+ addAtomsDerivatives( 0, dd*cangle*ddik + sangle*invdd*dir );
+ addAtomsDerivatives( 1, -dd*cangle*(ddik+ddij) - sangle*invdd*dir );
+ addAtomsDerivatives( 2, dd*cangle*ddij );
+ addBoxDerivatives( -dd*cangle*(Tensor(normal,ddij)+Tensor(dir,ddik)) - sangle*invdd*Tensor(dir,dir) );
+
+ return val;
+}
+
+Vector InPlaneDistances::getCentralAtom(){
+ addCentralAtomDerivatives( 0, Tensor::identity() );
+ return getPosition(0);
+}
+
+}
+}
diff --git a/src/reference/DotProductDistance.cpp b/src/reference/DotProductDistance.cpp
new file mode 100644
index 000000000..7dd4d8f86
--- /dev/null
+++ b/src/reference/DotProductDistance.cpp
@@ -0,0 +1,54 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-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/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "MetricRegister.h"
+#include "ArgumentOnlyDistance.h"
+#include "core/Value.h"
+
+namespace PLMD {
+
+class DotProductDistance : public ArgumentOnlyDistance {
+public:
+ DotProductDistance( const ReferenceConfigurationOptions& ro );
+ void read( const PDB& );
+ double calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared );
+};
+
+PLUMED_REGISTER_METRIC(DotProductDistance,"DOTPRODUCT")
+
+DotProductDistance::DotProductDistance( const ReferenceConfigurationOptions& ro ):
+ReferenceConfiguration(ro),
+ArgumentOnlyDistance(ro)
+{
+}
+
+void DotProductDistance::read( const PDB& pdb ){
+ readArgumentsFromPDB( pdb );
+}
+
+double DotProductDistance::calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
+ double dot=0.0;
+ for (unsigned long i=0; i<vals.size(); ++i) dot+=getReferenceArgument(i)*arg[i];
+ return -log(dot);
+}
+
+
+}
--
GitLab