1 /*
2 * (c) 2004 Iowa State University
3 * see the LICENSE file in the top level directory
4 */
5
6 /* Frame.cpp
7
8 Class to organize data which is specific to one geometry point
9
10 BMB 4/1998
11 Added Bond Order and Hydrogen Bond generation 12/1998
12 Fixed problem reading incomplete sets of UHF orbitals in ParseEigenVectors - 10/2000 BMB
13 Moved normal mode parsing to ParseNormalModes - 1/2001 BMB
14 Updated CI orbital parser for addition of sym - 8/2002 BMB
15 Adjusted the MP2 orbital dat file parser for reduced variational space - 9/2002 BMB
16 */
17
18 #include "Globals.h"
19 #include "GlobalExceptions.h"
20 #include "MyTypes.h"
21 #include "myFiles.h"
22 #include "Frame.h"
23 #include "Progress.h"
24 #include "Prefs.h"
25 #include "Gradient.h"
26 #include "AtomTypeList.h"
27 #include "wx/wx.h"
28 #include <string.h>
29 #include <new>
30 #include <ctype.h>
31 #include "MolDisplayWin.h"
32
33 #if defined(WIN32)
34 #undef AddAtom
35 #endif
36
37 extern WinPrefs *gPreferences;
38
Frame(MolDisplayWin * myMolWin)39 Frame::Frame(MolDisplayWin * myMolWin) {
40 MolWin = myMolWin;
41 Energy = 0.0;
42 time = 0.0;
43 IRCPt = 0;
44 Atoms = NULL;
45 Bonds = NULL;
46 NumAtoms = 0;
47 AtomAllocation = 0;
48 NumBonds = 0;
49 BondAllocation = 0;
50 SpecialAtoms = NULL;
51 Vibs = NULL;
52 SurfaceList = NULL;
53 Gradient = NULL;
54
55 NextFrame = NULL;
56 PreviousFrame = NULL;
57
58 natoms_selected = 0;
59 targeted_atom = -1;
60 }
61
62 /*
63 Frame& Frame::operator= (const Frame& f)
64 {
65 delete [] Atoms;
66 Atoms = NULL;
67 delete [] Bonds;
68 Bonds = NULL;
69
70 this->NumAtoms = 0;
71 this->NumBonds = 0;
72 this->AtomAllocation = 0;
73 this->BondAllocation = 0;
74
75 for ( int i = 0; i < f.NumAtoms; i++)
76 AddAtom(f.Atoms[i].Type, f.Atoms[i].Position);
77
78 for ( int i = 0; i < f.NumBonds; i++)
79 AddBond(f.Bonds[i].Atom1, f.Bonds[i].Atom2, f.Bonds[i].Order);
80
81 return *this;
82 }
83 */
84 //not needed currently
85
~Frame(void)86 Frame::~Frame(void) {
87 if (Atoms) delete [] Atoms;
88 if (Bonds) delete [] Bonds;
89 if (SpecialAtoms) delete [] SpecialAtoms;
90 if (Vibs) delete Vibs;
91 if (Orbs.size() > 0) {
92 std::vector<OrbitalRec *>::const_iterator OrbSet = Orbs.begin();
93 while (OrbSet != Orbs.end()) {
94 delete (*OrbSet);
95 OrbSet++;
96 }
97 }
98 while (SurfaceList) {
99 Surface * temp = SurfaceList->NextSurface;
100 delete SurfaceList;
101 SurfaceList = temp;
102 }
103 if (Gradient) delete Gradient;
104 if (NextFrame) {
105 if (PreviousFrame) {
106 PreviousFrame->SetNextFrame(NextFrame);
107 }
108 NextFrame->SetPreviousFrame(PreviousFrame);
109 } else if (PreviousFrame) {
110 PreviousFrame->SetNextFrame(NULL);
111 }
112 }
DeleteOrbitals(void)113 void Frame::DeleteOrbitals(void) {
114 if (Orbs.size() > 0) {
115 std::vector<OrbitalRec *>::const_iterator OrbSet = Orbs.begin();
116 while (OrbSet != Orbs.end()) {
117 delete (*OrbSet);
118 OrbSet++;
119 }
120 Orbs.clear();
121 }
122 }
123
SetNextFrame(Frame * next)124 void Frame::SetNextFrame(Frame * next) { NextFrame = next; }
125
SetPreviousFrame(Frame * previous)126 void Frame::SetPreviousFrame(Frame * previous) { PreviousFrame = previous; }
127
GetNextFrame(void)128 Frame * Frame::GetNextFrame(void) { return NextFrame; }
129
GetPreviousFrame(void)130 Frame * Frame::GetPreviousFrame(void) { return PreviousFrame; }
131
AddAtom(long AtomType,const CPoint3D & AtomPosition,long index)132 mpAtom *Frame::AddAtom(long AtomType, const CPoint3D & AtomPosition,
133 long index) {
134
135 /* AtomType is the atom's atomic number, starting with hydrogen at 1. It is used to index the
136 atom types array elsewhere so must be validated. */
137
138 mpAtom * result = NULL;
139
140 if (NumAtoms>=AtomAllocation) IncreaseAtomAllocation(MAX(NumAtoms,10));
141
142 if (NumAtoms<AtomAllocation) {
143 if ((index<=-1)||(index>=NumAtoms)) {//Add to the end of the list
144 index = NumAtoms;
145 } else { //insert the atom into the middle of the list
146 for (int i=NumAtoms; i>index; i--) {
147 Atoms[i] = Atoms[i-1];
148 }
149
150 // Adjust bonds that connect higher-numbered atoms.
151 for (int i = 0; i < NumBonds; ++i) {
152 if (Bonds[i].Atom1 >= index) Bonds[i].Atom1++;
153 if (Bonds[i].Atom2 >= index) Bonds[i].Atom2++;
154 }
155
156 if (targeted_atom >= index) {
157 targeted_atom++;
158 }
159 }
160 if ((AtomType>0)&&(AtomType<=kMaxAtomTypes))
161 Atoms[index].Type = AtomType;
162 else
163 Atoms[index].Type = 114; //Convert invalid types to a max type to clearly indicate a problem
164 Atoms[index].Position = AtomPosition;
165 Atoms[index].flags = 0;
166 Atoms[index].SetDefaultCoordinationNumber();
167 Atoms[index].SetDefaultLonePairCount();
168 result = &Atoms[index];
169 NumAtoms++;
170 }
171 //Delete any orbitals and normal modes
172 if (Vibs) {
173 delete Vibs;
174 Vibs = NULL;
175 }
176 DeleteOrbitals();
177 while (SurfaceList) DeleteSurface(0);
178 return result;
179
180 }
181
AddAtom(const mpAtom & atm,long index,const CPoint3D * pos)182 mpAtom * Frame::AddAtom(const mpAtom& atm, long index, const CPoint3D *pos) {
183
184 mpAtom * result = NULL;
185
186 if (NumAtoms>=AtomAllocation) IncreaseAtomAllocation(MAX(NumAtoms,10));
187
188 if (NumAtoms<AtomAllocation) {
189 if ((index<=-1)||(index>=NumAtoms)) {//Add to the end of the list
190 index = NumAtoms;
191 } else { //insert the atom into the middle of the list
192 for (int i=NumAtoms; i>index; i--) {
193 Atoms[i] = Atoms[i-1];
194 }
195
196 // Adjust bonds that connect higher-numbered atoms.
197 for (int i = 0; i < NumBonds; ++i) {
198 if (Bonds[i].Atom1 >= index) Bonds[i].Atom1++;
199 if (Bonds[i].Atom2 >= index) Bonds[i].Atom2++;
200 }
201
202 if (targeted_atom >= index) {
203 targeted_atom++;
204 }
205 }
206 Atoms[index] = atm;
207 result = &Atoms[index];
208 if (pos) {
209 SetAtomPosition(index, *pos);
210 }
211 if (atm.GetSelectState()) {
212 natoms_selected++;
213 }
214 NumAtoms++;
215 }
216 //Delete any orbitals and normal modes
217 if (Vibs) {
218 delete Vibs;
219 Vibs = NULL;
220 }
221 DeleteOrbitals();
222 while (SurfaceList) DeleteSurface(0);
223 return result;
224 }
225
IncreaseAtomAllocation(long NumAdditional)226 bool Frame::IncreaseAtomAllocation(long NumAdditional) {
227 if (AtomAllocation+NumAdditional < NumAtoms) return false;
228 mpAtom * temp = new mpAtom[AtomAllocation+NumAdditional];
229 if (temp) {
230 if (Atoms != NULL) {
231 memcpy(temp, Atoms, NumAtoms*sizeof(mpAtom));
232 delete [] Atoms;
233 }
234 Atoms = temp;
235 AtomAllocation += NumAdditional;
236 } else return false;
237 return true;
238 }
ReadGradient(BufferFile * Buffer,wxFileOffset SearchLength)239 void Frame::ReadGradient(BufferFile * Buffer, wxFileOffset SearchLength) {
240 wxFileOffset SavedPos = Buffer->GetFilePos(), npos;
241 bool found=false;
242 int Style = 1;
243
244 if (Buffer->LocateKeyWord(" NSERCH", 7, SearchLength)) { //Search for gradient data at the end of an optimization point
245 found = true;
246 } else if (Buffer->LocateKeyWord("UNITS ARE HARTREE/BOHR E'X", 29, SearchLength)) {
247 Buffer->BackupnLines(6);
248 Style = 2;
249 found = true;
250 }
251 if (found) {
252 Gradient = new GradientData;
253 if (Gradient) {
254 npos = Buffer->GetFilePos();
255 if (Buffer->LocateKeyWord("COORDINATES (BOHR) GRADIENT (HARTREE/BOHR)",
256 60, SearchLength)) {
257 Style=0;
258 Buffer->SetFilePos(npos);
259 }
260 if (!Gradient->ParseGAMESSGradient(Buffer, NumAtoms, SearchLength, Style)) {
261 delete Gradient;
262 Gradient = NULL;
263 }
264 }
265 }
266
267 Buffer->SetFilePos(SavedPos); //reset position so other code isn't broken
268 }
GradientVectorAvailable(void) const269 bool Frame::GradientVectorAvailable(void) const {
270 bool result = false;
271 if (Gradient) {
272 result = Gradient->GradientVectorAvailable(NumAtoms);
273 }
274 return result;
275 }
RetrieveAtomGradient(long theAtom,CPoint3D & gradVector) const276 bool Frame::RetrieveAtomGradient(long theAtom, CPoint3D & gradVector) const {
277 if (Gradient) return Gradient->RetrieveAtomGradient(theAtom, gradVector);
278 return false;
279 }
280
GetRMSGradient(void) const281 float Frame::GetRMSGradient(void) const {
282 float result=-1.0f;
283 if (Gradient) {
284 result = Gradient->GetRMS();
285 }
286 return result;
287 }
GetMaxGradient(void) const288 float Frame::GetMaxGradient(void) const {
289 float result=-1.0;
290 if (Gradient) {
291 result = Gradient->GetMaximum();
292 }
293 return result;
294 }
SetRMSGradient(float v)295 void Frame::SetRMSGradient(float v) {
296 if (!Gradient) Gradient = new GradientData;
297 if (Gradient) Gradient->SetRMS(v);
298 }
SetMaximumGradient(float v)299 void Frame::SetMaximumGradient(float v) {
300 if (!Gradient) Gradient = new GradientData;
301 if (Gradient) Gradient->SetMaximum(v);
302 }
AddSpecialAtom(CPoint3D Vector,long AtomNum)303 bool Frame::AddSpecialAtom(CPoint3D Vector, long AtomNum) {
304 if (!SpecialAtoms) SpecialAtoms = new CPoint3D[AtomAllocation];
305 if (!SpecialAtoms) return false;
306 if (AtomNum >= AtomAllocation) return false;
307 SpecialAtoms[AtomNum] = Vector;
308 return true;
309 }
AddBond(long Atom1,long Atom2,const BondOrder & b)310 bool Frame::AddBond(long Atom1, long Atom2, const BondOrder & b) {
311 bool result = false;
312 //Validate the pair of atom references
313 if ((Atom1>=0)&&(Atom2>=0)&&(Atom1!=Atom2)&&(Atom1<NumAtoms)&&(Atom2<NumAtoms)) {
314 if (NumBonds == BondAllocation) IncreaseBondAllocation(MAX(10, NumBonds/10));
315 if (NumBonds<BondAllocation) {
316 //To improve scanning speed always place the higher atom # first
317 if (Atom1 > Atom2) {
318 Bonds[NumBonds].Atom1 = Atom1;
319 Bonds[NumBonds].Atom2 = Atom2;
320 } else {
321 Bonds[NumBonds].Atom1 = Atom2;
322 Bonds[NumBonds].Atom2 = Atom1;
323 }
324 Bonds[NumBonds].Order = b;
325 Bonds[NumBonds].Highlite = 0;
326 NumBonds++;
327 result = true;
328 }
329 }
330 return result;
331 }
332
GetBondLength(long atom1,long atom2,float * length)333 bool Frame::GetBondLength(long atom1, long atom2, float * length) {
334
335 bool result = false;
336
337 if ((atom1 >= 0)&&(atom2>=0)&&(atom1<NumAtoms)&&(atom2<NumAtoms)&&(atom1!=atom2)) {
338 CPoint3D offset;
339 offset.x = Atoms[atom1].Position.x - Atoms[atom2].Position.x;
340 offset.y = Atoms[atom1].Position.y - Atoms[atom2].Position.y;
341 offset.z = Atoms[atom1].Position.z - Atoms[atom2].Position.z;
342 *length = offset.Magnitude();
343 result = true;
344 }
345
346 return result;
347 }
GetBondAngle(long atom1,long atom2,long atom3,float * result)348 bool Frame::GetBondAngle(long atom1, long atom2, long atom3, float * result) {
349 bool res = false;
350
351 if ((atom1 >= 0)&&(atom2>=0)&&(atom3>=0)&&(atom1<NumAtoms)&&(atom2<NumAtoms)&&
352 (atom3<NumAtoms)&&(atom1!=atom2)&&(atom1!=atom3)&&(atom2!=atom3)) {
353
354 CPoint3D offset;
355 float length1, length2, length3;
356 offset.x = Atoms[atom2].Position.x - Atoms[atom1].Position.x;
357 offset.y = Atoms[atom2].Position.y - Atoms[atom1].Position.y;
358 offset.z = Atoms[atom2].Position.z - Atoms[atom1].Position.z;
359 length1 = offset.Magnitude();
360 offset.x = Atoms[atom1].Position.x - Atoms[atom3].Position.x;
361 offset.y = Atoms[atom1].Position.y - Atoms[atom3].Position.y;
362 offset.z = Atoms[atom1].Position.z - Atoms[atom3].Position.z;
363 length3 = offset.Magnitude();
364 offset.x = Atoms[atom3].Position.x - Atoms[atom2].Position.x;
365 offset.y = Atoms[atom3].Position.y - Atoms[atom2].Position.y;
366 offset.z = Atoms[atom3].Position.z - Atoms[atom2].Position.z;
367 length2 = offset.Magnitude();
368 if ((fabs(length1)>0.0001)&&(fabs(length3)>0.0001)&&(fabs(length2)>0.0001)) {
369 *result = acos(((length1*length1+length2*length2-length3*length3)/
370 (2*length1*length2)));
371 *result *= kRadToDegree;
372 } else *result = 0.0;
373 res = true;
374 }
375 return res;
376 }
GetBondDihedral(long atom1,long bondAtom,long AngleAtom,long DihedralAtom,float * angle)377 bool Frame::GetBondDihedral(long atom1, long bondAtom, long AngleAtom, long DihedralAtom,
378 float * angle) {
379 bool result = false;
380
381 if ((atom1 >= 0)&&(bondAtom>=0)&&(AngleAtom>=0)&&(DihedralAtom>=0)&&(atom1<NumAtoms)&&
382 (bondAtom<NumAtoms)&&(AngleAtom<NumAtoms)&&(DihedralAtom<NumAtoms)&&
383 (atom1!=bondAtom)&&(atom1!=AngleAtom)&&(atom1!=DihedralAtom)&&
384 (bondAtom!=AngleAtom)&&(bondAtom!=DihedralAtom)&&(AngleAtom!=DihedralAtom)) {
385
386 CPoint3D BondVector, offset2, offset3;
387 //Bond Length
388 BondVector.x = Atoms[bondAtom].Position.x - Atoms[atom1].Position.x;
389 BondVector.y = Atoms[bondAtom].Position.y - Atoms[atom1].Position.y;
390 BondVector.z = Atoms[bondAtom].Position.z - Atoms[atom1].Position.z;
391 //Bond angle
392 offset2.x = Atoms[atom1].Position.x - Atoms[AngleAtom].Position.x;
393 offset2.y = Atoms[atom1].Position.y - Atoms[AngleAtom].Position.y;
394 offset2.z = Atoms[atom1].Position.z - Atoms[AngleAtom].Position.z;
395 offset2.x = Atoms[AngleAtom].Position.x - Atoms[bondAtom].Position.x;
396 offset2.y = Atoms[AngleAtom].Position.y - Atoms[bondAtom].Position.y;
397 offset2.z = Atoms[AngleAtom].Position.z - Atoms[bondAtom].Position.z;
398 //Dihedral angle
399 offset3.x = Atoms[DihedralAtom].Position.x - Atoms[AngleAtom].Position.x;
400 offset3.y = Atoms[DihedralAtom].Position.y - Atoms[AngleAtom].Position.y;
401 offset3.z = Atoms[DihedralAtom].Position.z - Atoms[AngleAtom].Position.z;
402 // float length3 = offset3.Magnitude();
403 CPoint3D UnitIJ = BondVector;
404 CPoint3D UnitJK = offset2;
405 CPoint3D UnitKL = offset3;
406 Normalize3D(&UnitIJ);
407 Normalize3D(&UnitJK);
408 Normalize3D(&UnitKL);
409 CPoint3D Normal1, Normal2;
410 CrossProduct3D(&UnitIJ, &UnitJK, &Normal1);
411 CrossProduct3D(&UnitJK, &UnitKL, &Normal2);
412 float DotPJ = DotProduct3D(&UnitIJ, &UnitJK);
413 float DotPK = DotProduct3D(&UnitJK, &UnitKL);
414 DotPJ = 1.0 - DotPJ*DotPJ;
415 DotPK = 1.0 - DotPK*DotPK;
416 if ((DotPJ > 0.0)||(DotPK > 0.0)) { //3 of the atom are linear, Bad!
417 float SinPJ = sqrt(DotPJ);
418 float SinPK = sqrt(DotPK);
419 float Dot = DotProduct3D(&Normal1, &Normal2)/(SinPJ*SinPK);
420 if (fabs(Dot) <= kCosErrorTolerance) { //Bad value for a cos
421 if (Dot > 0.9999) Dot = 1.0;
422 else if (Dot < -0.9999) Dot = -1.0;
423 *angle = acos(Dot);
424 float Pi = acos(-1.0);
425 if (fabs(*angle) < kZeroTolerance) *angle = 0.0;
426 else if (fabs((*angle)-Pi) < kZeroTolerance) *angle = Pi;
427 float Sense = DotProduct3D(&Normal2, &BondVector);
428 if (Sense < 0.0) *angle = -*angle;
429 *angle *= 180.0/Pi;
430 result = true;
431 }
432 }
433 }
434 return result;
435 }
GetBondLength(long ibond)436 float Frame::GetBondLength(long ibond) {
437 CPoint3D offset;
438 long atom1 = Bonds[ibond].Atom1;
439 long atom2 = Bonds[ibond].Atom2;
440
441 offset.x = Atoms[atom1].Position.x - Atoms[atom2].Position.x;
442 offset.y = Atoms[atom1].Position.y - Atoms[atom2].Position.y;
443 offset.z = Atoms[atom1].Position.z - Atoms[atom2].Position.z;
444
445 return offset.Magnitude();
446 }
ChangeBond(long theBond,short whichPart,long theAtom)447 void Frame::ChangeBond(long theBond, short whichPart, long theAtom) {
448 if ((theBond >= 0)&&(theBond < NumBonds)&&(theAtom>=0)&&(theAtom<NumAtoms)) {
449 if (whichPart == 1) Bonds[theBond].Atom1 = theAtom;
450 else if (whichPart == 2) Bonds[theBond].Atom2 = theAtom;
451 }
452 }
GetBondAtom(long theBond,short thePart)453 long Frame::GetBondAtom(long theBond, short thePart) {
454 if ((theBond >= 0)&&(theBond < NumBonds)) {
455 if (thePart == 1) return Bonds[theBond].Atom1;
456 else if (thePart == 2) return Bonds[theBond].Atom2;
457 }
458 return -1;
459 }
DeleteBond(long BondNum)460 void Frame::DeleteBond(long BondNum) {
461 if (BondNum < NumBonds) {
462 NumBonds--;
463 for (long i=BondNum; i<NumBonds; i++) {
464 Bonds[i] = Bonds[i+1];
465 }
466 }
467 }
IncreaseBondAllocation(long NumAdditional)468 bool Frame::IncreaseBondAllocation(long NumAdditional) {
469 if (BondAllocation+NumAdditional < NumBonds) return false;
470 Bond * temp = new Bond[BondAllocation+NumAdditional];
471 if (temp) {
472 if (Bonds != NULL) {
473 memcpy(temp, Bonds, NumBonds*sizeof(Bond));
474 delete [] Bonds;
475 }
476 Bonds = temp;
477 BondAllocation += NumAdditional;
478 } else return false;
479 return true;
480 }
DeleteAtom(long AtomNum)481 void Frame::DeleteAtom(long AtomNum) { //remove the atom and pull down any higher atoms
482 if ((AtomNum>=0)&&(AtomNum<NumAtoms)) {
483 if (GetAtomSelection(AtomNum)) {
484 natoms_selected--;
485 }
486 if ((AtomNum<(NumAtoms-1))&&(NumAtoms>1))
487 memcpy(&(Atoms[AtomNum]), &(Atoms[AtomNum+1]), (NumAtoms-AtomNum - 1)*sizeof(mpAtom));
488 NumAtoms--;
489
490 if (targeted_atom == AtomNum) {
491 targeted_atom = -1;
492 } else if (targeted_atom > AtomNum) {
493 targeted_atom--;
494 }
495
496 //remove this atom from the bond list
497 /* for (long ii=0; ii<NumBonds; ii++) { */
498 for (long ii = NumBonds - 1; ii >= 0; ii--) {
499 if (Bonds[ii].Atom1 == AtomNum || Bonds[ii].Atom2 == AtomNum) {
500 DeleteBond(ii);
501 } else {
502 if (Bonds[ii].Atom1 > AtomNum) {
503 Bonds[ii].Atom1--;
504 }
505 if (Bonds[ii].Atom2 > AtomNum) {
506 Bonds[ii].Atom2--;
507 }
508 }
509 }
510
511 // Delete any orbitals and normal modes
512 if (Vibs) {
513 delete Vibs;
514 Vibs = NULL;
515 }
516 DeleteOrbitals();
517 while (SurfaceList) DeleteSurface(0);
518 }
519 }
SetAtomType(long theAtom,short atmType)520 bool Frame::SetAtomType(long theAtom, short atmType) {
521 bool result = false;
522 if ((theAtom>=0)&&(theAtom<NumAtoms)) {
523 result = Atoms[theAtom].SetType(atmType);
524 Atoms[theAtom].SetDefaultCoordinationNumber();
525 Atoms[theAtom].SetDefaultLonePairCount();
526 }
527 return result;
528 }
GetAtomType(long iatom) const529 short Frame::GetAtomType(long iatom) const {
530 short result=0;
531 if (Atoms&&(iatom>=0)&&(iatom<NumAtoms)) result = Atoms[iatom].GetType();
532 return result;
533 }
GetAtomPosition(long iatom,CPoint3D & p) const534 bool Frame::GetAtomPosition(long iatom, CPoint3D & p) const {
535 bool result=false;
536 if (Atoms&&(iatom>=0)&&(iatom<NumAtoms)) {
537 result = true;
538 p = Atoms[iatom].Position;
539 }
540 return result;
541 }
SetAtomPosition(long theAtom,const CPoint3D & pos)542 bool Frame::SetAtomPosition(long theAtom, const CPoint3D & pos) {
543 bool result = false;
544
545 if ((theAtom>=0)&&(theAtom<NumAtoms)) {
546 Atoms[theAtom].Position = pos;
547 result = true;
548 }
549 return result;
550 }
551
SetAtomSelection(long atom_id,bool select_it)552 bool Frame::SetAtomSelection(long atom_id, bool select_it) {
553
554 // In order to properly maintain the number of atoms that are selected, all
555 // atom selection should be done by calling this function. Calculating the
556 // number of atoms selected on the fly is computationally wasteful and is
557 // needed frequently, so we want to keep track of the number here.
558
559 // If the atom_id is out of bounds, return false.
560 if (atom_id < 0 || atom_id >= NumAtoms) {
561 return false;
562 }
563
564 // If this atom is wanting to be selected, but it's also a
565 // symmetry dependent atom and we're in symmetry edit mode,
566 // then we do not let it get selected.
567 if (select_it && MolWin->InSymmetryEditMode() &&
568 !Atoms[atom_id].IsSymmetryUnique()) {
569 return false;
570 }
571
572 // Otherwise, select or deselect the atom and return true.
573 if (select_it != Atoms[atom_id].GetSelectState()) {
574 if (select_it) {
575 natoms_selected++;
576 } else {
577 natoms_selected--;
578 }
579 }
580 Atoms[atom_id].SetSelectState(select_it);
581 return true;
582
583 }
584
GetAtomSelection(long atom_id) const585 bool Frame::GetAtomSelection(long atom_id) const {
586 return atom_id >= 0 && atom_id < NumAtoms &&
587 Atoms[atom_id].GetSelectState();
588 }
589
GetNumAtomsSelected(void) const590 int Frame::GetNumAtomsSelected(void) const {
591 return natoms_selected;
592 }
593
GetNumElectrons(void) const594 long Frame::GetNumElectrons(void) const {
595 long result=0;
596 for (long i=0; i<NumAtoms; i++) result += Atoms[i].Type;
597 return result;
598 };
599
600 /**
601 * Reevaluates the atom's bonds in this frame. If KeepOldBonds is true and
602 * selectedOnly is true, then all bonds between selected atoms are retained and
603 * so are all bonds between unselected atoms; the only bonds that change are
604 * those between an unselected atom and a selected one. If KeepOldBonds is
605 * true and selectedOnly is false, all incoming bonds are retained and new ones
606 * can be added anywhere. If KeepOldBonds is false, all incoming bonds are
607 * deleted. In this case, selectedOnly being true will add new bonds only if
608 * they involve at least one selected atom, or if true, they may be added
609 * anywhere.
610 * @param Prefs Preferences used to determining bonding sensitivity.
611 * @param KeepOldBonds Flag indicating which bonds to leave alone.
612 * @param ProgressInd A progress window to indicate long operational status
613 * @param selectedOnly Flag indicating which atoms' bonds to consider.
614 */
SetBonds(WinPrefs * Prefs,bool KeepOldBonds,Progress * ProgressInd,bool selectedOnly)615 void Frame::SetBonds(WinPrefs *Prefs, bool KeepOldBonds, Progress * ProgressInd, bool selectedOnly) {
616
617 long iatm, jatm, maxbonds;
618 CPoint3D offset;
619 float distance, AutoDist=-1.0;
620 bool newBond=true;
621
622 if (ProgressInd) ProgressInd->ChangeText("Determining bonding...");
623 maxbonds = NumAtoms*8;
624 Bond * OldBonds = Bonds;
625 long NumOldBonds = NumBonds;
626 long OldBondAllocation = BondAllocation;
627 Bonds = new Bond[maxbonds];
628 if (Bonds == NULL) {
629 MessageAlert("Insufficient Memory for Set Bond\nLength operation. Old Bonds left untouched.");
630 Bonds = OldBonds;
631 return;
632 }
633 NumBonds = 0;
634 BondAllocation = maxbonds;
635 if (KeepOldBonds) {
636 //Here we copy the old bond array over preserving the existing bonds
637 NumBonds = 0;
638 for (long ibond=0; ibond<NumOldBonds; ibond++) {
639 // In selectedOnly mode we only copy bonds that connect atoms
640 // in the same selection state. In this way, we reevaluate bonds
641 // between the interface between selected and unselected.
642 if (selectedOnly) {
643 if (Atoms[OldBonds[ibond].Atom1].GetSelectState() ==
644 Atoms[OldBonds[ibond].Atom2].GetSelectState()) {
645 Bonds[NumBonds] = OldBonds[ibond];
646 NumBonds++;
647 }
648 } else {
649 Bonds[NumBonds] = OldBonds[ibond];
650 NumBonds++;
651 }
652 }
653 }
654 float AutoScale = Prefs->GetAutoBondScale();
655 bool AutoBond = Prefs->GetAutoBond();
656 bool HHBondFlag = Prefs->GetHHBondFlag();
657 bool GuessBondOrder = Prefs->DetermineBondOrder();
658 float MaxBondLength = Prefs->GetMaxBondLength();
659 MaxBondLength *= MaxBondLength;
660 BondOrder lOrder;
661 long iType, jType;
662 float workTotal = 100.0f/NumAtoms;
663 if (AutoBond && Prefs->AllowHydrogenBonds()) workTotal /= 3.0f;
664 for (iatm=0; iatm<NumAtoms; iatm++) {
665 if (ProgressInd) {
666 if (!ProgressInd->UpdateProgress((float) iatm*workTotal)) {
667 delete [] Bonds;
668 Bonds = OldBonds;
669 NumBonds = NumOldBonds;
670 BondAllocation = OldBondAllocation;
671 return;
672 }
673 }
674 iType = Atoms[iatm].GetType();
675 if (iType > 115) continue;
676 long iSize = Prefs->GetAtomSize(iType-1);
677 CPoint3D iPos = Atoms[iatm].Position;
678 bool iSelectState = !Atoms[iatm].GetSelectState();
679 for (jatm=iatm+1; jatm<NumAtoms; jatm++) {
680 if (selectedOnly && iSelectState && !Atoms[jatm].GetSelectState())
681 continue;
682 jType = Atoms[jatm].GetType();
683 if (HHBondFlag) /*if both atoms are H's don't allow bonds if HHBondFlag is set*/
684 if ((iType == 1)&&(jType == 1)) continue;
685 if (jType > 115) continue;
686 offset = iPos - Atoms[jatm].Position;
687 distance = offset.x * offset.x + offset.y * offset.y + offset.z * offset.z;
688 //We'll compare the squares of the distances to avoid the sqrt
689 // distance = offset.Magnitude();
690 newBond = false;
691 if (distance <= MaxBondLength) {
692 newBond = true;
693 lOrder = kSingleBond;
694 } else if (AutoBond && !newBond) {
695 AutoDist = AutoScale*((float) (iSize + Prefs->GetAtomSize(jType-1)));
696 AutoDist *= AutoDist;
697 if (distance <= AutoDist) {
698 newBond = true;
699 lOrder = kSingleBond;
700 }
701 if (newBond && GuessBondOrder) {//See if this might qualify as a multiple bond
702 if ((iType != 1) && (jType != 1)) {
703 if (distance <= (.725 * .725 * AutoDist))
704 lOrder = kTripleBond;
705 else if (distance <= (0.80 * 0.80 * AutoDist))
706 lOrder = kDoubleBond;
707 }
708 }
709 }
710 if (newBond) {
711 if (KeepOldBonds) {
712 if (BondExists(iatm, jatm) >= 0) newBond = false;
713 }
714 if (newBond) {
715 if (Atoms[iatm].IsSIMOMMAtom() != Atoms[jatm].IsSIMOMMAtom())
716 continue; //bonds are not allowed between SIMOMM and ab intio atoms
717 if (! AddBond(iatm, jatm, lOrder)) {
718 MessageAlert("Insufficient Memory for Set Bond\nLength operation. Old Bonds left untouched.");
719 delete [] Bonds;
720 Bonds = OldBonds;
721 NumBonds = NumOldBonds;
722 BondAllocation = OldBondAllocation;
723 return;
724 }
725 }
726 }
727 }
728 }
729 if (AutoBond && Prefs->AllowHydrogenBonds()) {
730 //Searching the bond list for existing bonds is killing this algorithm when there are large numbers
731 //of atoms. Try creating a vector of the existing bonds to iatm and jatm early then search that short
732 //list later.
733 std::vector<long> bondsToI, bondsToJ;
734 for (iatm=0; iatm<NumAtoms; iatm++) {
735 if (ProgressInd) {
736 if (!ProgressInd->UpdateProgress((float) (2*iatm+NumAtoms)*workTotal)) {
737 delete [] Bonds;
738 Bonds = OldBonds;
739 NumBonds = NumOldBonds;
740 BondAllocation = OldBondAllocation;
741 return;
742 }
743 }
744 iType = Atoms[iatm].GetType();
745 //only consider H bonds with N, O, F, P, S, Cl, Se, and Br
746 if (iType != 1) continue; //Loop over hydrogens
747 if (selectedOnly && !Atoms[iatm].GetSelectState()) continue;
748 float iSize = AutoScale * ((float) (Prefs->GetAtomSize(iType-1)));
749 CPoint3D iPos = Atoms[iatm].Position;
750 bondsToI.clear();
751 for (long i=0; i<NumBonds; i++) {
752 if ((Bonds[i].Atom1 == iatm)||(Bonds[i].Atom2 == iatm)) bondsToI.push_back(i);
753 }
754 for (jatm=0; jatm<NumAtoms; jatm++) {
755 if (selectedOnly && !Atoms[jatm].GetSelectState()) continue;
756 jType = Atoms[jatm].GetType();
757 if (!(((jType>=7)&&(jType<=9))||((jType>=15)&&(jType<=17))||
758 (jType==34)||(jType==35))) continue;
759
760 //At this point iatm is a Hydrogen and jatm is a N, O, F, P, S, Cl, Se, or Br
761 std::vector<long>::const_iterator it_i;
762 bool ijbond=false;
763 for (it_i = bondsToI.begin(); it_i < bondsToI.end(); it_i++) {
764 if (*it_i == jatm) {
765 ijbond = true;
766 break;
767 }
768 }
769 if (!ijbond) { //can't be an existing bond
770 newBond = false;
771 offset = iPos - Atoms[jatm].Position;
772 distance = offset.x*offset.x + offset.y*offset.y + offset.z*offset.z;
773 // distance = offset.Magnitude();
774 // AutoDist = AutoScale*((float) (Prefs->GetAtomSize(iType-1) +
775 // Prefs->GetAtomSize(jType-1)));
776 AutoDist = iSize + AutoScale*((float) (Prefs->GetAtomSize(jType-1)));
777 float testDistance = 1.6*1.6*AutoDist*AutoDist;
778 if (distance <= testDistance) newBond = true;
779
780 if (newBond) {
781 bondsToJ.clear();
782 for (long i=0; i<NumBonds; i++) {
783 if ((Bonds[i].Atom1 == jatm)||(Bonds[i].Atom2 == jatm)) bondsToJ.push_back(i);
784 }
785 //scan the bond list, prevent h-bonds between any pair of
786 //atoms which are bonded to the same atom
787 for (long i=0; i<NumAtoms; i++) {
788 long ib = -1;
789 std::vector<long>::const_iterator it_iBond;
790 for (it_iBond = bondsToI.begin(); it_iBond < bondsToI.end(); it_iBond++) {
791 if (*it_iBond == i) {
792 ib = *it_iBond;
793 break;
794 }
795 }
796 //Don't filter on preexisting hydrogen bonds
797 if (ib >= 0)
798 if (Bonds[ib].Order == kHydrogenBond) {
799 //already a H-bond to this H. This is ok if the H
800 //is between the heavy atoms. However, if one distance
801 //is much bigger than the other the longer one should
802 //be dropped.
803 float angle;
804 GetBondAngle(i, iatm, jatm, &angle);
805 if (angle < 70.0) {
806 GetBondLength(i, iatm, &testDistance);
807 float angle2;
808 GetBondAngle(i, jatm, iatm, &angle2);
809 if (angle2 > 95.0) {
810 if (distance > (testDistance*testDistance)) {
811 newBond = false;
812 break;
813 } else { //delete the other bond
814 DeleteBond(ib);
815 }
816 }
817 }
818 ib = -1;
819 }
820 long jb = -1;
821 for (it_i = bondsToJ.begin(); it_i < bondsToJ.end(); it_i++) {
822 if (*it_i == i) {
823 jb = *it_i;
824 break;
825 }
826 }
827 if (jb >= 0)
828 if (Bonds[jb].Order == kHydrogenBond) jb = -1;
829 if ((ib>=0) && (jb>=0)) {
830 //Both atoms have an existing single (or higher) bond
831 //and it thus doesn't make sense for them to have an H bond
832 newBond = false;
833 break;
834 }
835 if (!newBond) break;
836 }
837 if (newBond) {
838 if (! AddBond(iatm, jatm, kHydrogenBond)) {
839 MessageAlert("Insufficient Memory for Set Bond\nLength operation. Old Bonds left untouched.");
840 delete [] Bonds;
841 Bonds = OldBonds;
842 NumBonds = NumOldBonds;
843 BondAllocation = OldBondAllocation;
844 return;
845 }
846 bondsToI.push_back(NumBonds);
847 }
848 }
849 }
850 }
851 }
852 }
853 if (OldBonds) delete [] OldBonds;
854 if (BondAllocation > NumBonds+50) {
855 Bond * temp = new Bond[NumBonds+5];
856 maxbonds = NumBonds+5;
857 if (temp) {
858 memcpy(temp, Bonds, NumBonds*sizeof(Bond));
859 delete [] Bonds;
860 Bonds = temp;
861 BondAllocation = maxbonds;
862 }
863 }
864 } /* SetBonds */
865
GetMP2Energy(void) const866 double Frame::GetMP2Energy(void) const {
867 return GetEnergy(PT2Energy);
868 }
GetKineticEnergy(void) const869 double Frame::GetKineticEnergy(void) const {
870 return GetEnergy(KineticEnergy);
871 }
GetEnergy(TypeOfEnergy t) const872 double Frame::GetEnergy(TypeOfEnergy t) const {
873 double result = 0.0;
874 std::vector<EnergyValue>::const_iterator it = Energies.begin();
875 while (it < Energies.end()) {
876 if ((*it).type == t) {
877 result = (*it).value;
878 break;
879 }
880 ++it;
881 }
882
883 return result;
884 }
SetEnergy(const double & val,TypeOfEnergy t)885 void Frame::SetEnergy(const double & val, TypeOfEnergy t) {
886 bool found = false;
887
888 std::vector<EnergyValue>::iterator it = Energies.begin();
889 while (it < Energies.end()) {
890 if ((*it).type == t) {
891 (*it).value = val;
892 found = true;
893 break;
894 }
895 ++it;
896 }
897 if (!found) Energies.push_back(EnergyValue(val, t));
898 }
899
BondExists(long a1,long a2) const900 long Frame::BondExists(long a1, long a2) const {
901 long result = -1;
902 long higherAtom;
903 long lowerAtom;
904 if (a1 > a2) {
905 higherAtom = a1;
906 lowerAtom = a2;
907 } else {
908 higherAtom = a2;
909 lowerAtom = a1;
910 }
911
912 //We assume that the atom pair is always listed with the higher atom in the Atom1 spot
913 for (long i=0; i<NumBonds; i++) {
914 if (Bonds[i].Atom1 == higherAtom) {
915 if (Bonds[i].Atom2 == lowerAtom) {
916 // if ((Bonds[i].Atom1 == a1 && Bonds[i].Atom2 == a2) ||
917 // (Bonds[i].Atom1 == a2 && Bonds[i].Atom2 == a1)) {
918 result = i;
919 break;
920 }
921 }
922 }
923 return result;
924 }
925
926 void ReadGAMESSlogVectors(BufferFile * Buffer, float *Vectors, long MaxFuncs, long NumFuncs);
927
ParseGAMESSGuessVectors(BufferFile * Buffer,long NumFuncs,TypeOfWavefunction t,Progress * lProgress)928 void Frame::ParseGAMESSGuessVectors(BufferFile * Buffer, long NumFuncs, TypeOfWavefunction t,
929 Progress * lProgress) {
930 //first we need to parse off the separate list of occupancies, then try
931 //using the normal eigenvector routine to parse the rest.
932 //What is the UHF behavior? There are a couple of problems. First UHF needs to be handled
933 //also ROHF printout seems broke since it prints alpha and beta orbitals separately
934 //finally the MCSCF guess appears to have all zero occupation #'s.
935 float * Occupancies=NULL;
936 char Line[kMaxLineLength+1];
937 int NumOrbs = NumFuncs;
938 //occupancies are printed out when using HCORE or Huekel guess types
939 if (Buffer->LocateKeyWord("ASSIGNED OCCUPANCIES", 20)) {
940 Buffer->SkipnLines(2);
941 Occupancies = new float[NumFuncs];
942 for (int i=0; i<NumFuncs; i++) Occupancies[i] = 0.0;
943 int iorb = 0;
944 while (iorb < NumOrbs) {
945 int imaxorb = MIN(10, NumOrbs-iorb); //Max of 10 orbitals per line
946 Buffer->GetLine(Line);
947 int LinePos = 0;
948 for (int jorb=0; jorb<imaxorb; jorb++) {
949 int nChar;
950 int ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(Occupancies[iorb]), &nChar);
951 if (ScanErr==1) LinePos += nChar;
952 else {
953 imaxorb = jorb;
954 if (jorb==0) { //No more orbitals found
955 imaxorb = 0;
956 NumOrbs = iorb;
957 }
958 break;
959 }
960 iorb++;
961 }
962 if (imaxorb <= 0) break;
963 }
964 NumOrbs = iorb;
965 }
966 Buffer->BackupnLines(1);
967 Buffer->SetFilePos(Buffer->FindBlankLine());
968 Buffer->SkipnLines(1);
969 if (NumOrbs > 0) {
970 OrbitalRec * orbs = ParseGAMESSEigenVectors(Buffer, NumFuncs, NumOrbs,
971 /*long NumBetaOrbs*/ 0, NumOrbs, 0,
972 t, lProgress);
973 if (orbs) {
974 if (Occupancies) orbs->SetOccupancy(Occupancies, NumFuncs);
975 orbs->setOrbitalType(GuessOrbital);
976 } else if (Occupancies) {
977 delete [] Occupancies;
978 Occupancies = NULL;
979 }
980 }
981 }
982 //Handle MCSCF vectors (Natural and Optimized)
ParseGAMESSMCSCFVectors(BufferFile * Buffer,long NumFuncs,long NumOrbs,Progress * lProgress)983 void Frame::ParseGAMESSMCSCFVectors(BufferFile * Buffer, long NumFuncs,
984 long NumOrbs, Progress * lProgress) {
985 long iorb, imaxorb=0, NumNOrbs=0, NumOptOrbs=0, TestOrb,
986 ScanErr, LinePos, jorb;
987 int nChar;
988 float *Vectors, *lEnergy, *OccNums;
989 char *SymType, Line[kMaxLineLength+1];
990
991 long StartPos = Buffer->GetFilePos(), NOrbPos, OptOrbPos;
992 //Look for the MCSCF Natural Orbitals, really should only read them in
993 //for FORS type MCSCF functions
994 if (Buffer->LocateKeyWord("MCSCF NATURAL ORBITALS", 22)) {
995 NOrbPos = Buffer->GetFilePos();
996 imaxorb =1;
997 NumNOrbs = NumOrbs;
998 } else if (Buffer->LocateKeyWord("-MCHF- NATURAL ORBITALS", 23)) {
999 NOrbPos = Buffer->GetFilePos();
1000 imaxorb =1;
1001 NumNOrbs = NumOrbs;
1002 }
1003 Buffer->SetFilePos(StartPos);
1004 if (Buffer->LocateKeyWord("MCSCF OPTIMIZED ORBITALS", 24)) {
1005 OptOrbPos = Buffer->GetFilePos();
1006 imaxorb++;
1007 NumOptOrbs = NumFuncs; //GAMESS now seems to print all optimized orbs
1008 } else if (Buffer->LocateKeyWord("-MCHF- OPTIMIZED ORBITALS", 25)) {
1009 OptOrbPos = Buffer->GetFilePos();
1010 imaxorb++;
1011 NumOptOrbs = NumFuncs; //GAMESS now seems to print all optimized orbs
1012 }
1013 Buffer->SetFilePos(StartPos);
1014 if (imaxorb==0) return; //No MC orbital vectors found
1015
1016 long maxfuncs = NumFuncs; //Always true
1017 OrbitalRec * OrbSet = NULL;
1018
1019 try {
1020 if (NumNOrbs) { //read in the Natural Orbitals
1021 OrbSet = new OrbitalRec(NumNOrbs, 0, maxfuncs);
1022 OrbSet->setOrbitalWavefunctionType(MCSCF);
1023 OrbSet->setOrbitalType(NaturalOrbital);
1024 Vectors = OrbSet->Vectors;
1025 lEnergy = OrbSet->Energy;
1026 SymType = OrbSet->SymType;
1027 bool SymmetriesFound = true;
1028 OccNums = OrbSet->OrbOccupation = new float [NumNOrbs];
1029 if (!OccNums) throw MemoryError();
1030 //Clear the occupation numbers just in case a full set is not read in
1031 for (iorb=0; iorb<NumNOrbs; iorb++) OccNums[iorb] = 0.0;
1032
1033 Buffer->SetFilePos(NOrbPos);
1034 Buffer->SkipnLines(3);
1035 iorb=0;
1036 while (iorb<NumNOrbs) {
1037 // Allow a little backgrounding and user cancels
1038 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1039 delete OrbSet;
1040 return;
1041 }
1042 imaxorb = ((10) > (NumNOrbs-iorb)) ? (NumNOrbs-iorb) : (10); //Max of 10 orbitals per line
1043 Buffer->GetLine(Line);
1044 LinePos = 0;
1045 for (jorb=0; jorb<imaxorb; jorb++) {
1046 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1047 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1048 else {
1049 imaxorb = jorb;
1050 if (jorb==0) { //No more orbitals found
1051 imaxorb = 0;
1052 NumNOrbs = iorb;
1053 }
1054 break;
1055 }
1056 }
1057 if (imaxorb <= 0) break;
1058 Buffer->GetLine(Line);
1059 //Older versions had a blank line, skip it if present
1060 if (IsBlank(Line)) Buffer->GetLine(Line);
1061 //first read in the orbital energy/occupation number of each orbital
1062 LinePos = 0;
1063 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital energies
1064 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(lEnergy[iorb+jorb]),&nChar);
1065 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1066 if (lEnergy[iorb+jorb] >= 0.0) { //Only core orbital energies are printed
1067 if (OccNums) OccNums[iorb+jorb] = lEnergy[iorb+jorb];
1068 lEnergy[iorb+jorb] = 0.0;
1069 } else if (OccNums) OccNums[iorb+jorb] = 2.0;
1070 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1071 }
1072 //In old versions orbital symetries are not printed for Natural orbitals, but new versions have them
1073 Buffer->GetLine(Line);
1074 if (!IsBlank(Line)) {
1075 LinePos = 0;
1076 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1077 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1078 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1079 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1080 }
1081 } else SymmetriesFound = false;
1082 //read in the vector block
1083 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*maxfuncs]), maxfuncs, imaxorb);
1084 iorb += imaxorb;
1085 Buffer->SkipnLines(1); //Skip blank line between blocks
1086 }
1087 if (!SymmetriesFound && (OrbSet->SymType!=NULL)) {
1088 delete [] OrbSet->SymType;
1089 OrbSet->SymType = NULL;
1090 }
1091 OrbSet->ReSize(NumNOrbs, 0);
1092 }
1093 }
1094 catch (...) {
1095 if (OrbSet != NULL) {
1096 delete OrbSet;
1097 OrbSet = NULL;
1098 }
1099 }
1100 if (OrbSet) {
1101 OrbSet->SetOrbitalOccupancy(NumNOrbs, 0); //We should know the occupation # for all natural orbs.
1102 Orbs.push_back(OrbSet);
1103 OrbSet = NULL;
1104 }
1105 try {
1106 if (NumOptOrbs) {
1107 OrbSet = new OrbitalRec(NumOptOrbs, 0, maxfuncs);
1108 Vectors = OrbSet->Vectors;
1109 lEnergy = OrbSet->Energy;
1110 SymType = OrbSet->SymType;
1111 OrbSet->setOrbitalWavefunctionType(MCSCF);
1112 OrbSet->setOrbitalType(OptimizedOrbital);
1113
1114 Buffer->SetFilePos(OptOrbPos); //move to the pre-determined start of the optimized orbitals
1115 Buffer->SkipnLines(3);
1116 iorb=0;
1117 while (iorb<NumOptOrbs) {
1118 // Allow a little backgrounding and user cancels
1119 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1120 delete OrbSet;
1121 return;
1122 }
1123 // imaxorb = MIN(10, NumOptOrbs-iorb); //Max of 10 orbitals per line
1124 imaxorb = ((10) > (NumOptOrbs-iorb)) ? (NumOptOrbs-iorb) : (10); //Max of 10 orbitals per line
1125 Buffer->GetLine(Line);
1126 LinePos = 0;
1127 for (jorb=0; jorb<imaxorb; jorb++) {
1128 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1129 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1130 else {
1131 imaxorb = jorb;
1132 if (jorb==0) { //No more orbitals found
1133 imaxorb = 0;
1134 NumOptOrbs = iorb;
1135 }
1136 break;
1137 }
1138 }
1139 if (imaxorb <= 0) break;
1140 //first read in the orbital energy/occupation number of each orbital
1141 Buffer->GetLine(Line);
1142 LinePos = 0;
1143 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital energies
1144 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(lEnergy[iorb+jorb]),&nChar);
1145 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1146 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1147 }
1148 Buffer->GetLine(Line);
1149 LinePos = 0;
1150 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1151 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1152 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1153 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1154 }
1155 //read in the vector block
1156 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*maxfuncs]), maxfuncs, imaxorb);
1157 iorb += imaxorb;
1158 Buffer->SkipnLines(1); //Skip blank line between blocks
1159 }
1160 OrbSet->ReSize(NumOptOrbs, 0);
1161 }
1162 }
1163 catch (...) {
1164 if (OrbSet != NULL) {
1165 delete OrbSet;
1166 OrbSet = NULL;
1167 }
1168 }
1169 if (OrbSet) {
1170 Orbs.push_back(OrbSet);
1171 }
1172 }
1173 //There is a relatively rare runtype based on MCSCF that produces diabatic orbitals
1174 //Mike describes these as neither the final MOs or Naturals orbitals (though there are occupancies)
ParseGAMESSMCSCFDiabaticVectors(BufferFile * Buffer,long NumFuncs,long NumOrbs,Progress * lProgress)1175 void Frame::ParseGAMESSMCSCFDiabaticVectors(BufferFile * Buffer, long NumFuncs,
1176 long NumOrbs, Progress * lProgress) {
1177 long iorb, imaxorb=0, NumNOrbs=0, TestOrb,
1178 ScanErr, LinePos, jorb;
1179 int nChar;
1180 float *Vectors, *lEnergy, *OccNums;
1181 char *SymType, Line[kMaxLineLength+1];
1182
1183 long maxfuncs = NumFuncs; //Always true
1184 OrbitalRec * OrbSet = NULL;
1185 NumNOrbs = NumOrbs;
1186
1187 try {
1188 OrbSet = new OrbitalRec(NumNOrbs, 0, maxfuncs);
1189 OrbSet->setOrbitalWavefunctionType(MCSCF);
1190 OrbSet->setOrbitalType(DiabaticMolecularOrbital);
1191 Vectors = OrbSet->Vectors;
1192 lEnergy = OrbSet->Energy;
1193 SymType = OrbSet->SymType;
1194 bool SymmetriesFound = true;
1195 OccNums = OrbSet->OrbOccupation = new float [NumNOrbs];
1196 if (!OccNums) throw MemoryError();
1197 //Clear the occupation numbers just in case a full set is not read in
1198 for (iorb=0; iorb<NumNOrbs; iorb++) OccNums[iorb] = 0.0;
1199
1200 Buffer->SetFilePos(Buffer->FindBlankLine());
1201 Buffer->SkipnLines(1);
1202 iorb=0;
1203 while (iorb<NumNOrbs) {
1204 // Allow a little backgrounding and user cancels
1205 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1206 delete OrbSet;
1207 return;
1208 }
1209 imaxorb = ((10) > (NumNOrbs-iorb)) ? (NumNOrbs-iorb) : (10); //Max of 10 orbitals per line
1210 Buffer->GetLine(Line);
1211 LinePos = 0;
1212 for (jorb=0; jorb<imaxorb; jorb++) {
1213 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1214 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1215 else {
1216 imaxorb = jorb;
1217 if (jorb==0) { //No more orbitals found
1218 imaxorb = 0;
1219 NumNOrbs = iorb;
1220 }
1221 break;
1222 }
1223 }
1224 if (imaxorb <= 0) break;
1225 Buffer->GetLine(Line);
1226 //Older versions had a blank line, skip it if present
1227 if (IsBlank(Line)) Buffer->GetLine(Line);
1228 //first read in the orbital energy/occupation number of each orbital
1229 LinePos = 0;
1230 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital energies
1231 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(lEnergy[iorb+jorb]),&nChar);
1232 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1233 if (lEnergy[iorb+jorb] >= 0.0) { //Only core orbital energies are printed
1234 if (OccNums) OccNums[iorb+jorb] = lEnergy[iorb+jorb];
1235 lEnergy[iorb+jorb] = 0.0;
1236 } else if (OccNums) OccNums[iorb+jorb] = 2.0;
1237 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1238 }
1239 //In old versions orbital symetries are not printed for Natural orbitals, but new versions have them
1240 Buffer->GetLine(Line);
1241 if (!IsBlank(Line)) {
1242 LinePos = 0;
1243 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1244 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1245 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1246 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1247 }
1248 } else SymmetriesFound = false;
1249 //read in the vector block
1250 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*maxfuncs]), maxfuncs, imaxorb);
1251 iorb += imaxorb;
1252 Buffer->SkipnLines(1); //Skip blank line between blocks
1253 }
1254 if (!SymmetriesFound && (OrbSet->SymType!=NULL)) {
1255 delete [] OrbSet->SymType;
1256 OrbSet->SymType = NULL;
1257 }
1258 OrbSet->ReSize(NumNOrbs, 0);
1259 }
1260 catch (...) {
1261 if (OrbSet != NULL) {
1262 delete OrbSet;
1263 OrbSet = NULL;
1264 }
1265 }
1266 if (OrbSet) {
1267 OrbSet->SetOrbitalOccupancy(NumNOrbs, 0); //We should know the occupation # for all natural orbs.
1268 Orbs.push_back(OrbSet);
1269 }
1270 }
1271
1272 //Reads in a general set of eigenvectors from a GAMESS log file. The number of beta
1273 //orbitals is used only to indicate the need to read in both alpha and beta sets.
1274 //It is possible that there will be fewer orbitals than NumOrbs. If so then this
1275 //routine will (hopefully) recognize that and truncate the read.
ParseGAMESSEigenVectors(BufferFile * Buffer,long NumFuncs,long NumOrbs,long NumBetaOrbs,const long & NumOccAlpha,const long & NumOccBeta,const TypeOfWavefunction & method,Progress * lProgress)1276 OrbitalRec * Frame::ParseGAMESSEigenVectors(BufferFile * Buffer, long NumFuncs, long NumOrbs,
1277 long NumBetaOrbs, const long & NumOccAlpha, const long & NumOccBeta,
1278 const TypeOfWavefunction & method, Progress * lProgress) {
1279 long iorb=0, imaxorb, maxfuncs, TestOrb, LinePos, ScanErr, jorb, OrbitalOffset=0;
1280 int nChar;
1281 float *Vectors, *lEnergy;
1282 char *SymType, Line[kMaxLineLength+1];
1283 bool energyError=false;
1284
1285 maxfuncs = NumFuncs; //Always true
1286 OrbitalRec * OrbSet = NULL;
1287 try {
1288 OrbSet = new OrbitalRec(NumOrbs, (NumBetaOrbs?NumOrbs:0), NumFuncs);
1289 Vectors = OrbSet->Vectors;
1290 lEnergy = OrbSet->Energy;
1291 SymType = OrbSet->SymType;
1292 OrbSet->setOrbitalWavefunctionType(RHF);
1293 OrbSet->setOrbitalType(OptimizedOrbital);
1294
1295 //UHF orbs are not always organized in the same format
1296 //To simplify just search for the 1 denoting the first vector
1297 if (NumBetaOrbs) Buffer->LocateKeyWord("1", 1, -1);
1298 while (iorb<NumOrbs) {
1299 // Allow a little backgrounding and user cancels
1300 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1301 delete OrbSet;
1302 return NULL;
1303 }
1304 imaxorb = MIN(10, NumOrbs-iorb); //Max of 10 orbitals per line
1305 Buffer->GetLine(Line);
1306 LinePos = 0;
1307 for (jorb=0; jorb<imaxorb; jorb++) {
1308 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1309 //Dmitri's npreo option skips printing of the first n orbitals. Need to handle a case
1310 //where the orbital printout doesn't start at 1.
1311 if ((iorb+jorb+1)==1) {
1312 if (ScanErr && (TestOrb > 1)) {
1313 OrbitalOffset = TestOrb - 1;
1314 OrbSet->setStartingOrbitalOffset(OrbitalOffset);
1315 }
1316 }
1317 if (ScanErr && (TestOrb==jorb+iorb+OrbitalOffset+1)) LinePos += nChar;
1318 else {
1319 imaxorb = jorb;
1320 if (jorb==0) { //No more orbitals found
1321 imaxorb = 0;
1322 NumOrbs = iorb;
1323 }
1324 break;
1325 }
1326 }
1327 if (imaxorb > 0) {
1328 //read in energy and the symmetry of each orbital (ie. A1 A2 B1�)
1329 //First comes the orbital energy. Unfortunately for very large system the numbers
1330 //can run together (no space in the output format) or exceed the field width (***'s)
1331 //Since that will probably only happen for very high energy orbitals (that don't
1332 //matter chemically) just use the last "good" value and echo a warning.
1333 //This will flag the case with *'s , but not when the values run together.
1334 Buffer->GetLine(Line);
1335 if (strlen(Line) <= 0) Buffer->GetLine(Line);
1336 LinePos = 0;
1337 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital energies
1338 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(lEnergy[iorb+jorb]),&nChar);
1339 if (ScanErr <= 0) {
1340 wxLogWarning(_("Error parsing the energy for orbital %ld. Using last energy and continuing."), iorb+jorb);
1341 if ((iorb+jorb)>0) lEnergy[iorb+jorb] = lEnergy[iorb+jorb-1];
1342 jorb++;
1343 while (jorb <imaxorb) {
1344 lEnergy[iorb+jorb] = lEnergy[iorb+jorb-1];
1345 jorb++;
1346 }
1347 energyError = true;
1348 }
1349 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1350 }
1351 Buffer->GetLine(Line);
1352 if ((strlen(Line) > 0)&&SymType) {
1353 LinePos = 0;
1354 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1355 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1356 if (ScanErr <= 0) throw DataError(); //Looks like the MO's are not complete
1357 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1358 }
1359 } else if (SymType) {
1360 delete [] OrbSet->SymType;
1361 SymType = OrbSet->SymType = NULL;
1362 }
1363 //read in the vector block
1364 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*maxfuncs]), maxfuncs, imaxorb);
1365 iorb += imaxorb;
1366 Buffer->SkipnLines(1); //Skip blank line between blocks
1367 }
1368 }
1369 OrbSet->ReSize(NumOrbs, (NumBetaOrbs?NumOrbs:0));
1370 if (NumBetaOrbs) {
1371 OrbSet->setOrbitalWavefunctionType(UHF);
1372 Vectors = OrbSet->VectorsB;
1373 lEnergy = OrbSet->EnergyB;
1374 SymType = OrbSet->SymTypeB;
1375 Buffer->BackupnLines(2);
1376 if (!Buffer->LocateKeyWord("BETA SET",8)) throw DataError();
1377 //UHF orbs are not always organized in the same format
1378 //To simplify just search for the 1 denoting the first vector
1379 Buffer->LocateKeyWord("1", 1, -1);
1380 // Buffer->SkipnLines(6);
1381 iorb=0;
1382
1383 while (iorb<NumOrbs) {
1384 // Allow a little backgrounding and user cancels
1385 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1386 delete OrbSet;
1387 return NULL;
1388 }
1389 imaxorb = MIN(10, NumOrbs-iorb); //Max of 10 orbitals per line
1390 Buffer->GetLine(Line);
1391 LinePos = 0;
1392 for (jorb=0; jorb<imaxorb; jorb++) {
1393 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1394 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1395 else {
1396 imaxorb = jorb;
1397 if (jorb==0) { //No more orbitals found
1398 imaxorb = 0;
1399 NumOrbs = iorb;
1400 }
1401 break;
1402 }
1403 }
1404 if (imaxorb > 0) {
1405 //First comes the orbital energy. Unfortunately for very large system the numbers
1406 //can run together (no space in the output format) or exceed the field width (***'s)
1407 //Since that will probably only happen for very high energy orbitals (that don't
1408 //matter chemically) just use the last "good" value and echo a warning.
1409 Buffer->GetLine(Line);
1410 if (strlen(Line) <= 0) Buffer->GetLine(Line);
1411 LinePos = 0;
1412 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital energies
1413 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(lEnergy[iorb+jorb]),&nChar);
1414 if (ScanErr <= 0) {
1415 wxLogWarning(_("Error parsing the energy for orbital %d. Using last energy and continuing."), iorb+jorb);
1416 if ((iorb+jorb)>0) lEnergy[iorb+jorb] = lEnergy[iorb+jorb-1];
1417 jorb++;
1418 while (jorb <imaxorb) {
1419 lEnergy[iorb+jorb] = lEnergy[iorb+jorb-1];
1420 jorb++;
1421 }
1422 }
1423 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1424 }
1425 Buffer->GetLine(Line); // Read in the symmetry of each orbital (ie. A1 A2 B1�)
1426 if ((strlen(Line) > 0)&&SymType) {
1427 LinePos = 0;
1428 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1429 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1430 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1431 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1432 }
1433 } else if (SymType) {
1434 delete [] OrbSet->SymTypeB;
1435 SymType = OrbSet->SymTypeB = NULL;
1436 }
1437 //read in the vector block
1438 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*maxfuncs]), maxfuncs, imaxorb);
1439 iorb += imaxorb;
1440 Buffer->SkipnLines(1); //Skip blank line between blocks
1441 }
1442 }
1443 }
1444 }
1445 catch (...) {
1446 if (OrbSet) {
1447 delete OrbSet;
1448 OrbSet = NULL;
1449 }
1450 }
1451 if (OrbSet != NULL) {
1452 if ((OrbSet->getNumAlphaOrbitals() + OrbSet->getNumBetaOrbitals())<=0) {
1453 delete OrbSet;
1454 OrbSet = NULL;
1455 } else {
1456 OrbSet->SetOrbitalOccupancy(NumOccAlpha, NumOccBeta);
1457 OrbSet->setOrbitalWavefunctionType(method);
1458 Orbs.push_back(OrbSet);
1459 }
1460 }
1461 return OrbSet;
1462 }
ReadMolDenOrbitals(BufferFile * Buffer,long NumFuncs)1463 void Frame::ReadMolDenOrbitals(BufferFile * Buffer, long NumFuncs) {
1464 //We don't have a lot of information about the orbitals so assume the largest case
1465
1466 OrbitalRec * OrbSet = NULL;
1467 try {
1468 //There doesn't seem to be a good way to predetermine the number of orbitals so dimension a
1469 //worse case
1470 OrbSet = new OrbitalRec(NumFuncs, NumFuncs, NumFuncs);
1471 char Line[kMaxLineLength+1];
1472 bool done=false;
1473 long alphaCount=0, betaCount=0;
1474
1475 while (!done && (Buffer->GetFilePos()<Buffer->GetFileSize())) {
1476 float * vector = NULL, energy=0.0, occ=-1;
1477 bool header=true, good=true, alphaSpin=true;
1478 char symLabel[5]="";
1479 while (header) {
1480 Buffer->GetLine(Line);
1481 if (FindKeyWord(Line,"[",1)>=0) {
1482 good = false;
1483 done = true;
1484 break;
1485 }
1486 int p;
1487 if ((p=FindKeyWord(Line, "ENE=", 4))>=0) {
1488 ConvertExponentStyle(Line);
1489 sscanf(&(Line[p+4]),"%f", &energy);
1490 } else if (FindKeyWord(Line, "SPIN=", 5)>=0) {
1491 if (FindKeyWord(Line, "BETA", 4)>=0) alphaSpin = false;
1492 } else if ((p=FindKeyWord(Line, "OCCUP=", 6))>=0) {
1493 ConvertExponentStyle(Line);
1494 sscanf(&(Line[p+6]),"%f", &occ);
1495 } else if ((p=FindKeyWord(Line, "SYM=", 4))>=0) {
1496 sscanf(&(Line[p+4]),"%4s", symLabel);
1497 } else {
1498 int test;
1499 float junk;
1500 ConvertExponentStyle(Line);
1501 if (sscanf(Line, "%d %f", &test, &junk) == 2 ) {//start of the vector
1502 Buffer->BackupnLines(1);
1503 header = false;
1504 break;
1505 }
1506 }
1507 }
1508 if (!good) break;
1509 if (alphaSpin) {
1510 vector = &(OrbSet->Vectors[alphaCount*NumFuncs]);
1511 OrbSet->Energy[alphaCount] = energy;
1512 if (occ > -1.0) {
1513 if (!(OrbSet->OrbOccupation)) {
1514 OrbSet->OrbOccupation = new float[NumFuncs];
1515 }
1516 OrbSet->OrbOccupation[alphaCount] = occ;
1517 }
1518 strncpy(&(OrbSet->SymType[alphaCount*5]), symLabel, 4);
1519 } else {
1520 vector = &(OrbSet->VectorsB[betaCount*NumFuncs]);
1521 OrbSet->EnergyB[betaCount] = energy;
1522 if (occ > -1.0) {
1523 if (!(OrbSet->OrbOccupationB)) {
1524 OrbSet->OrbOccupationB = new float[NumFuncs];
1525 }
1526 OrbSet->OrbOccupationB[betaCount] = occ;
1527 }
1528 strncpy(&(OrbSet->SymTypeB[betaCount*5]), symLabel, 4);
1529 }
1530
1531 for (int i=0; i<NumFuncs; i++) { //loop over the vector itself
1532 Buffer->GetLine(Line);
1533 ConvertExponentStyle(Line);
1534 int iline;
1535 float vec;
1536 if (sscanf(Line, "%d %f", &iline, &vec) == 2) {
1537 if (iline != (i+1)) throw DataError();
1538 vector[i] = vec;
1539 } else throw DataError();
1540 }
1541 if (alphaSpin) alphaCount++;
1542 else betaCount++;
1543 }
1544
1545 OrbSet->NumAlphaOrbs = alphaCount;
1546 if (OrbSet->OrbOccupation) OrbSet->NumOccupiedAlphaOrbs = alphaCount;
1547 OrbSet->NumBetaOrbs = betaCount;
1548 if (OrbSet->OrbOccupationB) OrbSet->NumOccupiedBetaOrbs = betaCount;
1549 }
1550 catch (...) {
1551 if (OrbSet) {
1552 delete OrbSet;
1553 OrbSet=NULL;
1554 }
1555 }
1556 if (OrbSet) Orbs.push_back(OrbSet);
1557 }
1558 //Handle CI vectors
ParseGAMESSCIVectors(BufferFile * Buffer,long NumFuncs,Progress * lProgress)1559 void Frame::ParseGAMESSCIVectors(BufferFile * Buffer, long NumFuncs, Progress * lProgress) {
1560 long iorb, jorb, imaxorb=0, NumNOrbs=0, TestOrb,
1561 ScanErr, LinePos;
1562 int nChar;
1563 float *Vectors, *OccNums;
1564 char Line[kMaxLineLength+1], *SymType=NULL;
1565 OrbitalRec * OrbSet = NULL;
1566
1567 if (Buffer->LocateKeyWord("NATURAL ORBITALS", 16)) {
1568 try {
1569 Buffer->SkipnLines(3); //setup to beginning of the first block
1570 bool Style = false;
1571 long begPos = Buffer->GetFilePos();
1572 //Determine the style of the orbital output
1573 Buffer->SkipnLines(1);
1574 Buffer->GetLine(Line);
1575 int t = strlen(Line);
1576 if (t > 5) Style=true;
1577 Buffer->SetFilePos(begPos);
1578
1579 OrbSet = new OrbitalRec(NumFuncs, 0, NumFuncs);
1580 OrbSet->setOrbitalWavefunctionType(CI);
1581 OrbSet->setOrbitalType(NaturalOrbital);
1582
1583 Vectors = OrbSet->Vectors;
1584 if (OrbSet->Energy) {
1585 delete [] OrbSet->Energy;
1586 OrbSet->Energy = NULL;
1587 }
1588 if (Style) {
1589 if (OrbSet->SymType == NULL) {
1590 OrbSet->SymType = new char [NumFuncs*5];
1591 if (!OrbSet->SymType) throw MemoryError();
1592 }
1593 SymType = OrbSet->SymType;
1594 } else {
1595 if (OrbSet->SymType) { //symmetry labels are not printed for NO's
1596 delete [] OrbSet->SymType;
1597 OrbSet->SymType = NULL;
1598 }
1599 }
1600 OccNums = OrbSet->OrbOccupation = new float [NumFuncs];
1601 if (!OccNums) throw MemoryError();
1602
1603 iorb=0;
1604 NumNOrbs = NumFuncs;
1605 while (iorb<NumNOrbs) {
1606 // Allow a little backgrounding and user cancels
1607 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1608 delete OrbSet;
1609 return;
1610 }
1611 imaxorb = MIN(10, NumNOrbs-iorb); //Max of 10 orbitals per line
1612 Buffer->GetLine(Line);
1613 LinePos = 0;
1614 for (jorb=0; jorb<imaxorb; jorb++) {
1615 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1616 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1617 else {
1618 imaxorb = jorb;
1619 if (jorb==0) { //No more orbitals found
1620 imaxorb = 0;
1621 NumNOrbs = iorb;
1622 }
1623 break;
1624 }
1625 }
1626 if (imaxorb <= 0) break;
1627 if (!Style) Buffer->SkipnLines(1); //Skip blank line
1628 //first read in the orbital occupation number of each orbital
1629 Buffer->GetLine(Line);
1630 LinePos = 0;
1631 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital occupations
1632 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(OccNums[iorb+jorb]),&nChar);
1633 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1634 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1635 }
1636 if (Style) { //New style has orbital symmetries
1637 Buffer->GetLine(Line);
1638 if (strlen(Line) > 0) {
1639 LinePos = 0;
1640 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1641 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1642 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1643 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1644 }
1645 }
1646 } else //old style has blank line
1647 Buffer->SkipnLines(1); //skip blank line
1648 //orbital symetries are not printed for Natural orbitals
1649 //read in the vector block
1650 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*NumFuncs]), NumFuncs, imaxorb);
1651 iorb += imaxorb;
1652 Buffer->SkipnLines(1); //Skip blank line between blocks
1653 }
1654 }
1655 catch (...) {
1656 if (OrbSet) {
1657 delete OrbSet;
1658 OrbSet=NULL;
1659 }
1660 }
1661 OrbSet->NumAlphaOrbs = NumNOrbs;
1662 OrbSet->NumOccupiedAlphaOrbs = NumNOrbs;
1663 if (OrbSet) Orbs.push_back(OrbSet);
1664 }
1665 }
1666 //Handle EOM-CC natural orbitals
ParseGAMESSEOM_CC_Vectors(BufferFile * Buffer,long NumFuncs,Progress * lProgress)1667 void Frame::ParseGAMESSEOM_CC_Vectors(BufferFile * Buffer, long NumFuncs, Progress * lProgress) {
1668 long iorb, jorb, imaxorb=0, NumNOrbs=0, TestOrb,
1669 ScanErr, LinePos;
1670 int nChar;
1671 float *Vectors, *OccNums;
1672 char Line[kMaxLineLength+1], *SymType=NULL;
1673 OrbitalRec * OrbSet = NULL;
1674
1675 if (Buffer->LocateKeyWord("NATURAL ORBITALS", 16)) {
1676 try {
1677 Buffer->SetFilePos(Buffer->FindBlankLine());
1678 Buffer->SkipnLines(1);
1679 lProgress->ChangeText("Reading EOM-CC Natural Orbitals");
1680
1681 bool Style = false;
1682 long begPos = Buffer->GetFilePos();
1683 //Determine the style of the orbital output
1684 Buffer->SkipnLines(1);
1685 Buffer->GetLine(Line);
1686 int t = strlen(Line);
1687 if (t > 5) Style=true;
1688 Buffer->SetFilePos(begPos);
1689
1690 OrbSet = new OrbitalRec(NumFuncs, 0, NumFuncs);
1691 OrbSet->setOrbitalWavefunctionType(EOM_CC);
1692 OrbSet->setOrbitalType(NaturalOrbital);
1693
1694 Vectors = OrbSet->Vectors;
1695 if (OrbSet->Energy) {
1696 delete [] OrbSet->Energy;
1697 OrbSet->Energy = NULL;
1698 }
1699 if (OrbSet->SymType == NULL) {
1700 OrbSet->SymType = new char [NumFuncs*5];
1701 if (!OrbSet->SymType) throw MemoryError();
1702 }
1703 SymType = OrbSet->SymType;
1704 OccNums = OrbSet->OrbOccupation = new float [NumFuncs];
1705 if (!OccNums) throw MemoryError();
1706
1707 iorb=0;
1708 NumNOrbs = NumFuncs;
1709 while (iorb<NumNOrbs) {
1710 // Allow a little backgrounding and user cancels
1711 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1712 delete OrbSet;
1713 return;
1714 }
1715 imaxorb = MIN(10, NumNOrbs-iorb); //Max of 10 orbitals per line
1716 Buffer->GetLine(Line);
1717 LinePos = 0;
1718 for (jorb=0; jorb<imaxorb; jorb++) {
1719 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1720 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1721 else {
1722 imaxorb = jorb;
1723 if (jorb==0) { //No more orbitals found
1724 imaxorb = 0;
1725 NumNOrbs = iorb;
1726 }
1727 break;
1728 }
1729 }
1730 if (imaxorb <= 0) break;
1731 //first read in the orbital occupation number of each orbital
1732 Buffer->GetLine(Line);
1733 LinePos = 0;
1734 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital occupations
1735 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(OccNums[iorb+jorb]),&nChar);
1736 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1737 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1738 }
1739 Buffer->GetLine(Line);
1740 if (strlen(Line) > 0) {
1741 LinePos = 0;
1742 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1743 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1744 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1745 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1746 }
1747 }
1748 //read in the vector block
1749 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*NumFuncs]), NumFuncs, imaxorb);
1750 iorb += imaxorb;
1751 Buffer->SkipnLines(1); //Skip blank line between blocks
1752 }
1753 }
1754 catch (...) {
1755 if (OrbSet) {
1756 delete OrbSet;
1757 OrbSet=NULL;
1758 }
1759 }
1760 OrbSet->NumAlphaOrbs = NumNOrbs;
1761 OrbSet->NumOccupiedAlphaOrbs = NumNOrbs;
1762 if (OrbSet) Orbs.push_back(OrbSet);
1763 }
1764 }
1765 //Parse UHF Natural Orbitals
1766 //
ParseUHFNOs(BufferFile * Buffer,long NumFuncs,Progress * lProgress)1767 void Frame::ParseUHFNOs(BufferFile * Buffer, long NumFuncs, Progress * lProgress) {
1768 long iorb, jorb, imaxorb=0, NumNOrbs=0, TestOrb,
1769 ScanErr, LinePos;
1770 int nChar;
1771 float *Vectors, *OccNums;
1772 char Line[kMaxLineLength+1], *SymType=NULL;
1773 OrbitalRec * OrbSet = NULL;
1774
1775 try {
1776 Buffer->SetFilePos(Buffer->FindBlankLine());
1777 Buffer->SkipnLines(1);
1778
1779 OrbSet = new OrbitalRec(NumFuncs, 0, NumFuncs);
1780 OrbSet->setOrbitalWavefunctionType(UHF);
1781 OrbSet->setOrbitalType(NaturalOrbital);
1782
1783 Vectors = OrbSet->Vectors;
1784 if (OrbSet->Energy) {
1785 delete [] OrbSet->Energy;
1786 OrbSet->Energy = NULL;
1787 }
1788 if (OrbSet->SymType == NULL) {
1789 OrbSet->SymType = new char [NumFuncs*5];
1790 if (!OrbSet->SymType) throw MemoryError();
1791 }
1792 SymType = OrbSet->SymType;
1793 OccNums = OrbSet->OrbOccupation = new float [NumFuncs];
1794 if (!OccNums) throw MemoryError();
1795
1796 iorb=0;
1797 NumNOrbs = NumFuncs;
1798 while (iorb<NumNOrbs) {
1799 // Allow a little backgrounding and user cancels
1800 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1801 delete OrbSet;
1802 return;
1803 }
1804 imaxorb = MIN(10, NumNOrbs-iorb); //Max of 10 orbitals per line
1805 Buffer->GetLine(Line);
1806 LinePos = 0;
1807 for (jorb=0; jorb<imaxorb; jorb++) {
1808 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1809 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1810 else {
1811 imaxorb = jorb;
1812 if (jorb==0) { //No more orbitals found
1813 imaxorb = 0;
1814 NumNOrbs = iorb;
1815 }
1816 break;
1817 }
1818 }
1819 if (imaxorb <= 0) break;
1820 //first read in the orbital occupation number of each orbital
1821 Buffer->GetLine(Line);
1822 LinePos = 0;
1823 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital occupations
1824 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(OccNums[iorb+jorb]),&nChar);
1825 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1826 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1827 }
1828 Buffer->GetLine(Line);
1829 if (strlen(Line) > 0) {
1830 LinePos = 0;
1831 for (jorb=0; jorb<imaxorb; jorb++) { //Get the orbital symmetries
1832 ScanErr = sscanf(&(Line[LinePos]), "%4s%n", &(SymType[(iorb+jorb)*5]),&nChar);
1833 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1834 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1835 }
1836 }
1837 //read in the vector block
1838 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*NumFuncs]), NumFuncs, imaxorb);
1839 iorb += imaxorb;
1840 Buffer->SkipnLines(1); //Skip blank line between blocks
1841 }
1842 }
1843 catch (...) {
1844 if (OrbSet) {
1845 delete OrbSet;
1846 OrbSet=NULL;
1847 }
1848 }
1849 if (OrbSet) {
1850 OrbSet->NumAlphaOrbs = NumNOrbs;
1851 OrbSet->NumOccupiedAlphaOrbs = NumNOrbs;
1852 if ((OrbSet->getNumAlphaOrbitals() + OrbSet->getNumBetaOrbitals())<=0) {
1853 delete OrbSet;
1854 OrbSet = NULL;
1855 } else {
1856 Orbs.push_back(OrbSet);
1857 }
1858 }
1859 }
1860 //Parse TD-DFT Natural Orbitals
1861 //
ParseTDDFTNOs(BufferFile * Buffer,long NumFuncs,Progress * lProgress)1862 void Frame::ParseTDDFTNOs(BufferFile * Buffer, long NumFuncs, Progress * lProgress) {
1863 long iorb, jorb, imaxorb=0, NumNOrbs=0, TestOrb,
1864 ScanErr, LinePos;
1865 int nChar;
1866 float *Vectors, *OccNums;
1867 char Line[kMaxLineLength+1];
1868 OrbitalRec * OrbSet = NULL;
1869
1870 try {
1871 Buffer->SetFilePos(Buffer->FindBlankLine());
1872 Buffer->SkipnLines(1);
1873
1874 OrbSet = new OrbitalRec(NumFuncs, 0, NumFuncs);
1875 OrbSet->setOrbitalWavefunctionType(TDDFT);
1876 OrbSet->setOrbitalType(NaturalOrbital);
1877
1878 Vectors = OrbSet->Vectors;
1879 if (OrbSet->Energy) {
1880 delete [] OrbSet->Energy;
1881 OrbSet->Energy = NULL;
1882 }
1883 if (OrbSet->SymType) {
1884 delete [] OrbSet->SymType;
1885 OrbSet->SymType = NULL;
1886 }
1887 OccNums = OrbSet->OrbOccupation = new float [NumFuncs];
1888 if (!OccNums) throw MemoryError();
1889
1890 iorb=0;
1891 NumNOrbs = NumFuncs;
1892 while (iorb<NumNOrbs) {
1893 // Allow a little backgrounding and user cancels
1894 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1895 delete OrbSet;
1896 return;
1897 }
1898 imaxorb = MIN(10, NumNOrbs-iorb); //Max of 10 orbitals per line
1899 Buffer->GetLine(Line);
1900 LinePos = 0;
1901 for (jorb=0; jorb<imaxorb; jorb++) {
1902 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1903 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1904 else {
1905 imaxorb = jorb;
1906 if (jorb==0) { //No more orbitals found
1907 imaxorb = 0;
1908 NumNOrbs = iorb;
1909 }
1910 break;
1911 }
1912 }
1913 if (imaxorb <= 0) break;
1914 //first read in the orbital occupation number of each orbital
1915 Buffer->GetLine(Line); //blank
1916 Buffer->GetLine(Line);
1917 LinePos = 0;
1918 for (jorb=0; jorb<imaxorb; jorb++) {//Grab the orbital occupations
1919 ScanErr = sscanf(&(Line[LinePos]), "%f%n", &(OccNums[iorb+jorb]),&nChar);
1920 if (ScanErr<=0) throw DataError(); //Looks like the MO's are not complete
1921 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
1922 }
1923 Buffer->GetLine(Line);
1924 //read in the vector block
1925 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*NumFuncs]), NumFuncs, imaxorb);
1926 iorb += imaxorb;
1927 Buffer->SkipnLines(1); //Skip blank line between blocks
1928 }
1929 }
1930 catch (...) {
1931 if (OrbSet) {
1932 delete OrbSet;
1933 OrbSet=NULL;
1934 }
1935 }
1936 if (OrbSet) {
1937 OrbSet->NumAlphaOrbs = NumNOrbs;
1938 OrbSet->NumOccupiedAlphaOrbs = NumNOrbs;
1939 if ((OrbSet->getNumAlphaOrbitals() + OrbSet->getNumBetaOrbitals())<=0) {
1940 delete OrbSet;
1941 OrbSet = NULL;
1942 } else {
1943 Orbs.push_back(OrbSet);
1944 }
1945 }
1946 }
1947 //Parse a set of Localized MOs from a GAMESS log file. Note that
1948 //there is no symmetry or orbital energy information available.
ParseGAMESSLMOs(BufferFile * Buffer,long NumFuncs,long NumAlphaOrbs,long NumOccBeta,Progress * lProgress,bool OrientedSet)1949 OrbitalRec * Frame::ParseGAMESSLMOs(BufferFile * Buffer, long NumFuncs, long NumAlphaOrbs,
1950 long NumOccBeta, Progress * lProgress, bool OrientedSet) {
1951 long iorb=0, imaxorb, maxfuncs, TestOrb, LinePos, ScanErr;
1952 int nChar;
1953 float *Vectors;
1954 char Line[150];
1955
1956 maxfuncs = NumFuncs; //Always true
1957 if (NumAlphaOrbs <= 0) throw DataError();
1958 OrbitalRec * TargetSet = NULL;
1959 try {
1960 TargetSet = new OrbitalRec(NumAlphaOrbs, NumOccBeta, NumFuncs);
1961 if (!OrientedSet) {
1962 TargetSet->setOrbitalType(LocalizedOrbital);
1963 } else {
1964 TargetSet->setOrbitalType(OrientedLocalizedOrbital);
1965 }
1966 // TargetSet->setOrbitalWavefunctionType(RHF);
1967 Vectors = TargetSet->Vectors;
1968 if (TargetSet->Energy) delete [] TargetSet->Energy;
1969 TargetSet->Energy = NULL;
1970 if (TargetSet->SymType) delete [] TargetSet->SymType;
1971 TargetSet->SymType = NULL;
1972 if (NumOccBeta) {
1973 if (TargetSet->EnergyB) delete [] TargetSet->EnergyB;
1974 TargetSet->EnergyB = NULL;
1975 if (TargetSet->SymTypeB) delete [] TargetSet->SymTypeB;
1976 TargetSet->SymTypeB = NULL;
1977 }
1978
1979 while (iorb<NumAlphaOrbs) {
1980 // Allow a little backgrounding and user cancels
1981 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
1982 delete TargetSet;
1983 return NULL;
1984 }
1985 imaxorb = MIN(10, NumAlphaOrbs-iorb); //Max of 10 orbitals per line
1986 Buffer->GetLine(Line);
1987 LinePos = 0;
1988 for (long jorb=0; jorb<imaxorb; jorb++) {
1989 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
1990 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
1991 else {
1992 imaxorb = jorb;
1993 if (jorb==0) { //No more orbitals found
1994 imaxorb = 0;
1995 NumAlphaOrbs = iorb;
1996 }
1997 break;
1998 }
1999 }
2000 if (imaxorb > 0) {
2001 Buffer->SkipnLines(1); //Skip blank line
2002 //read in the vector block
2003 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*maxfuncs]), maxfuncs, imaxorb);
2004 iorb += imaxorb;
2005 Buffer->SkipnLines(1); //Skip blank line between blocks
2006 }
2007 }
2008 if (NumOccBeta) {
2009 TargetSet->setOrbitalWavefunctionType(UHF);
2010 Vectors = TargetSet->VectorsB;
2011 //GAMESS now punchs out the Fock operator for ruedenburg type localization
2012 //which I need to skip if it is present
2013 if (Buffer->LocateKeyWord("FOCK OPERATOR FOR THE LOCALIZED ORBITALS",40))
2014 Buffer->SkipnLines(1);
2015 if (!Buffer->LocateKeyWord("LOCALIZED ORBITALS",18)) throw DataError();
2016 Buffer->SkipnLines(2);
2017 iorb=0;
2018
2019 while (iorb<NumOccBeta) {
2020 // Allow a little backgrounding and user cancels
2021 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
2022 delete TargetSet;
2023 return NULL;
2024 }
2025 imaxorb = MIN(10, NumOccBeta-iorb); //Max of 10 orbitals per line
2026 Buffer->GetLine(Line);
2027 LinePos = 0;
2028 for (long jorb=0; jorb<imaxorb; jorb++) {
2029 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
2030 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
2031 else {
2032 imaxorb = jorb;
2033 if (jorb==0) { //No more orbitals found
2034 imaxorb = 0;
2035 NumOccBeta = iorb;
2036 }
2037 break;
2038 }
2039 }
2040 if (imaxorb > 0) {
2041 Buffer->SkipnLines(1); //Skip blank line
2042 //read in the vector block
2043 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*maxfuncs]), maxfuncs, imaxorb);
2044 iorb += imaxorb;
2045 Buffer->SkipnLines(1); //Skip blank line between blocks
2046 }
2047 }
2048 }
2049 TargetSet->ReSize(NumAlphaOrbs, NumOccBeta);
2050 }
2051 catch (...) {
2052 if (TargetSet) {
2053 delete TargetSet;
2054 TargetSet = NULL;
2055 }
2056 }
2057 if (TargetSet) Orbs.push_back(TargetSet);
2058 //Kludge for Jan
2059 // if (!Orbs->EigenVectors && NumAlphaOrbs == 1) {
2060 // Orbs->EigenVectors = Orbs->LocalOrbitals;
2061 // Orbs->LocalOrbitals = NULL;
2062 // Orbs->NumOccupiedAlphaOrbs = 1;
2063 // if (NumOccBeta == 1) Orbs->NumOccupiedBetaOrbs = 1;
2064 // }
2065 return TargetSet;
2066 }
ParseGVBGIOrbitals(BufferFile * Buffer,const long & NumFuncs,const long & NumGVBPairs,Progress * lProgress)2067 void Frame::ParseGVBGIOrbitals(BufferFile * Buffer, const long & NumFuncs, const long & NumGVBPairs, Progress * lProgress) {
2068 long iorb, imaxorb, TestOrb, LinePos, ScanErr, jorb, NumOrbs;
2069 int nChar;
2070 char Line[kMaxLineLength+1];
2071 NumOrbs = NumFuncs;
2072
2073 for (long iPair=0; iPair<NumGVBPairs; iPair++) {
2074 // the next non-blank line should be "PAIR n"
2075 if (!Buffer->LocateKeyWord("PAIR", 4, Buffer->GetFilePos() + 250)) {
2076 wxLogMessage(_("Unable to locate an expected set of GVB pair orbitals. Skipped."));
2077 break; //If PAIR does not appear in the next few lines something is incorrect
2078 }
2079 Buffer->SkipnLines(2);
2080
2081 OrbitalRec * OrbSet = NULL;
2082 iorb=0;
2083
2084 try {
2085 OrbSet = new OrbitalRec(NumFuncs, 0, NumFuncs);
2086 OrbSet->setOrbitalWavefunctionType(GVB);
2087 OrbSet->setOrbitalType(NaturalOrbital);
2088
2089 float * Vectors = OrbSet->Vectors;
2090 if (OrbSet->Energy) {
2091 delete [] OrbSet->Energy;
2092 OrbSet->Energy = NULL;
2093 }
2094 if (OrbSet->SymType) { //symmetry labels are not printed for NO's
2095 delete [] OrbSet->SymType;
2096 OrbSet->SymType = NULL;
2097 }
2098
2099 while (iorb<NumOrbs) {
2100 // Allow a little backgrounding and user cancels
2101 if (!lProgress->UpdateProgress(Buffer->GetPercentRead())) {
2102 delete OrbSet;
2103 return;
2104 }
2105 imaxorb = MIN(10, NumOrbs-iorb); //Max of 10 orbitals per line
2106 Buffer->GetLine(Line);
2107 LinePos = 0;
2108 for (jorb=0; jorb<imaxorb; jorb++) {
2109 ScanErr = sscanf(&(Line[LinePos]), "%ld%n", &TestOrb, &nChar);
2110 if (ScanErr && (TestOrb==jorb+iorb+1)) LinePos += nChar;
2111 else {
2112 imaxorb = jorb;
2113 if (jorb==0) { //No more orbitals found
2114 imaxorb = 0;
2115 NumOrbs = iorb;
2116 }
2117 break;
2118 }
2119 }
2120 Buffer->SkipnLines(1);
2121 if (imaxorb > 0) {
2122 //read in the vector block
2123 ReadGAMESSlogVectors(Buffer, &(Vectors[iorb*NumFuncs]), NumFuncs, imaxorb);
2124 iorb += imaxorb;
2125 Buffer->SkipnLines(1); //Skip blank line between blocks
2126 }
2127 }
2128 }
2129 catch (...) {
2130 if (OrbSet) {
2131 delete OrbSet;
2132 OrbSet = NULL;
2133 }
2134 }
2135 if (OrbSet != NULL) {
2136 OrbSet->ReSize(NumOrbs, 0);
2137 Orbs.push_back(OrbSet);
2138 }
2139 Buffer->BackupnLines(2);
2140 }
2141 }
2142
ReadMP2Vectors(BufferFile * Buffer,BufferFile * DatBuffer,long NumFuncs,Progress *,long * readflag)2143 void Frame::ReadMP2Vectors(BufferFile * Buffer, BufferFile * DatBuffer, long NumFuncs, Progress * /*lProgress*/,
2144 long * readflag) {
2145
2146 if (*readflag == 1) { //If first pass query the user about reading vectors
2147 // 1 = Yes
2148 // 0 = No
2149 #ifdef __wxBuild__
2150 int result = wxMessageBox(wxT("Do you wish to read the MP2 natural orbitals from the .dat file\?"),
2151 wxT(""), wxYES_NO | wxICON_QUESTION);
2152 if(result == wxYES)
2153 *readflag = 1;
2154 else
2155 *readflag = 0;
2156 #else
2157 *readflag = YesOrNoDialog(131, 2);
2158 #endif
2159 }
2160 if (!*readflag) return;
2161
2162 if (!DatBuffer)
2163 DatBuffer = OpenDatFile();
2164 if (DatBuffer) {
2165 DatBuffer->SetFilePos(0);
2166 *readflag = 2;
2167 bool found = false;
2168 char Line[kMaxLineLength];
2169 double testEnergy;
2170 while (!found && DatBuffer->LocateKeyWord("MP2 NATURAL ORBITALS, E(MP2)=", 29,-1)) {
2171 DatBuffer->GetLine(Line);
2172 sscanf(&(Line[30]), "%lf", &testEnergy);
2173 if (fabs(testEnergy - GetMP2Energy()) < 1.0e-8) found = true;
2174 }
2175 if (found) { //Found the correct MP2 vectors
2176 OrbitalRec * OrbSet=NULL;
2177 try {
2178 OrbSet = new OrbitalRec(NumFuncs, 0, NumFuncs);
2179
2180 // OrbSet->ReadVecGroup(DatBuffer, NumFuncs, 0);
2181 OrbSet->ReadVecGroup(DatBuffer, NumFuncs, RHFMP2);
2182
2183 long lNumOrbs = OrbSet->NumAlphaOrbs;
2184 float * OccNums = OrbSet->OrbOccupation = new float [lNumOrbs];
2185 if (!OccNums) throw MemoryError();
2186 Buffer->GetLine(Line);
2187 long scanerr, LinePos;
2188 int nchar;
2189 LinePos = 0;
2190 for (long ifunc=0; ifunc<lNumOrbs; ifunc++) {
2191 scanerr = sscanf(&(Line[LinePos]), "%f%n", &(OccNums[ifunc]), &nchar);
2192 if (scanerr != 1) {
2193 Buffer->GetLine(Line);
2194 LinePos = 0;
2195 scanerr = sscanf(Line, "%f%n", &(OccNums[ifunc]), &nchar);
2196 if (scanerr != 1) throw DataError();
2197 }
2198 LinePos += nchar;
2199 }
2200 OrbSet->NumOccupiedAlphaOrbs = lNumOrbs;
2201 }
2202 catch (...) {
2203 if (OrbSet) {
2204 delete OrbSet;
2205 OrbSet = NULL;
2206 }
2207 }
2208 if (OrbSet) {
2209 OrbSet->setOrbitalWavefunctionType(RHFMP2);
2210 OrbSet->setOrbitalType(NaturalOrbital);
2211 Orbs.push_back(OrbSet);
2212 }
2213 }
2214 }
2215 }
2216
GetNumMMAtoms(void)2217 long Frame::GetNumMMAtoms(void) {
2218 long result = 0;
2219 for (int i=0; i<NumAtoms; i++)
2220 if (Atoms[i].IsSIMOMMAtom()) result++;
2221 return result;
2222 }
2223
2224 //call to read in the coefficient vectors for orbitals. This routine reads one block
2225 //of vectors which is in the same format for all wavefunction types. You should call
2226 //with the Buffer positioned at the beginning of the line containing the first
2227 //coef. and the Vectors positioned to accept it.
ReadGAMESSlogVectors(BufferFile * Buffer,float * Vectors,long MaxFuncs,long NumFuncs)2228 void ReadGAMESSlogVectors(BufferFile * Buffer, float *Vectors, long MaxFuncs, long NumFuncs) {
2229 long LinePos, iorb;
2230 int nChar;
2231 char LineText[kMaxLineLength+1];
2232
2233 for (long ifunc=0; ifunc<MaxFuncs; ifunc++) {
2234 Buffer->GetLine(LineText);
2235 LinePos = 15; //skip past the junk at the beginning of each line
2236 if (isalpha(LineText[LinePos])) LinePos++;
2237 for (iorb=0; iorb<NumFuncs; iorb++) { //get the coef for each orbital/function
2238 sscanf(&(LineText[LinePos]), "%f%n", &(Vectors[ifunc+iorb*MaxFuncs]),&nChar);
2239 LinePos+=nChar; //nChar contains the number of char's read by sscanf including spaces
2240 }
2241 }
2242 }
GetAtomTypes(void)2243 AtomTypeList * Frame::GetAtomTypes(void) {
2244 AtomTypeList * list = new AtomTypeList;
2245 if (list) {
2246 for (long i=0; i<NumAtoms; i++) {
2247 list->AddAtomType(Atoms[i].Type);
2248 }
2249 }
2250 return list;
2251 }
2252 //Read in the full set of normal modes, if present from a GAMESS log file
ParseNormalModes(BufferFile * Buffer,Progress * ProgressInd,WinPrefs * Prefs)2253 void Frame::ParseNormalModes(BufferFile * Buffer, Progress * ProgressInd, WinPrefs * Prefs) {
2254
2255 try {//catch errors locally
2256 if (Buffer->LocateKeyWord("NORMAL COORDINATE ANALYSIS", 25)) {
2257 //If there are no normal modes we just fall out the bottom and return
2258 ProgressInd->ChangeText("Reading Normal Modes");
2259 if (!ProgressInd->UpdateProgress((100*Buffer->GetFilePos())/Buffer->GetFileLength()))
2260 {throw UserCancel();}
2261
2262 //If this is a SIMOMM run then only the MM atoms have modes and are assumed
2263 //to be first in the atom list.
2264 //Update! Apparently the current code does print out any vibration information
2265 //for any MM atoms.
2266 long NumSIMOMMAtoms = GetNumMMAtoms();
2267 long NumModes = 3*NumAtoms;
2268 long NumActiveAtoms = NumAtoms;
2269 if (NumSIMOMMAtoms > 0) {
2270 NumActiveAtoms = NumAtoms - NumSIMOMMAtoms;
2271 NumModes = (NumActiveAtoms)*3;
2272 }
2273 VibRec * lVibs = Vibs = new VibRec(NumModes, NumAtoms);
2274
2275 //insufficient memory available to store the normal modes, so abort reading them in...
2276 // should alert the user to the reason he didn't get normal modes...
2277 // if (!lVibs) {
2278 // break;
2279 // }
2280 long imode=0, ifreq, LastPass=0, test, cmode;
2281 int nchar;
2282 char LineText[kMaxLineLength], token[kMaxLineLength];
2283 while (imode<NumModes) {
2284 if (!Buffer->LocateKeyWord("FREQUENCY:", 10)) break;
2285 long NumVibs = MIN(9, NumModes-imode);
2286 // long NumVibs = ((9) > (NumModes-imode)) ? (NumModes-imode) : (9);
2287 Buffer->GetLine(LineText);
2288 long LinePos = 11;
2289 long LineLength = strlen(LineText);
2290 for (unsigned long icol=0; icol<NumVibs; icol++) {
2291 if (LinePos<LineLength) {
2292 token[0]='\0';
2293 int iscan = sscanf(&(LineText[LinePos]), "%s%n", token, &nchar);
2294 LinePos += nchar;
2295 long tokenlength = strlen(token);
2296 if ((iscan>0)&&(nchar>0)&&(tokenlength>0)) {
2297 if (LineText[LinePos+1] == 'I') {
2298 token[tokenlength] = 'i';
2299 tokenlength++;
2300 token[tokenlength]=0;
2301 }
2302 lVibs->Frequencies.push_back(std::string(token));
2303 LinePos+=2;
2304 } else NumVibs = icol;
2305 } else NumVibs = icol;
2306 }
2307 imode += NumVibs;
2308 //RAMAN type of runs now allow a symmetry
2309 if (Buffer->LocateKeyWord("SYMMETRY:", 9, Buffer->GetFilePos()+132)) {
2310 Buffer->GetLine(LineText);
2311 LinePos = 10;
2312 if ((imode == NumVibs)&&(lVibs->Symmetry.empty())) {
2313 lVibs->Symmetry.reserve(NumModes);
2314 }
2315 long LineLength = strlen(LineText);
2316 for (unsigned long icol=0; icol<NumVibs; icol++) {
2317 if (LinePos<LineLength) {
2318 token[0]='\0';
2319 int iscan = sscanf(&(LineText[LinePos]), "%s%n", token, &nchar);
2320 LinePos += nchar;
2321 long test = strlen(token);
2322 if ((iscan>0)&&(nchar>0)&&(test>0)) {
2323 lVibs->Symmetry.push_back(std::string(token));
2324 } else NumVibs = icol;
2325 } else NumVibs = icol;
2326 }
2327 }
2328 if (Buffer->LocateKeyWord("REDUCED MASS:", 13, Buffer->GetFilePos()+132)) {
2329 Buffer->GetLine(LineText);
2330 LinePos = 14;
2331 if ((imode == NumVibs)&&(lVibs->ReducedMass.empty())) {
2332 lVibs->ReducedMass.reserve(NumModes);
2333 }
2334 LineLength = strlen(LineText);
2335 long tVib = NumVibs;
2336 float rmass;
2337 for (long icol=0; icol<tVib; icol++) {
2338 if (LinePos<LineLength) {
2339 test = sscanf(&(LineText[LinePos]), "%255s%n", token, &nchar);
2340 LinePos += nchar;
2341 if (test) {
2342 if (token[0] != '*') {
2343 test = sscanf(token, "%f", &rmass);
2344 if (test)
2345 lVibs->ReducedMass.push_back(rmass);
2346 } else lVibs->ReducedMass.push_back(10000.0);
2347 } else tVib = icol;
2348 } else tVib = icol;
2349 }
2350 }
2351 if (Buffer->LocateKeyWord("INTENSITY:", 10, Buffer->GetFilePos()+132)) {
2352 Buffer->GetLine(LineText);
2353 LinePos = 11;
2354 LineLength = strlen(LineText);
2355 long tVib = NumVibs;
2356 float Inten;
2357 for (long icol=0; icol<tVib; icol++) {
2358 if (LinePos<LineLength) {
2359 test = sscanf(&(LineText[LinePos]), "%255s%n", token, &nchar);
2360 LinePos += nchar;
2361 if (test) {
2362 if (token[0] != '*') {
2363 test = sscanf(token, "%f", &Inten);
2364 if (test)
2365 lVibs->Intensities[icol+LastPass] = Inten;
2366 } else lVibs->Intensities[icol+LastPass] = 10000.0;
2367 } else tVib = icol;
2368 } else tVib = icol;
2369 }
2370 }
2371 //In the 2010 version Mike is changing Raman intensity to Raman Activity
2372 if ((Buffer->LocateKeyWord("RAMAN INTENSITY:", 16, Buffer->GetFilePos()+132))||
2373 (Buffer->LocateKeyWord("RAMAN ACTIVITY:", 15, Buffer->GetFilePos()+132))) {
2374 Buffer->GetLine(LineText);
2375 LinePos = 17;
2376 if ((imode == NumVibs)&&(lVibs->RamanIntensity.empty())) {
2377 lVibs->RamanIntensity.reserve(NumModes);
2378 }
2379 LineLength = strlen(LineText);
2380 long tVib = NumVibs;
2381 float raman;
2382 for (long icol=0; icol<tVib; icol++) {
2383 if (LinePos<LineLength) {
2384 test = sscanf(&(LineText[LinePos]), "%255s%n", token, &nchar);
2385 LinePos += nchar;
2386 if (test) {
2387 if (token[0] != '*') {
2388 test = sscanf(token, "%f", &raman);
2389 if (test)
2390 lVibs->RamanIntensity.push_back(raman);
2391 } else lVibs->RamanIntensity.push_back(10000.0);
2392 } else tVib = icol;
2393 } else tVib = icol;
2394 }
2395 }
2396 if (Buffer->LocateKeyWord("DEPOLARIZATION:", 15, Buffer->GetFilePos()+132)) {
2397 Buffer->GetLine(LineText);
2398 LinePos = 16;
2399 if ((imode == NumVibs)&&(lVibs->Depolarization.empty())) {
2400 lVibs->Depolarization.reserve(NumModes);
2401 }
2402 LineLength = strlen(LineText);
2403 long tVib = NumVibs;
2404 float depol;
2405 for (long icol=0; icol<tVib; icol++) {
2406 if (LinePos<LineLength) {
2407 test = sscanf(&(LineText[LinePos]), "%255s%n", token, &nchar);
2408 LinePos += nchar;
2409 if (test) {
2410 if (token[0] != '*') {
2411 test = sscanf(token, "%f", &depol);
2412 if (test)
2413 lVibs->Depolarization.push_back(depol);
2414 } else lVibs->Depolarization.push_back(10000.0);
2415 } else tVib = icol;
2416 } else tVib = icol;
2417 }
2418 }
2419 test = Buffer->FindBlankLine();
2420 Buffer->SetFilePos(test);
2421 Buffer->SkipnLines(1);
2422
2423 for (test=NumSIMOMMAtoms; test<NumAtoms; test++) {
2424 if (!ProgressInd->UpdateProgress((100*Buffer->GetFilePos())/Buffer->GetFileLength()))
2425 {throw UserCancel();}
2426
2427 Buffer->GetLine(LineText);
2428 LinePos = 20;
2429 float lmass = Prefs->GetSqrtAtomMass(Atoms[test].GetType()-1);
2430 for (ifreq=LastPass; ifreq<imode; ifreq++) {
2431 cmode = test + ifreq*NumAtoms;
2432 int iscan = sscanf(&(LineText[LinePos]), "%f%n", &((lVibs->NormMode[cmode]).x), &nchar);
2433 if (iscan <= 0) throw DataError();
2434 LinePos += nchar;
2435 (lVibs->NormMode[cmode]).x *= lmass;
2436 }
2437 Buffer->GetLine(LineText);
2438 LinePos = 20;
2439 for (ifreq=LastPass; ifreq<imode; ifreq++) {
2440 cmode = test + ifreq*NumAtoms;
2441 int iscan = sscanf(&(LineText[LinePos]), "%f%n", &((lVibs->NormMode[cmode]).y), &nchar);
2442 if (iscan <= 0) throw DataError();
2443 LinePos += nchar;
2444 (lVibs->NormMode[cmode]).y *= lmass;
2445 }
2446 Buffer->GetLine(LineText);
2447 LinePos = 20;
2448 for (ifreq=LastPass; ifreq<imode; ifreq++) {
2449 cmode = test + ifreq*NumAtoms;
2450 int iscan = sscanf(&(LineText[LinePos]), "%f%n", &((lVibs->NormMode[cmode]).z), &nchar);
2451 if (iscan <= 0) throw DataError();
2452 LinePos += nchar;
2453 (lVibs->NormMode[cmode]).z *= lmass;
2454 }
2455 }
2456 //ok apparently its the simomm atoms that are missing...
2457 // if (NumSIMOMMAtoms > 0) { //zero out remaining atoms for SIMOMM runs
2458 // for (int i=NumSIMOMMAtoms; i<NumAtoms; i++) {
2459 // cmode = i+ifreq*NumAtoms;
2460 // (lVibs->NormMode[cmode]).x = 0.0;
2461 // (lVibs->NormMode[cmode]).y = 0.0;
2462 // (lVibs->NormMode[cmode]).z = 0.0;
2463 // }
2464 // }
2465 LastPass = imode;
2466 sprintf(LineText, "Read in %ld normal modes", LastPass);
2467 ProgressInd->ChangeText(LineText);
2468 if (!ProgressInd->UpdateProgress((100*Buffer->GetFilePos())/Buffer->GetFileLength()))
2469 {throw UserCancel();}
2470 }
2471 lVibs->NumModes = LastPass;
2472 }
2473 } //trap errors here and delete the VibRec
2474 catch (std::bad_alloc) {//Memory error, cleanup and return.
2475 if (Vibs) { delete Vibs; Vibs = NULL; }
2476 MessageAlert("Insufficient memory to read normal modes. Normal Modes will be skipped!");
2477 }
2478 catch (UserCancel) {//We need to rethrow this one since the whole operation should be aborted
2479 if (Vibs) { delete Vibs; Vibs = NULL; }
2480 throw;
2481 }
2482 catch (...) {//File and data errors. Either way delete the vectors and return.
2483 if (Vibs) { delete Vibs; Vibs = NULL; }
2484 MessageAlert("Error parsing normal modes. Normal modes will be skipped.");
2485 }
2486 }
ParseMolDenFrequencies(BufferFile * Buffer,WinPrefs * Prefs)2487 void Frame::ParseMolDenFrequencies(BufferFile * Buffer, WinPrefs * Prefs) {
2488 try {//catch errors locally
2489 if (Buffer->LocateKeyWord("[FREQ]", 6)) {
2490 Buffer->SkipnLines(1);
2491 long NumModes = 3*NumAtoms;
2492 VibRec * lVibs = Vibs = new VibRec(NumModes, NumAtoms);
2493
2494 char LineText[kMaxLineLength], token[kMaxLineLength];
2495 for (int i=0; i<NumModes; i++) {
2496 Buffer->GetLine(LineText);
2497 if (LineText[0] == '[') {
2498 NumModes = i;
2499 break;
2500 }
2501 int c = sscanf(LineText, "%s", token);
2502 if (c == 1) {
2503 lVibs->Frequencies.push_back(std::string(token));
2504 } else {
2505 NumModes = i;
2506 break;
2507 }
2508 if (Buffer->GetFilePos() >= Buffer->GetFileSize()) {
2509 NumModes = i;
2510 break;
2511 }
2512 }
2513 Buffer->SetFilePos(0);
2514 if (Buffer->LocateKeyWord("[INT]", 5)) {
2515 Buffer->SkipnLines(1);
2516 char token2[kMaxLineLength];
2517 for (int i=0; i<NumModes; i++) {
2518 Buffer->GetLine(LineText);
2519 if (LineText[0] == '[') break;
2520 int c = sscanf(LineText, "%s %s", token, token2);
2521 float inten;
2522 if (c >= 1) {
2523 int check = sscanf(token, "%f", &inten);
2524 if (check != 1) inten = 10000.0;
2525 lVibs->Intensities[i] = inten;
2526 if (c == 2) {
2527 check = sscanf(token, "%f", &inten);
2528 if (check != 1) inten = 10000.0;
2529 lVibs->RamanIntensity.push_back(inten);
2530 }
2531 } else break;
2532 if (Buffer->GetFilePos() >= Buffer->GetFileSize()) break;
2533 }
2534 }
2535 Buffer->SetFilePos(0);
2536 if (Buffer->LocateKeyWord("[FR-NORM-COORD]", 14)) {
2537 Buffer->SkipnLines(1);
2538 for (int i=0; i<NumModes; i++) {
2539 Buffer->GetLine(LineText);
2540 int check;
2541 sscanf(LineText, "%s %d", token, &check);
2542 if (check != (i+1)) {
2543 NumModes = i;
2544 break;
2545 }
2546 for (int j=0; j<NumAtoms; j++) {
2547 Buffer->GetLine(LineText);
2548 int imode = j + i*NumAtoms;
2549 check = sscanf(LineText, "%f %f %f", &(lVibs->NormMode[imode].x), &(lVibs->NormMode[imode].y),
2550 &(lVibs->NormMode[imode].z));
2551 float lmass = Prefs->GetSqrtAtomMass(Atoms[j].GetType()-1);
2552 lVibs->NormMode[imode] *= lmass;
2553 }
2554 }
2555 }
2556 lVibs->NumModes = NumModes;
2557 }
2558 } //trap errors here and delete the VibRec
2559 catch (std::bad_alloc) {//Memory error, cleanup and return.
2560 if (Vibs) { delete Vibs; Vibs = NULL; }
2561 MessageAlert("Insufficient memory to read normal modes. Normal Modes will be skipped!");
2562 }
2563 catch (...) {//File and data errors. Either way delete the vectors and return.
2564 if (Vibs) { delete Vibs; Vibs = NULL; }
2565 MessageAlert("Error parsing normal modes. Normal modes will be skipped.");
2566 }
2567 }
2568
toggleMMAtomVisibility(void)2569 void Frame::toggleMMAtomVisibility(void) {
2570 for (int i=0; i<NumAtoms; i++) {
2571 if (Atoms[i].IsSIMOMMAtom())
2572 Atoms[i].SetInvisibility(true-Atoms[i].GetInvisibility());
2573 }
2574 }
2575
toggleAbInitioVisibility(void)2576 void Frame::toggleAbInitioVisibility(void) {
2577 for (int i=0; i<NumAtoms; i++) {
2578 if (! (Atoms[i].IsSIMOMMAtom()))
2579 Atoms[i].SetInvisibility(true-Atoms[i].GetInvisibility());
2580 }
2581 }
2582
toggleEFPVisibility(void)2583 void Frame::toggleEFPVisibility(void) {
2584 for (int i=0; i<NumAtoms; i++) {
2585 if (Atoms[i].IsEffectiveFragment())
2586 Atoms[i].SetInvisibility(true-Atoms[i].GetInvisibility());
2587 }
2588 }
2589
resetAllSelectState()2590 void Frame::resetAllSelectState() {
2591
2592 // Typically, atom selection should be done through SetAtomSelection, but
2593 // we can simulate the effect here.
2594 for (int i = 0; i < NumAtoms; i++)
2595 Atoms[i].SetSelectState(false);
2596 natoms_selected = 0;
2597
2598 for (int i = 0; i < NumBonds; i++)
2599 Bonds[i].SetSelectState(false);
2600 }
2601
GetAtomNumBonds(int atom_id) const2602 int Frame::GetAtomNumBonds(int atom_id) const {
2603
2604 int num_bonds = 0;
2605
2606 for (int i = 0; i < NumBonds; i++) {
2607 if (Bonds[i].Atom1 == atom_id || Bonds[i].Atom2 == atom_id) {
2608 if (Bonds[i].Order > kHydrogenBond) num_bonds++;
2609 }
2610 }
2611
2612 return num_bonds;
2613
2614 }
2615