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 : std::string fmt;
48 : std::string file;
49 : std::string description;
50 : std::vector<unsigned> pdb_resid_indices;
51 : std::vector<double> beta, occupancy;
52 : std::vector<std::string> argnames;
53 : std::vector<AtomNumber> pdb_atom_indices;
54 : void buildArgnames();
55 : void printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const ;
56 : public:
57 : static void registerKeywords( Keywords& keys );
58 : explicit DumpPDB(const ActionOptions&);
59 7 : void calculate() override {}
60 7 : void apply() override {}
61 : void update() override ;
62 : };
63 :
64 : PLUMED_REGISTER_ACTION(DumpPDB,"DUMPPDB")
65 :
66 23 : void DumpPDB::registerKeywords( Keywords& keys ) {
67 23 : Action::registerKeywords( keys );
68 23 : ActionPilot::registerKeywords( keys );
69 23 : ActionWithArguments::registerKeywords( keys );
70 23 : keys.use("ARG");
71 46 : keys.add("optional","ATOMS","value containing positions of atoms that should be output");
72 46 : keys.add("compulsory","STRIDE","0","the frequency with which the atoms should be output");
73 46 : keys.add("compulsory","FILE","the name of the file on which to output these quantities");
74 46 : keys.add("compulsory","FMT","%f","the format that should be used to output real numbers");
75 46 : keys.add("compulsory","OCCUPANCY","1.0","vector of values to output in the occupancy column of the pdb file");
76 46 : keys.add("compulsory","BETA","1.0","vector of values to output in the beta column of the pdb file");
77 46 : keys.add("optional","DESCRIPTION","the title to use for your PDB output");
78 46 : keys.add("optional","ATOM_INDICES","the indices of the atoms in your PDB output");
79 46 : keys.add("optional","RESIDUE_INDICES","the indices of the residues in your PDB output");
80 46 : keys.add("optional","ARG_NAMES","the names of the arguments that are being output");
81 23 : }
82 :
83 15 : DumpPDB::DumpPDB(const ActionOptions&ao):
84 : Action(ao),
85 : ActionWithArguments(ao),
86 15 : ActionPilot(ao) {
87 30 : parse("FILE",file);
88 15 : if(file.length()==0) {
89 0 : error("name out output file was not specified");
90 : }
91 : std::vector<std::string> atoms;
92 30 : parseVector("ATOMS",atoms);
93 15 : if( atoms.size()>0 ) {
94 : std::vector<Value*> atom_args;
95 8 : interpretArgumentList( atoms, plumed.getActionSet(), this, atom_args );
96 8 : if( atom_args.size()!=1 ) {
97 0 : error("only one action should appear in input to atom args");
98 : }
99 8 : std::vector<Value*> args( getArguments() );
100 8 : args.push_back( atom_args[0] );
101 8 : requestArguments( args );
102 : std::vector<std::string> indices;
103 16 : parseVector("ATOM_INDICES",indices);
104 : std::vector<Value*> xvals,yvals,zvals,masv,chargev;
105 8 : ActionAtomistic::getAtomValuesFromPlumedObject(plumed,xvals,yvals,zvals,masv,chargev);
106 8 : ActionAtomistic::interpretAtomList( indices, xvals, this, pdb_atom_indices );
107 8 : log.printf(" printing atoms : ");
108 153 : for(unsigned i=0; i<indices.size(); ++i) {
109 145 : log.printf("%d ", pdb_atom_indices[i].serial() );
110 : }
111 8 : log.printf("\n");
112 8 : if( pdb_atom_indices.size()!=atom_args[0]->getShape()[1]/3 ) {
113 0 : error("mismatch between size of matrix containing positions and vector of atom indices");
114 : }
115 8 : beta.resize( atom_args[0]->getShape()[1]/3 );
116 8 : occupancy.resize( atom_args[0]->getShape()[1]/3 );
117 8 : parseVector("OCCUPANCY", occupancy );
118 8 : parseVector("BETA", beta );
119 16 : parseVector("RESIDUE_INDICES",pdb_resid_indices);
120 8 : if( pdb_resid_indices.size()==0 ) {
121 5 : pdb_resid_indices.resize( pdb_atom_indices.size() );
122 111 : for(unsigned i=0; i<pdb_resid_indices.size(); ++i) {
123 106 : pdb_resid_indices[i] = pdb_atom_indices[i].serial();
124 : }
125 3 : } else if( pdb_resid_indices.size()!=pdb_atom_indices.size() ) {
126 0 : error("mismatch between number of atom indices provided and number of residue indices provided");
127 : }
128 8 : }
129 15 : log.printf(" printing configurations to PDB file to file named %s \n", file.c_str() );
130 30 : parseVector("ARG_NAMES",argnames);
131 15 : if( argnames.size()==0 ) {
132 14 : buildArgnames();
133 : }
134 15 : parse("FMT",fmt);
135 15 : fmt=" "+fmt;
136 15 : if( getStride()==0 ) {
137 8 : setStride(0);
138 8 : log.printf(" printing pdb at end of calculation \n");
139 : }
140 :
141 30 : parse("DESCRIPTION",description);
142 15 : if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->hasDerivatives() ) {
143 0 : error("argument for printing of PDB should be vector or matrix");
144 : }
145 15 : unsigned nv=getPntrToArgument(0)->getShape()[0];
146 44 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
147 29 : if( getPntrToArgument(i)->getRank()==0 || getPntrToArgument(i)->hasDerivatives() ) {
148 0 : error("argument for printing of PDB should be vector or matrix");
149 : }
150 29 : if( getPntrToArgument(i)->getShape()[0]!=nv ) {
151 0 : error("for printing to pdb files all arguments must have same number of values");
152 : }
153 : }
154 15 : }
155 :
156 2143 : void DumpPDB::printAtom( OFile& opdbf, const unsigned& anum, const unsigned& rnum, const Vector& pos, const double& m, const double& q ) const {
157 2143 : if( rnum>999 ) {
158 0 : plumed_merror("atom number too large to be used as residue number");
159 : }
160 : std::array<char,6> at;
161 2143 : const char* msg = h36::hy36encode(5,anum,&at[0]);
162 2143 : plumed_assert(msg==nullptr) << msg;
163 2143 : at[5]=0;
164 2143 : double lunits = plumed.getUnits().getLength()/0.1;
165 2143 : opdbf.printf("ATOM %s X RES %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
166 2143 : &at[0], rnum, lunits*pos[0], lunits*pos[1], lunits*pos[2], m, q );
167 2143 : }
168 :
169 23 : void DumpPDB::buildArgnames() {
170 23 : unsigned nargs = getNumberOfArguments();
171 23 : if( pdb_atom_indices.size()>0 ) {
172 17 : nargs = nargs - 1;
173 : }
174 : if( nargs<0 ) {
175 : return;
176 : }
177 :
178 23 : argnames.resize(0);
179 23 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
180 : if( getPntrToArgument(0)->getRank()==2 ) {
181 : nvals = getPntrToArgument(0)->getShape()[0];
182 : }
183 48 : for(unsigned i=0; i<nargs; ++i) {
184 25 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
185 0 : error("all arguments should have same number of values");
186 : }
187 25 : if( getPntrToArgument(i)->getRank()==1 ) {
188 17 : argnames.push_back( getPntrToArgument(i)->getName() );
189 8 : } else if( getPntrToArgument(i)->getRank()==2 ) {
190 8 : (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames );
191 : }
192 25 : getPntrToArgument(i)->buildDataStore();
193 : }
194 : }
195 :
196 19 : void DumpPDB::update() {
197 19 : OFile opdbf;
198 19 : opdbf.link(*this);
199 19 : opdbf.setBackupString("analysis");
200 19 : opdbf.open( file );
201 19 : std::size_t psign=fmt.find("%");
202 19 : plumed_assert( psign!=std::string::npos );
203 38 : std::string descr2="%s=%-" + fmt.substr(psign+1) + " ";
204 :
205 : unsigned totargs = 0;
206 19 : unsigned nvals = getPntrToArgument(0)->getShape()[0];
207 56 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
208 37 : if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
209 0 : error("all arguments should have same number of values");
210 : }
211 37 : if( getPntrToArgument(i)->getRank()==1 ) {
212 19 : totargs += 1;
213 18 : } else if( getPntrToArgument(i)->getRank()==2 ) {
214 18 : totargs += getPntrToArgument(i)->getShape()[1];
215 : }
216 : }
217 19 : if( totargs!=argnames.size() ) {
218 9 : buildArgnames();
219 : }
220 :
221 19 : if( description.length()>0 ) {
222 11 : opdbf.printf("# %s AT STEP %lld TIME %f \n", description.c_str(), getStep(), getTime() );
223 : }
224 19 : unsigned nargs = getNumberOfArguments(), atomarg=0;
225 19 : if( pdb_atom_indices.size()>0 ) {
226 9 : if( nargs>1 ) {
227 3 : atomarg = nargs - 1;
228 : nargs = nargs-1;
229 : } else {
230 : nargs = 0;
231 : }
232 : }
233 448 : for(unsigned i=0; i<nvals; ++i) {
234 : unsigned n=0;
235 1269 : for(unsigned j=0; j<nargs; j++) {
236 : Value* thisarg=getPntrToArgument(j);
237 840 : opdbf.printf("REMARK ");
238 840 : if( thisarg->getRank()==1 ) {
239 405 : opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i) );
240 405 : n++;
241 435 : } else if( thisarg->getRank()==2 ) {
242 435 : unsigned ncols = thisarg->getShape()[1];
243 1297 : for(unsigned k=0; k<ncols; ++k) {
244 862 : opdbf.printf( descr2.c_str(), argnames[n].c_str(), thisarg->get(i*ncols+k) );
245 862 : n++;
246 : }
247 : }
248 840 : opdbf.printf("\n");
249 : }
250 429 : if( pdb_atom_indices.size()>0 ) {
251 138 : unsigned npos = pdb_atom_indices.size();
252 138 : Vector pos;
253 2281 : for(unsigned k=0; k<npos; ++k) {
254 2143 : pos[0]=getPntrToArgument(atomarg)->get(npos*(3*i+0) + k);
255 2143 : pos[1]=getPntrToArgument(atomarg)->get(npos*(3*i+1) + k);
256 2143 : pos[2]=getPntrToArgument(atomarg)->get(npos*(3*i+2) + k);
257 2143 : printAtom( opdbf, pdb_atom_indices[k].serial(), pdb_resid_indices[k], pos, occupancy[k], beta[k] );
258 : }
259 : }
260 429 : opdbf.printf("END\n");
261 : }
262 19 : opdbf.close();
263 19 : }
264 :
265 : }
266 : }
|
