Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "GenericMolInfo.h"
23 : #include "ActionRegister.h"
24 : #include "ActionSet.h"
25 : #include "PlumedMain.h"
26 : #include "tools/MolDataClass.h"
27 : #include "tools/PDB.h"
28 : #include "tools/Communicator.h"
29 : #include "tools/Pbc.h"
30 : #include "tools/Tools.h"
31 : #include "config/Config.h"
32 : #include <cmath>
33 : #include <algorithm>
34 : #include <map>
35 :
36 : namespace PLMD {
37 :
38 :
39 : /*
40 : This action is defined in core/ as it is used by other actions.
41 : Anyway, it is registered in setup/, so that excluding that module from
42 : compilation will exclude it from plumed.
43 : */
44 :
45 :
46 168 : void GenericMolInfo::registerKeywords( Keywords& keys ) {
47 168 : ActionSetup::registerKeywords(keys);
48 168 : keys.add("compulsory","STRUCTURE","a file in pdb format containing a reference structure. "
49 : "This is used to defines the atoms in the various residues, chains, etc . "
50 : "For more details on the PDB file format visit http://www.wwpdb.org/docs.html");
51 168 : keys.add("compulsory","MOLTYPE","protein","what kind of molecule is contained in the pdb file - usually not needed since protein/RNA/DNA are compatible");
52 168 : keys.add("compulsory","PYTHON_BIN","default","python interpreter");
53 168 : keys.add("atoms","CHAIN","(for masochists ( mostly Davide Branduardi ) ) The atoms involved in each of the chains of interest in the structure.");
54 168 : keys.add("hidden","STRIDE","frequency for resetting the python interpreter. Should be 1.");
55 168 : keys.addFlag("WHOLE", false, "The reference structure is whole, i.e. not broken by PBC. If the pdb file has a CRYST1 record with space group P 1, this is checked along each chain, and a warning is written to the log if consecutive atoms of a residue, or of two consecutive protein or nucleic-acid residues, are on opposite sides of the periodic boundary");
56 168 : }
57 :
58 328 : GenericMolInfo::~GenericMolInfo() {
59 : // empty destructor to delete unique_ptr
60 492 : }
61 :
62 : namespace {
63 : // MOLINFO WHOLE takes the reference structure as whole. Check it along each chain in the box of the
64 : // CRYST1 record: consecutive atoms of a residue, or of two consecutive protein or nucleic-acid residues,
65 : // are close, so a step between them that is not its own minimum image crosses the periodic boundary.
66 : // Returns a warning, or an empty string.
67 16 : std::string checkWhole(const PDB & pdb,double scale) {
68 : // skip a cell that cannot be the periodic box: none or degenerate, the 1 A cube written when there
69 : // is no unit cell, or a crystallographic cell (simulation programs write the space group P 1)
70 16 : const Tensor & box(pdb.getBoxVec());
71 16 : const double det=box.determinant();
72 : bool unitcube=true;
73 64 : for(unsigned k=0; k<3; k++) {
74 48 : unitcube=unitcube && std::fabs(pdb.getBoxAxs()[k]-scale)<1e-3*scale && std::fabs(pdb.getBoxAng()[k]-90.0)<1e-3;
75 : }
76 16 : std::string group=pdb.getSpaceGroup();
77 16 : group.erase(std::remove(group.begin(),group.end(),' '),group.end());
78 16 : if(unitcube || (group.length()>0 && group!="P1") || !std::isfinite(det) || det<=0.0) {
79 7 : return "";
80 : }
81 9 : Pbc pbc;
82 9 : pbc.setBox(box);
83 268 : auto polymer=[](std::string r) {
84 268 : Tools::ltrim(r);
85 268 : Tools::trim(r);
86 380 : return MolDataClass::allowedResidue("protein",r) || MolDataClass::allowedResidue("dna",r) || MolDataClass::allowedResidue("rna",r);
87 : };
88 9 : const auto & at(pdb.getAtomNumbers());
89 9 : const auto & pos(pdb.getPositions());
90 : unsigned n=0;
91 : std::string list;
92 2240 : for(unsigned i=1; i<at.size(); i++) {
93 2231 : if(pdb.getChainID(at[i])!=pdb.getChainID(at[i-1])) {
94 2 : continue;
95 : }
96 2229 : if(pdb.getResidueNumber(at[i])!=pdb.getResidueNumber(at[i-1])
97 2497 : && !(polymer(pdb.getResidueName(at[i-1])) && polymer(pdb.getResidueName(at[i])))) {
98 34 : continue;
99 : }
100 : // the difference is a lattice vector, or zero up to rounding
101 2195 : if((delta(pos[i-1],pos[i])-pbc.distance(pos[i-1],pos[i])).modulo()>1e-3*std::cbrt(det)) {
102 15 : if(++n<=10) {
103 11 : std::string r=pdb.getResidueName(at[i]);
104 11 : Tools::ltrim(r);
105 11 : Tools::trim(r);
106 33 : list+=" "+std::to_string(at[i].serial())+" ("+r+std::to_string(pdb.getResidueNumber(at[i]))+")";
107 : }
108 : }
109 : }
110 9 : if(n==0) {
111 7 : return "";
112 : }
113 : return "the reference structure in the MOLINFO pdb file may not be whole: in the box of its CRYST1 record, "
114 6 : +std::to_string(n)+" atoms are across the periodic boundary from the atom before them in their chain:"+list
115 7 : +(n>10 ? " and "+std::to_string(n-10)+" more" : "")
116 2 : +". WHOLEMOLECULES will then split these molecules in every frame: make the pdb whole"
117 : +" (for example with gmx trjconv -pbc whole), or check that its CRYST1 record is the box the structure was written in";
118 : }
119 : }
120 :
121 164 : GenericMolInfo::GenericMolInfo( const ActionOptions&ao ):
122 : Action(ao),
123 : ActionAnyorder(ao),
124 : ActionPilot(ao),
125 : ActionAtomistic(ao),
126 164 : iswhole_(false) {
127 164 : plumed_assert(getStride()==1);
128 : // Read what is contained in the pdb file
129 164 : parse("MOLTYPE",mytype);
130 :
131 : // check if whole
132 164 : parseFlag("WHOLE", iswhole_);
133 :
134 164 : auto* latestmoldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
135 164 : if( latestmoldat ) {
136 4 : log<<" overriding last MOLINFO with label " << latestmoldat->getLabel()<<"\n";
137 : }
138 :
139 164 : std::vector<AtomNumber> backbone;
140 328 : parseAtomList("CHAIN",backbone);
141 164 : if( read_backbone.size()==0 ) {
142 0 : for(unsigned i=1;; ++i) {
143 328 : parseAtomList("CHAIN",i,backbone);
144 164 : if( backbone.size()==0 ) {
145 : break;
146 : }
147 0 : read_backbone.push_back(backbone);
148 0 : backbone.resize(0);
149 : }
150 : } else {
151 0 : read_backbone.push_back(backbone);
152 : }
153 164 : if( read_backbone.size()==0 ) {
154 164 : parse("STRUCTURE",reference);
155 :
156 164 : if( ! pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength())) {
157 0 : plumed_merror("missing input file " + reference );
158 : }
159 164 : if( iswhole_ ) {
160 16 : const std::string msg=checkWhole(pdb,usingNaturalUnits() ? 1.0 : 0.1/getUnits().getLength());
161 16 : if( msg.length()>0 ) {
162 2 : warning(msg);
163 : }
164 : }
165 :
166 164 : std::vector<std::string> chains;
167 164 : pdb.getChainNames( chains );
168 164 : log.printf(" pdb file named %s contains %u chains \n",reference.c_str(), static_cast<unsigned>(chains.size()));
169 478 : for(unsigned i=0; i<chains.size(); ++i) {
170 : int start,end;
171 : std::string errmsg;
172 314 : pdb.getResidueRange( chains[i], start, end, errmsg );
173 314 : if( errmsg.length()!=0 ) {
174 0 : error( errmsg );
175 : }
176 : AtomNumber astart,aend;
177 314 : pdb.getAtomRange( chains[i], astart, aend, errmsg );
178 314 : if( errmsg.length()!=0 ) {
179 0 : error( errmsg );
180 : }
181 314 : log.printf(" chain named %s contains residues %d to %d and atoms %u to %u \n",chains[i].c_str(),start,end,astart.serial(),aend.serial());
182 : }
183 :
184 : std::string python_bin;
185 164 : parse("PYTHON_BIN",python_bin);
186 164 : if(python_bin=="no") {
187 0 : log<<" python interpreter disabled\n";
188 : } else {
189 164 : pythonCmd=config::getEnvCommand();
190 164 : if(python_bin!="default") {
191 0 : log<<" forcing python interpreter: "<<python_bin<<"\n";
192 0 : pythonCmd+=" env PLUMED_PYTHON_BIN="+python_bin;
193 : }
194 : bool sorted=true;
195 164 : const auto & at=pdb.getAtomNumbers();
196 135762 : for(unsigned i=0; i<at.size(); i++) {
197 135598 : if(at[i].index()!=i) {
198 : sorted=false;
199 : }
200 : }
201 164 : if(!sorted) {
202 2 : log<<" PDB is not sorted, python interpreter will be disabled\n";
203 162 : } else if(!Subprocess::available()) {
204 0 : log<<" subprocess is not available, python interpreter will be disabled\n";
205 : } else {
206 162 : enablePythonInterpreter=true;
207 : }
208 : }
209 164 : }
210 164 : }
211 :
212 3 : std::map<std::string,std::string> GenericMolInfo::getSpecialKeywords() {
213 : std::map<std::string,std::string> mapkeys;
214 3 : mapkeys.insert( std::pair<std::string,std::string>("@mda:","atom selection built using syntax for <a href=\\\"https://docs.mdanalysis.org/stable/documentation_pages/selections.html\\\">MDAnalysis</a>") );
215 3 : mapkeys.insert( std::pair<std::string,std::string>("@mdt:","atom selection built using syntax for <a href=\\\"https://www.mdtraj.org/1.9.8.dev0/index.html\\\">mdtraj</a>") );
216 3 : mapkeys.insert( std::pair<std::string,std::string>("@vmdexec:","atom selection built using syntax for <a href=\\\"https://www.ks.uiuc.edu/Research/vmd/\\\">VMD</a>") );
217 3 : mapkeys.insert( std::pair<std::string,std::string>("@vmd:","atom selection built using syntax for <a href=\\\"https://www.ks.uiuc.edu/Research/vmd/\\\">VMD</a> python module") );
218 3 : mapkeys.insert( std::pair<std::string,std::string>("@nucleic","all atoms that are part of a DNA or RNA molecule") );
219 3 : mapkeys.insert( std::pair<std::string,std::string>("@protein","all atoms that are part of a protein") );
220 3 : mapkeys.insert( std::pair<std::string,std::string>("@water","all water molecules") );
221 3 : mapkeys.insert( std::pair<std::string,std::string>("@ions","all the ions") );
222 3 : mapkeys.insert( std::pair<std::string,std::string>("@hydrogens","all hydrogen atoms") );
223 3 : mapkeys.insert( std::pair<std::string,std::string>("@nonhydrogens","all non hydrogen atoms") );
224 3 : mapkeys.insert( std::pair<std::string,std::string>("@phi-","the four atoms that are required to calculate the phi dihedral for") );
225 3 : mapkeys.insert( std::pair<std::string,std::string>("@psi-","the four atoms that are required to calculate the psi dihedral for") );
226 3 : mapkeys.insert( std::pair<std::string,std::string>("@omega-","the four atoms that are required to calculate the omega dihedral for") );
227 3 : mapkeys.insert( std::pair<std::string,std::string>("@chi1-","the four atoms that are required to calculate the chi1 dihedral for") );
228 3 : mapkeys.insert( std::pair<std::string,std::string>("@chi2-","the four atoms that are required to calculate the chi2 dihedral for") );
229 3 : mapkeys.insert( std::pair<std::string,std::string>("@chi3-","the four atoms that are required to calculate the chi3 dihedral for") );
230 3 : mapkeys.insert( std::pair<std::string,std::string>("@chi4-","the four atoms that are required to calculate the chi4 dihedral for") );
231 3 : mapkeys.insert( std::pair<std::string,std::string>("@chi5-","the four atoms that are required to calculate the chi5 dihedral for") );
232 3 : mapkeys.insert( std::pair<std::string,std::string>("@sidechain-","the protein sidechain atoms in") );
233 3 : mapkeys.insert( std::pair<std::string,std::string>("@back-","the protein/dna/rna backbone atoms in") );
234 3 : mapkeys.insert( std::pair<std::string,std::string>("@alpha-","the four atoms that are required to calculate the alpha backbone dihedral for") );
235 3 : mapkeys.insert( std::pair<std::string,std::string>("@beta-","the four atoms that are required to calculate the beta backbone dihedral for") );
236 3 : mapkeys.insert( std::pair<std::string,std::string>("@gamma-","the four atoms that are required to calculate the gamma backbone dihedral for") );
237 3 : mapkeys.insert( std::pair<std::string,std::string>("@delta-","the four atoms that are required to calculate the delta backbone dihedral for") );
238 3 : mapkeys.insert( std::pair<std::string,std::string>("@epsilon-","the four atoms that are required to calculate the backbone epsilon dihedral for") );
239 3 : mapkeys.insert( std::pair<std::string,std::string>("@zeta-","the four atoms that are required to calculate the zeta backbone dihedral for") );
240 3 : mapkeys.insert( std::pair<std::string,std::string>("@v0-","the four atoms that are required to calculate the v0 sugar dihedral for") );
241 3 : mapkeys.insert( std::pair<std::string,std::string>("@v1-","the four atoms that are required to calculate the v1 sugar dihedral for") );
242 3 : mapkeys.insert( std::pair<std::string,std::string>("@v2-","the four atoms that are required to calculate the v2 sugar dihedral for") );
243 3 : mapkeys.insert( std::pair<std::string,std::string>("@v3-","the four atoms that are required to calculate the v3 sugar dihedral for") );
244 3 : mapkeys.insert( std::pair<std::string,std::string>("@v4-","the four atoms that are required to calculate the v4 sugar dihedral for") );
245 3 : mapkeys.insert( std::pair<std::string,std::string>("@chi-","the four atoms that are required to alcullate the chi dihedral for") );
246 3 : mapkeys.insert( std::pair<std::string,std::string>("@sugar-","the heavy atoms of the sugar in") );
247 3 : mapkeys.insert( std::pair<std::string,std::string>("@base-","the heavy atoms of the base in") );
248 3 : mapkeys.insert( std::pair<std::string,std::string>("@lcs-","an ordered triplet of atoms on the 6-membered ring of the nucleobase in") );
249 3 : return mapkeys;
250 : }
251 :
252 60 : void GenericMolInfo::getBackbone( std::vector<std::string>& restrings, const std::string& fortype, std::vector< std::vector<AtomNumber> >& backbone ) {
253 60 : if( fortype!=mytype ) {
254 0 : error("cannot calculate a variable designed for " + fortype + " molecules for molecule type " + mytype );
255 : }
256 60 : if( MolDataClass::numberOfAtomsPerResidueInBackbone( mytype )==0 ) {
257 0 : error("backbone is not defined for molecule type " + mytype );
258 : }
259 :
260 60 : if( read_backbone.size()!=0 ) {
261 0 : if( restrings.size()!=1 ) {
262 0 : error("cannot interpret anything other than all for residues when using CHAIN keywords");
263 : }
264 0 : if( restrings[0]!="all" ) {
265 0 : error("cannot interpret anything other than all for residues when using CHAIN keywords");
266 : }
267 0 : backbone.resize( read_backbone.size() );
268 0 : for(unsigned i=0; i<read_backbone.size(); ++i) {
269 0 : backbone[i].resize( read_backbone[i].size() );
270 0 : for(unsigned j=0; j<read_backbone[i].size(); ++j) {
271 0 : backbone[i][j]=read_backbone[i][j];
272 : }
273 : }
274 : } else {
275 : bool useter=false; // This is used to deal with terminal groups in WHOLEMOLECULES
276 : // residues are identified by (chain, residue number), so that chains that use
277 : // the same residue numbers (e.g. all starting from 1) are not mixed up
278 60 : std::vector< std::pair<std::string,int> > residues;
279 60 : std::vector<std::string> chains;
280 : {
281 60 : std::vector<std::string> allchains;
282 60 : pdb.getChainNames( allchains );
283 : // a chain can appear in more than one block of the pdb file: list it once
284 756 : for(const auto & c : allchains) {
285 696 : if( std::find(chains.begin(),chains.end(),c)==chains.end() ) {
286 696 : chains.push_back(c);
287 : }
288 : }
289 60 : }
290 60 : if( restrings.size()==1 ) {
291 58 : useter=( restrings[0].find("ter")!=std::string::npos );
292 : }
293 60 : if( restrings.size()==1 && restrings[0].find("all")!=std::string::npos ) {
294 : // list explicitly the residues that are present in each chain:
295 : // this works also with negative residue numbers and with gaps in the numbering
296 748 : for(const auto & c : chains) {
297 692 : const std::vector<int> res = pdb.getResiduesInChain( c );
298 692 : if( res.empty() ) {
299 : continue;
300 : }
301 : std::size_t first=0, last=res.size();
302 692 : if( !useter ) {
303 1384 : if( MolDataClass::isTerminalGroup( mytype, pdb.getResidueName( res.front(), c ) ) ) {
304 : first++;
305 : }
306 1384 : if( last>first && MolDataClass::isTerminalGroup( mytype, pdb.getResidueName( res.back(), c ) ) ) {
307 690 : last--;
308 : }
309 : }
310 6310 : for(std::size_t k=first; k<last; ++k) {
311 5618 : residues.emplace_back( c, res[k] );
312 : }
313 : }
314 : } else {
315 4 : Tools::interpretRanges(restrings);
316 : // chains containing each residue number
317 : std::map< int, std::vector<std::string> > chains_of_residue;
318 8 : for(const auto & c : chains) {
319 60 : for(const int r : pdb.getResiduesInChain( c )) {
320 56 : chains_of_residue[r].push_back( c );
321 : }
322 : }
323 : std::string ambiguous;
324 28 : for(const auto & rs : restrings) {
325 : int r;
326 24 : if( !Tools::convertNoexcept( rs, r ) ) {
327 0 : error("cannot interpret " + rs + " as a residue number");
328 : }
329 : const auto p=chains_of_residue.find( r );
330 24 : if( p==chains_of_residue.end() ) {
331 0 : error("residue " + rs + " is not in the reference pdb file");
332 : }
333 : // a residue number present in more than one chain is taken from the first chain, as in previous versions
334 24 : if( p->second.size()>1 ) {
335 0 : ambiguous += " " + rs;
336 : }
337 24 : residues.emplace_back( p->second[0], r );
338 : }
339 4 : if( ambiguous.length()>0 ) {
340 0 : warning("residues" + ambiguous + " are present in more than one chain of the reference pdb file: they are taken from the first chain that contains them. Use RESIDUES=all to include all the chains");
341 : }
342 : }
343 60 : if( residues.empty() ) {
344 0 : error("no residues found in the reference pdb file");
345 : }
346 :
347 : // Convert the list of involved residues into a list of segments of consecutive residues of the same chain
348 60 : std::vector< std::vector< std::pair<std::string,int> > > segments;
349 60 : std::vector< std::pair<std::string,int> > thissegment;
350 60 : thissegment.push_back( residues[0] );
351 5642 : for(unsigned i=1; i<residues.size(); ++i) {
352 5582 : if( residues[i-1].second+1!=residues[i].second || residues[i-1].first!=residues[i].first ) {
353 638 : segments.push_back(thissegment);
354 638 : thissegment.resize(0);
355 : }
356 5582 : thissegment.push_back( residues[i] );
357 : }
358 60 : segments.push_back( thissegment );
359 :
360 : // And now get the backbone atoms from each segment
361 60 : backbone.resize( segments.size() );
362 60 : std::vector<AtomNumber> atomnumbers;
363 758 : for(unsigned i=0; i<segments.size(); ++i) {
364 6340 : for(unsigned j=0; j<segments[i].size(); ++j) {
365 5642 : const std::string & chain=segments[i][j].first;
366 5642 : const int resnum=segments[i][j].second;
367 : std::string num;
368 5642 : Tools::convert( resnum, num );
369 5642 : num += " of chain " + chain;
370 5642 : std::string resname=pdb.getResidueName( resnum, chain );
371 5642 : if( !MolDataClass::allowedResidue(mytype, resname) ) {
372 0 : error("residue " + num + " is not recognized for moltype " + mytype );
373 : }
374 5642 : if( !useter && MolDataClass::isTerminalGroup( mytype, resname ) ) {
375 0 : error("residue " + num + " appears to be a terminal group");
376 : }
377 5642 : if( resname=="GLY" ) {
378 0 : warning("GLY residues are achiral - assuming HA1 atom is in CB position");
379 : }
380 5642 : MolDataClass::getBackboneForResidue( mytype, resnum, chain, pdb, atomnumbers );
381 5642 : if( atomnumbers.size()==0 ) {
382 0 : error("Could not find required backbone atom in residue number " + num );
383 : } else {
384 33852 : for(unsigned k=0; k<atomnumbers.size(); ++k) {
385 28210 : backbone[i].push_back( atomnumbers[k] );
386 : }
387 : }
388 5642 : atomnumbers.resize(0);
389 : }
390 : }
391 60 : }
392 60 : }
393 :
394 801 : void GenericMolInfo::interpretSymbol( const std::string& symbol, std::vector<AtomNumber>& atoms ) {
395 801 : if(Tools::startWith(symbol,"mdt:") || Tools::startWith(symbol,"mda:") || Tools::startWith(symbol,"vmd:") || Tools::startWith(symbol,"vmdexec:")) {
396 :
397 24 : plumed_assert(enablePythonInterpreter);
398 :
399 48 : log<<" symbol " + symbol + " will be sent to python interpreter\n";
400 24 : if(!selector_running) {
401 3 : log<<" MOLINFO "<<getLabel()<<": starting python interpreter\n";
402 3 : if(comm.Get_rank()==0) {
403 8 : selector=Tools::make_unique<Subprocess>(pythonCmd+" \""+config::getPlumedRoot()+"\"/scripts/selector.sh --pdb " + reference);
404 2 : selector->stop();
405 : }
406 3 : selector_running=true;
407 : }
408 :
409 24 : atoms.resize(0);
410 :
411 24 : if(comm.Get_rank()==0) {
412 16 : int ok=0;
413 : std::string error_msg;
414 : // this is a complicated way to have the exception propagated with MPI.
415 : // It is necessary since only one process calls the selector.
416 : // Probably I should recycle the exception propagation at library boundaries
417 : // to allow transferring the exception to other processes.
418 : try {
419 16 : plumed_assert(selector) << "Python interpreter is disabled, selection " + symbol + " cannot be interpreted";
420 : auto h=selector->contStop(); // stops again when it goes out of scope
421 16 : (*selector) << symbol << "\n";
422 16 : selector->flush();
423 : std::string res;
424 16 : std::vector<std::string> words;
425 : while(true) {
426 16 : selector->getline(res);
427 32 : words=Tools::getWords(res);
428 16 : if(!words.empty() && words[0]=="Error") {
429 0 : plumed_error()<<res;
430 : }
431 16 : if(!words.empty() && words[0]=="Selection:") {
432 : break;
433 : }
434 : }
435 : words.erase(words.begin());
436 712 : for(const auto & w : words) {
437 : int n;
438 696 : if(w.empty()) {
439 0 : continue;
440 : }
441 696 : Tools::convert(w,n);
442 1392 : atoms.push_back(AtomNumber::serial(n));
443 : }
444 16 : ok=1;
445 32 : } catch (const Exception & e) {
446 0 : error_msg=e.what();
447 0 : }
448 16 : comm.Bcast(ok,0);
449 16 : if(!ok) {
450 0 : size_t s=error_msg.length();
451 0 : comm.Bcast(s,0);
452 0 : comm.Bcast(error_msg,0);
453 0 : throw Exception()<<error_msg;
454 : }
455 16 : size_t nat=atoms.size();
456 16 : comm.Bcast(nat,0);
457 16 : comm.Bcast(atoms,0);
458 : } else {
459 8 : int ok=0;
460 : std::string error_msg;
461 8 : comm.Bcast(ok,0);
462 8 : if(!ok) {
463 : size_t s;
464 0 : comm.Bcast(s,0);
465 0 : error_msg.resize(s);
466 0 : comm.Bcast(error_msg,0);
467 0 : throw Exception()<<error_msg;
468 : }
469 8 : size_t nat=0;
470 8 : comm.Bcast(nat,0);
471 8 : atoms.resize(nat);
472 8 : comm.Bcast(atoms,0);
473 : }
474 24 : log<<" selection interpreted using ";
475 24 : if(Tools::startWith(symbol,"mdt:")) {
476 12 : log<<"mdtraj "<<cite("McGibbon et al, Biophys. J., 109, 1528 (2015)")<<"\n";
477 : }
478 24 : if(Tools::startWith(symbol,"mda:")) {
479 36 : log<<"MDAnalysis "<<cite("Gowers et al, Proceedings of the 15th Python in Science Conference, doi:10.25080/majora-629e541a-00e (2016)")<<"\n";
480 : }
481 24 : if(Tools::startWith(symbol,"vmdexec:")) {
482 0 : log<<"VMD "<<cite("Humphrey, Dalke, and Schulten, K., J. Molec. Graphics, 14, 33 (1996)")<<"\n";
483 : }
484 24 : if(Tools::startWith(symbol,"vmd:")) {
485 0 : log<<"VMD (github.com/Eigenstate/vmd-python) "<<cite("Humphrey, Dalke, and Schulten, K., J. Molec. Graphics, 14, 33 (1996)")<<"\n";
486 : }
487 24 : return;
488 : }
489 777 : MolDataClass::specialSymbol( mytype, symbol, pdb, atoms );
490 777 : if(atoms.empty()) {
491 0 : error(symbol + " not found in your MOLINFO structure");
492 : }
493 : }
494 :
495 72126 : std::string GenericMolInfo::getAtomName(AtomNumber a)const {
496 72126 : return pdb.getAtomName(a);
497 : }
498 :
499 3846 : bool GenericMolInfo::checkForAtom(AtomNumber a)const {
500 3846 : return pdb.checkForAtom(a);
501 : }
502 :
503 62194 : int GenericMolInfo::getResidueNumber(AtomNumber a)const {
504 62194 : return pdb.getResidueNumber(a);
505 : }
506 :
507 19512 : unsigned GenericMolInfo::getPDBsize()const {
508 19512 : return pdb.size();
509 : }
510 :
511 61190 : std::string GenericMolInfo::getResidueName(AtomNumber a)const {
512 61190 : return pdb.getResidueName(a);
513 : }
514 :
515 468 : std::string GenericMolInfo::getChainID(AtomNumber a)const {
516 468 : return pdb.getChainID(a);
517 : }
518 :
519 3789 : Vector GenericMolInfo::getPosition(AtomNumber a)const {
520 3789 : return pdb.getPosition(a);
521 : }
522 :
523 29503 : bool GenericMolInfo::isWhole() const {
524 29503 : return iswhole_;
525 : }
526 :
527 9653 : void GenericMolInfo::prepare() {
528 9653 : if(selector_running) {
529 3 : log<<" MOLINFO "<<getLabel()<<": killing python interpreter\n";
530 : selector.reset();
531 3 : selector_running=false;
532 : }
533 9653 : }
534 :
535 : }
|
