From fc6caf6a9ffdcb43c5047711bc0cee8a4b6606aa Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Tue, 26 Jun 2012 12:27:17 +0200
Subject: [PATCH] Added function that calculates the distance from a point in
cv space. The point in CV space can be specified using a pdb file, hence the
changes to Action, Value, etc - all these changes are all so that you can
calculate CVs from a pdb input file. Whilst I was messing about I noticed a
potential problem with natural units and input pdb files. As we are in
natural units there is no way of knowing how to convert the pdb from
angstroms (unit in pdb) to the natural length unit. For these reasons I
assume that when one is in natural units pdb input files are also in natural
units, this required changes in ColvarRMSD and my new ColvarTarget routine as
well as a small change in Atom.h.
---
src/Action.cpp | 12 +++++++-
src/Action.h | 9 ++++++
src/ActionAtomistic.cpp | 11 ++-----
src/ActionAtomistic.h | 2 +-
src/Atoms.h | 6 ++++
src/ColvarRMSD.cpp | 5 ++-
src/ColvarTarget.cpp | 68 +++++++++++++++++++++++++++++++++++++++++
src/PDB.cpp | 3 +-
src/PDB.h | 2 +-
src/TargetDist.cpp | 47 ++++++++++++++++++++++++++++
src/TargetDist.h | 29 ++++++++++++++++++
src/Value.cpp | 6 ++++
src/Value.h | 2 ++
13 files changed, 186 insertions(+), 16 deletions(-)
create mode 100644 src/ColvarTarget.cpp
create mode 100644 src/TargetDist.cpp
create mode 100644 src/TargetDist.h
diff --git a/src/Action.cpp b/src/Action.cpp
index df502c02b..bbd8e6791 100644
--- a/src/Action.cpp
+++ b/src/Action.cpp
@@ -185,7 +185,17 @@ void Action::warning( const std::string & msg ){
log.printf("WARNING for action %s with label %s : %s \n", name.c_str(), label.c_str(), msg.c_str() );
}
-
+void Action::calculateFromPDB( const PDB& pdb ){
+ activate();
+ for(Dependencies::iterator p=after.begin();p!=after.end();++p){
+ ActionWithValue*av=dynamic_cast<ActionWithValue*>(*p);
+ if(av){ av->clearInputForces(); av->clearDerivatives(); }
+ (*p)->readAtomsFromPDB( pdb );
+ (*p)->calculate();
+ }
+ readAtomsFromPDB( pdb );
+ calculate();
+}
diff --git a/src/Action.h b/src/Action.h
index 7919a2fd9..b5432059f 100644
--- a/src/Action.h
+++ b/src/Action.h
@@ -7,6 +7,7 @@
#include "Value.h"
#include "Tools.h"
#include "Log.h"
+#include "PDB.h"
namespace PLMD{
@@ -203,6 +204,14 @@ public:
FILE *fopen(const char *path, const char *mode);
int fclose(FILE*fp);
+
+/// Calculate the action given a pdb file as input. This is used to initialize
+/// things like distance from a point in CV map space given a pdb as an input file
+ void calculateFromPDB( const PDB& pdb );
+/// This is overwritten in ActionAtomistic so that we can read
+/// the atoms from the pdb input file rather than taking them from the
+/// MD code
+ virtual void readAtomsFromPDB( const PDB& pdb ){}
};
/////////////////////
diff --git a/src/ActionAtomistic.cpp b/src/ActionAtomistic.cpp
index ce8a12264..65d59ba47 100644
--- a/src/ActionAtomistic.cpp
+++ b/src/ActionAtomistic.cpp
@@ -174,19 +174,12 @@ void ActionAtomistic::applyForces(){
atoms.forceOnEnergy+=forceOnEnergy;
}
-void ActionAtomistic::readAndCalculate( const PDB& pdb ){
+void ActionAtomistic::readAtomsFromPDB( const PDB& pdb ){
for(unsigned j=0;j<indexes.size();j++){
if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
- if( pdb.getAtomNumbers()[j].index()!=j ) error("there are atoms missing in the pdb file");
+ if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");
positions[j]=pdb.getPositions()[indexes[j].index()];
}
for(unsigned j=0;j<indexes.size();j++) charges[j]=pdb.getBeta()[indexes[j].index()];
for(unsigned j=0;j<indexes.size();j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
- prepare(); calculate();
}
-
-
-
-
-
-
diff --git a/src/ActionAtomistic.h b/src/ActionAtomistic.h
index bf2ef85c2..cc64f87cf 100644
--- a/src/ActionAtomistic.h
+++ b/src/ActionAtomistic.h
@@ -100,7 +100,7 @@ public:
const std::set<AtomNumber> & getUnique()const;
/// Read in an input file containing atom positions and calculate the action for the atomic
/// configuration therin
- void readAndCalculate( const PDB& pdb );
+ void readAtomsFromPDB( const PDB& pdb );
};
inline
diff --git a/src/Atoms.h b/src/Atoms.h
index 1d1aa5a9b..22d1b7456 100644
--- a/src/Atoms.h
+++ b/src/Atoms.h
@@ -152,6 +152,7 @@ public:
void readBinary(std::istream&);
double getKBoltzmann()const;
double getMDKBoltzmann()const;
+ bool usingNaturalUnits()const;
void setNaturalUnits(bool n){naturalUnits=n;};
void setMDNaturalUnits(bool n){MDnaturalUnits=n;};
};
@@ -171,6 +172,11 @@ ActionWithVirtualAtom* Atoms::getVirtualAtomsAction(AtomNumber i)const{
return virtualAtomsActions[i.index()-getNatoms()];
}
+inline
+bool Atoms::usingNaturalUnits() const {
+ return naturalUnits;
+}
+
}
diff --git a/src/ColvarRMSD.cpp b/src/ColvarRMSD.cpp
index 4dc139bd4..cbccf5972 100644
--- a/src/ColvarRMSD.cpp
+++ b/src/ColvarRMSD.cpp
@@ -72,9 +72,8 @@ PLUMED_COLVAR_INIT(ao),rmsd(log)
addValueWithDerivatives(); setNotPeriodic();
PDB pdb;
- // read everything in ang and transform to nm
-
- pdb.read(reference,0.1/atoms.getUnits().length);
+ // read everything in ang and transform to nm if we are not in natural units
+ pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().length);
rmsd.set(pdb,type);
diff --git a/src/ColvarTarget.cpp b/src/ColvarTarget.cpp
new file mode 100644
index 000000000..e417138fc
--- /dev/null
+++ b/src/ColvarTarget.cpp
@@ -0,0 +1,68 @@
+#include "Function.h"
+#include "PlumedMain.h"
+#include "ActionRegister.h"
+#include "PDB.h"
+#include "TargetDist.h"
+#include "Atoms.h"
+
+using namespace std;
+
+namespace PLMD {
+
+//+PLUMEDOC FUNCTION TARGET
+/**
+
+The pythagorean distance from a particular structure measured in the space defined by some
+set of collective variables
+
+*/
+//+ENDPLUMEDOC
+
+class TargetFrame : public Function {
+private:
+ TargetDist target;
+ std::vector<double> derivs;
+public:
+ TargetFrame(const ActionOptions&);
+ virtual void calculate();
+ static void registerKeywords(Keywords& keys );
+};
+
+PLUMED_REGISTER_ACTION(TargetFrame,"TARGET")
+
+void TargetFrame::registerKeywords(Keywords& keys){
+ Function::registerKeywords(keys);
+ keys.use("ARG");
+ keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure");
+ keys.add("optional","REFERENCE_VEC","the vector of values for the CVs at the reference point (if you use this you don't need REFERENCE)");
+}
+
+TargetFrame::TargetFrame(const ActionOptions&ao):
+Action(ao),
+Function(ao),
+target(log)
+{
+ std::vector<double> targ;
+ parseVector("REFERENCE_VEC",targ);
+ if( targ.size()!=0 ){
+ target.read( targ, getArguments() );
+ } else {
+ string reference;
+ parse("REFERENCE",reference);
+ PDB pdb;
+ pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().length);
+ printf("Read pdb file with %d atoms inside\n",pdb.size());
+ target.read( pdb, getArguments() );
+ }
+ checkRead();
+ derivs.resize( getNumberOfArguments() );
+ addValueWithDerivatives(); setNotPeriodic();
+}
+
+void TargetFrame::calculate(){
+ double r=target.calculate( derivs );
+ setValue(r);
+ for(unsigned i=0;i<derivs.size();i++) setDerivative(i,derivs[i]);
+}
+
+}
diff --git a/src/PDB.cpp b/src/PDB.cpp
index 083654c58..a15707a68 100644
--- a/src/PDB.cpp
+++ b/src/PDB.cpp
@@ -27,7 +27,8 @@ unsigned PDB::size()const{
return positions.size();
}
-void PDB::read(const std::string&file,double scale){
+void PDB::read(const std::string&file,bool naturalUnits,double scale){
+ if(naturalUnits) scale=1.0;
FILE* fp=fopen(file.c_str(),"r");
//cerr<<file<<endl;
string line;
diff --git a/src/PDB.h b/src/PDB.h
index 09c652bf4..9bf54cba2 100644
--- a/src/PDB.h
+++ b/src/PDB.h
@@ -18,7 +18,7 @@ class PDB{
std::vector<AtomNumber> numbers;
public:
/// Read the pdb from a file, scaling positions by a factor scale
- void read(const std::string&file,double scale);
+ void read(const std::string&file,bool naturalUnits,double scale);
/// Access to the position array
const std::vector<Vector> & getPositions()const;
/// Access to the occupancy array
diff --git a/src/TargetDist.cpp b/src/TargetDist.cpp
new file mode 100644
index 000000000..cea6f5c37
--- /dev/null
+++ b/src/TargetDist.cpp
@@ -0,0 +1,47 @@
+#include "TargetDist.h"
+
+namespace PLMD {
+
+void TargetDist::read( const PDB& pdb, std::vector<Value*> ar ){
+ // Clear values in target actions
+ for(unsigned i=0;i<ar.size();++i){
+ (ar[i]->getPntrToAction())->clearInputForces();
+ (ar[i]->getPntrToAction())->clearDerivatives();
+ }
+
+ // Caclulate target actions from input in PDB file
+ std::vector<double> targ( ar.size() );
+ for(unsigned i=0;i<ar.size();++i){
+ if( ar[i]->valueHasBeenSet() ){
+ targ[i]=ar[i]->get();
+ } else {
+ ActionWithValue* vv=ar[i]->getPntrToAction();
+ (ar[i]->getPntrToAction())->calculateFromPDB( pdb );
+ targ[i]=ar[i]->get();
+ }
+ }
+ read( targ, ar );
+}
+
+void TargetDist::read( const std::vector<double>& targ, std::vector<Value*> ar ){
+ plumed_assert( targ.size()==ar.size() );
+
+ target.resize( ar.size() ); args.resize( ar.size() );
+ log.printf(" distance from this point in cv space : ");
+ for(unsigned i=0;i<target.size();++i){ log.printf("%f ", targ[i]); target[i]=targ[i]; args[i]=ar[i]; }
+ log.printf("\n");
+}
+
+double TargetDist::calculate( std::vector<double>& derivs ){
+ plumed_assert( derivs.size()==args.size() );
+ double dist=0, tmp;
+ for(unsigned i=0;i<args.size();++i){
+ tmp=args[i]->difference( target[i], args[i]->get() );
+ derivs[i]=tmp; dist+=tmp*tmp;
+ }
+ dist=sqrt(dist);
+ for(unsigned i=0;i<args.size();++i) derivs[i]/=dist;
+ return dist;
+}
+
+}
diff --git a/src/TargetDist.h b/src/TargetDist.h
new file mode 100644
index 000000000..23a931f03
--- /dev/null
+++ b/src/TargetDist.h
@@ -0,0 +1,29 @@
+#ifndef __PLUMED_TargetDist_h
+#define __PLUMED_TargetDist_h
+
+#include "Value.h"
+#include "ActionWithValue.h"
+#include "PDB.h"
+#include <vector>
+#include <string>
+
+namespace PLMD{
+
+class Log;
+class PDB;
+
+class TargetDist {
+private:
+ std::vector<Value*> args;
+ std::vector<double> target;
+ Log &log;
+public:
+ TargetDist(Log& log) : log(log) {};
+ void read( const PDB& pdb, std::vector<Value*> args );
+ void read( const std::vector<double>& targ, std::vector<Value*> ar );
+ double calculate( std::vector<double>& derivs );
+};
+
+}
+
+#endif
diff --git a/src/Value.cpp b/src/Value.cpp
index 9c12ce19e..09fd5fb0c 100644
--- a/src/Value.cpp
+++ b/src/Value.cpp
@@ -103,6 +103,12 @@ double Value::projection(const Value& v1,const Value&v2){
return proj;
}
+ActionWithValue* Value::getPntrToAction(){
+ plumed_assert( action!=NULL );
+ return action;
+}
+
+
diff --git a/src/Value.h b/src/Value.h
index 6dafffffc..55b5e4b46 100644
--- a/src/Value.h
+++ b/src/Value.h
@@ -109,6 +109,8 @@ public:
double difference(double)const;
/// Calculate the difference between two values of this function
double difference(double,double)const;
+/// This returns the pointer to the action where this value is calculated
+ ActionWithValue* getPntrToAction();
/// This sets up the gradients
void setGradients();
static double projection(const Value&,const Value&);
--
GitLab