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 :
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 56 : void SetupMolInfo::registerKeywords( Keywords& keys ) {
42 56 : ActionSetup::registerKeywords(keys);
43 224 : 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 : "For more details on the PDB file format visit http://www.wwpdb.org/docs.html");
46 280 : 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 224 : keys.add("atoms","CHAIN","(for masochists ( mostly Davide Branduardi ) ) The atoms involved in each of the chains of interest in the structure.");
48 56 : }
49 :
50 275 : SetupMolInfo::~SetupMolInfo() {
51 55 : delete &pdb;
52 110 : }
53 :
54 55 : SetupMolInfo::SetupMolInfo( const ActionOptions&ao ):
55 : Action(ao),
56 : ActionSetup(ao),
57 : ActionAtomistic(ao),
58 110 : pdb(*new(PDB))
59 : {
60 : // Read what is contained in the pdb file
61 110 : parse("MOLTYPE",mytype);
62 :
63 110 : std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
64 55 : if( moldat.size()!=0 ) error("cannot use more than one MOLINFO action in input");
65 :
66 : std::vector<AtomNumber> backbone;
67 110 : parseAtomList("CHAIN",backbone);
68 55 : if( read_backbone.size()==0 ) {
69 0 : for(unsigned i=1;; ++i) {
70 110 : parseAtomList("CHAIN",i,backbone);
71 55 : if( backbone.size()==0 ) break;
72 0 : read_backbone.push_back(backbone);
73 0 : backbone.resize(0);
74 : }
75 : } else {
76 0 : read_backbone.push_back(backbone);
77 : }
78 55 : if( read_backbone.size()==0 ) {
79 110 : std::string reference; parse("STRUCTURE",reference);
80 :
81 165 : if( ! pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()))plumed_merror("missing input file " + reference );
82 :
83 110 : std::vector<std::string> chains; pdb.getChainNames( chains );
84 165 : log.printf(" pdb file named %s contains %u chains \n",reference.c_str(), static_cast<unsigned>(chains.size()));
85 434 : for(unsigned i=0; i<chains.size(); ++i) {
86 : unsigned start,end; std::string errmsg;
87 216 : pdb.getResidueRange( chains[i], start, end, errmsg );
88 108 : if( errmsg.length()!=0 ) error( errmsg );
89 : AtomNumber astart,aend;
90 216 : pdb.getAtomRange( chains[i], astart, aend, errmsg );
91 108 : if( errmsg.length()!=0 ) error( errmsg );
92 324 : 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());
93 : }
94 : }
95 55 : }
96 :
97 25 : void SetupMolInfo::getBackbone( std::vector<std::string>& restrings, const std::string& fortype, std::vector< std::vector<AtomNumber> >& backbone ) {
98 50 : if( fortype!=mytype ) error("cannot calculate a variable designed for " + fortype + " molecules for molecule type " + mytype );
99 25 : if( MolDataClass::numberOfAtomsPerResidueInBackbone( mytype )==0 ) error("backbone is not defined for molecule type " + mytype );
100 :
101 25 : if( read_backbone.size()!=0 ) {
102 0 : if( restrings.size()!=1 ) error("cannot interpret anything other than all for residues when using CHAIN keywords");
103 0 : if( restrings[0]!="all" ) error("cannot interpret anything other than all for residues when using CHAIN keywords");
104 0 : backbone.resize( read_backbone.size() );
105 0 : for(unsigned i=0; i<read_backbone.size(); ++i) {
106 0 : backbone[i].resize( read_backbone[i].size() );
107 0 : for(unsigned j=0; j<read_backbone[i].size(); ++j) backbone[i][j]=read_backbone[i][j];
108 : }
109 : } else {
110 : bool useter=false; // This is used to deal with terminal groups in WHOLEMOLECULES
111 25 : if( restrings.size()==1 ) {
112 23 : useter=( restrings[0].find("ter")!=std::string::npos );
113 23 : if( restrings[0].find("all")!=std::string::npos ) {
114 42 : std::vector<std::string> chains; pdb.getChainNames( chains );
115 1023 : for(unsigned i=0; i<chains.size(); ++i) {
116 : unsigned r_start, r_end; std::string errmsg, mm, nn;
117 654 : pdb.getResidueRange( chains[i], r_start, r_end, errmsg );
118 327 : if( !useter ) {
119 327 : std::string resname = pdb.getResidueName( r_start );
120 327 : if( MolDataClass::isTerminalGroup( mytype, resname ) ) r_start++;
121 654 : resname = pdb.getResidueName( r_end );
122 327 : if( MolDataClass::isTerminalGroup( mytype, resname ) ) r_end--;
123 : }
124 327 : Tools::convert(r_start,mm); Tools::convert(r_end,nn);
125 390 : if(i==0) restrings[0] = mm + "-" + nn;
126 918 : else restrings.push_back( mm + "-" + nn );
127 : }
128 : }
129 : }
130 25 : Tools::interpretRanges(restrings);
131 :
132 : // Convert the list of involved residues into a list of segments of chains
133 25 : int nk, nj; std::vector< std::vector<unsigned> > segments;
134 : std::vector<unsigned> thissegment;
135 50 : Tools::convert(restrings[0],nk); thissegment.push_back(nk);
136 7931 : for(unsigned i=1; i<restrings.size(); ++i) {
137 5254 : Tools::convert(restrings[i-1],nk);
138 2627 : Tools::convert(restrings[i],nj);
139 10200 : if( (nk+1)!=nj || pdb.getChainID(nk)!=pdb.getChainID(nj) ) {
140 308 : segments.push_back(thissegment);
141 308 : thissegment.resize(0);
142 : }
143 5254 : thissegment.push_back(nj);
144 : }
145 25 : segments.push_back( thissegment );
146 :
147 : // And now get the backbone atoms from each segment
148 25 : backbone.resize( segments.size() );
149 : std::vector<AtomNumber> atomnumbers;
150 1049 : for(unsigned i=0; i<segments.size(); ++i) {
151 8622 : for(unsigned j=0; j<segments[i].size(); ++j) {
152 5304 : std::string resname=pdb.getResidueName( segments[i][j] );
153 2652 : if( !MolDataClass::allowedResidue(mytype, resname) ) {
154 0 : std::string num; Tools::convert( segments[i][j], num );
155 0 : error("residue " + num + " is not recognized for moltype " + mytype );
156 : }
157 2652 : if( !useter && MolDataClass::isTerminalGroup( mytype, resname ) ) {
158 0 : std::string num; Tools::convert( segments[i][j], num );
159 0 : error("residue " + num + " appears to be a terminal group");
160 : }
161 2652 : if( resname=="GLY" ) warning("GLY residues are achiral - assuming HA1 atom is in CB position");
162 5304 : MolDataClass::getBackboneForResidue( mytype, segments[i][j], pdb, atomnumbers );
163 2652 : if( atomnumbers.size()==0 ) {
164 0 : std::string num; Tools::convert( segments[i][j], num );
165 0 : error("Could not find required backbone atom in residue number " + num );
166 : } else {
167 45084 : for(unsigned k=0; k<atomnumbers.size(); ++k) backbone[i].push_back( atomnumbers[k] );
168 : }
169 2652 : atomnumbers.resize(0);
170 : }
171 : }
172 : }
173 25 : }
174 :
175 502 : void SetupMolInfo::interpretSymbol( const std::string& symbol, std::vector<AtomNumber>& atoms )const {
176 502 : MolDataClass::specialSymbol( mytype, symbol, pdb, atoms );
177 502 : }
178 :
179 9012 : std::string SetupMolInfo::getAtomName(AtomNumber a)const {
180 9012 : return pdb.getAtomName(a);
181 : }
182 :
183 3284 : unsigned SetupMolInfo::getResidueNumber(AtomNumber a)const {
184 3284 : return pdb.getResidueNumber(a);
185 : }
186 :
187 6244 : std::string SetupMolInfo::getResidueName(AtomNumber a)const {
188 6244 : return pdb.getResidueName(a);
189 : }
190 :
191 4839 : }
|
