1 // $Id: rigidbody.cpp 699 2009-05-24 21:30:19Z asaladin $
2 /****************************************************************************
3 * Copyright (C) 2006-2008 Adrien Saladin *
4 * adrien.saladin@gmail.com *
5 * Copyright (C) 2008 Pierre Poulain *
6 * Copyright (C) 2008 Sebastien Fiorucci *
7 * Copyright (C) 2008 Chantal Prevost *
8 * Copyright (C) 2008 Martin Zacharias *
9 * *
10 * This program is free software: you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation, either version 3 of the License, or *
13 * (at your option) any later version. *
14 * *
15 * This program is distributed in the hope that it will be useful, *
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18 * GNU General Public License for more details. *
19 * *
20 * You should have received a copy of the GNU General Public License *
21 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
22 * *
23 ***************************************************************************/
24
25
26 #include "rigidbody.h"
27 #include "atomselection.h"
28 #include "geometry.h"
29 #include "pdbio.h"
30
31
32 namespace PTools{
33
34
Rigidbody()35 Rigidbody::Rigidbody()
36 {
37 ResetMatrix();
38 }
39
40
41
Rigidbody(std::string filename)42 Rigidbody::Rigidbody(std::string filename)
43 {
44 ReadPDB(filename,*this);
45 ResetMatrix();
46 }
47
48
Rigidbody(const Rigidbody & model)49 Rigidbody::Rigidbody(const Rigidbody& model)
50 : CoordsArray(model)
51 {
52 //this copy constructor is needed because dbl[4][4] is not
53 // automatically copied with the default copy constructor
54 //TODO: verifier si c'est toujours le cas ...
55
56 this->mForces = model.mForces;
57 this->mAtomProp = model.mAtomProp;
58 this-> _description = model._description;
59
60 }
61
62
AddAtom(const Atomproperty & at,Coord3D co)63 void Rigidbody::AddAtom(const Atomproperty& at, Coord3D co)
64 {
65 mAtomProp.push_back(at);
66 AddCoord(co);
67 }
68
69
CopyAtom(uint i) const70 Atom Rigidbody::CopyAtom(uint i) const
71 {
72 Atom at(mAtomProp[i],GetCoords(i));
73 return at;
74 }
75
SetAtom(uint pos,const Atom & atom)76 void Rigidbody::SetAtom(uint pos, const Atom& atom)
77 {
78
79 if (/*pos<0 always true||*/ pos >= this->Size())
80 {
81 std::string message = "SetAtom: position ";
82 message += pos;
83 message += " is out of range";
84 throw std::out_of_range(message);
85 }
86 Atomproperty atp(atom);
87 Coord3D co(atom.GetCoords());
88 SetAtomProperty(pos, atp);
89 SetCoords(pos,co);
90 }
91
92
93
AddAtom(const Atom & at)94 void Rigidbody::AddAtom(const Atom& at)
95 {
96 Atomproperty atp(at);
97 Coord3D co = at.GetCoords();
98 AddAtom(atp,co);
99 }
100
101
FindCenter() const102 Coord3D Rigidbody::FindCenter() const
103 {
104 Coord3D center(0.0,0.0,0.0);
105 for (uint i=0; i< this->Size() ; i++)
106 {
107 center = center + GetCoords(i);
108 }
109 return ( (1.0/(dbl)this->Size())*center);
110 }
111
112
CenterToOrigin()113 void Rigidbody::CenterToOrigin()
114 {
115 Coord3D c = FindCenter();
116 Translate(Coord3D()-c);
117 }
118
RadiusGyration()119 dbl Rigidbody::RadiusGyration()
120 {
121 Coord3D c = this->FindCenter();
122 dbl r=0.0;
123 for (uint i=0; i< this->Size(); i++)
124 {
125 r += Norm2( c - this->GetCoords(i) );
126 }
127
128 dbl result = sqrt( r/ (double) this->Size() );
129 return result;
130 }
131
Radius()132 dbl Rigidbody::Radius()
133 {
134 Coord3D center = this->FindCenter();
135 uint size = this->Size();
136 dbl radius = 0.0;
137 for (uint i=0; i < size; i++)
138 {
139 dbl rad=Norm(center - this->GetCoords(i));
140 if (radius < rad) {radius=rad;}
141 }
142 return radius;
143 }
144
145
Translate(const Coord3D & tr)146 void Rigidbody::Translate(const Coord3D& tr)
147 {
148 CoordsArray::Translate(tr);
149 }
150
AttractEulerRotate(dbl phi,dbl ssi,dbl rot)151 void Rigidbody::AttractEulerRotate(dbl phi, dbl ssi, dbl rot)
152 {
153 CoordsArray::AttractEulerRotate(phi, ssi, rot);
154 }
155
156
157
158
159
SelectAllAtoms() const160 AtomSelection Rigidbody::SelectAllAtoms() const
161 {
162 AtomSelection newsel;
163 newsel.SetRigid(*this);
164 for (uint i=0; i < Size(); i++)
165 {
166 newsel.AddAtomIndex(i);
167 }
168
169
170 return newsel;
171
172 }
173
174
SelectAtomType(std::string atomtype)175 AtomSelection Rigidbody::SelectAtomType(std::string atomtype)
176 {
177 AtomSelection newsel;
178 newsel.SetRigid(*this);
179
180 for (uint i=0; i<Size(); i++)
181 {
182 if ( mAtomProp[i].GetType()==atomtype)
183 newsel.AddAtomIndex(i);
184 }
185
186 return newsel;
187 }
188
189
SelectResidType(std::string residtype)190 AtomSelection Rigidbody::SelectResidType(std::string residtype)
191 {
192 AtomSelection newsel;
193 newsel.SetRigid(*this);
194
195 for (uint i=0; i<Size(); i++)
196 {
197 if (mAtomProp[i].GetResidType()==residtype)
198 newsel.AddAtomIndex(i);
199 }
200 return newsel;
201 }
202
203
SelectChainId(std::string chainId)204 AtomSelection Rigidbody::SelectChainId(std::string chainId) {
205 AtomSelection newsel;
206 newsel.SetRigid(*this);
207 for (uint i=0; i<Size(); i++)
208 {
209 if (mAtomProp[i].GetChainId()==chainId)
210 newsel.AddAtomIndex(i);
211 }
212 return newsel;
213 }
214
SelectResRange(uint start,uint stop)215 AtomSelection Rigidbody::SelectResRange(uint start, uint stop)
216 {
217 AtomSelection newsel;
218 newsel.SetRigid(*this);
219
220 for (uint i=0; i < Size(); i++)
221 {
222 Atomproperty& atp ( mAtomProp[i] );
223 if (atp.GetResidId() >=start && atp.GetResidId() <= stop) newsel.AddAtomIndex(i);
224 }
225 return newsel;
226 }
227
228
CA()229 AtomSelection Rigidbody::CA() {
230 return SelectAtomType("CA");
231 }
232
isBackbone(const std::string & atomtype)233 bool isBackbone(const std::string & atomtype)
234 {
235
236 const std::string bbtypes[] = {"N", "CA", "C", "O"};
237 int const bbsize = sizeof(bbtypes)/sizeof(std::string);
238
239 for (int i =0; i<bbsize; i++)
240 {
241 if (atomtype == bbtypes[i] ) return true;
242 }
243
244 return false;
245 }
246
247
Backbone()248 AtomSelection Rigidbody::Backbone()
249 {
250 AtomSelection newsel;
251 newsel.SetRigid(*this);
252
253 for (uint i=0; i<this->Size(); i++)
254 {
255 if (isBackbone(CopyAtom(i).GetType()) )
256 {
257 newsel.AddAtomIndex(i);
258 }
259
260 }
261 return newsel;
262 }
263
264
265 /// operator +
operator +(const Rigidbody & rig)266 Rigidbody Rigidbody::operator+(const Rigidbody& rig) {
267 Rigidbody rigFinal(*this);
268 for (uint i=0; i< rig.Size() ; i++) {
269 rigFinal.AddCoord(rig.GetCoords(i));
270 rigFinal.mAtomProp.push_back(rig.mAtomProp[i]);
271 }
272 return rigFinal;
273 }
274
275
ABrotate(const Coord3D & A,const Coord3D & B,dbl theta)276 void Rigidbody::ABrotate(const Coord3D& A, const Coord3D& B, dbl theta)
277 {
278 PTools::ABrotate(A,B, *this, theta);
279 }
280
281
282
283 // unused by UGENE
284 /*std::string Rigidbody::PrintPDB() const
285 {
286 uint size=this->Size();
287
288 std::string output;
289 for (uint i=0; i < size ; i++)
290 {
291 Atom at(mAtomProp[i], this->GetCoords(i));
292 output = output + at.ToPdbString();
293 }
294 return output;
295 }*/
296
297
ApplyMatrix(const Matrix & mat)298 void Rigidbody::ApplyMatrix(const Matrix& mat)
299 {
300
301 dbl mat44[4][4];
302 for(uint i=0; i<4;i++)
303 for(uint j=0;j<4;j++)
304 mat44[i][j]=mat(i,j);
305 CoordsArray::MatrixMultiply(mat44);
306 }
307
308
309
310
311 } //namespace PTools
312