Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "core/ActionWithArguments.h"
      23             : #include "core/ActionWithValue.h"
      24             : #include "core/ActionAtomistic.h"
      25             : #include "core/ActionPilot.h"
      26             : #include "core/ActionRegister.h"
      27             : #include "core/PlumedMain.h"
      28             : #include "tools/OFile.h"
      29             : #include "tools/h36.h"
      30             : 
      31             : namespace PLMD {
      32             : namespace generic {
      33             : 
      34             : //+PLUMEDOC PRINTANALYSIS DUMPPDB
      35             : /*
      36             : Output PDB file.
      37             : 
      38             : 
      39             : \par Examples
      40             : 
      41             : */
      42             : //+ENDPLUMEDOC
      43             : 
      44             : class DumpPDB :
      45             :   public ActionWithArguments,
      46             :   public ActionPilot
      47             : {
      48             :   std::string fmt;
      49             :   std::string file;
      50             :   std::string description;
      51             :   std::vector<unsigned> pdb_resid_indices;
      52             :   std::vector<double> beta, occupancy;
      53             :   std::vector<std::string> argnames;
      54             :   std::vector<AtomNumber> pdb_atom_indices;
      55             :   void buildArgnames();
      56             :   void printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const ;
      57             : public:
      58             :   static void registerKeywords( Keywords& keys );
      59             :   explicit DumpPDB(const ActionOptions&);
      60           7 :   void calculate() override {}
      61           7 :   void apply() override {}
      62             :   void update() override ;
      63             : };
      64             : 
      65             : PLUMED_REGISTER_ACTION(DumpPDB,"DUMPPDB")
      66             : 
      67          21 : void DumpPDB::registerKeywords( Keywords& keys ) {
      68          21 :   Action::registerKeywords( keys );
      69          21 :   ActionPilot::registerKeywords( keys );
      70          21 :   ActionWithArguments::registerKeywords( keys );
      71          42 :   keys.addInputKeyword("optional","ARG","vector/matrix","the values that are being output in the PDB file");
      72          42 :   keys.add("optional","ATOMS","value containing positions of atoms that should be output");
      73          42 :   keys.add("compulsory","STRIDE","0","the frequency with which the atoms should be output");
      74          42 :   keys.add("compulsory","FILE","the name of the file on which to output these quantities");
      75          42 :   keys.add("compulsory","FMT","%f","the format that should be used to output real numbers");
      76          42 :   keys.add("compulsory","OCCUPANCY","1.0","vector of values to output in the occupancy column of the pdb file");
      77          42 :   keys.add("compulsory","BETA","1.0","vector of values to output in the beta column of the pdb file");
      78          42 :   keys.add("optional","DESCRIPTION","the title to use for your PDB output");
      79          42 :   keys.add("optional","ATOM_INDICES","the indices of the atoms in your PDB output");
      80          42 :   keys.add("optional","RESIDUE_INDICES","the indices of the residues in your PDB output");
      81          42 :   keys.add("optional","ARG_NAMES","the names of the arguments that are being output");
      82          21 : }
      83             : 
      84          15 : DumpPDB::DumpPDB(const ActionOptions&ao):
      85             :   Action(ao),
      86             :   ActionWithArguments(ao),
      87          15 :   ActionPilot(ao)
      88             : {
      89          30 :   parse("FILE",file);
      90          15 :   if(file.length()==0) error("name out output file was not specified");
      91          30 :   std::vector<std::string> atoms; parseVector("ATOMS",atoms);
      92          15 :   if( atoms.size()>0 ) {
      93           8 :     std::vector<Value*> atom_args; interpretArgumentList( atoms, plumed.getActionSet(), this, atom_args );
      94           8 :     if( atom_args.size()!=1 ) error("only one action should appear in input to atom args");
      95           8 :     std::vector<Value*> args( getArguments() ); args.push_back( atom_args[0] ); requestArguments( args );
      96          16 :     std::vector<std::string> indices; parseVector("ATOM_INDICES",indices); std::vector<Value*> xvals,yvals,zvals,masv,chargev;
      97           8 :     ActionAtomistic::getAtomValuesFromPlumedObject(plumed,xvals,yvals,zvals,masv,chargev);
      98           8 :     ActionAtomistic::interpretAtomList( indices, xvals, this, pdb_atom_indices );
      99           8 :     log.printf("  printing atoms : ");
     100         153 :     for(unsigned i=0; i<indices.size(); ++i) log.printf("%d ", pdb_atom_indices[i].serial() );
     101           8 :     log.printf("\n");
     102           8 :     if( pdb_atom_indices.size()!=atom_args[0]->getShape()[1]/3 ) error("mismatch between size of matrix containing positions and vector of atom indices");
     103           8 :     beta.resize( atom_args[0]->getShape()[1]/3 ); occupancy.resize( atom_args[0]->getShape()[1]/3 );
     104          16 :     parseVector("OCCUPANCY", occupancy ); parseVector("BETA", beta );
     105          16 :     parseVector("RESIDUE_INDICES",pdb_resid_indices);
     106           8 :     if( pdb_resid_indices.size()==0 ) {
     107           6 :       pdb_resid_indices.resize( pdb_atom_indices.size() );
     108         125 :       for(unsigned i=0; i<pdb_resid_indices.size(); ++i) pdb_resid_indices[i] = pdb_atom_indices[i].serial();
     109           2 :     } else if( pdb_resid_indices.size()!=pdb_atom_indices.size() ) error("mismatch between number of atom indices provided and number of residue indices provided");
     110           8 :   }
     111          15 :   log.printf("  printing configurations to PDB file to file named %s \n", file.c_str() );
     112          30 :   parseVector("ARG_NAMES",argnames); if( argnames.size()==0 ) buildArgnames();
     113          30 :   parse("FMT",fmt); fmt=" "+fmt;
     114          15 :   if( getStride()==0 ) { setStride(0); log.printf("  printing pdb at end of calculation \n"); }
     115             : 
     116          30 :   parse("DESCRIPTION",description);
     117          15 :   if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->hasDerivatives() ) error("argument for printing of PDB should be vector or matrix");
     118          15 :   unsigned nv=getPntrToArgument(0)->getShape()[0];
     119          44 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     120          29 :     if( getPntrToArgument(i)->getRank()==0 || getPntrToArgument(i)->hasDerivatives() ) error("argument for printing of PDB should be vector or matrix");
     121          29 :     if( getPntrToArgument(i)->getShape()[0]!=nv ) error("for printing to pdb files all arguments must have same number of values");
     122             :   }
     123          15 : }
     124             : 
     125        2143 : void DumpPDB::printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const {
     126        2143 :   if( rnum>999 ) plumed_merror("atom number too large to be used as residue number");
     127        2143 :   std::array<char,6> at; const char* msg = h36::hy36encode(5,anum,&at[0]);
     128        2143 :   plumed_assert(msg==nullptr) << msg; at[5]=0; double lunits = plumed.getUnits().getLength()/0.1;
     129        2143 :   opdbf.printf("ATOM  %s  X   RES  %4u    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
     130        2143 :                &at[0], rnum, lunits*pos[0], lunits*pos[1], lunits*pos[2], m, q );
     131        2143 : }
     132             : 
     133          23 : void DumpPDB::buildArgnames() {
     134          23 :   unsigned nargs = getNumberOfArguments(); if( pdb_atom_indices.size()>0 ) nargs = nargs - 1;
     135             :   if( nargs<0 ) return;
     136             : 
     137          23 :   argnames.resize(0); unsigned nvals = getPntrToArgument(0)->getShape()[0];
     138             :   if( getPntrToArgument(0)->getRank()==2 ) nvals = getPntrToArgument(0)->getShape()[0];
     139          48 :   for(unsigned i=0; i<nargs; ++i) {
     140          25 :     if( getPntrToArgument(i)->getShape()[0]!=nvals ) error("all arguments should have same number of values");
     141          25 :     if( getPntrToArgument(i)->getRank()==1 ) {
     142          17 :       argnames.push_back( getPntrToArgument(i)->getName() );
     143           8 :     } else if( getPntrToArgument(i)->getRank()==2 ) {
     144           8 :       (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames );
     145             :     }
     146          25 :     getPntrToArgument(i)->buildDataStore();
     147             :   }
     148             : }
     149             : 
     150          19 : void DumpPDB::update() {
     151          19 :   OFile opdbf; opdbf.link(*this);
     152          19 :   opdbf.setBackupString("analysis");
     153          19 :   opdbf.open( file );
     154          19 :   std::size_t psign=fmt.find("%"); plumed_assert( psign!=std::string::npos );
     155          38 :   std::string descr2="%s=%-" + fmt.substr(psign+1) + " ";
     156             : 
     157          19 :   unsigned totargs = 0; unsigned nvals = getPntrToArgument(0)->getShape()[0];
     158          56 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     159          37 :     if( getPntrToArgument(i)->getShape()[0]!=nvals ) error("all arguments should have same number of values");
     160          37 :     if( getPntrToArgument(i)->getRank()==1 ) totargs += 1;
     161          18 :     else if( getPntrToArgument(i)->getRank()==2 ) totargs += getPntrToArgument(i)->getShape()[1];
     162             :   }
     163          19 :   if( totargs!=argnames.size() ) buildArgnames();
     164             : 
     165          19 :   if( description.length()>0 ) opdbf.printf("# %s AT STEP %lld TIME %f \n", description.c_str(), getStep(), getTime() );
     166          19 :   unsigned nargs = getNumberOfArguments(), atomarg=0;
     167          19 :   if( pdb_atom_indices.size()>0 ) {
     168           9 :     if( nargs>1 ) { atomarg = nargs - 1; nargs = nargs-1; }
     169             :     else nargs = 0;
     170             :   }
     171         448 :   for(unsigned i=0; i<nvals; ++i) {
     172             :     unsigned n=0;
     173        1269 :     for(unsigned j=0; j<nargs; j++) {
     174         840 :       Value* thisarg=getPntrToArgument(j); opdbf.printf("REMARK ");
     175         840 :       if( thisarg->getRank()==1 ) {
     176         405 :         opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i) ); n++;
     177         435 :       } else if( thisarg->getRank()==2 ) {
     178         435 :         unsigned ncols = thisarg->getShape()[1];
     179        1297 :         for(unsigned k=0; k<ncols; ++k) { opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i*ncols+k) ); n++; }
     180             :       }
     181         840 :       opdbf.printf("\n");
     182             :     }
     183         429 :     if( pdb_atom_indices.size()>0 ) {
     184         138 :       unsigned npos = pdb_atom_indices.size(); Vector pos;
     185        2281 :       for(unsigned k=0; k<npos; ++k) {
     186        2143 :         pos[0]=getPntrToArgument(atomarg)->get(npos*(3*i+0) + k);
     187        2143 :         pos[1]=getPntrToArgument(atomarg)->get(npos*(3*i+1) + k);
     188        2143 :         pos[2]=getPntrToArgument(atomarg)->get(npos*(3*i+2) + k);
     189        2143 :         printAtom( opdbf, pdb_atom_indices[k].serial(), pdb_resid_indices[k], pos, occupancy[k], beta[k] );
     190             :       }
     191             :     }
     192         429 :     opdbf.printf("END\n");
     193             :   }
     194          19 :   opdbf.close();
     195          19 : }
     196             : 
     197             : }
     198             : }