1 /*
2 * (c) 2004 Iowa State University
3 * see the LICENSE file in the top level directory
4 */
5
6 /* MEP.cpp
7 Molecular Electrostatic Potential calculations
8 Adapted from GAMESS by BMB 10/1997
9 bug fixes and MEP surface code moved here by BMB 1/1998
10 Removed Mac window dependancies BMB 4/1998
11 Fixed MP stack size problem in MEP3D surfaces BMB 9/2000
12 */
13
14 #include "Globals.h"
15 #include "GlobalExceptions.h"
16 #include "Progress.h"
17 #include "MoleculeData.h"
18 #include "Frame.h"
19 #include "BasisSet.h"
20 #include "RysPolynomial.h"
21 #include "GaussHermite.h"
22 #include "SurfaceTypes.h"
23 #if defined(powerc) && defined(MacintoshBuild)
24 #include <Multiprocessing.h>
25 #endif
26 #ifdef __wxBuild__
27 #include <wx/stopwatch.h>
28 #endif
29 #include <iostream>
30
31 extern Boolean gMPAware;
32 extern long gNumProcessors;
33
34 //Clearly there is some work needed to support H and I functions
35 //Some is dimensions and powers. Looks like Root6 in RysPol has been completely redone as well, though
36 //that might not be necessary. Some expansion for the number of roots is needed.
37
38 #define kMaxShellType 6 // 0 index, 6 gives I shell
39 #define kMaxAngularMomentumSum 84 //s-g functions - 1+3+6+10+15+21(h)+(28(I)
40 class AngularIndeces {
41 public:
42 short x[kMaxAngularMomentumSum], y[kMaxAngularMomentumSum], z[kMaxAngularMomentumSum];
43 AngularIndeces();
44 };
45
46 //setup arrays of x, y, z powers for each angular momenta (through g)
AngularIndeces()47 AngularIndeces::AngularIndeces() {
48 for (long i=0; i<kMaxAngularMomentumSum; i++)
49 x[i] = y[i] = z[i] = 0;
50 x[1] = y[2] = z[3] = 1; //P x, y, z
51 x[4] = y[5] = z[6] = 2; //D x2, y2, z2
52 x[7] = x[8] = y[7] = y[9] = z[8] = z[9] = 1;//D xy,xz,yz
53 x[10] = y[11] = z[12] = 3; //F x3, y3, z3
54 x[13] = x[14] = y[15] = y[16] = z[17] = z[18] = 2;//F x2y, x2z, xy2, y2z, xz2, yz2, xyz
55 x[15] = x[17] = x[19] = y[13] = y[18] = y[19] = z[14] = z[16] = z[19] = 1;
56 x[20] = y[21] = z[22] = 4; //G x4, y4, z4
57 x[23] = x[24] = y[25] = y[26] = z[27] = z[28] = 3; //G x3y, x3z, xy3, y3z, xz3, yz3
58 x[29] = x[30] = x[32] = //G x2y2, x2z2, y2z2, x2yz, xy2z, xyz2
59 y[29] = y[31] = y[33] = z[30] = z[31] = z[34] = 2;
60 x[25] = x[27] = x[33] = x[34] = y[23] = y[28] = y[32] = y[34] =
61 z[24] = z[26] = z[32] = z[33] = 1;
62 x[35] = y[36] = z[37] = 5; //H x5, y5, z5
63 x[38] = x[39] = y[40] = y[41] = z[42] = z[43] = 4; //H x4y, x4z, xy4, y4z, xz4, yz4
64 y[38] = z[39] = x[40] = z[41] = x[42] = y[43] = 1;
65 x[44] = x[45] = y[46] = y[47] = z[48] = z[49] = 3; //H x3y2, x3z2, x2y3, y3z2, x2z3, y2z3
66 y[44] = z[45] = x[46] = z[47] = x[48] = y[49] = 2;
67 x[50] = y[51] = z[52] = 3; //H x3yz, xy3z, xyz3
68 y[50] = z[50] = x[51] = z[51] = x[52] = y[52] = 1;
69 x[53] = y[53] = x[54] = z[54] = y[55] = z[55] = 2; //H x2y2z, x2yz2, xy2z2
70 z[53] = y[54] = x[55] = 1;
71 x[56] = y[57] = z[58] = 6; //I x6, y6, z6
72 x[59] = x[60] = y[61] = y[62] = z[63] = z[64] = 5; //I x5y, x5z, xy5, y5z, xz5, yz5
73 y[59] = z[60] = x[61] = z[62] = x[63] = y[64] = 1;
74 x[65] = x[66] = y[67] = y[68] = z[69] = z[70] = 4; //I x4y2, x4z2, x2y4, y4z2, x2z4, y2z4
75 y[65] = z[66] = x[67] = z[68] = x[69] = y[70] = 2;
76 x[71] = y[72] = z[73] = 4; //I x4yz, xy4z, xyz4
77 y[71] = z[71] = x[72] = z[72] = x[73] = y[73] = 1;
78 x[74] = y[74] = x[75] = z[75] = y[76] = z[76] = 3; //I x3y3, x3z3, y3z3
79 x[77] = x[78] = y[79] = y[80] = z[81] = z[82] = 3; //I x3y2z, x3yz2, x2y3z, xy3z2, x2yz3, xy2z3
80 y[77] = z[78] = x[79] = z[80] = x[81] = y[82] = 2;
81 z[77] = y[78] = z[79] = x[80] = y[81] = x[82] = 1;
82 x[83] = y[83] = z[83] = 2; //I x2y2z2
83 }
84 //NOTE: This routine assumes that the atomic coordinates are in Angstroms and are thus
85 //converted to Bohr when used.
CalculateMOGrid(MoleculeData * lData,Progress * lProgress)86 void MEP2DSurface::CalculateMOGrid(MoleculeData *lData, Progress * lProgress) {
87 long NumPoints = NumGridPoints*NumGridPoints;
88 GaussHermiteData GHData;
89
90 //ProfilerClear();
91 //ProfilerSetStatus(1);
92 InitGaussHermite(&GHData); //setup the Gauss-Hermite roots and weights data
93
94 Frame * lFrame = lData->GetCurrentFramePtr();
95 BasisSet * Basis = lData->GetBasisSet();
96 AODensity * TotalAODensity = lFrame->GetAODensity(Basis, OrbSet);
97
98 if ((Grid && (GridAllocation != NumPoints))||!TotalAODensity) {
99 FreeGrid();
100 }
101
102 if (!Grid) {
103 //Attempt to allocate memory for the 2D grid
104 AllocateGrid(NumPoints);
105 }
106 //If the memory allocation failed the MO calculation can not be done
107 if ((Grid == NULL)||(TotalAODensity==NULL)) return;
108 //If sufficient memory is available then setup pointers for fast local use
109 float * lGrid;
110 lGrid = Grid;
111 //Store the Grid mins and incs at the beginning of the grid
112 GridMax = -1.0e20;
113 GridMin = 1.0e20;
114 //loop over the plotting grid in the x, y, and z directions
115 long n=0;
116 float XGridValue, YGridValue, ZGridValue, junk1, junk2;
117 CPoint3D lOrigin, lXInc, lYInc;
118
119 lOrigin = Origin * kAng2BohrConversion;
120 lXInc = XInc * kAng2BohrConversion;
121 lYInc = YInc * kAng2BohrConversion;
122 for (long iXPt=0; iXPt<NumGridPoints; iXPt++) {
123 //Give up the CPU and check for cancels
124 if (!lProgress->UpdateProgress(100*iXPt/NumGridPoints)) { //User canceled so clean things up and abort
125 FreeGrid();
126 return;
127 }
128
129 XGridValue = lOrigin.x + iXPt*lXInc.x;
130 YGridValue = lOrigin.y + iXPt*lXInc.y;
131 ZGridValue = lOrigin.z + iXPt*lXInc.z;
132 for (long iYPt=0; iYPt<NumGridPoints; iYPt++) {
133 lGrid[n] = lFrame->CalculateMEP(XGridValue, YGridValue, ZGridValue,
134 Basis, TotalAODensity, &GHData, &junk1, &junk2);
135
136 GridMax = MAX(GridMax, lGrid[n]);
137 GridMin = MIN(GridMin, lGrid[n]);
138 n++;
139 XGridValue += lYInc.x;
140 YGridValue += lYInc.y;
141 ZGridValue += lYInc.z;
142 }
143 }
144 //ProfilerSetStatus(0);
145 //ProfilerDump("\pMEP profile");
146 }
147 #ifdef powerc
148 typedef struct MEP3DGridData {
149 MEP3DSurface * Surf;
150 long xStart;
151 long xEnd;
152 Frame * lFrame;
153 BasisSet * Basis;
154 AODensity * TotalAODensity;
155 GaussHermiteData * GHData;
156 float GridMax;
157 long PercentDone;
158 MPQueueID TerminationQueueID;
159 } MEP3DGridData;
Calc3DMEPGrid(MEP3DGridData * Data)160 static OSErr Calc3DMEPGrid(MEP3DGridData * Data) {
161 Data->GridMax = Data->Surf->CalculateGrid(Data->xStart, Data->xEnd, Data->lFrame, Data->Basis,
162 Data->TotalAODensity, Data->GHData,
163 NULL, &(Data->PercentDone), true);
164 return noErr;
165 }
166 #endif
167 #ifdef __wxBuild__
168 typedef struct MEP3DGridData {
169 MEP3DSurface * Surf;
170 long xStart;
171 long xEnd;
172 Frame * lFrame;
173 BasisSet * Basis;
174 AODensity * TotalAODensity;
175 GaussHermiteData * GHData;
176 long PercentDone;
177 float GridMax;
178 } MEP3DGridData;
179 //derive a class from wxThread that we will use to create the child threads
180 class MEP3DThread : public wxThread {
181 public:
182 MEP3DThread(MEP3DGridData * data);
183 virtual void * Entry();
184
GetPercentDone(void) const185 long GetPercentDone(void) const {return myData->PercentDone;};
GetGridMax(void) const186 float GetGridMax(void) const {return myData->GridMax;};
187
188 MEP3DGridData * myData;
189 };
MEP3DThread(MEP3DGridData * data)190 MEP3DThread::MEP3DThread(MEP3DGridData * data) : wxThread(wxTHREAD_JOINABLE) {
191 myData = data;
192 }
Entry()193 void * MEP3DThread::Entry() {
194 myData->GridMax = myData->Surf->CalculateGrid(myData->xStart, myData->xEnd, myData->lFrame, myData->Basis,
195 myData->TotalAODensity, myData->GHData,
196 NULL, &(myData->PercentDone), true);
197 myData->PercentDone = 100;
198 return NULL;
199 }
200 typedef MEP3DThread * MEP3DThreadPtr;
201
202 #endif
CalculateGrid(long xStart,long xEnd,Frame * lFrame,BasisSet * Basis,AODensity * TotalAODensity,GaussHermiteData * GHData,Progress * lProgress,long * PercentDone,bool MPTask)203 float MEP3DSurface::CalculateGrid(long xStart, long xEnd, Frame * lFrame, BasisSet * Basis,
204 AODensity * TotalAODensity, GaussHermiteData * GHData,
205 Progress * lProgress, long * PercentDone, bool MPTask) {
206 float lGridMax = -1.0e20, junk1, junk2;
207 float lGridMin = 1.0e20;
208 float XGridValue, YGridValue, ZGridValue;
209 long iXPt, iYPt, iZPt;
210 long n=xStart*(NumYGridPoints*NumZGridPoints);
211 float * lGrid = Grid;
212
213 XGridValue = Origin.x + xStart * XGridInc;
214 for (iXPt=xStart; iXPt<xEnd; iXPt++) {
215 //Note the percent done for MP tasks such that the progress dlog can be updated
216 if (MPTask) *PercentDone = 100*(iXPt-xStart)/(xEnd-xStart);
217 YGridValue = Origin.y;
218 for (iYPt=0; iYPt<NumYGridPoints; iYPt++) {
219 if (!MPTask) {
220 //Give up the CPU and check for cancels
221 if (!lProgress->UpdateProgress(100*iXPt/xEnd)) { //User canceled so clean things up and abort
222 FreeGrid();
223 return 0.0;
224 }
225 }
226 ZGridValue = Origin.z;
227 for (iZPt=0; iZPt<NumZGridPoints; iZPt++) {
228 lGrid[n] = lFrame->CalculateMEP(XGridValue, YGridValue, ZGridValue,
229 Basis, TotalAODensity, GHData, &junk1, &junk2);
230 lGridMax = MAX(lGridMax, lGrid[n]);
231 lGridMin = MIN(lGridMin, lGrid[n]);
232 n++;
233 ZGridValue += ZGridInc;
234 #ifdef __wxBuild__
235 if (MPTask) {
236 if ((wxThread::This())->TestDestroy()) {
237 //Getting here means the caller has requested a cancel so exit
238 return 0.0;
239 // Exit(); //Exit gracefully exits the thread
240 }
241 }
242 #endif
243 }
244 YGridValue += YGridInc;
245 }
246 XGridValue += XGridInc;
247 }
248 lGridMax = MAX(lGridMax, fabs(lGridMin));
249 return lGridMax;
250 }
251 //Generates a 3D grid over the specified MO in the plane specified.
252 //NOTE: This routine assumes that the atomic coordinates are in Angstroms and are thus
253 //converted to Bohr when used.
CalculateMEPGrid(MoleculeData * lData,Progress * lProgress)254 void MEP3DSurface::CalculateMEPGrid(MoleculeData *lData, Progress * lProgress) {
255 GaussHermiteData GHData;
256
257 InitGaussHermite(&GHData); //setup the Gauss-Hermite roots and weights data
258
259 Frame * lFrame = lData->GetCurrentFramePtr();
260 BasisSet * Basis = lData->GetBasisSet();
261 AODensity * TotalAODensity = lFrame->GetAODensity(Basis, OrbSet);
262 if (!TotalAODensity) throw MemoryError();
263
264 SetupGridParameters(lFrame); //Define the 3D volume
265
266 //allocate the memory needed for the grid
267 long NumPoints = NumXGridPoints*NumYGridPoints*NumZGridPoints;
268
269 if (Grid && (GridAllocation != NumPoints)) FreeGrid();
270 if (!Grid) {
271 //Attempt to allocate memory for the 2D grid
272 AllocateGrid(NumPoints);
273 }
274 //If the memory allocation failed the MO calculation can not be done
275 if (Grid == NULL) return;
276 //If sufficient memory is available then setup pointers for fast local use
277 float * lGrid;
278 lGrid = Grid;
279 //Get the appropriate MO vector
280 //Store the Grid mins and incs at the beginning of the grid
281 GridMax = -1.0e20;
282 #if defined(powerc) && defined(MacintoshBuild)
283 if (gMPAware && (gNumProcessors>=1)) {
284 OSErr status;
285 long i, start=0;
286 MEP3DGridData * DataPtrs = new MEP3DGridData[gNumProcessors];
287 if (!DataPtrs) throw MemoryError();
288 //Test allocations
289 MPQueueID * TerminationQueues = new MPQueueID[gNumProcessors];
290 MPTaskID * TaskIDs = new MPTaskID[gNumProcessors];
291 if (!TerminationQueues || !TaskIDs) {
292 delete [] DataPtrs;
293 throw MemoryError();
294 }
295 for (i=0; i<gNumProcessors; i++) {
296 status = MPCreateQueue(&(TerminationQueues[i])); /* Create the queue which will report the completion of the task. */
297 if (status != noErr) throw DataError();
298 DataPtrs[i].Surf = this;
299 DataPtrs[i].TerminationQueueID = TerminationQueues[i];
300 DataPtrs[i].xStart = start;
301 start += NumXGridPoints/gNumProcessors;
302 if (i+1 == gNumProcessors) start = NumXGridPoints;
303 DataPtrs[i].xEnd = start;
304 DataPtrs[i].lFrame = lFrame;
305 DataPtrs[i].Basis = Basis;
306 DataPtrs[i].TotalAODensity = TotalAODensity;
307 DataPtrs[i].GHData = &GHData;
308 DataPtrs[i].PercentDone = 0;
309 status = MPCreateTask((TaskProc)Calc3DMEPGrid, /* This is the task function. */
310 &(DataPtrs[i]), /* This is the parameter to the task function. */
311 16*1024, /* stack size, 0 will use system default of 4KB which is too small for MEPs */
312 TerminationQueues[i], /* We'll use this to sense task completion. */
313 0, /* We won't use the first part of the termination message. */
314 0, /* We won't use the second part of the termination message. */
315 0, /* Use the normal task options. (Currently this *must* be zero!) */
316 &(TaskIDs[i])); /* Here's where the ID of the new task will go. */
317 if (status != noErr) throw DataError();
318 }
319 long NumTasksRunning=gNumProcessors;
320 long TotalPercent;
321 lProgress->SetRunTime(1);
322 lProgress->SetSleepTime(20);
323 while (NumTasksRunning > 0) {
324 TotalPercent = 0;
325 for (i=0; i<gNumProcessors; i++) {
326 TotalPercent += DataPtrs[i].PercentDone;
327 if (TerminationQueues[i] != 0) {
328 long temp;
329 status = MPWaitOnQueue(TerminationQueues[i], 0, 0, (void **) &temp, kDurationImmediate);
330 if (status == noErr) {
331 GridMax = MAX(GridMax, DataPtrs[i].GridMax);
332 status = MPDeleteQueue(TerminationQueues[i]);
333 TerminationQueues[i] = 0;
334 NumTasksRunning --;
335 }
336 }
337 }
338 TotalPercent /= gNumProcessors;
339 //Give up the CPU and check for cancels
340 if (!lProgress->UpdateProgress(TotalPercent)) { //User canceled so clean things up and abort
341 for (i=0; i<gNumProcessors; i++) {
342 if (TerminationQueues[i] != 0) {
343 MPTerminateTask(TaskIDs[i], noErr);
344 long temp;
345 status = MPWaitOnQueue(TerminationQueues[i], 0, 0, (void **) &temp, kDurationForever);
346 if (status == noErr) {
347 status = MPDeleteQueue(TerminationQueues[i]);
348 TerminationQueues[i] = 0;
349 NumTasksRunning --;
350 }
351 }
352 }
353 FreeGrid();
354 }
355 }
356 delete [] TerminationQueues;
357 delete [] TaskIDs;
358 delete [] DataPtrs;
359 lProgress->ResetTimes();
360 } else
361 #endif
362 #ifdef __wxBuild__
363 if (wxThread::GetCPUCount() > 1) {
364 int myCPUCount = wxThread::GetCPUCount();
365 //Ok we have multiple CPUs so create a worker thread for each CPU and split up the grid calculation
366 //not sure if this is needed. It appears the default is 1?
367 // wxThread::SetConcurrency(myCPUCount);
368
369 MEP3DGridData * DataPtrs = new MEP3DGridData[myCPUCount];
370 //Clear pointers in case allocation fails and I need to delete part of them
371 int i;
372 MEP3DThreadPtr * myThreads = new MEP3DThreadPtr[myCPUCount];
373 long start=0;
374 for (i=0; i<myCPUCount; i++) {
375 DataPtrs[i].Surf = this;
376 DataPtrs[i].xStart = start;
377 start += NumXGridPoints/myCPUCount;
378 if (i+1 == myCPUCount) start = NumXGridPoints;
379 DataPtrs[i].xEnd = start;
380 DataPtrs[i].lFrame = lFrame;
381 DataPtrs[i].Basis = Basis;
382 DataPtrs[i].TotalAODensity = TotalAODensity;
383 DataPtrs[i].GHData = &GHData;
384 DataPtrs[i].PercentDone = 0;
385
386 //Create the actual thread
387 myThreads[i] = new MEP3DThread(&(DataPtrs[i]));
388 myThreads[i]->Create(128*1024);
389 //and fire it off, note my class is always a joinable thread so we will wait on it eventually
390 myThreads[i]->Run();
391 }
392
393 long NumTasksRunning=myCPUCount;
394 long TotalPercent;
395 while (NumTasksRunning > 0) {
396 TotalPercent = 0;
397 for (i=0; i<myCPUCount; i++) {
398 TotalPercent += DataPtrs[i].PercentDone;
399 if (myThreads[i] != NULL) {
400 if (! myThreads[i]->IsAlive()) {
401 //The thread is terminating
402 myThreads[i]->Wait();
403 GridMax = MAX(GridMax, DataPtrs[i].GridMax);
404 delete myThreads[i];
405 myThreads[i] = NULL;
406 NumTasksRunning --;
407 }
408 }
409 }
410 TotalPercent /= myCPUCount;
411 //Give up the CPU and check for cancels
412 if (!lProgress->UpdateProgress(TotalPercent)) { //User canceled so clean things up and abort
413 for (i=0; i<myCPUCount; i++) {
414 if (myThreads[i] != NULL) {
415 //The docs seem to indicate that Wait causes a thread to terminate,
416 //however this does not seem to be true. Wait appears to have the
417 //expected semantics that it blocks until the thread exits. Delete
418 //appears to be what is needed.
419 myThreads[i]->Delete();
420 delete myThreads[i];
421 myThreads[i] = NULL;
422 NumTasksRunning --;
423 }
424 }
425 GridMax = 1.0;
426 FreeGrid();
427 } else {
428 wxMilliSleep(100);
429 }
430 }
431
432 delete [] DataPtrs;
433 delete [] myThreads;
434 } else
435 #endif
436 {
437 long junk; //just a placeholder when cooperative task
438 GridMax = CalculateGrid(0,NumXGridPoints,lFrame, Basis,
439 TotalAODensity, &GHData, lProgress, &junk, false);
440 }
441 //Unlock the grid handle and return
442 Origin *= kBohr2AngConversion;
443 XGridInc *= kBohr2AngConversion;
444 YGridInc *= kBohr2AngConversion;
445 ZGridInc *= kBohr2AngConversion;
446 }
447 #ifdef __wxBuild__
448 typedef struct TED3DSurfGridData {
449 Surf3DBase * Surf;
450 MoleculeData * MolData;
451 AODensity * TotalDensity;
452 long xStart;
453 long xEnd;
454 long PercentDone;
455 GaussHermiteData * GHData;
456 float * grid;
457 CPoint3D * Contour;
458 } TED3DSurfGridData;
459 //derive a class from wxThread that we will use to create the child threads
460 class TEDMEP3DThread : public wxThread {
461 public:
462 TEDMEP3DThread(TED3DSurfGridData * data);
463 virtual void * Entry();
464
GetPercentDone(void) const465 long GetPercentDone(void) const {return myData->PercentDone;};
466
467 TED3DSurfGridData * myData;
468 };
TEDMEP3DThread(TED3DSurfGridData * data)469 TEDMEP3DThread::TEDMEP3DThread(TED3DSurfGridData * data) : wxThread(wxTHREAD_JOINABLE) {
470 myData = data;
471 }
Entry()472 void * TEDMEP3DThread::Entry() {
473 myData->Surf->CalculateSurfaceValuesGrid(myData->MolData, NULL, myData->TotalDensity,
474 myData->xStart, myData->xEnd, myData->GHData,
475 myData->grid, myData->Contour,
476 &(myData->PercentDone), true);
477 myData->PercentDone = 100;
478 return NULL;
479 }
480 typedef TEDMEP3DThread * TEDMEP3DThreadPtr;
481
482 #endif
CalculateSurfaceValues(MoleculeData * lData,Progress * lProgress)483 void Surf3DBase::CalculateSurfaceValues(MoleculeData *lData, Progress * lProgress) {
484 GaussHermiteData GHData;
485
486 if (!ContourHndl) return;
487 InitGaussHermite(&GHData); //setup the Gauss-Hermite roots and weights data
488
489 Frame * lFrame = lData->GetCurrentFramePtr();
490 BasisSet * Basis = lData->GetBasisSet();
491 AODensity * TotalAODensity = lFrame->GetAODensity(Basis, OrbSet);
492
493 if ((List && (ListAllocation != NumVertices))||!TotalAODensity)
494 FreeList();
495
496 if (!List) {
497 //Attempt to allocate memory for the psuedo 2D grid
498 AllocateList(NumVertices);
499 }
500 //If the memory allocation failed the MO calculation can not be done
501 if (List == NULL) return;
502 //If sufficient memory is available then setup pointers for fast local use
503 float * lGrid;
504 CPoint3D * Contour;
505 lGrid = List;
506 Contour = ContourHndl;
507
508 #ifdef __wxBuild__
509 if (wxThread::GetCPUCount() > 1) {
510 int myCPUCount = wxThread::GetCPUCount();
511 //Ok we have multiple CPUs so create a worker thread for each CPU and split up the grid calculation
512 //not sure if this is needed. It appears the default is 1?
513 // wxThread::SetConcurrency(myCPUCount);
514
515 TED3DSurfGridData * DataPtrs = new TED3DSurfGridData[myCPUCount];
516 int i;
517 TEDMEP3DThreadPtr * myThreads = new TEDMEP3DThreadPtr[myCPUCount];
518 long start=0;
519 for (i=0; i<myCPUCount; i++) {
520 DataPtrs[i].Surf = this;
521 DataPtrs[i].MolData = lData;
522 DataPtrs[i].TotalDensity = TotalAODensity;
523 DataPtrs[i].xStart = start;
524 start += NumVertices/myCPUCount;
525 if (i+1 == myCPUCount) start = NumVertices;
526 DataPtrs[i].xEnd = start;
527 DataPtrs[i].GHData = &GHData;
528 DataPtrs[i].grid = lGrid;
529 DataPtrs[i].Contour = Contour;
530 DataPtrs[i].PercentDone = 0;
531
532 //Create the actual thread
533 myThreads[i] = new TEDMEP3DThread(&(DataPtrs[i]));
534 myThreads[i]->Create(128*1024);
535 //and fire it off, note my class is always a joinable thread so we will wait on it eventually
536 myThreads[i]->Run();
537 }
538
539 long NumTasksRunning=myCPUCount;
540 long TotalPercent;
541 while (NumTasksRunning > 0) {
542 TotalPercent = 0;
543 for (i=0; i<myCPUCount; i++) {
544 TotalPercent += DataPtrs[i].PercentDone;
545 if (myThreads[i] != NULL) {
546 if (! myThreads[i]->IsAlive()) {
547 //The thread is terminating
548 myThreads[i]->Wait();
549 delete myThreads[i];
550 myThreads[i] = NULL;
551 NumTasksRunning --;
552 }
553 }
554 }
555 TotalPercent /= myCPUCount;
556 //Give up the CPU and check for cancels
557 if (!lProgress->UpdateProgress(TotalPercent)) { //User canceled so clean things up and abort
558 for (i=0; i<myCPUCount; i++) {
559 if (myThreads[i] != NULL) {
560 //The docs seem to indicate that Wait causes a thread to terminate,
561 //however this does not seem to be true. Wait appears to have the
562 //expected semantics that it blocks until the thread exits. Delete
563 //appears to be what is needed.
564 myThreads[i]->Delete();
565 delete myThreads[i];
566 myThreads[i] = NULL;
567 NumTasksRunning --;
568 }
569 }
570 FreeList();
571 } else {
572 wxMilliSleep(100);
573 }
574 }
575
576 delete [] DataPtrs;
577 delete[] myThreads;
578 } else
579 #endif
580 {
581 long PercentDone;
582 CalculateSurfaceValuesGrid(lData, lProgress, TotalAODensity, 0, NumVertices,
583 &GHData, lGrid, Contour, &PercentDone, false);
584 }
585 }
CalculateSurfaceValuesGrid(MoleculeData * lData,Progress * lProgress,AODensity * TotalAODensity,long start,long end,GaussHermiteData * GHData,float * lGrid,CPoint3D * Contour,long * PercentDone,bool MPTask)586 void Surf3DBase::CalculateSurfaceValuesGrid(MoleculeData * lData, Progress * lProgress,
587 AODensity * TotalAODensity, long start, long end,
588 GaussHermiteData * GHData, float * lGrid,
589 CPoint3D * Contour, long * PercentDone, bool MPTask) {
590 Frame * lFrame = lData->GetCurrentFramePtr();
591 BasisSet * Basis = lData->GetBasisSet();
592 #ifdef __wxBuild__
593 wxStopWatch timer;
594 long CheckTime = timer.Time();
595 #endif
596 //Store the Grid mins and incs at the beginning of the grid
597 //loop over the plotting grid in the x, y, and z directions
598 float XGridValue, YGridValue, ZGridValue, junk1, junk2;
599 long backpoint = 0;
600 for (long iPt=start; iPt<end; iPt++) {
601 *PercentDone = 100*(iPt-start)/(end-start);
602 if ((backpoint > 30)&&(lProgress != NULL)) {
603 //Give up the CPU and check for cancels
604 if (!lProgress->UpdateProgress(*PercentDone)) { //User canceled so clean things up and abort
605 FreeList();
606 return;
607 }
608 backpoint = 0;
609 }
610
611 XGridValue = Contour[iPt].x*kAng2BohrConversion; //Contour values must be converted to bohrs
612 YGridValue = Contour[iPt].y*kAng2BohrConversion;
613 ZGridValue = Contour[iPt].z*kAng2BohrConversion;
614
615 lGrid[iPt] = lFrame->CalculateMEP(XGridValue, YGridValue, ZGridValue,
616 Basis, TotalAODensity, GHData, &junk1, &junk2);
617 backpoint++;
618 #ifdef __wxBuild__
619 if (MPTask) {
620 long tempTime = timer.Time();
621 if ((tempTime - CheckTime) > 50) {
622 CheckTime = tempTime;
623 if ((wxThread::This())->TestDestroy()) {
624 //Getting here means the caller has requested a cancel so exit
625 return;
626 }
627 }
628 }
629 #endif
630 }
631 }
632
633 #define kMaxShellTypeCubed ((kMaxShellType+1)*(kMaxShellType+1)*(kMaxShellType+1))
634 #define kMaxAngFuncSquare (28*28) //Arrays dimensioned as the square of the max number of angular functions
635 // 28 is for max of I functions
636 //Calculates the MEP value at the specified x,y,z coordinate
CalculateMEP(float X_value,float Y_value,float Z_value,BasisSet * Basis,AODensity * TotalAODensity,GaussHermiteData * GHData,float * ElectronicMEP,float * NuclearMEP)637 float Frame::CalculateMEP(float X_value, float Y_value, float Z_value,
638 BasisSet *Basis, AODensity * TotalAODensity, GaussHermiteData * GHData,
639 float * ElectronicMEP, float * NuclearMEP) {
640 double ElectronicPotential=0.0, NuclearPotential=0.0; //results of the MEP calculation
641 long iatom, jatom, iprim, jprim, iroot, ijx[kMaxAngFuncSquare], ijy[kMaxAngFuncSquare], ijz[kMaxAngFuncSquare],
642 iang, jang, jmax, ij,nn, mm, NumJPrims, JPrimMax, TotalIShells, TotalJShells,
643 nx, ny, nz, in, jn, li, lj;
644 double xi, yi, zi, xj, yj, zj, r2, temp;
645 double xin[kMaxShellTypeCubed], yin[kMaxShellTypeCubed], zin[kMaxShellTypeCubed],
646 wint[kMaxAngFuncSquare], dij[kMaxAngFuncSquare];
647 Root RysRoots;
648 GaussHermiteIntegral GHInt;
649 AngularIndeces AngIndex;
650 std::vector<BasisShell> & Shells = Basis->Shells;
651 float * Density = TotalAODensity->GetDensityArray();
652 long * DensityIndex = TotalAODensity->GetDensityIndex();
653 short * DensityCheck = TotalAODensity->GetDensityCheck();
654 std::vector<long> & NuclearCharge = Basis->NuclearCharge;
655
656 double Pi212= 1.1283791670955;
657
658 TotalIShells = -1;
659 long CheckIndex = -1;
660 for (iatom=0; iatom < NumAtoms; iatom++) {
661 if (iatom > Basis->MapLength) continue;
662 xi = Atoms[iatom].Position.x*kAng2BohrConversion;
663 yi = Atoms[iatom].Position.y*kAng2BohrConversion;
664 zi = Atoms[iatom].Position.z*kAng2BohrConversion;
665 GHInt.xi = xi;
666 GHInt.yi = yi;
667 GHInt.zi = zi;
668 for (long ishell=Basis->BasisMap[2*iatom]; ishell<=Basis->BasisMap[2*iatom+1]; ishell++) {
669 TotalIShells ++;
670 long Lishell = Shells[ishell].GetShellType();
671 long NumIFuncs = Shells[ishell].GetNumFuncs(false);
672 long IAngMin = Shells[ishell].GetAngularStart(false);
673 long IAngMax = IAngMin + NumIFuncs;
674 long iprimMax = Shells[ishell].GetNumPrimitives();
675 //IAngMin is subtracted because it is added back later when the index is used
676 long LocIindex = DensityIndex[TotalIShells] - IAngMin;
677
678 TotalJShells = -1;
679 for (jatom=0; jatom <= iatom; jatom++) {
680 xj = Atoms[jatom].Position.x*kAng2BohrConversion;
681 yj = Atoms[jatom].Position.y*kAng2BohrConversion;
682 zj = Atoms[jatom].Position.z*kAng2BohrConversion;
683 GHInt.xj = xj;
684 GHInt.yj = yj;
685 GHInt.zj = zj;
686 r2 = (xi-xj)*(xi-xj) + (yi-yj)*(yi-yj) + (zi-zj)*(zi-zj);
687 long jShellMax = Basis->BasisMap[2*jatom+1];
688 if (iatom == jatom) jShellMax = ishell;
689 for (long jshell=Basis->BasisMap[2*jatom]; jshell<=jShellMax; jshell++) {
690 TotalJShells ++;
691 bool IequalsJ = (iatom == jatom) && (ishell == jshell);
692 long Ljshell = Shells[jshell].GetShellType();
693 long NumJFuncs = Shells[jshell].GetNumFuncs(false);
694 long JAngMin = Shells[jshell].GetAngularStart(false);
695 long JAngMax = JAngMin + NumJFuncs;
696 //JAngMin is subtracted because it is added back later when the index is used
697 long LocJindex = DensityIndex[TotalJShells]-JAngMin;
698 //Attempt to screen based on the density
699 CheckIndex ++;
700 if ( ! DensityCheck[CheckIndex]) continue;
701
702 NumJPrims = JPrimMax = Shells[jshell].GetNumPrimitives();
703
704 RysRoots.nRoots = (Lishell + Ljshell)/2 + 1;
705
706 ij=0;
707 jmax = JAngMax;
708 for (iang=IAngMin; iang<IAngMax; iang++) {
709 nx = (kMaxShellType+1)*AngIndex.x[iang];
710 ny = (kMaxShellType+1)*AngIndex.y[iang];
711 nz = (kMaxShellType+1)*AngIndex.z[iang];
712
713 if (IequalsJ) jmax = iang+1;
714 for (jang=JAngMin; jang<jmax; jang++) {
715 ijx[ij] = nx + AngIndex.x[jang];
716 ijy[ij] = ny + AngIndex.y[jang];
717 ijz[ij] = nz + AngIndex.z[jang];
718 ij++;
719 }
720 }
721 for (iprim=0; iprim<kMaxAngFuncSquare; iprim++) wint[iprim]=0.0;
722
723 for (iprim=0; iprim<iprimMax; iprim++) {
724 float Ai = Shells[ishell].Exponent[iprim];
725
726 if (IequalsJ) JPrimMax = iprim+1;
727 for (jprim=0; jprim<JPrimMax; jprim++) {
728 float Aj = Shells[jshell].Exponent[jprim];
729
730 double AiAj = Ai + Aj;
731 double OneDAiAj = 1.0/AiAj;
732 double Fi = Pi212 * OneDAiAj;
733
734 double AAx = (Ai*xi + Aj*xj);
735 double AAy = (Ai*yi + Aj*yj);
736 double AAz = (Ai*zi + Aj*zj);
737
738 double Ax = AAx*OneDAiAj;
739 double Ay = AAy*OneDAiAj;
740 double Az = AAz*OneDAiAj;
741
742 temp = Ai*Aj*r2*OneDAiAj;
743 //skip products with exponentials smaller than exp^-46 (GAMESS default)
744 if (temp > 46.0) continue;
745 double ExpFactor = Fi*exp(-temp);
746 bool Double = IequalsJ && (iprim != jprim);
747 nn=0;
748 jmax = JAngMax;
749 double DensityFactorI, DensityFactorJ;
750 //Compute the primitive coefficient factors
751 for (iang=IAngMin; iang<IAngMax; iang++) {
752 if (IAngMin == iang)
753 DensityFactorI = ExpFactor * Shells[ishell].NormCoef[iprim];
754 else if ((iang == 1) && (NumIFuncs == 4))
755 DensityFactorI = ExpFactor * Shells[ishell].NormCoef[iprim+iprimMax];
756 //Angular normalization factors are applied already in the AO density
757 if (IequalsJ) jmax = iang+1;
758 for (jang=JAngMin; jang<jmax; jang++) {
759 if (JAngMin == jang) {
760 DensityFactorJ = DensityFactorI * Shells[jshell].NormCoef[jprim];
761 if (Double) {
762 if (jang != 0 || iang == 0) DensityFactorJ += DensityFactorJ;
763 else DensityFactorJ += Shells[ishell].NormCoef[iprim]*
764 Shells[jshell].NormCoef[jprim+NumJPrims]*ExpFactor;
765 }
766 } else if ((jang == 1) && (NumJFuncs == 4)) {
767 DensityFactorJ = DensityFactorI * Shells[jshell].NormCoef[jprim+NumJPrims];
768 if (Double) DensityFactorJ += DensityFactorJ;
769 }
770
771 dij[nn] = DensityFactorJ;
772 nn++;
773 }
774 }
775
776 RysRoots.x = AiAj *((Ax-X_value)*(Ax-X_value) + (Ay-Y_value)*(Ay-Y_value) + (Az-Z_value)*(Az-Z_value));
777 if (RysRoots.nRoots <= 3) Root123(&RysRoots);
778 else if (RysRoots.nRoots == 4) Root4(&RysRoots);
779 else if (RysRoots.nRoots == 5) Root5(&RysRoots);
780 else Root6(&RysRoots);
781
782 mm=0;
783 for (iroot=0; iroot<RysRoots.nRoots; iroot++) {
784 double UU = AiAj*RysRoots.root[iroot];
785 double WW = RysRoots.weight[iroot];
786 double tt = 1.0/(AiAj+UU);
787
788 GHInt.T = sqrt(tt);
789 GHInt.x0 = (AAx + UU*X_value)*tt;
790 GHInt.y0 = (AAy + UU*Y_value)*tt;
791 GHInt.z0 = (AAz + UU*Z_value)*tt;
792
793 in = mm;
794 for (iang=0; iang<=Lishell; iang++) {
795 GHInt.NumI = iang;
796 for (jang=0; jang<=Ljshell; jang++) {
797 jn = in+jang;
798 GHInt.NumJ = jang;
799 DoGaussHermite(GHData, &GHInt);
800
801 xin[jn] = GHInt.xIntegral;
802 yin[jn] = GHInt.yIntegral;
803 zin[jn] = GHInt.zIntegral * WW;
804 }
805 in += (kMaxShellType+1);
806 }
807 mm += ((kMaxShellType+1)*(kMaxShellType+1));
808 } //Loop over orbital products
809 for (iang=0; iang<ij; iang++) {
810 nx = ijx[iang];
811 ny = ijy[iang];
812 nz = ijz[iang];
813 double sum = 0.0;
814 mm = 0;
815 for (iroot=0; iroot<RysRoots.nRoots; iroot++) {
816 sum += xin[nx+mm] * yin[ny+mm] * zin[nz+mm];
817 mm += ((kMaxShellType+1)*(kMaxShellType+1));
818 }
819 wint[iang] += sum*dij[iang];
820 }
821 }
822 } //End of primitive loops
823 long index = 0;
824 double dw;
825 jmax = JAngMax;
826 temp = 0.0;
827 for (iang=IAngMin; iang<IAngMax; iang++) {
828 li = LocIindex + iang;
829 in = li*(li+1)/2;
830 if (IequalsJ) jmax = iang+1;
831 for (jang=JAngMin; jang<jmax; jang++) {
832 lj = LocJindex + jang;
833 dw = Density[lj + in]*wint[index];
834 if (li != lj) dw += dw;
835 temp += dw;
836 index ++;
837 }
838 }
839 ElectronicPotential -= temp;
840 }
841 }
842 } //now for the nuclear contribution
843 xi -= X_value;
844 yi -= Y_value;
845 zi -= Z_value;
846 r2 = sqrt(xi*xi + yi*yi + zi*zi);
847 if (r2 > 0.001) NuclearPotential += NuclearCharge[iatom] / r2;
848 }
849 *ElectronicMEP = ElectronicPotential;
850 *NuclearMEP = NuclearPotential;
851 return (ElectronicPotential + NuclearPotential);
852 }
AODensity(void)853 AODensity::AODensity(void) {
854 DensityArray = NULL;
855 IndexArray = NULL;
856 DensityCheckArray=NULL;
857 Dimension = 0;
858 }
~AODensity(void)859 AODensity::~AODensity(void) {
860 if (DensityArray) delete[] DensityArray;
861 if (IndexArray) delete[] IndexArray;
862 if (DensityCheckArray) delete[] DensityCheckArray;
863 }
SetupMemory(long NumBasisFunctions)864 void AODensity::SetupMemory(long NumBasisFunctions) {
865 if (DensityArray) {delete[] DensityArray; DensityArray = NULL;}
866 if (IndexArray) {delete[] IndexArray; IndexArray = NULL;}
867 if (DensityCheckArray) {delete[] DensityCheckArray; DensityCheckArray=NULL;}
868 if (NumBasisFunctions > 0) {
869 DensityArray = new float[(NumBasisFunctions*(NumBasisFunctions+1))/2];
870 IndexArray = new long[NumBasisFunctions];
871 DensityCheckArray = new short[(NumBasisFunctions*(NumBasisFunctions+1))/2];
872 if (!DensityArray || !IndexArray || !DensityCheckArray) throw MemoryError();
873 Dimension = NumBasisFunctions;
874 }
875 }
GetAODensity(BasisSet * lBasis,long NumAtoms)876 AODensity * OrbitalRec::GetAODensity(BasisSet * lBasis, long NumAtoms) {
877 if (!TotalAODensity) {
878 long NumBasisFunctions = lBasis->GetNumBasisFuncs(false);
879 TotalAODensity = new AODensity;
880 if (TotalAODensity) {
881 TotalAODensity->SetupMemory(NumBasisFunctions);
882 GetAODensityMatrix(TotalAODensity->GetDensityArray(), NumOccupiedAlphaOrbs,
883 NumOccupiedBetaOrbs, NumBasisFunctions);
884 lBasis->GetShellIndexArray(TotalAODensity->GetDensityIndex());
885
886 std::vector<BasisShell> & Shells = lBasis->Shells;
887 float * Density = TotalAODensity->GetDensityArray();
888 long * DensityIndex = TotalAODensity->GetDensityIndex();
889 short * DensityCheck = TotalAODensity->GetDensityCheck();
890 long iatom, jatom, TotalJShells, iang, jang, li, in, lj, jmax;
891 double DensityFactorI, DensityFactorJ;
892 long CheckIndex=0;
893 long TotalIShells = -1;
894 double sqrt3 = sqrt(3.0);
895 double onedsqrt3 = 1/sqrt3;
896 double sqrt5 = sqrt(5.0);
897 double sqrt7 = sqrt(7.0);
898 double sqrt11 = sqrt(11.0);
899 for (iatom=0; iatom < NumAtoms; iatom++) {
900 if (iatom > lBasis->MapLength) continue;
901 for (long ishell=lBasis->BasisMap[2*iatom]; ishell<=lBasis->BasisMap[2*iatom+1]; ishell++) {
902 TotalIShells ++;
903 // long Lishell = Shells[ishell].GetShellType();
904 long NumIFuncs = Shells[ishell].GetNumFuncs(false);
905 long IAngMin = Shells[ishell].GetAngularStart(false);
906 long IAngMax = IAngMin + NumIFuncs;
907 //IAngMin is subtracted because it is added back later when the index is used
908 long LocIindex = DensityIndex[TotalIShells] - IAngMin;
909
910 TotalJShells = -1;
911 for (jatom=0; jatom <= iatom; jatom++) {
912 long jShellMax = lBasis->BasisMap[2*jatom+1];
913 if (iatom == jatom) jShellMax = ishell;
914 for (long jshell=lBasis->BasisMap[2*jatom]; jshell<=jShellMax; jshell++) {
915 TotalJShells ++;
916 bool IequalsJ = (iatom == jatom) && (ishell == jshell);
917 // long Ljshell = Shells[jshell].GetShellType();
918 long NumJFuncs = Shells[jshell].GetNumFuncs(false);
919 long JAngMin = Shells[jshell].GetAngularStart(false);
920 long JAngMax = JAngMin + NumJFuncs;
921 //JAngMin is subtracted because it is added back later when the index is used
922 long LocJindex = DensityIndex[TotalJShells]-JAngMin;
923 //Attempt to screen based on the density
924 short NonZeroShells=0;
925 jmax = JAngMax;
926 for (iang=IAngMin; iang<IAngMax; iang++) {
927 if (IAngMin == iang)
928 DensityFactorI = 1.0;
929 if (iang >= (IAngMin+3)) {
930 switch (iang) { //properly normalize d, f, and g's
931 case 7:
932 case 19:
933 case 32:
934 case 50: //H shell func 15 - x3yz
935 case 65: //I shell func 10 - x4y2
936 case 71: //I shell func 16 - x4yz
937 DensityFactorI *= sqrt3;
938 break;
939 case 13:
940 case 77: //I shell func 22 - x3y2z
941 DensityFactorI *= sqrt5;
942 break;
943 case 23:
944 DensityFactorI *= sqrt7;
945 break;
946 case 29:
947 case 53: //H shell func 18 - x2y2z
948 case 83: //I shell func 27 - x2y2z2
949 DensityFactorI *= sqrt5*onedsqrt3;
950 break;
951 case 38:
952 DensityFactorI *= 3.0; //sqrt(9) - H shell func 4, x4y
953 break;
954 case 44:
955 DensityFactorI *= sqrt7*onedsqrt3; //H shell 10, x3y2
956 break;
957 case 59: //I shell func 4 - x5y
958 DensityFactorI *= sqrt11;
959 break;
960 case 74: //I shell func 18 - x3y3
961 DensityFactorI *= sqrt7/(sqrt3*sqrt5);
962 break;
963 }
964 }
965 li = LocIindex + iang;
966 in = li*(li+1)/2;
967 if (IequalsJ) jmax = iang+1;
968 for (jang=JAngMin; jang<jmax; jang++) {
969 lj = LocJindex + jang;
970 if (JAngMin == jang)
971 DensityFactorJ = DensityFactorI;
972 if (jang >= (JAngMin+3)) {
973 switch (jang) { //properly normalize d, f, and g's
974 case 7:
975 case 19:
976 case 32:
977 case 50: //H shell func 15 - x3yz
978 case 65: //I shell func 10 - x4y2
979 case 71: //I shell func 16 - x4yz
980 DensityFactorJ *= sqrt3;
981 break;
982 case 13:
983 case 77: //I shell func 22 - x3y2z
984 DensityFactorJ *= sqrt5;
985 break;
986 case 23:
987 DensityFactorJ *= sqrt7;
988 break;
989 case 29:
990 case 53: //H shell func 18 - x2y2z
991 case 83: //I shell func 27 - x2y2z2
992 DensityFactorJ *= sqrt5*onedsqrt3;
993 break;
994 case 38:
995 DensityFactorJ *= 3.0; //sqrt(9) - H shell func 4, x4y
996 break;
997 case 44:
998 DensityFactorJ *= sqrt7*onedsqrt3; //H shell 10, x3y2
999 break;
1000 case 59: //I shell func 4 - x5y
1001 DensityFactorJ *= sqrt11;
1002 break;
1003 case 74: //I shell func 18 - x3y3
1004 DensityFactorJ *= sqrt7/(sqrt3*sqrt5);
1005 break;
1006 }
1007 }
1008 Density[lj+in] *= DensityFactorJ;
1009 if (fabs(Density[lj+in]) > 0.0001) NonZeroShells++;
1010 }
1011 }
1012 DensityCheck[CheckIndex] = NonZeroShells;
1013 CheckIndex++;
1014 }
1015 }
1016 }
1017 }
1018
1019 }
1020 if (!TotalAODensity) throw MemoryError();
1021 }
1022 return TotalAODensity;
1023 }
1024
GetAODensityMatrix(float * AODensityArray,long NumOccAlpha,long NumOccBeta,long NumBasisFuncs) const1025 void OrbitalRec::GetAODensityMatrix(float * AODensityArray, long NumOccAlpha, long NumOccBeta,
1026 long NumBasisFuncs) const {
1027 long n, i, j, m;
1028
1029 if (!AODensityArray) throw MemoryError(); //Memory must be allocated by the caller
1030 float *OccupancyA=NULL, *OccupancyB=NULL;
1031 bool OccMemIsLocal = false;
1032 OccupancyA = OrbOccupation;
1033 OccupancyB = OrbOccupationB;
1034 if (!OrbOccupation) {
1035 if (BaseWavefunction == RHF) {
1036 OccupancyA = new float[NumOccAlpha];
1037 if (OccupancyA) for (i=0; i<NumOccAlpha; i++) OccupancyA[i]=2.0;
1038 OccMemIsLocal = true;
1039 } else if (BaseWavefunction == ROHF) {
1040 OccupancyA = new float[NumOccAlpha];
1041 if (OccupancyA) for (i=0; i<NumOccAlpha; i++) {
1042 if (i<NumOccBeta) OccupancyA[i]=2.0;
1043 else OccupancyA[i]=1.0;
1044 }
1045 OccMemIsLocal = true;
1046 } else if ((BaseWavefunction == UHF)&&(OrbitalType != NaturalOrbital)) {
1047 OccupancyA = new float[NumOccAlpha];
1048 if (OccupancyA) for (i=0; i<NumOccAlpha; i++) OccupancyA[i]=1.0;
1049 OccupancyB = new float[NumOccBeta];
1050 if (OccupancyB) for (i=0; i<NumOccBeta; i++) OccupancyB[i]=1.0;
1051 OccMemIsLocal = true;
1052 }
1053 }
1054 if (!OccupancyA) throw MemoryError();
1055 //clear the AO density array
1056 n=0;
1057 for (i=0; i<NumBasisFuncs; i++) {
1058 for (j=0; j<=i; j++) {
1059 AODensityArray[n]=0.0;
1060 n++;
1061 }
1062 }
1063
1064 for (m=0; m<NumOccAlpha; m++) {
1065 n=0;
1066 for (i=0; i<NumBasisFuncs; i++) {
1067 for (j=0; j<=i; j++) {
1068 AODensityArray[n] += Vectors[i+m*NumBasisFuncs]*OccupancyA[m]*Vectors[j+m*NumBasisFuncs];
1069 n++;
1070 }
1071 }
1072 }
1073 if (((BaseWavefunction == UHF)&&(OrbitalType != NaturalOrbital))||OccupancyB) {
1074 if (!OccupancyB) throw MemoryError();
1075 for (m=0; m<NumOccBeta; m++) {
1076 n=0;
1077 for (i=0; i<NumBasisFuncs; i++) {
1078 for (j=0; j<=i; j++) {
1079 AODensityArray[n] += VectorsB[i+m*NumBasisFuncs]*OccupancyB[m]*VectorsB[j+m*NumBasisFuncs];
1080 n++;
1081 }
1082 }
1083 }
1084 }
1085 //delete memory if it was allocated above
1086 if (OccMemIsLocal)
1087 delete[] OccupancyA;
1088 if (OccMemIsLocal && OccupancyB) delete[] OccupancyB;
1089 }
1090
1091