Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 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 "SetupMolInfo.h"
23 : #include "Atoms.h"
24 : #include "ActionRegister.h"
25 : #include "ActionSet.h"
26 : #include "PlumedMain.h"
27 : #include "tools/MolDataClass.h"
28 : #include "tools/PDB.h"
29 : #include "config/Config.h"
30 :
31 : namespace PLMD {
32 :
33 :
34 : /*
35 : This action is defined in core/ as it is used by other actions.
36 : Anyway, it is registered in setup/, so that excluding that module from
37 : compilation will exclude it from plumed.
38 : */
39 :
40 :
41 75 : void SetupMolInfo::registerKeywords( Keywords& keys ) {
42 75 : ActionSetup::registerKeywords(keys);
43 : keys.add("compulsory","STRUCTURE","a file in pdb format containing a reference structure. "
44 : "This is used to defines the atoms in the various residues, chains, etc . "
45 300 : "For more details on the PDB file format visit http://www.wwpdb.org/docs.html");
46 375 : keys.add("compulsory","MOLTYPE","protein","what kind of molecule is contained in the pdb file - usually not needed since protein/RNA/DNA are compatible");
47 375 : keys.add("compulsory","PYTHON_BIN","default","python interpreter");
48 300 : keys.add("atoms","CHAIN","(for masochists ( mostly Davide Branduardi ) ) The atoms involved in each of the chains of interest in the structure.");
49 300 : keys.add("hidden","STRIDE","frequency for resetting the python interpreter. Should be 1.");
50 75 : }
51 :
52 370 : SetupMolInfo::~SetupMolInfo() {
53 : // empty destructor to delete unique_ptr
54 148 : }
55 :
56 74 : SetupMolInfo::SetupMolInfo( const ActionOptions&ao ):
57 : Action(ao),
58 : ActionSetup(ao),
59 : ActionPilot(ao),
60 296 : ActionAtomistic(ao)
61 : {
62 74 : plumed_assert(getStride()==1);
63 : // Read what is contained in the pdb file
64 148 : parse("MOLTYPE",mytype);
65 :
66 148 : std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
67 74 : if( moldat.size()!=0 ) error("cannot use more than one MOLINFO action in input");
68 :
69 : std::vector<AtomNumber> backbone;
70 148 : parseAtomList("CHAIN",backbone);
71 74 : if( read_backbone.size()==0 ) {
72 0 : for(unsigned i=1;; ++i) {
73 148 : parseAtomList("CHAIN",i,backbone);
74 74 : if( backbone.size()==0 ) break;
75 0 : read_backbone.push_back(backbone);
76 0 : backbone.resize(0);
77 0 : }
78 : } else {
79 0 : read_backbone.push_back(backbone);
80 : }
81 74 : if( read_backbone.size()==0 ) {
82 148 : parse("STRUCTURE",reference);
83 :
84 222 : if( ! pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()))plumed_merror("missing input file " + reference );
85 :
86 74 : std::vector<std::string> chains; pdb.getChainNames( chains );
87 148 : log.printf(" pdb file named %s contains %u chains \n",reference.c_str(), static_cast<unsigned>(chains.size()));
88 338 : for(unsigned i=0; i<chains.size(); ++i) {
89 : unsigned start,end; std::string errmsg;
90 132 : pdb.getResidueRange( chains[i], start, end, errmsg );
91 132 : if( errmsg.length()!=0 ) error( errmsg );
92 : AtomNumber astart,aend;
93 132 : pdb.getAtomRange( chains[i], astart, aend, errmsg );
94 132 : if( errmsg.length()!=0 ) error( errmsg );
95 264 : log.printf(" chain named %s contains residues %u to %u and atoms %u to %u \n",chains[i].c_str(),start,end,astart.serial(),aend.serial());
96 : }
97 :
98 : std::string python_bin;
99 148 : parse("PYTHON_BIN",python_bin);
100 74 : if(python_bin=="no") {
101 0 : log<<" python interpreter disabled\n";
102 : } else {
103 148 : pythonCmd=config::getEnvCommand();
104 74 : if(python_bin!="default") {
105 0 : log<<" forcing python interpreter: "<<python_bin<<"\n";
106 0 : pythonCmd+=" env PLUMED_PYTHON_BIN="+python_bin;
107 : }
108 : bool sorted=true;
109 74 : const auto & at=pdb.getAtomNumbers();
110 77392 : for(unsigned i=0; i<at.size(); i++) {
111 38659 : if(at[i].index()!=i) sorted=false;
112 : }
113 74 : if(!sorted) {
114 1 : log<<" PDB is not sorted, python interpreter will be disabled\n";
115 73 : } else if(!Subprocess::available()) {
116 0 : log<<" subprocess is not available, python interpreter will be disabled\n";
117 : } else {
118 73 : enablePythonInterpreter=true;
119 : }
120 74 : }
121 : }
122 74 : }
123 :
124 25 : void SetupMolInfo::getBackbone( std::vector<std::string>& restrings, const std::string& fortype, std::vector< std::vector<AtomNumber> >& backbone ) {
125 50 : if( fortype!=mytype ) error("cannot calculate a variable designed for " + fortype + " molecules for molecule type " + mytype );
126 25 : if( MolDataClass::numberOfAtomsPerResidueInBackbone( mytype )==0 ) error("backbone is not defined for molecule type " + mytype );
127 :
128 25 : if( read_backbone.size()!=0 ) {
129 0 : if( restrings.size()!=1 ) error("cannot interpret anything other than all for residues when using CHAIN keywords");
130 0 : if( restrings[0]!="all" ) error("cannot interpret anything other than all for residues when using CHAIN keywords");
131 0 : backbone.resize( read_backbone.size() );
132 0 : for(unsigned i=0; i<read_backbone.size(); ++i) {
133 0 : backbone[i].resize( read_backbone[i].size() );
134 0 : for(unsigned j=0; j<read_backbone[i].size(); ++j) backbone[i][j]=read_backbone[i][j];
135 : }
136 : } else {
137 : bool useter=false; // This is used to deal with terminal groups in WHOLEMOLECULES
138 25 : if( restrings.size()==1 ) {
139 23 : useter=( restrings[0].find("ter")!=std::string::npos );
140 23 : if( restrings[0].find("all")!=std::string::npos ) {
141 21 : std::vector<std::string> chains; pdb.getChainNames( chains );
142 675 : for(unsigned i=0; i<chains.size(); ++i) {
143 : unsigned r_start, r_end; std::string errmsg, mm, nn;
144 327 : pdb.getResidueRange( chains[i], r_start, r_end, errmsg );
145 327 : if( !useter ) {
146 327 : std::string resname = pdb.getResidueName( r_start );
147 327 : if( MolDataClass::isTerminalGroup( mytype, resname ) ) r_start++;
148 654 : resname = pdb.getResidueName( r_end );
149 327 : if( MolDataClass::isTerminalGroup( mytype, resname ) ) r_end--;
150 : }
151 327 : Tools::convert(r_start,mm); Tools::convert(r_end,nn);
152 369 : if(i==0) restrings[0] = mm + "-" + nn;
153 918 : else restrings.push_back( mm + "-" + nn );
154 21 : }
155 : }
156 : }
157 25 : Tools::interpretRanges(restrings);
158 :
159 : // Convert the list of involved residues into a list of segments of chains
160 : int nk, nj; std::vector< std::vector<unsigned> > segments;
161 : std::vector<unsigned> thissegment;
162 50 : Tools::convert(restrings[0],nk); thissegment.push_back(nk);
163 5304 : for(unsigned i=1; i<restrings.size(); ++i) {
164 5254 : Tools::convert(restrings[i-1],nk);
165 2627 : Tools::convert(restrings[i],nj);
166 9584 : if( (nk+1)!=nj || pdb.getChainID(nk)!=pdb.getChainID(nj) ) {
167 308 : segments.push_back(thissegment);
168 308 : thissegment.resize(0);
169 : }
170 5254 : thissegment.push_back(nj);
171 : }
172 25 : segments.push_back( thissegment );
173 :
174 : // And now get the backbone atoms from each segment
175 25 : backbone.resize( segments.size() );
176 : std::vector<AtomNumber> atomnumbers;
177 716 : for(unsigned i=0; i<segments.size(); ++i) {
178 5637 : for(unsigned j=0; j<segments[i].size(); ++j) {
179 2652 : std::string resname=pdb.getResidueName( segments[i][j] );
180 2652 : if( !MolDataClass::allowedResidue(mytype, resname) ) {
181 0 : std::string num; Tools::convert( segments[i][j], num );
182 0 : error("residue " + num + " is not recognized for moltype " + mytype );
183 : }
184 2652 : if( !useter && MolDataClass::isTerminalGroup( mytype, resname ) ) {
185 0 : std::string num; Tools::convert( segments[i][j], num );
186 0 : error("residue " + num + " appears to be a terminal group");
187 : }
188 2652 : if( resname=="GLY" ) warning("GLY residues are achiral - assuming HA1 atom is in CB position");
189 5304 : MolDataClass::getBackboneForResidue( mytype, segments[i][j], pdb, atomnumbers );
190 2652 : if( atomnumbers.size()==0 ) {
191 0 : std::string num; Tools::convert( segments[i][j], num );
192 0 : error("Could not find required backbone atom in residue number " + num );
193 : } else {
194 29172 : for(unsigned k=0; k<atomnumbers.size(); ++k) backbone[i].push_back( atomnumbers[k] );
195 : }
196 2652 : atomnumbers.resize(0);
197 : }
198 25 : }
199 : }
200 25 : }
201 :
202 563 : void SetupMolInfo::interpretSymbol( const std::string& symbol, std::vector<AtomNumber>& atoms ) {
203 2797 : if(Tools::startWith(symbol,"mdt:") || Tools::startWith(symbol,"mda:") || Tools::startWith(symbol,"vmd:") || Tools::startWith(symbol,"vmdexec:")) {
204 :
205 8 : plumed_assert(enablePythonInterpreter);
206 :
207 32 : log<<" symbol " + symbol + " will be sent to python interpreter\n";
208 8 : if(!selector) {
209 1 : log<<" MOLINFO "<<getLabel()<<": starting python interpreter\n";
210 1 : if(comm.Get_rank()==0) {
211 7 : selector.reset(new Subprocess(pythonCmd+" \""+config::getPlumedRoot()+"\"/scripts/selector.sh --pdb " + reference));
212 1 : selector->stop();
213 : }
214 : }
215 :
216 8 : if(comm.Get_rank()==0) {
217 8 : int ok=0;
218 : std::string error_msg;
219 : // this is a complicated way to have the exception propagated with MPI.
220 : // It is necessary since only one process calls the selector.
221 : // Probably I should recycle the exception propagation at library boundaries
222 : // to allow transfering the exception to other processes.
223 : try {
224 8 : plumed_assert(selector) << "Python interpreter is disabled, selection " + symbol + " cannot be interpreted";
225 : auto h=selector->contStop(); // stops again when it goes out of scope
226 8 : (*selector) << symbol << "\n";
227 8 : selector->flush();
228 : std::string res;
229 8 : std::vector<std::string> words;
230 : while(true) {
231 8 : selector->getline(res);
232 16 : words=Tools::getWords(res);
233 16 : if(!words.empty() && words[0]=="Error") plumed_error()<<res;
234 16 : if(!words.empty() && words[0]=="Selection:") break;
235 : }
236 8 : words.erase(words.begin());
237 8 : atoms.resize(0);
238 364 : for(auto & w : words) {
239 : int n;
240 348 : if(w.empty()) continue;
241 348 : Tools::convert(w,n);
242 696 : atoms.push_back(AtomNumber::serial(n));
243 : }
244 16 : ok=1;
245 0 : } catch (Exception & e) {
246 0 : error_msg=e.what();
247 : }
248 8 : comm.Bcast(ok,0);
249 8 : if(!ok) {
250 0 : size_t s=error_msg.length();
251 0 : comm.Bcast(s,0);
252 0 : comm.Bcast(error_msg,0);
253 0 : throw Exception()<<error_msg;
254 : }
255 8 : size_t nat=atoms.size();
256 8 : comm.Bcast(nat,0);
257 8 : comm.Bcast(atoms,0);
258 : } else {
259 0 : int ok=0;
260 : std::string error_msg;
261 0 : comm.Bcast(ok,0);
262 0 : if(!ok) {
263 : size_t s;
264 0 : comm.Bcast(s,0);
265 0 : error_msg.resize(s);
266 0 : comm.Bcast(error_msg,0);
267 0 : throw Exception()<<error_msg;
268 : }
269 0 : size_t nat=0;
270 0 : comm.Bcast(nat,0);
271 0 : comm.Bcast(atoms,0);
272 : }
273 8 : log<<" selection interpreted using ";
274 20 : if(Tools::startWith(symbol,"mdt:")) log<<"mdtraj "<<cite("McGibbon et al, Biophys. J., 109, 1528 (2015)")<<"\n";
275 28 : if(Tools::startWith(symbol,"mda:")) log<<"MDAnalysis "<<cite("Gowers et al, Proceedings of the 15th Python in Science Conference, doi:10.25080/majora-629e541a-00e (2016)")<<"\n";
276 16 : if(Tools::startWith(symbol,"vmdexec:")) log<<"VMD "<<cite("Humphrey, Dalke, and Schulten, K., J. Molec. Graphics, 14, 33 (1996)")<<"\n";
277 16 : if(Tools::startWith(symbol,"vmd:")) log<<"VMD (github.com/Eigenstate/vmd-python) "<<cite("Humphrey, Dalke, and Schulten, K., J. Molec. Graphics, 14, 33 (1996)")<<"\n";
278 563 : return;
279 : }
280 555 : MolDataClass::specialSymbol( mytype, symbol, pdb, atoms );
281 555 : if(atoms.empty()) error(symbol + " is not a label plumed knows");
282 : }
283 :
284 27456 : std::string SetupMolInfo::getAtomName(AtomNumber a)const {
285 27456 : return pdb.getAtomName(a);
286 : }
287 :
288 4818 : unsigned SetupMolInfo::getResidueNumber(AtomNumber a)const {
289 4818 : return pdb.getResidueNumber(a);
290 : }
291 :
292 9258 : std::string SetupMolInfo::getResidueName(AtomNumber a)const {
293 9258 : return pdb.getResidueName(a);
294 : }
295 :
296 5029 : void SetupMolInfo::prepare() {
297 5029 : if(selector) {
298 1 : log<<" MOLINFO "<<getLabel()<<": killing python interpreter\n";
299 1 : selector.reset();
300 : }
301 5029 : }
302 :
303 5874 : }
|
