1 /*
2 * (c) 2004 Iowa State University
3 * see the LICENSE file in the top level directory
4 */
5
6 /* Internals.cpp
7
8 Classes to handle internal coordinates
9
10 BMB 4/1998
11 */
12
13 #include "Globals.h"
14 #include "GlobalExceptions.h"
15 #include "Internals.h"
16 #include "MoleculeData.h"
17 #include "Frame.h"
18 #include "Math3D.h"
19 #include <iostream>
20
Internals(void)21 Internals::Internals(void) {
22 MOPacStyle = NULL;
23 GeneralStyle = NULL;
24 }
Internals(Internals * copy)25 Internals::Internals(Internals * copy) {
26 MOPacStyle = NULL;
27 GeneralStyle = NULL;
28 if (copy->MOPacStyle)
29 MOPacStyle = new MOPacInternals(copy->MOPacStyle);
30 }
~Internals(void)31 Internals::~Internals(void) {
32 if (MOPacStyle) delete MOPacStyle;
33 if (GeneralStyle) delete GeneralStyle;
34 }
CreateMOPacInternals(long Num)35 void Internals::CreateMOPacInternals(long Num) {
36 if (!MOPacStyle) {
37 MOPacStyle = new MOPacInternals(Num);
38 }
39 }
MOPacInternals(long Num)40 MOPacInternals::MOPacInternals(long Num) {
41 ConnectionAtoms = new long[Num];
42 Values = new float[Num];
43 Type = new char[Num];
44 if (!ConnectionAtoms || !Values || !Type) throw MemoryError();
45 Allocation = Num;
46 Count = 0;
47 }
MOPacInternals(MOPacInternals * copy)48 MOPacInternals::MOPacInternals(MOPacInternals * copy) {
49 ConnectionAtoms = new long[copy->Count];
50 Values = new float[copy->Count];
51 Type = new char[copy->Count];
52 if (!ConnectionAtoms || !Values || !Type) throw MemoryError();
53 Allocation = copy->Count;
54 Count = copy->Count;
55 for (long i=0; i<Count; i++) {
56 ConnectionAtoms[i] = copy->ConnectionAtoms[i];
57 Values[i] = copy->Values[i];
58 Type[i] = copy->Type[i];
59 }
60 }
~MOPacInternals(void)61 MOPacInternals::~MOPacInternals(void) {
62 if (ConnectionAtoms) delete [] ConnectionAtoms;
63 if (Values) delete [] Values;
64 if (Type) delete [] Type;
65 }
WriteCoordinatesToFile(BufferFile * File,MoleculeData * MainData,WinPrefs * Prefs)66 void Internals::WriteCoordinatesToFile(BufferFile * File, MoleculeData * MainData, WinPrefs * Prefs) {
67 if (!MOPacStyle) {
68 CreateMOPacInternals(3*MainData->GetMaximumAtomCount());
69 if (MOPacStyle) MOPacStyle->GuessInit(MainData);
70 }
71 if (MOPacStyle) MOPacStyle->WriteCoordinatesToFile(File, MainData, Prefs);
72 }
WriteMPCZMatCoordinatesToFile(BufferFile * File,MoleculeData * MainData,WinPrefs * Prefs)73 void Internals::WriteMPCZMatCoordinatesToFile(BufferFile * File, MoleculeData * MainData, WinPrefs * Prefs) {
74 if (!MOPacStyle) {
75 CreateMOPacInternals(3*MainData->GetMaximumAtomCount());
76 if (MOPacStyle) MOPacStyle->GuessInit(MainData);
77 }
78 if (MOPacStyle) MOPacStyle->WriteMPCZMatCoordinatesToFile(File, MainData, Prefs);
79 }
80 //Guess the connection atom set, then setup the values
GuessInit(MoleculeData * MainData,long theAtom,bool keepOld)81 void MOPacInternals::GuessInit(MoleculeData * MainData, long theAtom, bool keepOld) {
82 Frame * lFrame = MainData->GetCurrentFramePtr();
83 if (3*lFrame->NumAtoms > Allocation) return;
84 if (lFrame->NumAtoms < 2) return;
85 long StartAtom=1, EndAtom=lFrame->NumAtoms;
86 //Nothing to do for the first atom
87 if (theAtom == 0) return;
88 if (theAtom>0) {
89 StartAtom=theAtom;
90 EndAtom=StartAtom+1;
91 }
92 // First pair of atoms
93 CPoint3D BondVector, testVector;
94 long BondedAtom, AngleAtom, DihedralAtom, i;
95 float BondLength, SecondLength, ThirdLength, testLength;
96 for (long iatom=StartAtom; iatom<EndAtom; iatom++) {
97 //Find the closest atom to choose as the bonded atom
98 BondedAtom = 0;
99 BondVector.x = lFrame->Atoms[iatom].Position.x - lFrame->Atoms[0].Position.x;
100 BondVector.y = lFrame->Atoms[iatom].Position.y - lFrame->Atoms[0].Position.y;
101 BondVector.z = lFrame->Atoms[iatom].Position.z - lFrame->Atoms[0].Position.z;
102 BondLength = BondVector.Magnitude();
103 SecondLength = ThirdLength = 1.0E10;
104 AngleAtom = DihedralAtom = 0;
105 for (i=1; i<iatom; i++) {
106 testVector.x = lFrame->Atoms[iatom].Position.x - lFrame->Atoms[i].Position.x;
107 testVector.y = lFrame->Atoms[iatom].Position.y - lFrame->Atoms[i].Position.y;
108 testVector.z = lFrame->Atoms[iatom].Position.z - lFrame->Atoms[i].Position.z;
109 testLength = testVector.Magnitude();
110 if (testLength<BondLength) {
111 ThirdLength = SecondLength;
112 DihedralAtom = AngleAtom;
113 SecondLength = BondLength;
114 AngleAtom = BondedAtom;
115 BondedAtom = i;
116 BondLength = testLength;
117 BondVector = testVector;
118 } else if (testLength<SecondLength) {
119 ThirdLength = SecondLength;
120 DihedralAtom = AngleAtom;
121 SecondLength = testLength;
122 AngleAtom = i;
123 } else if (testLength<ThirdLength) {
124 ThirdLength = testLength;
125 DihedralAtom = i;
126 }
127 } //the three closest atoms are now known
128 /* Originally I simply used the closest 3 atoms as the connection atoms.
129 However, this apparently does not work so well for large molecules that
130 fold back on themselves, especially those with hydrogen bonds.
131 So instead try using the connection list for the atom defining the bond.
132 */
133 if (keepOld) {
134 if ((ConnectionAtoms[3*iatom]>=0)&&(ConnectionAtoms[3*iatom]<iatom)) {
135 if (AngleAtom == ConnectionAtoms[3*iatom]) AngleAtom = BondedAtom;
136 if (DihedralAtom == ConnectionAtoms[3*iatom]) DihedralAtom = BondedAtom;
137 BondedAtom = ConnectionAtoms[3*iatom];
138 }
139 }
140 //I ran into an issue in an Oh struture where the three defining atoms are linear.
141 //Needless to say this doesn't result in a good definition so test for that here.
142 //Get the angle between the three defining atoms, if its 0 or 180 they are linear.
143 float angle = 0.0f;
144 bool valid = lFrame->GetBondAngle(BondedAtom, AngleAtom, DihedralAtom, &angle);
145 if (valid && ((fabs(angle)>179.9)||(fabs(angle)<0.1))) {
146 float newLength = 1.0E10;
147 //find the closest atom that isn't linear with the other two
148 for (i=1; i<iatom; i++) {
149 if ((i!=BondedAtom)&&(i!=AngleAtom)&&(i!=DihedralAtom)) {
150 testVector.x = lFrame->Atoms[iatom].Position.x - lFrame->Atoms[i].Position.x;
151 testVector.y = lFrame->Atoms[iatom].Position.y - lFrame->Atoms[i].Position.y;
152 testVector.z = lFrame->Atoms[iatom].Position.z - lFrame->Atoms[i].Position.z;
153 testLength = testVector.Magnitude();
154 if (testLength < newLength) {
155 lFrame->GetBondAngle(BondedAtom, AngleAtom, i, &angle);
156 if (!((fabs(angle)>179.9)||(fabs(angle)<0.1))) {
157 DihedralAtom = i;
158 newLength = testLength;
159 }
160 }
161 }
162 }
163 }
164 ConnectionAtoms[3*iatom] = BondedAtom;
165 Values[3*iatom] = BondLength;
166 Type[3*iatom] = 0;
167 if (iatom >= 2) {
168 if (keepOld &&(ConnectionAtoms[3*iatom+1]>=0)&&(ConnectionAtoms[3*iatom+1]<iatom)&&
169 (ConnectionAtoms[3*iatom+1]!=BondedAtom)) {
170 if (DihedralAtom == ConnectionAtoms[3*iatom+1]) DihedralAtom = BondedAtom;
171 AngleAtom = ConnectionAtoms[3*iatom+1];
172 } else {
173 //Bond angle - atom -> BondedAtom -> AngleAtom
174 if ((BondedAtom > 0)&&(ConnectionAtoms[3*BondedAtom]!=BondedAtom)&&(lFrame->Atoms[iatom].Type>1)&&
175 (lFrame->Atoms[ConnectionAtoms[3*BondedAtom]].Type>1)) {
176 ConnectionAtoms[3*iatom+1] = ConnectionAtoms[3*BondedAtom];
177 if (DihedralAtom == ConnectionAtoms[3*BondedAtom]) DihedralAtom = AngleAtom;
178 AngleAtom = ConnectionAtoms[3*BondedAtom];
179 } else
180 ConnectionAtoms[3*iatom+1] = AngleAtom;
181 }
182 Type[3*iatom+1] = 1;
183 if (iatom >= 3) {
184 if (keepOld &&(ConnectionAtoms[3*iatom+2]>=0)&&(ConnectionAtoms[3*iatom+2]<iatom)&&
185 (ConnectionAtoms[3*iatom+2]!=BondedAtom)&&(ConnectionAtoms[3*iatom+2]!=AngleAtom)) {
186 if (DihedralAtom == ConnectionAtoms[3*iatom+1]) DihedralAtom = BondedAtom;
187 DihedralAtom = ConnectionAtoms[3*iatom+2];
188 } else {
189 //Dihedral Angle
190 if ((BondedAtom > 2)&&(ConnectionAtoms[3*BondedAtom+1]!=AngleAtom)&&(lFrame->Atoms[iatom].Type>1)&&
191 (lFrame->Atoms[ConnectionAtoms[3*BondedAtom+1]].Type>1))
192 ConnectionAtoms[3*iatom+2] = ConnectionAtoms[3*BondedAtom+1];
193 else {
194 if (DihedralAtom != ConnectionAtoms[3*iatom+1])
195 ConnectionAtoms[3*iatom+2] = DihedralAtom;
196 else
197 ConnectionAtoms[3*iatom+2] = AngleAtom;
198 }
199 }
200 Type[3*iatom+2] = 2;
201 }
202 }
203 }
204 if (theAtom > 0)
205 UpdateInternalValuesAtom(MainData, theAtom);
206 else { //update all the atoms
207 Count = 3*lFrame->NumAtoms;
208 CartesiansToInternals(MainData);
209 }
210 }
AddInternalCoordinate(long whichAtom,long connectedAtom,short type,float value)211 void MOPacInternals::AddInternalCoordinate(long whichAtom, long connectedAtom, short type, float value) {
212 if (3*whichAtom > Allocation) return;
213 if (connectedAtom > whichAtom) return; //forward references are not allowed
214 if ((type < 0)||(type>2)) return; //invalid type
215 if (3*whichAtom > Count) Count = 3*whichAtom;
216 ConnectionAtoms[3*whichAtom+type] = connectedAtom;
217 Values[3*whichAtom+type] = value;
218 }
219 //Take the given cartesians and connection list and recalculate
220 //the values of each internal coordinate
CartesiansToInternals(MoleculeData * MainData)221 void MOPacInternals::CartesiansToInternals(MoleculeData * MainData) {
222 Frame * lFrame = MainData->GetCurrentFramePtr();
223 if (3*lFrame->NumAtoms > Count) return;
224 if (lFrame->NumAtoms < 2) return;
225
226 for (long i=1; i<lFrame->NumAtoms; i++) {
227 UpdateInternalValuesAtom(MainData, i);
228 }
229 }
UpdateInternalValuesAtom(MoleculeData * MainData,long theAtom)230 void MOPacInternals::UpdateInternalValuesAtom(MoleculeData * MainData, long theAtom) {
231 Frame * lFrame = MainData->GetCurrentFramePtr();
232 if ((theAtom <= 0)||(theAtom*3 > Count)) return; //sanity check
233 long BondedAtom, AngleAtom, DihedralAtom;
234 float BondLength, BondAngle, Dihedral=0.0;
235 CPoint3D BondVector, offset2, offset3;
236 //Bond Length
237 BondedAtom = ConnectionAtoms[3*theAtom];
238 BondVector.x = lFrame->Atoms[BondedAtom].Position.x - lFrame->Atoms[theAtom].Position.x;
239 BondVector.y = lFrame->Atoms[BondedAtom].Position.y - lFrame->Atoms[theAtom].Position.y;
240 BondVector.z = lFrame->Atoms[BondedAtom].Position.z - lFrame->Atoms[theAtom].Position.z;
241 BondLength = BondVector.Magnitude();
242 Values[3*theAtom] = BondLength;
243 if (theAtom <= 1) return;
244 //Bond angle
245 AngleAtom = ConnectionAtoms[3*theAtom+1];
246 offset2.x = lFrame->Atoms[theAtom].Position.x - lFrame->Atoms[AngleAtom].Position.x;
247 offset2.y = lFrame->Atoms[theAtom].Position.y - lFrame->Atoms[AngleAtom].Position.y;
248 offset2.z = lFrame->Atoms[theAtom].Position.z - lFrame->Atoms[AngleAtom].Position.z;
249 float length3 = offset2.Magnitude();
250 offset2.x = lFrame->Atoms[AngleAtom].Position.x - lFrame->Atoms[BondedAtom].Position.x;
251 offset2.y = lFrame->Atoms[AngleAtom].Position.y - lFrame->Atoms[BondedAtom].Position.y;
252 offset2.z = lFrame->Atoms[AngleAtom].Position.z - lFrame->Atoms[BondedAtom].Position.z;
253 float length2 = offset2.Magnitude();
254 if ((fabs(BondLength)>0.0001)&&(fabs(length3)>0.0001)&&(fabs(length2)>0.0001)) {
255 float Radians = (BondLength*BondLength+length2*length2-length3*length3)/
256 (2*BondLength*length2);
257 //Make sure the angle is within the allowable values for acos (-1 to 1)
258 if (Radians > 1.0) Radians = 1.0;
259 else if (Radians < -1.0) Radians = -1.0;
260 BondAngle = acos(Radians);
261 BondAngle *= kRadToDegree;
262 } else BondAngle = 0.0;
263 Values[3*theAtom+1] = BondAngle;
264 //Dihedral angle
265 if (theAtom > 2) {
266 DihedralAtom = ConnectionAtoms[3*theAtom+2];
267 offset3.x = lFrame->Atoms[DihedralAtom].Position.x - lFrame->Atoms[AngleAtom].Position.x;
268 offset3.y = lFrame->Atoms[DihedralAtom].Position.y - lFrame->Atoms[AngleAtom].Position.y;
269 offset3.z = lFrame->Atoms[DihedralAtom].Position.z - lFrame->Atoms[AngleAtom].Position.z;
270 // float length3 = offset3.Magnitude();
271 CPoint3D UnitIJ = BondVector;
272 CPoint3D UnitJK = offset2;
273 CPoint3D UnitKL = offset3;
274 Normalize3D(&UnitIJ);
275 Normalize3D(&UnitJK);
276 Normalize3D(&UnitKL);
277 CPoint3D Normal1, Normal2;
278 CrossProduct3D(&UnitIJ, &UnitJK, &Normal1);
279 CrossProduct3D(&UnitJK, &UnitKL, &Normal2);
280 float DotPJ = DotProduct3D(&UnitIJ, &UnitJK);
281 float DotPK = DotProduct3D(&UnitJK, &UnitKL);
282 DotPJ = 1.0 - DotPJ*DotPJ;
283 DotPK = 1.0 - DotPK*DotPK;
284 if ((DotPJ > 0.0)||(DotPK > 0.0)) { //3 of the atom are linear, Bad!
285 float SinPJ = sqrt(DotPJ);
286 float SinPK = sqrt(DotPK);
287 float Dot = DotProduct3D(&Normal1, &Normal2)/(SinPJ*SinPK);
288 if (fabs(Dot) <= kCosErrorTolerance) { //Bad value for a cos
289 if (Dot > 0.9999) Dot = 1.0;
290 else if (Dot < -0.9999) Dot = -1.0;
291 Dihedral = acos(Dot);
292 float Pi = acos(-1.0);
293 if (fabs(Dihedral) < kZeroTolerance) Dihedral = 0.0;
294 else if (fabs(Dihedral-Pi) < kZeroTolerance) Dihedral = Pi;
295 float Sense = DotProduct3D(&Normal2, &BondVector);
296 if (Sense < 0.0) Dihedral = -Dihedral;
297 Dihedral *= 180.0/Pi;
298 }
299 Values[3*theAtom+2] = Dihedral;
300 }
301 }
302 }
303 //Take the connections list and values to create a set of cartesian coordinates
InternalsToCartesians(MoleculeData * MainData,WinPrefs * Prefs,long ChangedAtom)304 void MOPacInternals::InternalsToCartesians(MoleculeData * MainData, WinPrefs * Prefs,
305 long ChangedAtom) {
306 Frame * lFrame = MainData->GetCurrentFramePtr();
307 double DegToRad = acos(-1.0)/180.0;
308
309 if (lFrame->NumAtoms>1) {
310 mpAtom * tempAtoms = new mpAtom[lFrame->NumAtoms];
311 if (tempAtoms) {
312 long i;
313 for (i=0; i<lFrame->NumAtoms; i++) tempAtoms[i].Type = lFrame->Atoms[i].Type;
314 if (ChangedAtom < 3) {
315 tempAtoms[0].Position.x = tempAtoms[0].Position.y = tempAtoms[0].Position.z = 0;
316 tempAtoms[1].Position.x = tempAtoms[1].Position.y = 0.0;
317 tempAtoms[1].Position.z = Values[3];
318 } else {
319 tempAtoms[0].Position.x = lFrame->Atoms[0].Position.x;
320 tempAtoms[0].Position.y = lFrame->Atoms[0].Position.y;
321 tempAtoms[0].Position.z = lFrame->Atoms[0].Position.z;
322 tempAtoms[1].Position.x = lFrame->Atoms[1].Position.x;
323 tempAtoms[1].Position.y = lFrame->Atoms[1].Position.y;
324 tempAtoms[1].Position.z = lFrame->Atoms[1].Position.z;
325 }
326 if (lFrame->NumAtoms > 2) {
327 if (ChangedAtom < 3) {
328 tempAtoms[2].Position.x = Values[6]*sin(Values[7]*DegToRad);
329 tempAtoms[2].Position.y = 0.0;
330 if (ConnectionAtoms[6]==0)
331 tempAtoms[2].Position.z = Values[6]*cos(Values[7]*DegToRad);
332 else tempAtoms[2].Position.z = tempAtoms[1].Position.z -
333 Values[6]*cos(Values[7]*DegToRad);
334 i=3; //add atoms using simple formula as long as molecule is linear
335 while ((tempAtoms[i-1].Position.x <= 1.0e-5)&&(i<lFrame->NumAtoms)) {
336 tempAtoms[i].Position.x = Values[3*i]*sin(Values[3*i+1]*DegToRad);
337 tempAtoms[i].Position.y = 0.0;
338 float test = tempAtoms[ConnectionAtoms[3*i]].Position.z -
339 tempAtoms[ConnectionAtoms[3*i+1]].Position.z;
340 tempAtoms[i].Position.z = tempAtoms[ConnectionAtoms[3*i]].Position.z -
341 Values[3*i]*cos(Values[3*i+1]*DegToRad)*((test>=0.0)?1.0:0.0);
342 i++;
343 }
344 } else {
345 for (long j=2; j<ChangedAtom; j++) {
346 tempAtoms[j].Position.x = lFrame->Atoms[j].Position.x;
347 tempAtoms[j].Position.y = lFrame->Atoms[j].Position.y;
348 tempAtoms[j].Position.z = lFrame->Atoms[j].Position.z;
349 }
350 i = ChangedAtom;
351 } //Now for general 3D molecules
352 for (; i<lFrame->NumAtoms; i++) {
353 long j, k, l;
354 j = ConnectionAtoms[3*i];
355 k = ConnectionAtoms[3*i+1];
356 l = ConnectionAtoms[3*i+2];
357 CPoint3D Vkl, Vjk;
358 Vkl.x = tempAtoms[k].Position.x - tempAtoms[l].Position.x;
359 Vkl.y = tempAtoms[k].Position.y - tempAtoms[l].Position.y;
360 Vkl.z = tempAtoms[k].Position.z - tempAtoms[l].Position.z;
361 CPoint3D Unitkl = Vkl;
362 Normalize3D(&Unitkl);
363 Vjk.x = tempAtoms[j].Position.x - tempAtoms[k].Position.x;
364 Vjk.y = tempAtoms[j].Position.y - tempAtoms[k].Position.y;
365 Vjk.z = tempAtoms[j].Position.z - tempAtoms[k].Position.z;
366 CPoint3D Unitjk = Vjk;
367 Normalize3D(&Unitjk);
368 CPoint3D KLxJK;
369 CrossProduct3D(&Unitkl, &Unitjk, &KLxJK);
370 CPoint3D A3;
371 double klDotjk = DotProduct3D(&Unitkl, &Unitjk);
372 klDotjk = 1.0 - klDotjk*klDotjk;
373 double SinA3;
374 if (klDotjk > 0.0001)
375 SinA3 = sqrt(klDotjk);
376 else
377 SinA3 = 0.001;
378 A3.x = KLxJK.x/SinA3;
379 A3.y = KLxJK.y/SinA3;
380 A3.z = KLxJK.z/SinA3;
381 CPoint3D A4;
382 CrossProduct3D(&A3, &Unitjk, &A4);
383 CPoint3D AtomVector;
384 AtomVector.x = Values[3*i]*((-Unitjk.x*cos(Values[3*i+1]*DegToRad)) +
385 A4.x*sin(Values[3*i+1]*DegToRad)*cos(Values[3*i+2]*DegToRad) +
386 A3.x*sin(Values[3*i+1]*DegToRad)*sin(Values[3*i+2]*DegToRad));
387 AtomVector.y = Values[3*i]*((-Unitjk.y*cos(Values[3*i+1]*DegToRad)) +
388 A4.y*sin(Values[3*i+1]*DegToRad)*cos(Values[3*i+2]*DegToRad) +
389 A3.y*sin(Values[3*i+1]*DegToRad)*sin(Values[3*i+2]*DegToRad));
390 AtomVector.z = Values[3*i]*((-Unitjk.z*cos(Values[3*i+1]*DegToRad)) +
391 A4.z*sin(Values[3*i+1]*DegToRad)*cos(Values[3*i+2]*DegToRad) +
392 A3.z*sin(Values[3*i+1]*DegToRad)*sin(Values[3*i+2]*DegToRad));
393 tempAtoms[i].Position.x = tempAtoms[j].Position.x + AtomVector.x;
394 tempAtoms[i].Position.y = tempAtoms[j].Position.y + AtomVector.y;
395 tempAtoms[i].Position.z = tempAtoms[j].Position.z + AtomVector.z;
396 }
397 } //rotate/translate the new coordinates to match the existing coords
398 //note this only attempts to match the first four atoms, the others
399 //will be allowed to wander
400 // MinimizeDifferences(lFrame->Atoms, tempAtoms, lFrame->NumAtoms, Prefs,
401 // MAX(lFrame->NumAtoms, 4));
402 if (ChangedAtom < 4)
403 MinimizeDifferences(lFrame->Atoms, tempAtoms, lFrame->NumAtoms, Prefs,
404 lFrame->NumAtoms);
405 for (i=0; i<lFrame->NumAtoms; i++)
406 lFrame->Atoms[i].Position = tempAtoms[i].Position;
407 delete [] tempAtoms;
408 }
409 }
410 }
411
AppendAtom(MoleculeData * MainData)412 void MOPacInternals::AppendAtom(MoleculeData * MainData) {
413 Frame * lFrame = MainData->GetCurrentFramePtr();
414 if (3*lFrame->NumAtoms > Allocation) {
415 long * templ = new long[3*lFrame->NumAtoms+12];
416 if (templ) {
417 memcpy(templ, ConnectionAtoms, Count*sizeof(long));
418 delete [] ConnectionAtoms;
419 ConnectionAtoms = templ;
420 } else throw MemoryError();
421 float * tempVal = new float[3*lFrame->NumAtoms+12];
422 if (tempVal) {
423 memcpy(tempVal, Values, Count*sizeof(float));
424 delete [] Values;
425 Values = tempVal;
426 } else throw MemoryError();
427 char * tempType = new char [3*lFrame->NumAtoms+12];
428 if (tempType) {
429 memcpy(tempType, Type, Count*sizeof(char));
430 delete [] Type;
431 Type = tempType;
432 } else throw MemoryError();
433 Allocation = 3*lFrame->NumAtoms + 12;
434 }
435 //Now add the new atoms connections to the end of the list - valid since
436 //new atoms are always appended to the end of the atom list
437 long newAtom = Count/3;
438 Count += 3;
439 GuessInit(MainData, newAtom, false);
440 }
DeleteAtom(MoleculeData * MainData,long WhichAtom)441 void MOPacInternals::DeleteAtom(MoleculeData * MainData, long WhichAtom) {
442 Frame * lFrame = MainData->GetCurrentFramePtr();
443 for (long i=WhichAtom; i<lFrame->NumAtoms; i++) { //pull down the atoms higher in the list
444 bool update = false;
445 for (long j=0; j<3; j++) {
446 ConnectionAtoms[3*i+j] = ConnectionAtoms[3*i+j+3];
447 Values[3*i+j] = Values[3*i+j+3];
448 Type[3*i+j] = Type[3*i+j+3];
449 if (ConnectionAtoms[3*i+j] == WhichAtom) {
450 //Have to change the reference to a different atom
451 update = true;
452 ConnectionAtoms[3*i+j] = -1;
453 } else if (ConnectionAtoms[3*i+j] > WhichAtom) ConnectionAtoms[3*i+j]--; //reduce connection atom count
454 }
455 if (update)
456 GuessInit(MainData, i, true);
457 }
458 if (Count > 0) Count -= 3; //reduce the number of internal coordinates
459 }
UpdateAtoms(MoleculeData * MainData)460 void MOPacInternals::UpdateAtoms(MoleculeData * MainData) {
461 Frame * lFrame = MainData->GetCurrentFramePtr();
462 while (lFrame->NumAtoms*3 > Count) { //More atoms, just append atoms to the end of the list
463 AppendAtom(MainData);
464 }
465 while (lFrame->NumAtoms*3 < Count) {
466 DeleteAtom(MainData, lFrame->NumAtoms+1);
467 }
468 }
ChangeAtomIndex(MoleculeData * MainData,long OldPosition,long NewPosition)469 void MOPacInternals::ChangeAtomIndex(MoleculeData * MainData, long OldPosition, long NewPosition) {
470 //rearrange an atom in the list. There are two cases, moving an atom up in the list
471 //and pulling it down. Moving up is always safe, just need to update all indexes and possibly the atom
472 //that is being moved.
473 //Moving down (later) in the list may force some atom references to change since the reference
474 //must be to an atom earlier in the list. (noop if the positions are the same)
475 if (NewPosition < OldPosition) {
476 long i0 = ConnectionAtoms[OldPosition*3];
477 long i1 = ConnectionAtoms[OldPosition*3 + 1];
478 long i2 = ConnectionAtoms[OldPosition*3 + 2];
479 float v0 = Values[OldPosition*3];
480 float v1 = Values[OldPosition*3 + 1];
481 float v2 = Values[OldPosition*3 + 2];
482 char t0 = Type[OldPosition*3];
483 char t1 = Type[OldPosition*3 + 1];
484 char t2 = Type[OldPosition*3 + 2];
485 for (long i=OldPosition; i>NewPosition; i--) {
486 ConnectionAtoms[i*3] = ConnectionAtoms[(i-1)*3];
487 //update the index to the new index for the same atom.
488 if (ConnectionAtoms[i*3] == OldPosition) ConnectionAtoms[3*i] = NewPosition;
489 else if ((ConnectionAtoms[i*3] >= NewPosition)&&(ConnectionAtoms[i*3] < OldPosition)) ConnectionAtoms[i*3] ++;
490 ConnectionAtoms[i*3+1] = ConnectionAtoms[(i-1)*3+1];
491 if (ConnectionAtoms[i*3+1] == OldPosition) ConnectionAtoms[3*i+1] = NewPosition;
492 else if ((ConnectionAtoms[i*3+1] >= NewPosition)&&(ConnectionAtoms[i*3+1] < OldPosition)) ConnectionAtoms[i*3+1] ++;
493 ConnectionAtoms[i*3+2] = ConnectionAtoms[(i-1)*3+2];
494 if (ConnectionAtoms[i*3+2] == OldPosition) ConnectionAtoms[3*i+2] = NewPosition;
495 else if ((ConnectionAtoms[i*3+2] >= NewPosition)&&(ConnectionAtoms[i*3+2] < OldPosition)) ConnectionAtoms[i*3+2] ++;
496 Values[i*3] = Values[(i-1)*3];
497 Values[i*3+1] = Values[(i-1)*3+1];
498 Values[i*3+2] = Values[(i-1)*3+2];
499 Type[i*3] = Type[(i-1)*3];
500 Type[i*3+1] = Type[(i-1)*3+1];
501 Type[i*3+2] = Type[(i-1)*3+2];
502 if ((i-1)<3) //The first three atoms do not have fully validated lists yet
503 GuessInit(MainData, i, true);
504 }
505 ConnectionAtoms[NewPosition*3] = i0;
506 ConnectionAtoms[NewPosition*3+1] = i1;
507 ConnectionAtoms[NewPosition*3+2] = i2;
508 Values[NewPosition*3] = v0;
509 Values[NewPosition*3+1] = v1;
510 Values[NewPosition*3+2] = v2;
511 Type[NewPosition*3] = t0;
512 Type[NewPosition*3+1] = t1;
513 Type[NewPosition*3+2] = t2;
514 bool update=false;
515 //make sure the moved atoms references point to atoms prior to it in the list
516 if (ConnectionAtoms[3*NewPosition] >= NewPosition) {
517 ConnectionAtoms[3*NewPosition] = -1;
518 update = true;
519 }
520 if ((ConnectionAtoms[3*NewPosition+1] >= NewPosition)||(ConnectionAtoms[3*NewPosition]==ConnectionAtoms[3*NewPosition+1])) {
521 ConnectionAtoms[3*NewPosition+1] = -1;
522 update = true;
523 }
524 if ((ConnectionAtoms[3*NewPosition+2] >= NewPosition)||
525 (ConnectionAtoms[3*NewPosition]==ConnectionAtoms[3*NewPosition+2])||
526 (ConnectionAtoms[3*NewPosition+1]==ConnectionAtoms[3*NewPosition+2])) {
527 ConnectionAtoms[3*NewPosition+2] = -1;
528 update = true;
529 }
530 if (update) GuessInit(MainData, NewPosition, true);
531 //Now update references for atoms higher up in the list
532 for (long i=(OldPosition+1); (i*3)<Count; i++) {
533 for (long j=0; j<3; j++) {
534 if (ConnectionAtoms[3*i+j] == OldPosition)
535 ConnectionAtoms[3*i+j] = NewPosition;
536 else if ((ConnectionAtoms[3*i+j] >= NewPosition) &&
537 (ConnectionAtoms[3*i+j] < OldPosition))
538 ConnectionAtoms[3*i+j]++;
539 }
540 }
541 } else if (NewPosition > OldPosition) {
542 //When moving down the list the atom being moved does not change, but any atom that references it
543 //before the NewPosition will have to be changed.
544 long i0 = ConnectionAtoms[OldPosition*3];
545 long i1 = ConnectionAtoms[OldPosition*3 + 1];
546 long i2 = ConnectionAtoms[OldPosition*3 + 2];
547 float v0 = Values[OldPosition*3];
548 float v1 = Values[OldPosition*3 + 1];
549 float v2 = Values[OldPosition*3 + 2];
550 char t0 = Type[OldPosition*3];
551 char t1 = Type[OldPosition*3 + 1];
552 char t2 = Type[OldPosition*3 + 2];
553 for (long i=OldPosition; i<NewPosition; i++) {
554 bool update = false;
555 ConnectionAtoms[i*3] = ConnectionAtoms[(i+1)*3];
556 ConnectionAtoms[i*3+1] = ConnectionAtoms[(i+1)*3+1];
557 ConnectionAtoms[i*3+2] = ConnectionAtoms[(i+1)*3+2];
558 Values[i*3] = Values[(i+1)*3];
559 Values[i*3+1] = Values[(i+1)*3+1];
560 Values[i*3+2] = Values[(i+1)*3+2];
561 Type[i*3] = Type[(i+1)*3];
562 Type[i*3+1] = Type[(i+1)*3+1];
563 Type[i*3+2] = Type[(i+1)*3+2];
564 //update the index to the new index for the same atom.
565 if (ConnectionAtoms[i*3] == OldPosition) {
566 //Since the newposition is later in the list this reference must change
567 ConnectionAtoms[3*i] = -1;
568 update = true;
569 } else if ((ConnectionAtoms[i*3] > OldPosition)&&(ConnectionAtoms[i*3] <= NewPosition))
570 ConnectionAtoms[i*3] --;
571 if (ConnectionAtoms[i*3+1] == OldPosition) {
572 ConnectionAtoms[3*i+1] = -1;
573 update = true;
574 } else if ((ConnectionAtoms[i*3+1] > OldPosition)&&(ConnectionAtoms[i*3+1] <= NewPosition))
575 ConnectionAtoms[i*3+1] --;
576 if (ConnectionAtoms[i*3+2] == OldPosition) {
577 ConnectionAtoms[3*i+2] = -1;
578 update = true;
579 } else if ((ConnectionAtoms[i*3+2] > OldPosition)&&(ConnectionAtoms[i*3+2] <= NewPosition))
580 ConnectionAtoms[i*3+2] --;
581
582 if (update) {
583 //A reference has changed. Confirm that the three index references are unique
584 //then compute the new values for this atom.
585 GuessInit(MainData, i, true);
586 }
587 }
588 ConnectionAtoms[NewPosition*3] = i0;
589 ConnectionAtoms[NewPosition*3+1] = i1;
590 ConnectionAtoms[NewPosition*3+2] = i2;
591 Values[NewPosition*3] = v0;
592 Values[NewPosition*3+1] = v1;
593 Values[NewPosition*3+2] = v2;
594 Type[NewPosition*3] = t0;
595 Type[NewPosition*3+1] = t1;
596 Type[NewPosition*3+2] = t2;
597 //Now update references for atoms higher up in the list
598 for (long i=(NewPosition+1); (i*3)<Count; i++) {
599 for (long j=0; j<3; j++) {
600 if (ConnectionAtoms[3*i+j] == OldPosition)
601 ConnectionAtoms[3*i+j] = NewPosition;
602 else if ((ConnectionAtoms[3*i+j] <= NewPosition) &&
603 (ConnectionAtoms[3*i+j] > OldPosition))
604 ConnectionAtoms[3*i+j]--;
605 }
606 }
607 }
608 }
609