From d6a289c071d938399ffe959f15839e5ff6ea43ce Mon Sep 17 00:00:00 2001
From: Michele Ceriotti <michele.ceriott@epfl.ch>
Date: Thu, 6 Aug 2015 18:32:59 +0200
Subject: [PATCH] Used MOLINFO to pass atom names to output files in OUTPUT_PDB

---
 src/analysis/OutputPDBFile.cpp           | 15 +++++++++++++--
 src/analysis/ReadDissimilarityMatrix.cpp |  6 ++++--
 src/reference/ReferenceAtoms.cpp         | 23 +++++++++++++++++------
 src/reference/ReferenceAtoms.h           |  2 +-
 src/reference/ReferenceConfiguration.cpp |  7 ++++---
 src/reference/ReferenceConfiguration.h   |  3 ++-
 6 files changed, 41 insertions(+), 15 deletions(-)

diff --git a/src/analysis/OutputPDBFile.cpp b/src/analysis/OutputPDBFile.cpp
index 0783b0fc2..073735835 100644
--- a/src/analysis/OutputPDBFile.cpp
+++ b/src/analysis/OutputPDBFile.cpp
@@ -24,7 +24,9 @@
 #include "reference/ReferenceArguments.h"
 #include "core/ActionRegister.h"
 #include "core/PlumedMain.h"
+#include "core/ActionSet.h"
 #include "core/Atoms.h"
+#include "core/SetupMolInfo.h"
 
 namespace PLMD {
 namespace analysis {
@@ -69,20 +71,29 @@ Action(ao),
 AnalysisWithDataCollection(ao),
 fmt("%f")
 {
+  // Find a moldata object
+  std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
+  if( moldat.empty() ) warning("PDB output files do not have atom types unless you use MOLDATA");
+
   parse("FILE",filename); parse("FMT",fmt);
   if( !getRestart() ){ OFile ofile; ofile.link(*this); ofile.setBackupString("analysis"); ofile.backupAllFiles(filename); }
   log.printf("  printing data to file named %s \n",filename.c_str() );
 }
 
 void OutputPDBFile::performAnalysis(){
+  // Find a moldata object
+  std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
+  if( moldat.size()>1 ) error("you should only have one MOLINFO action in your input file"); 
+  SetupMolInfo* mymoldat=NULL; if( moldat.size()==1 ) mymoldat=moldat[0];
+
   // Output the embedding in plumed pdb format
   OFile afile; afile.link(*this); afile.setBackupString("analysis"); std::size_t psign=fmt.find("%");
   afile.open( filename.c_str() ); std::string descr="REMARK WEIGHT=%-" + fmt.substr(psign+1) + " TYPE=" + getMetricName() + "\n";
   for(unsigned j=0;j<getNumberOfDataPoints();++j){
       afile.printf("DESCRIPTION: analysis data from calculation done at time %f \n",getLabel().c_str(),getTime() );
       afile.printf(descr.c_str(),getWeight(j) ); 
-      if( plumed.getAtoms().usingNaturalUnits() ) getReferenceConfiguration(j,false)->print( 1.0, afile, fmt );
-      else getReferenceConfiguration(j,false)->print( plumed.getAtoms().getUnits().getLength()/0.1, afile, fmt );
+      if( plumed.getAtoms().usingNaturalUnits() ) getReferenceConfiguration(j,false)->print( 1.0, mymoldat, afile, fmt );
+      else getReferenceConfiguration(j,false)->print( plumed.getAtoms().getUnits().getLength()/0.1, mymoldat, afile, fmt );
   }
   afile.close();
 }
diff --git a/src/analysis/ReadDissimilarityMatrix.cpp b/src/analysis/ReadDissimilarityMatrix.cpp
index 393167fc9..717803ac1 100644
--- a/src/analysis/ReadDissimilarityMatrix.cpp
+++ b/src/analysis/ReadDissimilarityMatrix.cpp
@@ -26,6 +26,7 @@
 #include "core/PlumedMain.h"
 #include "core/ActionSet.h"
 #include "core/ActionRegister.h"
+#include "core/ActionSetup.h"
 #include "tools/IFile.h"
 
 //+PLUMEDOC ANALYSIS READ_DISSIMILARITY_MATRIX
@@ -95,8 +96,9 @@ fake_data(NULL)
   } else {
      fake_data=metricRegister().create<ReferenceConfiguration>( "OPTIMAL" );
   }
-
-  if( mytraj.length()>0 && plumed.getActionSet().size()!=1 ) error("should only be this action and the READ_ANALYSIS_FRAMES command in the input file");
+  
+  std::vector<ActionSetup*> setupActions=plumed.getActionSet().select<ActionSetup*>();
+  if( mytraj.length()>0 && (plumed.getActionSet().size()-setupActions.size())!=1 ) error("should only be this action and the READ_ANALYSIS_FRAMES command in the input file");
   if( mytraj.length()==0 && plumed.getActionSet().size()!=0 ) error("read dissimilarity matrix command must be at top of input file");
 
   parse("FILE",fname);
diff --git a/src/reference/ReferenceAtoms.cpp b/src/reference/ReferenceAtoms.cpp
index 48c4fbbad..6a6748222 100644
--- a/src/reference/ReferenceAtoms.cpp
+++ b/src/reference/ReferenceAtoms.cpp
@@ -20,6 +20,7 @@
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 #include "ReferenceAtoms.h"
+#include "core/SetupMolInfo.h"
 #include "tools/OFile.h"
 #include "tools/PDB.h"
 
@@ -41,12 +42,22 @@ void ReferenceAtoms::readAtomsFromPDB( const PDB& pdb ){
   }
 }
 
-void ReferenceAtoms::printAtoms( const double& lunits, OFile& ofile ) const {
-  for(unsigned i=0;i<reference_atoms.size();++i){
-      ofile.printf("ATOM  %4d X    RES   %4u %8.3f %8.3f %8.3f %6.2f %6.2f\n",
-        indices[i].serial(), i, 
-        lunits*reference_atoms[i][0], lunits*reference_atoms[i][1], lunits*reference_atoms[i][2],
-        align[i], displace[i] );
+void ReferenceAtoms::printAtoms( const double& lunits, SetupMolInfo* mymoldat, OFile& ofile ) const {
+  if( !mymoldat ){
+      for(unsigned i=0;i<reference_atoms.size();++i){
+          ofile.printf("ATOM  %4d X    RES   %4u %8.3f %8.3f %8.3f %6.2f %6.2f\n",
+            indices[i].serial(), i, 
+            lunits*reference_atoms[i][0], lunits*reference_atoms[i][1], lunits*reference_atoms[i][2],
+            align[i], displace[i] );
+      }
+  } else {
+      for(unsigned i=0;i<reference_atoms.size();++i){
+          ofile.printf("ATOM  %5d %-4s %3s  %4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
+            indices[i].serial(), mymoldat->getAtomName(indices[i]).c_str(), 
+            mymoldat->getResidueName(indices[i]).c_str(), mymoldat->getResidueNumber(indices[i]),
+            lunits*reference_atoms[i][0], lunits*reference_atoms[i][1], lunits*reference_atoms[i][2],
+            align[i], displace[i] );
+      }
   }
 }
 
diff --git a/src/reference/ReferenceAtoms.h b/src/reference/ReferenceAtoms.h
index 6fc17ced1..a3dde109a 100644
--- a/src/reference/ReferenceAtoms.h
+++ b/src/reference/ReferenceAtoms.h
@@ -88,7 +88,7 @@ public:
 /// Set the positions of the reference atoms
   virtual void setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in )=0;
 /// Print the atomic positions
-  void printAtoms( const double& lunits, OFile& ofile ) const ;
+  void printAtoms( const double& lunits, SetupMolInfo* mymoldat, OFile& ofile ) const ;
 /// Return all atom indexes
   const std::vector<AtomNumber>& getAbsoluteIndexes();
 /// This returns how many atoms there should be
diff --git a/src/reference/ReferenceConfiguration.cpp b/src/reference/ReferenceConfiguration.cpp
index 288aafe15..496580e76 100644
--- a/src/reference/ReferenceConfiguration.cpp
+++ b/src/reference/ReferenceConfiguration.cpp
@@ -25,6 +25,7 @@
 #include "core/Value.h"
 #include "tools/OFile.h"
 #include "tools/PDB.h"
+#include "core/SetupMolInfo.h"
 
 namespace PLMD{
 
@@ -101,10 +102,10 @@ double ReferenceConfiguration::calculate( const std::vector<Vector>& pos, const
 
 void ReferenceConfiguration::print( const double& lunits, OFile& ofile, const double& time, const double& weight, const double& old_norm ){
   ofile.printf("REMARK TIME=%f LOG_WEIGHT=%f OLD_NORM=%f\n",time, weight, old_norm );
-  print( lunits, ofile, "%f" );  // HARD CODED FORMAT HERE AS THIS IS FOR CHECKPOINT FILE
+  print( lunits, NULL, ofile, "%f" );  // HARD CODED FORMAT HERE AS THIS IS FOR CHECKPOINT FILE
 }
 
-void ReferenceConfiguration::print( const double& lunits, OFile& ofile, const std::string& fmt ){
+void ReferenceConfiguration::print( const double& lunits, SetupMolInfo* mymoldat, OFile& ofile, const std::string& fmt ){
   ReferenceArguments* args=dynamic_cast<ReferenceArguments*>(this);
   if( property_names.size()>0 ){
       ofile.printf("REMARK PROPERTIES=%s", property_names[0].c_str() ); 
@@ -123,7 +124,7 @@ void ReferenceConfiguration::print( const double& lunits, OFile& ofile, const st
   }
   if(args) args->printArguments( ofile, fmt );
   ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>(this);
-  if(atoms) atoms->printAtoms( lunits, ofile );
+  if(atoms) atoms->printAtoms( lunits, mymoldat, ofile );
   ofile.printf("END\n");
 }
 
diff --git a/src/reference/ReferenceConfiguration.h b/src/reference/ReferenceConfiguration.h
index 081d38e94..d122f3e64 100644
--- a/src/reference/ReferenceConfiguration.h
+++ b/src/reference/ReferenceConfiguration.h
@@ -37,6 +37,7 @@ class Value;
 class Pbc;
 class OFile;
 class PDB;
+class SetupMolInfo;
 
 /// \ingroup TOOLBOX
 /// Abstract base class for calculating the distance from a reference configuration.
@@ -120,7 +121,7 @@ public:
   void copyDerivatives( const ReferenceConfiguration* );
 /// Print a pdb file containing the reference configuration
   void print( const double& lunits, OFile& ofile, const double& time, const double& weight, const double& old_norm );
-  void print( const double& lunits, OFile& ofile, const std::string& fmt );
+  void print( const double& lunits, SetupMolInfo* mymoldat, OFile& ofile, const std::string& fmt );
 /// Get one of the referene arguments
   virtual double getReferenceArgument( const unsigned& i ) const { plumed_error(); return 0.0; }
 /// These are overwritten in ReferenceArguments and ReferenceAtoms but are required here 
-- 
GitLab