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 ); keys.use("ARG");
71 42 : keys.add("optional","ATOMS","value containing positions of atoms that should be output");
72 42 : keys.add("compulsory","STRIDE","0","the frequency with which the atoms should be output");
73 42 : keys.add("compulsory","FILE","the name of the file on which to output these quantities");
74 42 : keys.add("compulsory","FMT","%f","the format that should be used to output real numbers");
75 42 : keys.add("compulsory","OCCUPANCY","1.0","vector of values to output in the occupancy column of the pdb file");
76 42 : keys.add("compulsory","BETA","1.0","vector of values to output in the beta column of the pdb file");
77 42 : keys.add("optional","DESCRIPTION","the title to use for your PDB output");
78 42 : keys.add("optional","ATOM_INDICES","the indices of the atoms in your PDB output");
79 42 : keys.add("optional","RESIDUE_INDICES","the indices of the residues in your PDB output");
80 42 : keys.add("optional","ARG_NAMES","the names of the arguments that are being output");
81 21 : }
82 :
83 15 : DumpPDB::DumpPDB(const ActionOptions&ao):
84 : Action(ao),
85 : ActionWithArguments(ao),
86 15 : ActionPilot(ao)
87 : {
88 30 : parse("FILE",file);
89 15 : if(file.length()==0) error("name out output file was not specified");
90 30 : std::vector<std::string> atoms; parseVector("ATOMS",atoms);
91 15 : if( atoms.size()>0 ) {
92 8 : std::vector<Value*> atom_args; interpretArgumentList( atoms, plumed.getActionSet(), this, atom_args );
93 8 : if( atom_args.size()!=1 ) error("only one action should appear in input to atom args");
94 8 : std::vector<Value*> args( getArguments() ); args.push_back( atom_args[0] ); requestArguments( args );
95 16 : std::vector<std::string> indices; parseVector("ATOM_INDICES",indices); std::vector<Value*> xvals,yvals,zvals,masv,chargev;
96 8 : ActionAtomistic::getAtomValuesFromPlumedObject(plumed,xvals,yvals,zvals,masv,chargev);
97 8 : ActionAtomistic::interpretAtomList( indices, xvals, this, pdb_atom_indices );
98 8 : log.printf(" printing atoms : ");
99 153 : for(unsigned i=0; i<indices.size(); ++i) log.printf("%d ", pdb_atom_indices[i].serial() );
100 8 : log.printf("\n");
101 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");
102 8 : beta.resize( atom_args[0]->getShape()[1]/3 ); occupancy.resize( atom_args[0]->getShape()[1]/3 );
103 16 : parseVector("OCCUPANCY", occupancy ); parseVector("BETA", beta );
104 16 : parseVector("RESIDUE_INDICES",pdb_resid_indices);
105 8 : if( pdb_resid_indices.size()==0 ) {
106 6 : pdb_resid_indices.resize( pdb_atom_indices.size() );
107 125 : for(unsigned i=0; i<pdb_resid_indices.size(); ++i) pdb_resid_indices[i] = pdb_atom_indices[i].serial();
108 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");
109 8 : }
110 15 : log.printf(" printing configurations to PDB file to file named %s \n", file.c_str() );
111 30 : parseVector("ARG_NAMES",argnames); if( argnames.size()==0 ) buildArgnames();
112 30 : parse("FMT",fmt); fmt=" "+fmt;
113 15 : if( getStride()==0 ) { setStride(0); log.printf(" printing pdb at end of calculation \n"); }
114 :
115 30 : parse("DESCRIPTION",description);
116 15 : if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->hasDerivatives() ) error("argument for printing of PDB should be vector or matrix");
117 15 : unsigned nv=getPntrToArgument(0)->getShape()[0];
118 44 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
119 29 : if( getPntrToArgument(i)->getRank()==0 || getPntrToArgument(i)->hasDerivatives() ) error("argument for printing of PDB should be vector or matrix");
120 29 : if( getPntrToArgument(i)->getShape()[0]!=nv ) error("for printing to pdb files all arguments must have same number of values");
121 : }
122 15 : }
123 :
124 2143 : void DumpPDB::printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const {
125 2143 : if( rnum>999 ) plumed_merror("atom number too large to be used as residue number");
126 2143 : std::array<char,6> at; const char* msg = h36::hy36encode(5,anum,&at[0]);
127 2143 : plumed_assert(msg==nullptr) << msg; at[5]=0; double lunits = plumed.getUnits().getLength()/0.1;
128 2143 : opdbf.printf("ATOM %s X RES %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
129 2143 : &at[0], rnum, lunits*pos[0], lunits*pos[1], lunits*pos[2], m, q );
130 2143 : }
131 :
132 23 : void DumpPDB::buildArgnames() {
133 23 : unsigned nargs = getNumberOfArguments(); if( pdb_atom_indices.size()>0 ) nargs = nargs - 1;
134 : if( nargs<0 ) return;
135 :
136 23 : argnames.resize(0); unsigned nvals = getPntrToArgument(0)->getShape()[0];
137 : if( getPntrToArgument(0)->getRank()==2 ) nvals = getPntrToArgument(0)->getShape()[0];
138 48 : for(unsigned i=0; i<nargs; ++i) {
139 25 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) error("all arguments should have same number of values");
140 25 : if( getPntrToArgument(i)->getRank()==1 ) {
141 17 : argnames.push_back( getPntrToArgument(i)->getName() );
142 8 : } else if( getPntrToArgument(i)->getRank()==2 ) {
143 8 : (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames );
144 : }
145 25 : getPntrToArgument(i)->buildDataStore();
146 : }
147 : }
148 :
149 19 : void DumpPDB::update() {
150 19 : OFile opdbf; opdbf.link(*this);
151 19 : opdbf.setBackupString("analysis");
152 19 : opdbf.open( file );
153 19 : std::size_t psign=fmt.find("%"); plumed_assert( psign!=std::string::npos );
154 38 : std::string descr2="%s=%-" + fmt.substr(psign+1) + " ";
155 :
156 19 : unsigned totargs = 0; unsigned nvals = getPntrToArgument(0)->getShape()[0];
157 56 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
158 37 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) error("all arguments should have same number of values");
159 37 : if( getPntrToArgument(i)->getRank()==1 ) totargs += 1;
160 18 : else if( getPntrToArgument(i)->getRank()==2 ) totargs += getPntrToArgument(i)->getShape()[1];
161 : }
162 19 : if( totargs!=argnames.size() ) buildArgnames();
163 :
164 19 : if( description.length()>0 ) opdbf.printf("# %s AT STEP %lld TIME %f \n", description.c_str(), getStep(), getTime() );
165 19 : unsigned nargs = getNumberOfArguments(), atomarg=0;
166 19 : if( pdb_atom_indices.size()>0 ) {
167 9 : if( nargs>1 ) { atomarg = nargs - 1; nargs = nargs-1; }
168 : else nargs = 0;
169 : }
170 448 : for(unsigned i=0; i<nvals; ++i) {
171 : unsigned n=0;
172 1269 : for(unsigned j=0; j<nargs; j++) {
173 840 : Value* thisarg=getPntrToArgument(j); opdbf.printf("REMARK ");
174 840 : if( thisarg->getRank()==1 ) {
175 405 : opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i) ); n++;
176 435 : } else if( thisarg->getRank()==2 ) {
177 435 : unsigned ncols = thisarg->getShape()[1];
178 1297 : for(unsigned k=0; k<ncols; ++k) { opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i*ncols+k) ); n++; }
179 : }
180 840 : opdbf.printf("\n");
181 : }
182 429 : if( pdb_atom_indices.size()>0 ) {
183 138 : unsigned npos = pdb_atom_indices.size(); Vector pos;
184 2281 : for(unsigned k=0; k<npos; ++k) {
185 2143 : pos[0]=getPntrToArgument(atomarg)->get(npos*(3*i+0) + k);
186 2143 : pos[1]=getPntrToArgument(atomarg)->get(npos*(3*i+1) + k);
187 2143 : pos[2]=getPntrToArgument(atomarg)->get(npos*(3*i+2) + k);
188 2143 : printAtom( opdbf, pdb_atom_indices[k].serial(), pdb_resid_indices[k], pos, occupancy[k], beta[k] );
189 : }
190 : }
191 429 : opdbf.printf("END\n");
192 : }
193 19 : opdbf.close();
194 19 : }
195 :
196 : }
197 : }
|
