1 /******************************************************************************
2
3 This source file is part of the Avogadro project.
4
5 Copyright (C) 2010 Eric C. Brown
6
7 This source code is released under the New BSD License, (the "License").
8
9 Unless required by applicable law or agreed to in writing, software
10 distributed under the License is distributed on an "AS IS" BASIS,
11 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 See the License for the specific language governing permissions and
13 limitations under the License.
14
15 ******************************************************************************/
16
17 #include <cmath>
18
19 #include "qtaimwavefunctionevaluator.h"
20
21 namespace Avogadro {
22 namespace QtPlugins {
23
QTAIMWavefunctionEvaluator(QTAIMWavefunction & wfn)24 QTAIMWavefunctionEvaluator::QTAIMWavefunctionEvaluator(QTAIMWavefunction& wfn)
25 {
26
27 m_nmo = wfn.numberOfMolecularOrbitals();
28 m_nprim = wfn.numberOfGaussianPrimitives();
29 m_nnuc = wfn.numberOfNuclei();
30
31 m_nucxcoord =
32 Map<const Matrix<qreal, Dynamic, 1>>(wfn.xNuclearCoordinates(), m_nnuc);
33 m_nucycoord =
34 Map<const Matrix<qreal, Dynamic, 1>>(wfn.yNuclearCoordinates(), m_nnuc);
35 m_nuczcoord =
36 Map<const Matrix<qreal, Dynamic, 1>>(wfn.zNuclearCoordinates(), m_nnuc);
37 m_nucz = Map<const Matrix<qint64, Dynamic, 1>>(wfn.nuclearCharges(), m_nnuc);
38 m_X0 = Map<const Matrix<qreal, Dynamic, 1>>(
39 wfn.xGaussianPrimitiveCenterCoordinates(), m_nprim, 1);
40 m_Y0 = Map<const Matrix<qreal, Dynamic, 1>>(
41 wfn.yGaussianPrimitiveCenterCoordinates(), m_nprim, 1);
42 m_Z0 = Map<const Matrix<qreal, Dynamic, 1>>(
43 wfn.zGaussianPrimitiveCenterCoordinates(), m_nprim, 1);
44 m_xamom = Map<const Matrix<qint64, Dynamic, 1>>(
45 wfn.xGaussianPrimitiveAngularMomenta(), m_nprim, 1);
46 m_yamom = Map<const Matrix<qint64, Dynamic, 1>>(
47 wfn.yGaussianPrimitiveAngularMomenta(), m_nprim, 1);
48 m_zamom = Map<const Matrix<qint64, Dynamic, 1>>(
49 wfn.zGaussianPrimitiveAngularMomenta(), m_nprim, 1);
50 m_alpha = Map<const Matrix<qreal, Dynamic, 1>>(
51 wfn.gaussianPrimitiveExponentCoefficients(), m_nprim, 1);
52 // TODO Implement screening for unoccupied molecular orbitals.
53 m_occno = Map<const Matrix<qreal, Dynamic, 1>>(
54 wfn.molecularOrbitalOccupationNumbers(), m_nmo, 1);
55 m_orbe = Map<const Matrix<qreal, Dynamic, 1>>(
56 wfn.molecularOrbitalEigenvalues(), m_nmo, 1);
57 m_coef = Map<const Matrix<qreal, Dynamic, Dynamic, RowMajor>>(
58 wfn.molecularOrbitalCoefficients(), m_nmo, m_nprim);
59 m_totalEnergy = wfn.totalEnergy();
60 m_virialRatio = wfn.virialRatio();
61
62 m_cutoff = log(1.e-15);
63
64 m_cdg000.resize(m_nmo);
65 m_cdg100.resize(m_nmo);
66 m_cdg010.resize(m_nmo);
67 m_cdg001.resize(m_nmo);
68 m_cdg200.resize(m_nmo);
69 m_cdg110.resize(m_nmo);
70 m_cdg101.resize(m_nmo);
71 m_cdg020.resize(m_nmo);
72 m_cdg011.resize(m_nmo);
73 m_cdg002.resize(m_nmo);
74 m_cdg300.resize(m_nmo);
75 m_cdg120.resize(m_nmo);
76 m_cdg102.resize(m_nmo);
77 m_cdg210.resize(m_nmo);
78 m_cdg030.resize(m_nmo);
79 m_cdg012.resize(m_nmo);
80 m_cdg201.resize(m_nmo);
81 m_cdg021.resize(m_nmo);
82 m_cdg003.resize(m_nmo);
83 m_cdg111.resize(m_nmo);
84 m_cdg400.resize(m_nmo);
85 m_cdg220.resize(m_nmo);
86 m_cdg202.resize(m_nmo);
87 m_cdg310.resize(m_nmo);
88 m_cdg130.resize(m_nmo);
89 m_cdg112.resize(m_nmo);
90 m_cdg301.resize(m_nmo);
91 m_cdg121.resize(m_nmo);
92 m_cdg103.resize(m_nmo);
93 m_cdg040.resize(m_nmo);
94 m_cdg022.resize(m_nmo);
95 m_cdg211.resize(m_nmo);
96 m_cdg031.resize(m_nmo);
97 m_cdg013.resize(m_nmo);
98 m_cdg004.resize(m_nmo);
99 }
100
molecularOrbital(const qint64 mo,const Matrix<qreal,3,1> xyz)101 qreal QTAIMWavefunctionEvaluator::molecularOrbital(
102 const qint64 mo, const Matrix<qreal, 3, 1> xyz)
103 {
104
105 qreal value = 0.0;
106
107 for (qint64 p = 0; p < m_nprim; ++p) {
108 qreal xx0 = xyz(0) - m_X0(p);
109 qreal yy0 = xyz(1) - m_Y0(p);
110 qreal zz0 = xyz(2) - m_Z0(p);
111
112 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
113
114 if (b0arg > m_cutoff) {
115 qreal ax0 = ipow(xx0, m_xamom(p));
116 qreal ay0 = ipow(yy0, m_yamom(p));
117 qreal az0 = ipow(zz0, m_zamom(p));
118
119 qreal b0 = exp(b0arg);
120
121 qreal dg000 = ax0 * ay0 * az0 * b0;
122
123 value += m_coef(mo, p) * dg000;
124 }
125 }
126
127 return value;
128 }
129
electronDensity(const Matrix<qreal,3,1> xyz)130 qreal QTAIMWavefunctionEvaluator::electronDensity(const Matrix<qreal, 3, 1> xyz)
131 {
132
133 qreal value;
134
135 m_cdg000.setZero();
136 for (qint64 p = 0; p < m_nprim; ++p) {
137 qreal xx0 = xyz(0) - m_X0(p);
138 qreal yy0 = xyz(1) - m_Y0(p);
139 qreal zz0 = xyz(2) - m_Z0(p);
140
141 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
142
143 if (b0arg > m_cutoff) {
144 qreal ax0 = ipow(xx0, m_xamom(p));
145 qreal ay0 = ipow(yy0, m_yamom(p));
146 qreal az0 = ipow(zz0, m_zamom(p));
147
148 qreal b0 = exp(b0arg);
149
150 qreal dg000 = ax0 * ay0 * az0 * b0;
151
152 for (qint64 m = 0; m < m_nmo; ++m) {
153 m_cdg000(m) += m_coef(m, p) * dg000;
154 }
155 }
156 }
157
158 value = 0.0;
159 for (qint64 m = 0; m < m_nmo; ++m) {
160 value += m_occno(m) * ipow(m_cdg000(m), 2);
161 }
162
163 return value;
164 }
165
gradientOfElectronDensity(Matrix<qreal,3,1> xyz)166 const Matrix<qreal, 3, 1> QTAIMWavefunctionEvaluator::gradientOfElectronDensity(
167 Matrix<qreal, 3, 1> xyz)
168 {
169
170 Matrix<qreal, 3, 1> value;
171
172 const qreal zero = 0.0;
173 const qreal one = 1.0;
174
175 m_cdg000.setZero();
176 m_cdg100.setZero();
177 m_cdg010.setZero();
178 m_cdg001.setZero();
179 for (qint64 p = 0; p < m_nprim; ++p) {
180 qreal xx0 = xyz(0) - m_X0(p);
181 qreal yy0 = xyz(1) - m_Y0(p);
182 qreal zz0 = xyz(2) - m_Z0(p);
183
184 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
185
186 if (b0arg > m_cutoff) {
187 qint64 aax0 = 1;
188 qint64 aay0 = 1;
189 qint64 aaz0 = 1;
190 qint64 aax1 = m_xamom(p);
191 qint64 aay1 = m_yamom(p);
192 qint64 aaz1 = m_zamom(p);
193
194 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
195 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
196 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
197
198 qreal ax1;
199 qreal ay1;
200 qreal az1;
201 if (m_xamom(p) < 1) {
202 ax1 = zero;
203 } else if (m_xamom(p) == 1) {
204 ax1 = one;
205 } else {
206 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
207 }
208
209 if (m_yamom(p) < 1) {
210 ay1 = zero;
211 } else if (m_yamom(p) == 1) {
212 ay1 = one;
213 } else {
214 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
215 }
216
217 if (m_zamom(p) < 1) {
218 az1 = zero;
219 } else if (m_zamom(p) == 1) {
220 az1 = one;
221 } else {
222 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
223 }
224
225 qreal b0 = exp(b0arg);
226
227 qreal bx1 = -2 * m_alpha(p) * xx0;
228 qreal by1 = -2 * m_alpha(p) * yy0;
229 qreal bz1 = -2 * m_alpha(p) * zz0;
230
231 qreal dg000 = ax0 * ay0 * az0 * b0;
232 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
233 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
234 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
235
236 for (qint64 m = 0; m < m_nmo; ++m) {
237 m_cdg000(m) += m_coef(m, p) * dg000;
238 m_cdg100(m) += m_coef(m, p) * dg100;
239 m_cdg010(m) += m_coef(m, p) * dg010;
240 m_cdg001(m) += m_coef(m, p) * dg001;
241 }
242 }
243 }
244
245 value.setZero();
246 for (qint64 m = 0; m < m_nmo; ++m) {
247 value(0) += m_occno(m) * m_cdg100(m) * m_cdg000(m);
248 value(1) += m_occno(m) * m_cdg010(m) * m_cdg000(m);
249 value(2) += m_occno(m) * m_cdg001(m) * m_cdg000(m);
250 }
251
252 return value;
253 }
254
hessianOfElectronDensity(const Matrix<qreal,3,1> xyz)255 const Matrix<qreal, 3, 3> QTAIMWavefunctionEvaluator::hessianOfElectronDensity(
256 const Matrix<qreal, 3, 1> xyz)
257 {
258
259 Matrix<qreal, 3, 3> value;
260
261 const qreal zero = 0.0;
262 const qreal one = 1.0;
263
264 m_cdg000.setZero();
265 m_cdg100.setZero();
266 m_cdg010.setZero();
267 m_cdg001.setZero();
268 m_cdg200.setZero();
269 m_cdg020.setZero();
270 m_cdg002.setZero();
271 m_cdg110.setZero();
272 m_cdg101.setZero();
273 m_cdg011.setZero();
274 for (qint64 p = 0; p < m_nprim; ++p) {
275 qreal xx0 = xyz(0) - m_X0(p);
276 qreal yy0 = xyz(1) - m_Y0(p);
277 qreal zz0 = xyz(2) - m_Z0(p);
278
279 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
280
281 if (b0arg > m_cutoff) {
282 qint64 aax0 = 1;
283 qint64 aay0 = 1;
284 qint64 aaz0 = 1;
285 qint64 aax1 = m_xamom(p);
286 qint64 aay1 = m_yamom(p);
287 qint64 aaz1 = m_zamom(p);
288 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
289 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
290 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
291
292 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
293 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
294 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
295
296 qreal ax1;
297 qreal ay1;
298 qreal az1;
299 if (m_xamom(p) < 1) {
300 ax1 = zero;
301 } else if (m_xamom(p) == 1) {
302 ax1 = one;
303 } else {
304 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
305 }
306
307 if (m_yamom(p) < 1) {
308 ay1 = zero;
309 } else if (m_yamom(p) == 1) {
310 ay1 = one;
311 } else {
312 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
313 }
314
315 if (m_zamom(p) < 1) {
316 az1 = zero;
317 } else if (m_zamom(p) == 1) {
318 az1 = one;
319 } else {
320 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
321 }
322
323 qreal ax2;
324 qreal ay2;
325 qreal az2;
326 if (m_xamom(p) < 2) {
327 ax2 = zero;
328 } else if (m_xamom(p) == 2) {
329 ax2 = one;
330 } else {
331 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
332 }
333
334 if (m_yamom(p) < 2) {
335 ay2 = zero;
336 } else if (m_yamom(p) == 2) {
337 ay2 = one;
338 } else {
339 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
340 }
341
342 if (m_zamom(p) < 2) {
343 az2 = zero;
344 } else if (m_zamom(p) == 2) {
345 az2 = one;
346 } else {
347 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
348 }
349
350 qreal b0 = exp(b0arg);
351
352 qreal bx1 = -2 * m_alpha(p) * xx0;
353 qreal by1 = -2 * m_alpha(p) * yy0;
354 qreal bz1 = -2 * m_alpha(p) * zz0;
355 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
356 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
357 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
358
359 qreal dg000 = ax0 * ay0 * az0 * b0;
360 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
361 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
362 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
363 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
364 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
365 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
366 qreal dg110 = az0 * b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1);
367 qreal dg101 = ay0 * b0 * (ax1 + ax0 * bx1) * (az1 + az0 * bz1);
368 qreal dg011 = ax0 * b0 * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
369
370 for (qint64 m = 0; m < m_nmo; ++m) {
371 m_cdg000(m) += m_coef(m, p) * dg000;
372 m_cdg100(m) += m_coef(m, p) * dg100;
373 m_cdg010(m) += m_coef(m, p) * dg010;
374 m_cdg001(m) += m_coef(m, p) * dg001;
375 m_cdg200(m) += m_coef(m, p) * dg200;
376 m_cdg020(m) += m_coef(m, p) * dg020;
377 m_cdg002(m) += m_coef(m, p) * dg002;
378 m_cdg110(m) += m_coef(m, p) * dg110;
379 m_cdg101(m) += m_coef(m, p) * dg101;
380 m_cdg011(m) += m_coef(m, p) * dg011;
381 }
382 }
383 }
384
385 value.setZero();
386 for (qint64 m = 0; m < m_nmo; ++m) {
387 value(0, 0) +=
388 2 * m_occno(m) * (ipow(m_cdg100(m), 2) + m_cdg000(m) * m_cdg200(m));
389 value(1, 1) +=
390 2 * m_occno(m) * (ipow(m_cdg010(m), 2) + m_cdg000(m) * m_cdg020(m));
391 value(2, 2) +=
392 2 * m_occno(m) * (ipow(m_cdg001(m), 2) + m_cdg000(m) * m_cdg002(m));
393 value(0, 1) +=
394 2 * m_occno(m) * (m_cdg100(m) * m_cdg010(m) + m_cdg000(m) * m_cdg110(m));
395 value(0, 2) +=
396 2 * m_occno(m) * (m_cdg100(m) * m_cdg001(m) + m_cdg000(m) * m_cdg101(m));
397 value(1, 2) +=
398 2 * m_occno(m) * (m_cdg010(m) * m_cdg001(m) + m_cdg000(m) * m_cdg011(m));
399 }
400 value(1, 0) = value(0, 1);
401 value(2, 0) = value(0, 2);
402 value(2, 1) = value(1, 2);
403
404 return value;
405 }
406
407 const Matrix<qreal, 3, 4>
gradientAndHessianOfElectronDensity(const Matrix<qreal,3,1> xyz)408 QTAIMWavefunctionEvaluator::gradientAndHessianOfElectronDensity(
409 const Matrix<qreal, 3, 1> xyz)
410 {
411
412 Matrix<qreal, 3, 1> gValue;
413 Matrix<qreal, 3, 3> hValue;
414 Matrix<qreal, 3, 4> value;
415
416 const qreal zero = 0.0;
417 const qreal one = 1.0;
418
419 m_cdg000.setZero();
420 m_cdg100.setZero();
421 m_cdg010.setZero();
422 m_cdg001.setZero();
423 m_cdg200.setZero();
424 m_cdg020.setZero();
425 m_cdg002.setZero();
426 m_cdg110.setZero();
427 m_cdg101.setZero();
428 m_cdg011.setZero();
429 for (qint64 p = 0; p < m_nprim; ++p) {
430 qreal xx0 = xyz(0) - m_X0(p);
431 qreal yy0 = xyz(1) - m_Y0(p);
432 qreal zz0 = xyz(2) - m_Z0(p);
433
434 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
435
436 if (b0arg > m_cutoff) {
437 qint64 aax0 = 1;
438 qint64 aay0 = 1;
439 qint64 aaz0 = 1;
440 qint64 aax1 = m_xamom(p);
441 qint64 aay1 = m_yamom(p);
442 qint64 aaz1 = m_zamom(p);
443 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
444 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
445 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
446
447 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
448 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
449 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
450
451 qreal ax1;
452 qreal ay1;
453 qreal az1;
454 if (m_xamom(p) < 1) {
455 ax1 = zero;
456 } else if (m_xamom(p) == 1) {
457 ax1 = one;
458 } else {
459 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
460 }
461
462 if (m_yamom(p) < 1) {
463 ay1 = zero;
464 } else if (m_yamom(p) == 1) {
465 ay1 = one;
466 } else {
467 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
468 }
469
470 if (m_zamom(p) < 1) {
471 az1 = zero;
472 } else if (m_zamom(p) == 1) {
473 az1 = one;
474 } else {
475 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
476 }
477
478 qreal ax2;
479 qreal ay2;
480 qreal az2;
481 if (m_xamom(p) < 2) {
482 ax2 = zero;
483 } else if (m_xamom(p) == 2) {
484 ax2 = one;
485 } else {
486 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
487 }
488
489 if (m_yamom(p) < 2) {
490 ay2 = zero;
491 } else if (m_yamom(p) == 2) {
492 ay2 = one;
493 } else {
494 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
495 }
496
497 if (m_zamom(p) < 2) {
498 az2 = zero;
499 } else if (m_zamom(p) == 2) {
500 az2 = one;
501 } else {
502 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
503 }
504
505 qreal b0 = exp(b0arg);
506
507 qreal bx1 = -2 * m_alpha(p) * xx0;
508 qreal by1 = -2 * m_alpha(p) * yy0;
509 qreal bz1 = -2 * m_alpha(p) * zz0;
510 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
511 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
512 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
513
514 qreal dg000 = ax0 * ay0 * az0 * b0;
515 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
516 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
517 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
518 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
519 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
520 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
521 qreal dg110 = az0 * b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1);
522 qreal dg101 = ay0 * b0 * (ax1 + ax0 * bx1) * (az1 + az0 * bz1);
523 qreal dg011 = ax0 * b0 * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
524
525 for (qint64 m = 0; m < m_nmo; ++m) {
526 m_cdg000(m) += m_coef(m, p) * dg000;
527 m_cdg100(m) += m_coef(m, p) * dg100;
528 m_cdg010(m) += m_coef(m, p) * dg010;
529 m_cdg001(m) += m_coef(m, p) * dg001;
530 m_cdg200(m) += m_coef(m, p) * dg200;
531 m_cdg020(m) += m_coef(m, p) * dg020;
532 m_cdg002(m) += m_coef(m, p) * dg002;
533 m_cdg110(m) += m_coef(m, p) * dg110;
534 m_cdg101(m) += m_coef(m, p) * dg101;
535 m_cdg011(m) += m_coef(m, p) * dg011;
536 }
537 }
538 }
539
540 gValue.setZero();
541 for (qint64 m = 0; m < m_nmo; ++m) {
542 gValue(0) += m_occno(m) * m_cdg100(m) * m_cdg000(m);
543 gValue(1) += m_occno(m) * m_cdg010(m) * m_cdg000(m);
544 gValue(2) += m_occno(m) * m_cdg001(m) * m_cdg000(m);
545 }
546
547 hValue.setZero();
548 for (qint64 m = 0; m < m_nmo; ++m) {
549 hValue(0, 0) +=
550 2 * m_occno(m) * (ipow(m_cdg100(m), 2) + m_cdg000(m) * m_cdg200(m));
551 hValue(1, 1) +=
552 2 * m_occno(m) * (ipow(m_cdg010(m), 2) + m_cdg000(m) * m_cdg020(m));
553 hValue(2, 2) +=
554 2 * m_occno(m) * (ipow(m_cdg001(m), 2) + m_cdg000(m) * m_cdg002(m));
555 hValue(0, 1) +=
556 2 * m_occno(m) * (m_cdg100(m) * m_cdg010(m) + m_cdg000(m) * m_cdg110(m));
557 hValue(0, 2) +=
558 2 * m_occno(m) * (m_cdg100(m) * m_cdg001(m) + m_cdg000(m) * m_cdg101(m));
559 hValue(1, 2) +=
560 2 * m_occno(m) * (m_cdg010(m) * m_cdg001(m) + m_cdg000(m) * m_cdg011(m));
561 }
562 hValue(1, 0) = hValue(0, 1);
563 hValue(2, 0) = hValue(0, 2);
564 hValue(2, 1) = hValue(1, 2);
565
566 value(0, 0) = gValue(0);
567 value(1, 0) = gValue(1);
568 value(2, 0) = gValue(2);
569 value(0, 1) = hValue(0, 0);
570 value(1, 1) = hValue(1, 0);
571 value(2, 1) = hValue(2, 0);
572 value(0, 2) = hValue(0, 1);
573 value(1, 2) = hValue(1, 1);
574 value(2, 2) = hValue(2, 1);
575 value(0, 3) = hValue(0, 2);
576 value(1, 3) = hValue(1, 2);
577 value(2, 3) = hValue(2, 2);
578
579 return value;
580 }
581
laplacianOfElectronDensity(const Matrix<qreal,3,1> xyz)582 qreal QTAIMWavefunctionEvaluator::laplacianOfElectronDensity(
583 const Matrix<qreal, 3, 1> xyz)
584 {
585
586 qreal value;
587
588 const qreal zero = 0.0;
589 const qreal one = 1.0;
590
591 m_cdg000.setZero();
592 m_cdg100.setZero();
593 m_cdg010.setZero();
594 m_cdg001.setZero();
595 m_cdg200.setZero();
596 m_cdg020.setZero();
597 m_cdg002.setZero();
598 for (qint64 p = 0; p < m_nprim; ++p) {
599 qreal xx0 = xyz(0) - m_X0(p);
600 qreal yy0 = xyz(1) - m_Y0(p);
601 qreal zz0 = xyz(2) - m_Z0(p);
602
603 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
604
605 if (b0arg > m_cutoff) {
606 qint64 aax0 = 1;
607 qint64 aay0 = 1;
608 qint64 aaz0 = 1;
609 qint64 aax1 = m_xamom(p);
610 qint64 aay1 = m_yamom(p);
611 qint64 aaz1 = m_zamom(p);
612 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
613 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
614 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
615
616 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
617 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
618 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
619
620 qreal ax1;
621 qreal ay1;
622 qreal az1;
623 if (m_xamom(p) < 1) {
624 ax1 = zero;
625 } else if (m_xamom(p) == 1) {
626 ax1 = one;
627 } else {
628 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
629 }
630
631 if (m_yamom(p) < 1) {
632 ay1 = zero;
633 } else if (m_yamom(p) == 1) {
634 ay1 = one;
635 } else {
636 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
637 }
638
639 if (m_zamom(p) < 1) {
640 az1 = zero;
641 } else if (m_zamom(p) == 1) {
642 az1 = one;
643 } else {
644 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
645 }
646
647 qreal ax2;
648 qreal ay2;
649 qreal az2;
650 if (m_xamom(p) < 2) {
651 ax2 = zero;
652 } else if (m_xamom(p) == 2) {
653 ax2 = one;
654 } else {
655 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
656 }
657
658 if (m_yamom(p) < 2) {
659 ay2 = zero;
660 } else if (m_yamom(p) == 2) {
661 ay2 = one;
662 } else {
663 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
664 }
665
666 if (m_zamom(p) < 2) {
667 az2 = zero;
668 } else if (m_zamom(p) == 2) {
669 az2 = one;
670 } else {
671 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
672 }
673
674 qreal b0 = exp(b0arg);
675
676 qreal bx1 = -2 * m_alpha(p) * xx0;
677 qreal by1 = -2 * m_alpha(p) * yy0;
678 qreal bz1 = -2 * m_alpha(p) * zz0;
679 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
680 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
681 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
682
683 qreal dg000 = ax0 * ay0 * az0 * b0;
684 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
685 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
686 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
687 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
688 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
689 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
690
691 for (qint64 m = 0; m < m_nmo; ++m) {
692 m_cdg000(m) += m_coef(m, p) * dg000;
693 m_cdg100(m) += m_coef(m, p) * dg100;
694 m_cdg010(m) += m_coef(m, p) * dg010;
695 m_cdg001(m) += m_coef(m, p) * dg001;
696 m_cdg200(m) += m_coef(m, p) * dg200;
697 m_cdg020(m) += m_coef(m, p) * dg020;
698 m_cdg002(m) += m_coef(m, p) * dg002;
699 }
700 }
701 }
702
703 value = 0.0;
704 for (qint64 m = 0; m < m_nmo; ++m) {
705 value +=
706 2 * m_occno(m) * (ipow(m_cdg100(m), 2) + m_cdg000(m) * m_cdg200(m)) +
707 2 * m_occno(m) * (ipow(m_cdg010(m), 2) + m_cdg000(m) * m_cdg020(m)) +
708 2 * m_occno(m) * (ipow(m_cdg001(m), 2) + m_cdg000(m) * m_cdg002(m));
709 }
710
711 return value;
712 }
713
714 const Matrix<qreal, 3, 1>
gradientOfElectronDensityLaplacian(const Matrix<qreal,3,1> xyz)715 QTAIMWavefunctionEvaluator::gradientOfElectronDensityLaplacian(
716 const Matrix<qreal, 3, 1> xyz)
717 {
718
719 Matrix<qreal, 3, 1> value;
720
721 const qreal zero = 0.0;
722 const qreal one = 1.0;
723
724 m_cdg000.setZero();
725 m_cdg100.setZero();
726 m_cdg010.setZero();
727 m_cdg001.setZero();
728 m_cdg200.setZero();
729 m_cdg020.setZero();
730 m_cdg002.setZero();
731 m_cdg110.setZero();
732 m_cdg101.setZero();
733 m_cdg011.setZero();
734 m_cdg300.setZero();
735 m_cdg120.setZero();
736 m_cdg102.setZero();
737 m_cdg210.setZero();
738 m_cdg030.setZero();
739 m_cdg012.setZero();
740 m_cdg201.setZero();
741 m_cdg021.setZero();
742 m_cdg003.setZero();
743 // m_cdg111.setZero();
744 for (qint64 p = 0; p < m_nprim; ++p) {
745 qreal xx0 = xyz(0) - m_X0(p);
746 qreal yy0 = xyz(1) - m_Y0(p);
747 qreal zz0 = xyz(2) - m_Z0(p);
748
749 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
750
751 if (b0arg > m_cutoff) {
752 qint64 aax0 = 1;
753 qint64 aay0 = 1;
754 qint64 aaz0 = 1;
755 qint64 aax1 = m_xamom(p);
756 qint64 aay1 = m_yamom(p);
757 qint64 aaz1 = m_zamom(p);
758 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
759 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
760 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
761 qint64 aax3 = m_xamom(p) * (m_xamom(p) - 1) * (m_xamom(p) - 2);
762 qint64 aay3 = m_yamom(p) * (m_yamom(p) - 1) * (m_yamom(p) - 2);
763 qint64 aaz3 = m_zamom(p) * (m_zamom(p) - 1) * (m_zamom(p) - 2);
764
765 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
766 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
767 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
768
769 qreal ax1;
770 qreal ay1;
771 qreal az1;
772 if (m_xamom(p) < 1) {
773 ax1 = zero;
774 } else if (m_xamom(p) == 1) {
775 ax1 = one;
776 } else {
777 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
778 }
779
780 if (m_yamom(p) < 1) {
781 ay1 = zero;
782 } else if (m_yamom(p) == 1) {
783 ay1 = one;
784 } else {
785 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
786 }
787
788 if (m_zamom(p) < 1) {
789 az1 = zero;
790 } else if (m_zamom(p) == 1) {
791 az1 = one;
792 } else {
793 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
794 }
795
796 qreal ax2;
797 qreal ay2;
798 qreal az2;
799 if (m_xamom(p) < 2) {
800 ax2 = zero;
801 } else if (m_xamom(p) == 2) {
802 ax2 = one;
803 } else {
804 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
805 }
806
807 if (m_yamom(p) < 2) {
808 ay2 = zero;
809 } else if (m_yamom(p) == 2) {
810 ay2 = one;
811 } else {
812 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
813 }
814
815 if (m_zamom(p) < 2) {
816 az2 = zero;
817 } else if (m_zamom(p) == 2) {
818 az2 = one;
819 } else {
820 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
821 }
822
823 qreal ax3;
824 qreal ay3;
825 qreal az3;
826 if (m_xamom(p) < 3) {
827 ax3 = zero;
828 } else if (m_xamom(p) == 3) {
829 ax3 = one;
830 } else {
831 ax3 = aax3 * ipow(xx0, m_xamom(p) - 3);
832 }
833
834 if (m_yamom(p) < 3) {
835 ay3 = zero;
836 } else if (m_yamom(p) == 3) {
837 ay3 = one;
838 } else {
839 ay3 = aay3 * ipow(yy0, m_yamom(p) - 3);
840 }
841
842 if (m_zamom(p) < 3) {
843 az3 = zero;
844 } else if (m_zamom(p) == 3) {
845 az3 = one;
846 } else {
847 az3 = aaz3 * ipow(zz0, m_zamom(p) - 3);
848 }
849
850 qreal b0 = exp(b0arg);
851
852 qreal bx1 = -2 * m_alpha(p) * xx0;
853 qreal by1 = -2 * m_alpha(p) * yy0;
854 qreal bz1 = -2 * m_alpha(p) * zz0;
855 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
856 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
857 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
858 qreal bx3 = (12 * ipow(m_alpha(p), 2) * xx0) -
859 (8 * ipow(m_alpha(p), 3) * ipow(xx0, 3));
860 qreal by3 = (12 * ipow(m_alpha(p), 2) * yy0) -
861 (8 * ipow(m_alpha(p), 3) * ipow(yy0, 3));
862 qreal bz3 = (12 * ipow(m_alpha(p), 2) * zz0) -
863 (8 * ipow(m_alpha(p), 3) * ipow(zz0, 3));
864
865 qreal dg000 = ax0 * ay0 * az0 * b0;
866 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
867 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
868 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
869 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
870 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
871 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
872 qreal dg110 = az0 * b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1);
873 qreal dg101 = ay0 * b0 * (ax1 + ax0 * bx1) * (az1 + az0 * bz1);
874 qreal dg011 = ax0 * b0 * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
875 qreal dg300 =
876 ay0 * az0 * b0 * (ax3 + 3 * ax2 * bx1 + 3 * ax1 * bx2 + ax0 * bx3);
877 qreal dg030 =
878 ax0 * az0 * b0 * (ay3 + 3 * ay2 * by1 + 3 * ay1 * by2 + ay0 * by3);
879 qreal dg003 =
880 ax0 * ay0 * b0 * (az3 + 3 * az2 * bz1 + 3 * az1 * bz2 + az0 * bz3);
881 qreal dg210 =
882 az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (ay1 + ay0 * by1);
883 qreal dg201 =
884 ay0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (az1 + az0 * bz1);
885 qreal dg120 =
886 az0 * b0 * (ax1 + ax0 * bx1) * (ay2 + 2 * ay1 * by1 + ay0 * by2);
887 qreal dg021 =
888 ax0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2) * (az1 + az0 * bz1);
889 qreal dg102 =
890 ay0 * b0 * (ax1 + ax0 * bx1) * (az2 + 2 * az1 * bz1 + az0 * bz2);
891 qreal dg012 =
892 ax0 * b0 * (ay1 + ay0 * by1) * (az2 + 2 * az1 * bz1 + az0 * bz2);
893 // qreal dg111 = b0*(ax1+ax0*bx1)*(ay1+ay0*by1)*(az1+az0*bz1);
894
895 for (qint64 m = 0; m < m_nmo; ++m) {
896 m_cdg000(m) += m_coef(m, p) * dg000;
897 m_cdg100(m) += m_coef(m, p) * dg100;
898 m_cdg010(m) += m_coef(m, p) * dg010;
899 m_cdg001(m) += m_coef(m, p) * dg001;
900 m_cdg200(m) += m_coef(m, p) * dg200;
901 m_cdg020(m) += m_coef(m, p) * dg020;
902 m_cdg002(m) += m_coef(m, p) * dg002;
903 m_cdg110(m) += m_coef(m, p) * dg110;
904 m_cdg101(m) += m_coef(m, p) * dg101;
905 m_cdg011(m) += m_coef(m, p) * dg011;
906 m_cdg300(m) += m_coef(m, p) * dg300;
907 m_cdg030(m) += m_coef(m, p) * dg030;
908 m_cdg003(m) += m_coef(m, p) * dg003;
909 m_cdg210(m) += m_coef(m, p) * dg210;
910 m_cdg201(m) += m_coef(m, p) * dg201;
911 m_cdg120(m) += m_coef(m, p) * dg120;
912 m_cdg021(m) += m_coef(m, p) * dg021;
913 m_cdg102(m) += m_coef(m, p) * dg102;
914 m_cdg012(m) += m_coef(m, p) * dg012;
915 // m_cdg111(m) += m_coef(m,p) * dg111;
916 }
917 }
918 }
919
920 qreal deriv300 = zero;
921 qreal deriv030 = zero;
922 qreal deriv003 = zero;
923 qreal deriv210 = zero;
924 qreal deriv201 = zero;
925 qreal deriv120 = zero;
926 qreal deriv021 = zero;
927 qreal deriv102 = zero;
928 qreal deriv012 = zero;
929 // qreal deriv111=zero;
930 for (qint64 m = 0; m < m_nmo; ++m) {
931 deriv300 += (m_occno(m) * (6 * m_cdg100(m) * m_cdg200(m) +
932 2 * m_cdg000(m) * m_cdg300(m)));
933 deriv030 += (m_occno(m) * (6 * m_cdg010(m) * m_cdg020(m) +
934 2 * m_cdg000(m) * m_cdg030(m)));
935 deriv003 += (m_occno(m) * (6 * m_cdg001(m) * m_cdg002(m) +
936 2 * m_cdg000(m) * m_cdg003(m)));
937 deriv210 += (m_occno(m) *
938 (2 * (2 * m_cdg100(m) * m_cdg110(m) +
939 m_cdg010(m) * m_cdg200(m) + m_cdg000(m) * m_cdg210(m))));
940 deriv201 += (m_occno(m) *
941 (2 * (2 * m_cdg100(m) * m_cdg101(m) +
942 m_cdg001(m) * m_cdg200(m) + m_cdg000(m) * m_cdg201(m))));
943 deriv120 += (m_occno(m) * (2 * (m_cdg020(m) * m_cdg100(m) +
944 2 * m_cdg010(m) * m_cdg110(m) +
945 m_cdg000(m) * m_cdg120(m))));
946 deriv021 += (m_occno(m) *
947 (2 * (2 * m_cdg010(m) * m_cdg011(m) +
948 m_cdg001(m) * m_cdg020(m) + m_cdg000(m) * m_cdg021(m))));
949 deriv102 += (m_occno(m) * (2 * (m_cdg002(m) * m_cdg100(m) +
950 2 * m_cdg001(m) * m_cdg101(m) +
951 m_cdg000(m) * m_cdg102(m))));
952 deriv012 += (m_occno(m) * (2 * (m_cdg002(m) * m_cdg010(m) +
953 2 * m_cdg001(m) * m_cdg011(m) +
954 m_cdg000(m) * m_cdg012(m))));
955 // deriv111+=(m_occno(m)*(
956 // 2*(m_cdg011(m)*m_cdg100(m)+m_cdg010(m)*m_cdg101(m)+m_cdg001(m)*m_cdg110(m)+m_cdg000(m)*m_cdg111(m))
957 // ));
958 }
959
960 value(0) = deriv300 + deriv120 + deriv102;
961 value(1) = deriv210 + deriv030 + deriv012;
962 value(2) = deriv201 + deriv021 + deriv003;
963
964 return value;
965 }
966
967 const Matrix<qreal, 3, 3>
hessianOfElectronDensityLaplacian(const Matrix<qreal,3,1> xyz)968 QTAIMWavefunctionEvaluator::hessianOfElectronDensityLaplacian(
969 const Matrix<qreal, 3, 1> xyz)
970 {
971
972 Matrix<qreal, 3, 3> value;
973
974 const qreal zero = 0.0;
975 const qreal one = 1.0;
976
977 m_cdg000.setZero();
978 m_cdg100.setZero();
979 m_cdg010.setZero();
980 m_cdg001.setZero();
981 m_cdg200.setZero();
982 m_cdg020.setZero();
983 m_cdg002.setZero();
984 m_cdg110.setZero();
985 m_cdg101.setZero();
986 m_cdg011.setZero();
987 m_cdg300.setZero();
988 m_cdg120.setZero();
989 m_cdg102.setZero();
990 m_cdg210.setZero();
991 m_cdg030.setZero();
992 m_cdg012.setZero();
993 m_cdg201.setZero();
994 m_cdg021.setZero();
995 m_cdg003.setZero();
996 m_cdg111.setZero();
997 m_cdg400.setZero();
998 m_cdg040.setZero();
999 m_cdg004.setZero();
1000 m_cdg310.setZero();
1001 m_cdg301.setZero();
1002 m_cdg130.setZero();
1003 m_cdg031.setZero();
1004 m_cdg103.setZero();
1005 m_cdg013.setZero();
1006 m_cdg220.setZero();
1007 m_cdg202.setZero();
1008 m_cdg022.setZero();
1009 m_cdg211.setZero();
1010 m_cdg121.setZero();
1011 m_cdg112.setZero();
1012
1013 for (qint64 p = 0; p < m_nprim; ++p) {
1014 qreal xx0 = xyz(0) - m_X0(p);
1015 qreal yy0 = xyz(1) - m_Y0(p);
1016 qreal zz0 = xyz(2) - m_Z0(p);
1017
1018 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
1019
1020 if (b0arg > m_cutoff) {
1021 qint64 aax0 = 1;
1022 qint64 aay0 = 1;
1023 qint64 aaz0 = 1;
1024 qint64 aax1 = m_xamom(p);
1025 qint64 aay1 = m_yamom(p);
1026 qint64 aaz1 = m_zamom(p);
1027 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
1028 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
1029 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
1030 qint64 aax3 = m_xamom(p) * (m_xamom(p) - 1) * (m_xamom(p) - 2);
1031 qint64 aay3 = m_yamom(p) * (m_yamom(p) - 1) * (m_yamom(p) - 2);
1032 qint64 aaz3 = m_zamom(p) * (m_zamom(p) - 1) * (m_zamom(p) - 2);
1033 qint64 aax4 =
1034 m_xamom(p) * (m_xamom(p) - 1) * (m_xamom(p) - 2) * (m_xamom(p) - 3);
1035 qint64 aay4 =
1036 m_yamom(p) * (m_yamom(p) - 1) * (m_yamom(p) - 2) * (m_xamom(p) - 3);
1037 qint64 aaz4 =
1038 m_zamom(p) * (m_zamom(p) - 1) * (m_zamom(p) - 2) * (m_xamom(p) - 3);
1039
1040 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
1041 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
1042 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
1043
1044 qreal ax1;
1045 qreal ay1;
1046 qreal az1;
1047 if (m_xamom(p) < 1) {
1048 ax1 = zero;
1049 } else if (m_xamom(p) == 1) {
1050 ax1 = one;
1051 } else {
1052 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
1053 }
1054
1055 if (m_yamom(p) < 1) {
1056 ay1 = zero;
1057 } else if (m_yamom(p) == 1) {
1058 ay1 = one;
1059 } else {
1060 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
1061 }
1062
1063 if (m_zamom(p) < 1) {
1064 az1 = zero;
1065 } else if (m_zamom(p) == 1) {
1066 az1 = one;
1067 } else {
1068 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
1069 }
1070
1071 qreal ax2;
1072 qreal ay2;
1073 qreal az2;
1074 if (m_xamom(p) < 2) {
1075 ax2 = zero;
1076 } else if (m_xamom(p) == 2) {
1077 ax2 = one;
1078 } else {
1079 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
1080 }
1081
1082 if (m_yamom(p) < 2) {
1083 ay2 = zero;
1084 } else if (m_yamom(p) == 2) {
1085 ay2 = one;
1086 } else {
1087 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
1088 }
1089
1090 if (m_zamom(p) < 2) {
1091 az2 = zero;
1092 } else if (m_zamom(p) == 2) {
1093 az2 = one;
1094 } else {
1095 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
1096 }
1097
1098 qreal ax3;
1099 qreal ay3;
1100 qreal az3;
1101 if (m_xamom(p) < 3) {
1102 ax3 = zero;
1103 } else if (m_xamom(p) == 3) {
1104 ax3 = one;
1105 } else {
1106 ax3 = aax3 * ipow(xx0, m_xamom(p) - 3);
1107 }
1108
1109 if (m_yamom(p) < 3) {
1110 ay3 = zero;
1111 } else if (m_yamom(p) == 3) {
1112 ay3 = one;
1113 } else {
1114 ay3 = aay3 * ipow(yy0, m_yamom(p) - 3);
1115 }
1116
1117 if (m_zamom(p) < 3) {
1118 az3 = zero;
1119 } else if (m_zamom(p) == 3) {
1120 az3 = one;
1121 } else {
1122 az3 = aaz3 * ipow(zz0, m_zamom(p) - 3);
1123 }
1124
1125 qreal ax4;
1126 qreal ay4;
1127 qreal az4;
1128 if (m_xamom(p) < 4) {
1129 ax4 = zero;
1130 } else if (m_xamom(p) == 4) {
1131 ax4 = one;
1132 } else {
1133 ax4 = aax4 * ipow(xx0, m_xamom(p) - 4);
1134 }
1135
1136 if (m_yamom(p) < 4) {
1137 ay4 = zero;
1138 } else if (m_yamom(p) == 4) {
1139 ay4 = one;
1140 } else {
1141 ay4 = aay4 * ipow(yy0, m_yamom(p) - 4);
1142 }
1143
1144 if (m_zamom(p) < 4) {
1145 az4 = zero;
1146 } else if (m_zamom(p) == 4) {
1147 az4 = one;
1148 } else {
1149 az4 = aaz4 * ipow(zz0, m_zamom(p) - 4);
1150 }
1151
1152 qreal b0 = exp(b0arg);
1153
1154 qreal bx1 = -2 * m_alpha(p) * xx0;
1155 qreal by1 = -2 * m_alpha(p) * yy0;
1156 qreal bz1 = -2 * m_alpha(p) * zz0;
1157 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
1158 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
1159 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
1160 qreal bx3 = (12 * ipow(m_alpha(p), 2) * xx0) -
1161 (8 * ipow(m_alpha(p), 3) * ipow(xx0, 3));
1162 qreal by3 = (12 * ipow(m_alpha(p), 2) * yy0) -
1163 (8 * ipow(m_alpha(p), 3) * ipow(yy0, 3));
1164 qreal bz3 = (12 * ipow(m_alpha(p), 2) * zz0) -
1165 (8 * ipow(m_alpha(p), 3) * ipow(zz0, 3));
1166 qreal bx4 = (12 * ipow(m_alpha(p), 2)) -
1167 (48 * ipow(m_alpha(p), 3) * ipow(xx0, 2)) +
1168 (16 * ipow(m_alpha(p), 4) * ipow(xx0, 4));
1169 qreal by4 = (12 * ipow(m_alpha(p), 2)) -
1170 (48 * ipow(m_alpha(p), 3) * ipow(yy0, 2)) +
1171 (16 * ipow(m_alpha(p), 4) * ipow(yy0, 4));
1172 qreal bz4 = (12 * ipow(m_alpha(p), 2)) -
1173 (48 * ipow(m_alpha(p), 3) * ipow(zz0, 2)) +
1174 (16 * ipow(m_alpha(p), 4) * ipow(zz0, 4));
1175
1176 qreal dg000 = ax0 * ay0 * az0 * b0;
1177 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
1178 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
1179 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
1180 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
1181 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
1182 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
1183 qreal dg110 = az0 * b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1);
1184 qreal dg101 = ay0 * b0 * (ax1 + ax0 * bx1) * (az1 + az0 * bz1);
1185 qreal dg011 = ax0 * b0 * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
1186 qreal dg300 =
1187 ay0 * az0 * b0 * (ax3 + 3 * ax2 * bx1 + 3 * ax1 * bx2 + ax0 * bx3);
1188 qreal dg030 =
1189 ax0 * az0 * b0 * (ay3 + 3 * ay2 * by1 + 3 * ay1 * by2 + ay0 * by3);
1190 qreal dg003 =
1191 ax0 * ay0 * b0 * (az3 + 3 * az2 * bz1 + 3 * az1 * bz2 + az0 * bz3);
1192 qreal dg210 =
1193 az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (ay1 + ay0 * by1);
1194 qreal dg201 =
1195 ay0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (az1 + az0 * bz1);
1196 qreal dg120 =
1197 az0 * b0 * (ax1 + ax0 * bx1) * (ay2 + 2 * ay1 * by1 + ay0 * by2);
1198 qreal dg021 =
1199 ax0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2) * (az1 + az0 * bz1);
1200 qreal dg102 =
1201 ay0 * b0 * (ax1 + ax0 * bx1) * (az2 + 2 * az1 * bz1 + az0 * bz2);
1202 qreal dg012 =
1203 ax0 * b0 * (ay1 + ay0 * by1) * (az2 + 2 * az1 * bz1 + az0 * bz2);
1204 qreal dg111 =
1205 b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
1206 qreal dg400 = ay0 * az0 * b0 * (ax4 + 4 * ax3 * bx1 + 6 * ax2 * bx2 +
1207 4 * ax1 * bx3 + ax0 * bx4);
1208 qreal dg040 = ax0 * az0 * b0 * (ay4 + 4 * ay3 * by1 + 6 * ay2 * by2 +
1209 4 * ay1 * by3 + ay0 * by4);
1210 qreal dg004 = ax0 * ay0 * b0 * (az4 + 4 * az3 * bz1 + 6 * az2 * bz2 +
1211 4 * az1 * bz3 + az0 * bz4);
1212 qreal dg310 = az0 * b0 *
1213 (ax3 + 3 * ax2 * bx1 + 3 * ax1 * bx2 + ax0 * bx3) *
1214 (ay1 + ay0 * by1);
1215 qreal dg301 = ay0 * b0 *
1216 (ax3 + 3 * ax2 * bx1 + 3 * ax1 * bx2 + ax0 * bx3) *
1217 (az1 + az0 * bz1);
1218 qreal dg130 = az0 * b0 * (ax1 + ax0 * bx1) *
1219 (ay3 + 3 * ay2 * by1 + 3 * ay1 * by2 + ay0 * by3);
1220 qreal dg031 = ax0 * b0 *
1221 (ay3 + 3 * ay2 * by1 + 3 * ay1 * by2 + ay0 * by3) *
1222 (az1 + az0 * bz1);
1223 qreal dg103 = ay0 * b0 * (ax1 + ax0 * bx1) *
1224 (az3 + 3 * az2 * bz1 + 3 * az1 * bz2 + az0 * bz3);
1225 qreal dg013 = ax0 * b0 * (ay1 + ay0 * by1) *
1226 (az3 + 3 * az2 * bz1 + 3 * az1 * bz2 + az0 * bz3);
1227 qreal dg220 = az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) *
1228 (ay2 + 2 * ay1 * by1 + ay0 * by2);
1229 qreal dg202 = ay0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) *
1230 (az2 + 2 * az1 * bz1 + az0 * bz2);
1231 qreal dg022 = ax0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2) *
1232 (az2 + 2 * az1 * bz1 + az0 * bz2);
1233 qreal dg211 = b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (ay1 + ay0 * by1) *
1234 (az1 + az0 * bz1);
1235 qreal dg121 = b0 * (ax1 + ax0 * bx1) * (ay2 + 2 * ay1 * by1 + ay0 * by2) *
1236 (az1 + az0 * bz1);
1237 qreal dg112 = b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1) *
1238 (az2 + 2 * az1 * bz1 + az0 * bz2);
1239
1240 for (qint64 m = 0; m < m_nmo; ++m) {
1241 m_cdg000(m) += m_coef(m, p) * dg000;
1242 m_cdg100(m) += m_coef(m, p) * dg100;
1243 m_cdg010(m) += m_coef(m, p) * dg010;
1244 m_cdg001(m) += m_coef(m, p) * dg001;
1245 m_cdg200(m) += m_coef(m, p) * dg200;
1246 m_cdg020(m) += m_coef(m, p) * dg020;
1247 m_cdg002(m) += m_coef(m, p) * dg002;
1248 m_cdg110(m) += m_coef(m, p) * dg110;
1249 m_cdg101(m) += m_coef(m, p) * dg101;
1250 m_cdg011(m) += m_coef(m, p) * dg011;
1251 m_cdg300(m) += m_coef(m, p) * dg300;
1252 m_cdg030(m) += m_coef(m, p) * dg030;
1253 m_cdg003(m) += m_coef(m, p) * dg003;
1254 m_cdg210(m) += m_coef(m, p) * dg210;
1255 m_cdg201(m) += m_coef(m, p) * dg201;
1256 m_cdg120(m) += m_coef(m, p) * dg120;
1257 m_cdg021(m) += m_coef(m, p) * dg021;
1258 m_cdg102(m) += m_coef(m, p) * dg102;
1259 m_cdg012(m) += m_coef(m, p) * dg012;
1260 m_cdg111(m) += m_coef(m, p) * dg111;
1261 m_cdg400(m) += m_coef(m, p) * dg400;
1262 m_cdg040(m) += m_coef(m, p) * dg040;
1263 m_cdg004(m) += m_coef(m, p) * dg004;
1264 m_cdg310(m) += m_coef(m, p) * dg310;
1265 m_cdg301(m) += m_coef(m, p) * dg301;
1266 m_cdg130(m) += m_coef(m, p) * dg130;
1267 m_cdg031(m) += m_coef(m, p) * dg031;
1268 m_cdg103(m) += m_coef(m, p) * dg103;
1269 m_cdg013(m) += m_coef(m, p) * dg013;
1270 m_cdg220(m) += m_coef(m, p) * dg220;
1271 m_cdg202(m) += m_coef(m, p) * dg202;
1272 m_cdg022(m) += m_coef(m, p) * dg022;
1273 m_cdg211(m) += m_coef(m, p) * dg211;
1274 m_cdg121(m) += m_coef(m, p) * dg121;
1275 m_cdg112(m) += m_coef(m, p) * dg112;
1276 }
1277 }
1278 }
1279
1280 qreal deriv400 = zero;
1281 qreal deriv040 = zero;
1282 qreal deriv004 = zero;
1283 qreal deriv310 = zero;
1284 qreal deriv301 = zero;
1285 qreal deriv130 = zero;
1286 qreal deriv031 = zero;
1287 qreal deriv103 = zero;
1288 qreal deriv013 = zero;
1289 qreal deriv220 = zero;
1290 qreal deriv202 = zero;
1291 qreal deriv022 = zero;
1292 qreal deriv211 = zero;
1293 qreal deriv121 = zero;
1294 qreal deriv112 = zero;
1295 for (qint64 m = 0; m < m_nmo; ++m) {
1296 deriv400 +=
1297 (m_occno(m) * (6 * ipow(m_cdg200(m), 2) + 8 * m_cdg100(m) * m_cdg300(m) +
1298 2 * m_cdg000(m) * m_cdg400(m)));
1299 deriv040 +=
1300 (m_occno(m) * (6 * ipow(m_cdg020(m), 2) + 8 * m_cdg010(m) * m_cdg030(m) +
1301 2 * m_cdg000(m) * m_cdg040(m)));
1302 deriv004 +=
1303 (m_occno(m) * (6 * ipow(m_cdg002(m), 2) + 8 * m_cdg001(m) * m_cdg003(m) +
1304 2 * m_cdg000(m) * m_cdg004(m)));
1305 deriv310 +=
1306 (m_occno(m) *
1307 (2 * (3 * m_cdg110(m) * m_cdg200(m) + 3 * m_cdg100(m) * m_cdg210(m) +
1308 m_cdg010(m) * m_cdg300(m) + m_cdg000(m) * m_cdg310(m))));
1309 deriv301 +=
1310 (m_occno(m) *
1311 (2 * (3 * m_cdg101(m) * m_cdg200(m) + 3 * m_cdg100(m) * m_cdg201(m) +
1312 m_cdg001(m) * m_cdg300(m) + m_cdg000(m) * m_cdg301(m))));
1313 deriv130 +=
1314 (m_occno(m) *
1315 (2 * (m_cdg030(m) * m_cdg100(m) + 3 * m_cdg020(m) * m_cdg110(m) +
1316 3 * m_cdg010(m) * m_cdg120(m) + m_cdg000(m) * m_cdg130(m))));
1317 deriv031 +=
1318 (m_occno(m) *
1319 (2 * (3 * m_cdg011(m) * m_cdg020(m) + 3 * m_cdg010(m) * m_cdg021(m) +
1320 m_cdg001(m) * m_cdg030(m) + m_cdg000(m) * m_cdg031(m))));
1321 deriv103 +=
1322 (m_occno(m) *
1323 (2 * (m_cdg003(m) * m_cdg100(m) + 3 * m_cdg002(m) * m_cdg101(m) +
1324 3 * m_cdg001(m) * m_cdg102(m) + m_cdg000(m) * m_cdg103(m))));
1325 deriv013 +=
1326 (m_occno(m) *
1327 (2 * (m_cdg003(m) * m_cdg010(m) + 3 * m_cdg002(m) * m_cdg011(m) +
1328 3 * m_cdg001(m) * m_cdg012(m) + m_cdg000(m) * m_cdg013(m))));
1329 deriv220 +=
1330 (m_occno(m) *
1331 (2 * (2 * ipow(m_cdg110(m), 2) + 2 * m_cdg100(m) * m_cdg120(m) +
1332 m_cdg020(m) * m_cdg200(m) + 2 * m_cdg010(m) * m_cdg210(m) +
1333 m_cdg000(m) * m_cdg220(m))));
1334 deriv202 +=
1335 (m_occno(m) *
1336 (2 * (2 * ipow(m_cdg101(m), 2) + 2 * m_cdg100(m) * m_cdg102(m) +
1337 m_cdg002(m) * m_cdg200(m) + 2 * m_cdg001(m) * m_cdg201(m) +
1338 m_cdg000(m) * m_cdg202(m))));
1339 deriv022 +=
1340 (m_occno(m) *
1341 (2 * (2 * ipow(m_cdg011(m), 2) + 2 * m_cdg010(m) * m_cdg012(m) +
1342 m_cdg002(m) * m_cdg020(m) + 2 * m_cdg001(m) * m_cdg021(m) +
1343 m_cdg000(m) * m_cdg022(m))));
1344 deriv211 +=
1345 (m_occno(m) *
1346 (2 * (2 * m_cdg101(m) * m_cdg110(m) + 2 * m_cdg100(m) * m_cdg111(m) +
1347 m_cdg011(m) * m_cdg200(m) + m_cdg010(m) * m_cdg201(m) +
1348 m_cdg001(m) * m_cdg210(m) + m_cdg000(m) * m_cdg211(m))));
1349 deriv121 +=
1350 (m_occno(m) *
1351 (2 * (m_cdg021(m) * m_cdg100(m) + m_cdg020(m) * m_cdg101(m) +
1352 2 * m_cdg011(m) * m_cdg110(m) + 2 * m_cdg010(m) * m_cdg111(m) +
1353 m_cdg001(m) * m_cdg120(m) + m_cdg000(m) * m_cdg121(m))));
1354 deriv112 +=
1355 (m_occno(m) *
1356 (2 * (m_cdg012(m) * m_cdg100(m) + 2 * m_cdg011(m) * m_cdg101(m) +
1357 m_cdg010(m) * m_cdg102(m) + m_cdg002(m) * m_cdg110(m) +
1358 2 * m_cdg001(m) * m_cdg111(m) + m_cdg000(m) * m_cdg112(m))));
1359 }
1360
1361 value(0, 0) = deriv400 + deriv220 + deriv202;
1362 value(1, 1) = deriv220 + deriv040 + deriv022;
1363 value(2, 2) = deriv202 + deriv022 + deriv004;
1364 value(0, 1) = deriv310 + deriv130 + deriv112;
1365 value(0, 2) = deriv301 + deriv121 + deriv103;
1366 value(1, 2) = deriv211 + deriv031 + deriv013;
1367 value(1, 0) = value(0, 1);
1368 value(2, 0) = value(0, 2);
1369 value(2, 1) = value(1, 2);
1370
1371 return value;
1372 }
1373
1374 const Matrix<qreal, 3, 4>
gradientAndHessianOfElectronDensityLaplacian(const Matrix<qreal,3,1> xyz)1375 QTAIMWavefunctionEvaluator::gradientAndHessianOfElectronDensityLaplacian(
1376 const Matrix<qreal, 3, 1> xyz)
1377 {
1378
1379 Matrix<qreal, 3, 1> gValue;
1380 Matrix<qreal, 3, 3> hValue;
1381 Matrix<qreal, 3, 4> value;
1382
1383 const qreal zero = 0.0;
1384 const qreal one = 1.0;
1385
1386 m_cdg000.setZero();
1387 m_cdg100.setZero();
1388 m_cdg010.setZero();
1389 m_cdg001.setZero();
1390 m_cdg200.setZero();
1391 m_cdg020.setZero();
1392 m_cdg002.setZero();
1393 m_cdg110.setZero();
1394 m_cdg101.setZero();
1395 m_cdg011.setZero();
1396 m_cdg300.setZero();
1397 m_cdg120.setZero();
1398 m_cdg102.setZero();
1399 m_cdg210.setZero();
1400 m_cdg030.setZero();
1401 m_cdg012.setZero();
1402 m_cdg201.setZero();
1403 m_cdg021.setZero();
1404 m_cdg003.setZero();
1405 m_cdg111.setZero();
1406 m_cdg400.setZero();
1407 m_cdg040.setZero();
1408 m_cdg004.setZero();
1409 m_cdg310.setZero();
1410 m_cdg301.setZero();
1411 m_cdg130.setZero();
1412 m_cdg031.setZero();
1413 m_cdg103.setZero();
1414 m_cdg013.setZero();
1415 m_cdg220.setZero();
1416 m_cdg202.setZero();
1417 m_cdg022.setZero();
1418 m_cdg211.setZero();
1419 m_cdg121.setZero();
1420 m_cdg112.setZero();
1421
1422 for (qint64 p = 0; p < m_nprim; ++p) {
1423 qreal xx0 = xyz(0) - m_X0(p);
1424 qreal yy0 = xyz(1) - m_Y0(p);
1425 qreal zz0 = xyz(2) - m_Z0(p);
1426
1427 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
1428
1429 if (b0arg > m_cutoff) {
1430 qint64 aax0 = 1;
1431 qint64 aay0 = 1;
1432 qint64 aaz0 = 1;
1433 qint64 aax1 = m_xamom(p);
1434 qint64 aay1 = m_yamom(p);
1435 qint64 aaz1 = m_zamom(p);
1436 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
1437 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
1438 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
1439 qint64 aax3 = m_xamom(p) * (m_xamom(p) - 1) * (m_xamom(p) - 2);
1440 qint64 aay3 = m_yamom(p) * (m_yamom(p) - 1) * (m_yamom(p) - 2);
1441 qint64 aaz3 = m_zamom(p) * (m_zamom(p) - 1) * (m_zamom(p) - 2);
1442 qint64 aax4 =
1443 m_xamom(p) * (m_xamom(p) - 1) * (m_xamom(p) - 2) * (m_xamom(p) - 3);
1444 qint64 aay4 =
1445 m_yamom(p) * (m_yamom(p) - 1) * (m_yamom(p) - 2) * (m_xamom(p) - 3);
1446 qint64 aaz4 =
1447 m_zamom(p) * (m_zamom(p) - 1) * (m_zamom(p) - 2) * (m_xamom(p) - 3);
1448
1449 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
1450 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
1451 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
1452
1453 qreal ax1;
1454 qreal ay1;
1455 qreal az1;
1456 if (m_xamom(p) < 1) {
1457 ax1 = zero;
1458 } else if (m_xamom(p) == 1) {
1459 ax1 = one;
1460 } else {
1461 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
1462 }
1463
1464 if (m_yamom(p) < 1) {
1465 ay1 = zero;
1466 } else if (m_yamom(p) == 1) {
1467 ay1 = one;
1468 } else {
1469 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
1470 }
1471
1472 if (m_zamom(p) < 1) {
1473 az1 = zero;
1474 } else if (m_zamom(p) == 1) {
1475 az1 = one;
1476 } else {
1477 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
1478 }
1479
1480 qreal ax2;
1481 qreal ay2;
1482 qreal az2;
1483 if (m_xamom(p) < 2) {
1484 ax2 = zero;
1485 } else if (m_xamom(p) == 2) {
1486 ax2 = one;
1487 } else {
1488 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
1489 }
1490
1491 if (m_yamom(p) < 2) {
1492 ay2 = zero;
1493 } else if (m_yamom(p) == 2) {
1494 ay2 = one;
1495 } else {
1496 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
1497 }
1498
1499 if (m_zamom(p) < 2) {
1500 az2 = zero;
1501 } else if (m_zamom(p) == 2) {
1502 az2 = one;
1503 } else {
1504 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
1505 }
1506
1507 qreal ax3;
1508 qreal ay3;
1509 qreal az3;
1510 if (m_xamom(p) < 3) {
1511 ax3 = zero;
1512 } else if (m_xamom(p) == 3) {
1513 ax3 = one;
1514 } else {
1515 ax3 = aax3 * ipow(xx0, m_xamom(p) - 3);
1516 }
1517
1518 if (m_yamom(p) < 3) {
1519 ay3 = zero;
1520 } else if (m_yamom(p) == 3) {
1521 ay3 = one;
1522 } else {
1523 ay3 = aay3 * ipow(yy0, m_yamom(p) - 3);
1524 }
1525
1526 if (m_zamom(p) < 3) {
1527 az3 = zero;
1528 } else if (m_zamom(p) == 3) {
1529 az3 = one;
1530 } else {
1531 az3 = aaz3 * ipow(zz0, m_zamom(p) - 3);
1532 }
1533
1534 qreal ax4;
1535 qreal ay4;
1536 qreal az4;
1537 if (m_xamom(p) < 4) {
1538 ax4 = zero;
1539 } else if (m_xamom(p) == 4) {
1540 ax4 = one;
1541 } else {
1542 ax4 = aax4 * ipow(xx0, m_xamom(p) - 4);
1543 }
1544
1545 if (m_yamom(p) < 4) {
1546 ay4 = zero;
1547 } else if (m_yamom(p) == 4) {
1548 ay4 = one;
1549 } else {
1550 ay4 = aay4 * ipow(yy0, m_yamom(p) - 4);
1551 }
1552
1553 if (m_zamom(p) < 4) {
1554 az4 = zero;
1555 } else if (m_zamom(p) == 4) {
1556 az4 = one;
1557 } else {
1558 az4 = aaz4 * ipow(zz0, m_zamom(p) - 4);
1559 }
1560
1561 qreal b0 = exp(b0arg);
1562
1563 qreal bx1 = -2 * m_alpha(p) * xx0;
1564 qreal by1 = -2 * m_alpha(p) * yy0;
1565 qreal bz1 = -2 * m_alpha(p) * zz0;
1566 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
1567 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
1568 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
1569 qreal bx3 = (12 * ipow(m_alpha(p), 2) * xx0) -
1570 (8 * ipow(m_alpha(p), 3) * ipow(xx0, 3));
1571 qreal by3 = (12 * ipow(m_alpha(p), 2) * yy0) -
1572 (8 * ipow(m_alpha(p), 3) * ipow(yy0, 3));
1573 qreal bz3 = (12 * ipow(m_alpha(p), 2) * zz0) -
1574 (8 * ipow(m_alpha(p), 3) * ipow(zz0, 3));
1575 qreal bx4 = (12 * ipow(m_alpha(p), 2)) -
1576 (48 * ipow(m_alpha(p), 3) * ipow(xx0, 2)) +
1577 (16 * ipow(m_alpha(p), 4) * ipow(xx0, 4));
1578 qreal by4 = (12 * ipow(m_alpha(p), 2)) -
1579 (48 * ipow(m_alpha(p), 3) * ipow(yy0, 2)) +
1580 (16 * ipow(m_alpha(p), 4) * ipow(yy0, 4));
1581 qreal bz4 = (12 * ipow(m_alpha(p), 2)) -
1582 (48 * ipow(m_alpha(p), 3) * ipow(zz0, 2)) +
1583 (16 * ipow(m_alpha(p), 4) * ipow(zz0, 4));
1584
1585 qreal dg000 = ax0 * ay0 * az0 * b0;
1586 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
1587 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
1588 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
1589 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
1590 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
1591 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
1592 qreal dg110 = az0 * b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1);
1593 qreal dg101 = ay0 * b0 * (ax1 + ax0 * bx1) * (az1 + az0 * bz1);
1594 qreal dg011 = ax0 * b0 * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
1595 qreal dg300 =
1596 ay0 * az0 * b0 * (ax3 + 3 * ax2 * bx1 + 3 * ax1 * bx2 + ax0 * bx3);
1597 qreal dg030 =
1598 ax0 * az0 * b0 * (ay3 + 3 * ay2 * by1 + 3 * ay1 * by2 + ay0 * by3);
1599 qreal dg003 =
1600 ax0 * ay0 * b0 * (az3 + 3 * az2 * bz1 + 3 * az1 * bz2 + az0 * bz3);
1601 qreal dg210 =
1602 az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (ay1 + ay0 * by1);
1603 qreal dg201 =
1604 ay0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (az1 + az0 * bz1);
1605 qreal dg120 =
1606 az0 * b0 * (ax1 + ax0 * bx1) * (ay2 + 2 * ay1 * by1 + ay0 * by2);
1607 qreal dg021 =
1608 ax0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2) * (az1 + az0 * bz1);
1609 qreal dg102 =
1610 ay0 * b0 * (ax1 + ax0 * bx1) * (az2 + 2 * az1 * bz1 + az0 * bz2);
1611 qreal dg012 =
1612 ax0 * b0 * (ay1 + ay0 * by1) * (az2 + 2 * az1 * bz1 + az0 * bz2);
1613 qreal dg111 =
1614 b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
1615 qreal dg400 = ay0 * az0 * b0 * (ax4 + 4 * ax3 * bx1 + 6 * ax2 * bx2 +
1616 4 * ax1 * bx3 + ax0 * bx4);
1617 qreal dg040 = ax0 * az0 * b0 * (ay4 + 4 * ay3 * by1 + 6 * ay2 * by2 +
1618 4 * ay1 * by3 + ay0 * by4);
1619 qreal dg004 = ax0 * ay0 * b0 * (az4 + 4 * az3 * bz1 + 6 * az2 * bz2 +
1620 4 * az1 * bz3 + az0 * bz4);
1621 qreal dg310 = az0 * b0 *
1622 (ax3 + 3 * ax2 * bx1 + 3 * ax1 * bx2 + ax0 * bx3) *
1623 (ay1 + ay0 * by1);
1624 qreal dg301 = ay0 * b0 *
1625 (ax3 + 3 * ax2 * bx1 + 3 * ax1 * bx2 + ax0 * bx3) *
1626 (az1 + az0 * bz1);
1627 qreal dg130 = az0 * b0 * (ax1 + ax0 * bx1) *
1628 (ay3 + 3 * ay2 * by1 + 3 * ay1 * by2 + ay0 * by3);
1629 qreal dg031 = ax0 * b0 *
1630 (ay3 + 3 * ay2 * by1 + 3 * ay1 * by2 + ay0 * by3) *
1631 (az1 + az0 * bz1);
1632 qreal dg103 = ay0 * b0 * (ax1 + ax0 * bx1) *
1633 (az3 + 3 * az2 * bz1 + 3 * az1 * bz2 + az0 * bz3);
1634 qreal dg013 = ax0 * b0 * (ay1 + ay0 * by1) *
1635 (az3 + 3 * az2 * bz1 + 3 * az1 * bz2 + az0 * bz3);
1636 qreal dg220 = az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) *
1637 (ay2 + 2 * ay1 * by1 + ay0 * by2);
1638 qreal dg202 = ay0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) *
1639 (az2 + 2 * az1 * bz1 + az0 * bz2);
1640 qreal dg022 = ax0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2) *
1641 (az2 + 2 * az1 * bz1 + az0 * bz2);
1642 qreal dg211 = b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2) * (ay1 + ay0 * by1) *
1643 (az1 + az0 * bz1);
1644 qreal dg121 = b0 * (ax1 + ax0 * bx1) * (ay2 + 2 * ay1 * by1 + ay0 * by2) *
1645 (az1 + az0 * bz1);
1646 qreal dg112 = b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1) *
1647 (az2 + 2 * az1 * bz1 + az0 * bz2);
1648
1649 for (qint64 m = 0; m < m_nmo; ++m) {
1650 m_cdg000(m) += m_coef(m, p) * dg000;
1651 m_cdg100(m) += m_coef(m, p) * dg100;
1652 m_cdg010(m) += m_coef(m, p) * dg010;
1653 m_cdg001(m) += m_coef(m, p) * dg001;
1654 m_cdg200(m) += m_coef(m, p) * dg200;
1655 m_cdg020(m) += m_coef(m, p) * dg020;
1656 m_cdg002(m) += m_coef(m, p) * dg002;
1657 m_cdg110(m) += m_coef(m, p) * dg110;
1658 m_cdg101(m) += m_coef(m, p) * dg101;
1659 m_cdg011(m) += m_coef(m, p) * dg011;
1660 m_cdg300(m) += m_coef(m, p) * dg300;
1661 m_cdg030(m) += m_coef(m, p) * dg030;
1662 m_cdg003(m) += m_coef(m, p) * dg003;
1663 m_cdg210(m) += m_coef(m, p) * dg210;
1664 m_cdg201(m) += m_coef(m, p) * dg201;
1665 m_cdg120(m) += m_coef(m, p) * dg120;
1666 m_cdg021(m) += m_coef(m, p) * dg021;
1667 m_cdg102(m) += m_coef(m, p) * dg102;
1668 m_cdg012(m) += m_coef(m, p) * dg012;
1669 m_cdg111(m) += m_coef(m, p) * dg111;
1670 m_cdg400(m) += m_coef(m, p) * dg400;
1671 m_cdg040(m) += m_coef(m, p) * dg040;
1672 m_cdg004(m) += m_coef(m, p) * dg004;
1673 m_cdg310(m) += m_coef(m, p) * dg310;
1674 m_cdg301(m) += m_coef(m, p) * dg301;
1675 m_cdg130(m) += m_coef(m, p) * dg130;
1676 m_cdg031(m) += m_coef(m, p) * dg031;
1677 m_cdg103(m) += m_coef(m, p) * dg103;
1678 m_cdg013(m) += m_coef(m, p) * dg013;
1679 m_cdg220(m) += m_coef(m, p) * dg220;
1680 m_cdg202(m) += m_coef(m, p) * dg202;
1681 m_cdg022(m) += m_coef(m, p) * dg022;
1682 m_cdg211(m) += m_coef(m, p) * dg211;
1683 m_cdg121(m) += m_coef(m, p) * dg121;
1684 m_cdg112(m) += m_coef(m, p) * dg112;
1685 }
1686 }
1687 }
1688
1689 qreal deriv300 = zero;
1690 qreal deriv030 = zero;
1691 qreal deriv003 = zero;
1692 qreal deriv210 = zero;
1693 qreal deriv201 = zero;
1694 qreal deriv120 = zero;
1695 qreal deriv021 = zero;
1696 qreal deriv102 = zero;
1697 qreal deriv012 = zero;
1698 qreal deriv400 = zero;
1699 qreal deriv040 = zero;
1700 qreal deriv004 = zero;
1701 qreal deriv310 = zero;
1702 qreal deriv301 = zero;
1703 qreal deriv130 = zero;
1704 qreal deriv031 = zero;
1705 qreal deriv103 = zero;
1706 qreal deriv013 = zero;
1707 qreal deriv220 = zero;
1708 qreal deriv202 = zero;
1709 qreal deriv022 = zero;
1710 qreal deriv211 = zero;
1711 qreal deriv121 = zero;
1712 qreal deriv112 = zero;
1713 for (qint64 m = 0; m < m_nmo; ++m) {
1714 deriv300 += (m_occno(m) * (6 * m_cdg100(m) * m_cdg200(m) +
1715 2 * m_cdg000(m) * m_cdg300(m)));
1716 deriv030 += (m_occno(m) * (6 * m_cdg010(m) * m_cdg020(m) +
1717 2 * m_cdg000(m) * m_cdg030(m)));
1718 deriv003 += (m_occno(m) * (6 * m_cdg001(m) * m_cdg002(m) +
1719 2 * m_cdg000(m) * m_cdg003(m)));
1720 deriv210 += (m_occno(m) *
1721 (2 * (2 * m_cdg100(m) * m_cdg110(m) +
1722 m_cdg010(m) * m_cdg200(m) + m_cdg000(m) * m_cdg210(m))));
1723 deriv201 += (m_occno(m) *
1724 (2 * (2 * m_cdg100(m) * m_cdg101(m) +
1725 m_cdg001(m) * m_cdg200(m) + m_cdg000(m) * m_cdg201(m))));
1726 deriv120 += (m_occno(m) * (2 * (m_cdg020(m) * m_cdg100(m) +
1727 2 * m_cdg010(m) * m_cdg110(m) +
1728 m_cdg000(m) * m_cdg120(m))));
1729 deriv021 += (m_occno(m) *
1730 (2 * (2 * m_cdg010(m) * m_cdg011(m) +
1731 m_cdg001(m) * m_cdg020(m) + m_cdg000(m) * m_cdg021(m))));
1732 deriv102 += (m_occno(m) * (2 * (m_cdg002(m) * m_cdg100(m) +
1733 2 * m_cdg001(m) * m_cdg101(m) +
1734 m_cdg000(m) * m_cdg102(m))));
1735 deriv012 += (m_occno(m) * (2 * (m_cdg002(m) * m_cdg010(m) +
1736 2 * m_cdg001(m) * m_cdg011(m) +
1737 m_cdg000(m) * m_cdg012(m))));
1738 // deriv111+=(m_occno(m)*(
1739 // 2*(m_cdg011(m)*m_cdg100(m)+m_cdg010(m)*m_cdg101(m)+m_cdg001(m)*m_cdg110(m)+m_cdg000(m)*m_cdg111(m))
1740 // ));
1741 deriv400 +=
1742 (m_occno(m) * (6 * ipow(m_cdg200(m), 2) + 8 * m_cdg100(m) * m_cdg300(m) +
1743 2 * m_cdg000(m) * m_cdg400(m)));
1744 deriv040 +=
1745 (m_occno(m) * (6 * ipow(m_cdg020(m), 2) + 8 * m_cdg010(m) * m_cdg030(m) +
1746 2 * m_cdg000(m) * m_cdg040(m)));
1747 deriv004 +=
1748 (m_occno(m) * (6 * ipow(m_cdg002(m), 2) + 8 * m_cdg001(m) * m_cdg003(m) +
1749 2 * m_cdg000(m) * m_cdg004(m)));
1750 deriv310 +=
1751 (m_occno(m) *
1752 (2 * (3 * m_cdg110(m) * m_cdg200(m) + 3 * m_cdg100(m) * m_cdg210(m) +
1753 m_cdg010(m) * m_cdg300(m) + m_cdg000(m) * m_cdg310(m))));
1754 deriv301 +=
1755 (m_occno(m) *
1756 (2 * (3 * m_cdg101(m) * m_cdg200(m) + 3 * m_cdg100(m) * m_cdg201(m) +
1757 m_cdg001(m) * m_cdg300(m) + m_cdg000(m) * m_cdg301(m))));
1758 deriv130 +=
1759 (m_occno(m) *
1760 (2 * (m_cdg030(m) * m_cdg100(m) + 3 * m_cdg020(m) * m_cdg110(m) +
1761 3 * m_cdg010(m) * m_cdg120(m) + m_cdg000(m) * m_cdg130(m))));
1762 deriv031 +=
1763 (m_occno(m) *
1764 (2 * (3 * m_cdg011(m) * m_cdg020(m) + 3 * m_cdg010(m) * m_cdg021(m) +
1765 m_cdg001(m) * m_cdg030(m) + m_cdg000(m) * m_cdg031(m))));
1766 deriv103 +=
1767 (m_occno(m) *
1768 (2 * (m_cdg003(m) * m_cdg100(m) + 3 * m_cdg002(m) * m_cdg101(m) +
1769 3 * m_cdg001(m) * m_cdg102(m) + m_cdg000(m) * m_cdg103(m))));
1770 deriv013 +=
1771 (m_occno(m) *
1772 (2 * (m_cdg003(m) * m_cdg010(m) + 3 * m_cdg002(m) * m_cdg011(m) +
1773 3 * m_cdg001(m) * m_cdg012(m) + m_cdg000(m) * m_cdg013(m))));
1774 deriv220 +=
1775 (m_occno(m) *
1776 (2 * (2 * ipow(m_cdg110(m), 2) + 2 * m_cdg100(m) * m_cdg120(m) +
1777 m_cdg020(m) * m_cdg200(m) + 2 * m_cdg010(m) * m_cdg210(m) +
1778 m_cdg000(m) * m_cdg220(m))));
1779 deriv202 +=
1780 (m_occno(m) *
1781 (2 * (2 * ipow(m_cdg101(m), 2) + 2 * m_cdg100(m) * m_cdg102(m) +
1782 m_cdg002(m) * m_cdg200(m) + 2 * m_cdg001(m) * m_cdg201(m) +
1783 m_cdg000(m) * m_cdg202(m))));
1784 deriv022 +=
1785 (m_occno(m) *
1786 (2 * (2 * ipow(m_cdg011(m), 2) + 2 * m_cdg010(m) * m_cdg012(m) +
1787 m_cdg002(m) * m_cdg020(m) + 2 * m_cdg001(m) * m_cdg021(m) +
1788 m_cdg000(m) * m_cdg022(m))));
1789 deriv211 +=
1790 (m_occno(m) *
1791 (2 * (2 * m_cdg101(m) * m_cdg110(m) + 2 * m_cdg100(m) * m_cdg111(m) +
1792 m_cdg011(m) * m_cdg200(m) + m_cdg010(m) * m_cdg201(m) +
1793 m_cdg001(m) * m_cdg210(m) + m_cdg000(m) * m_cdg211(m))));
1794 deriv121 +=
1795 (m_occno(m) *
1796 (2 * (m_cdg021(m) * m_cdg100(m) + m_cdg020(m) * m_cdg101(m) +
1797 2 * m_cdg011(m) * m_cdg110(m) + 2 * m_cdg010(m) * m_cdg111(m) +
1798 m_cdg001(m) * m_cdg120(m) + m_cdg000(m) * m_cdg121(m))));
1799 deriv112 +=
1800 (m_occno(m) *
1801 (2 * (m_cdg012(m) * m_cdg100(m) + 2 * m_cdg011(m) * m_cdg101(m) +
1802 m_cdg010(m) * m_cdg102(m) + m_cdg002(m) * m_cdg110(m) +
1803 2 * m_cdg001(m) * m_cdg111(m) + m_cdg000(m) * m_cdg112(m))));
1804 }
1805
1806 gValue(0) = deriv300 + deriv120 + deriv102;
1807 gValue(1) = deriv210 + deriv030 + deriv012;
1808 gValue(2) = deriv201 + deriv021 + deriv003;
1809
1810 hValue(0, 0) = deriv400 + deriv220 + deriv202;
1811 hValue(1, 1) = deriv220 + deriv040 + deriv022;
1812 hValue(2, 2) = deriv202 + deriv022 + deriv004;
1813 hValue(0, 1) = deriv310 + deriv130 + deriv112;
1814 hValue(0, 2) = deriv301 + deriv121 + deriv103;
1815 hValue(1, 2) = deriv211 + deriv031 + deriv013;
1816 hValue(1, 0) = hValue(0, 1);
1817 hValue(2, 0) = hValue(0, 2);
1818 hValue(2, 1) = hValue(1, 2);
1819
1820 value(0, 0) = gValue(0);
1821 value(1, 0) = gValue(1);
1822 value(2, 0) = gValue(2);
1823 value(0, 1) = hValue(0, 0);
1824 value(1, 1) = hValue(1, 0);
1825 value(2, 1) = hValue(2, 0);
1826 value(0, 2) = hValue(0, 1);
1827 value(1, 2) = hValue(1, 1);
1828 value(2, 2) = hValue(2, 1);
1829 value(0, 3) = hValue(0, 2);
1830 value(1, 3) = hValue(1, 2);
1831 value(2, 3) = hValue(2, 2);
1832
1833 return value;
1834 }
1835
kineticEnergyDensityG(Matrix<qreal,3,1> xyz)1836 qreal QTAIMWavefunctionEvaluator::kineticEnergyDensityG(Matrix<qreal, 3, 1> xyz)
1837 {
1838
1839 qreal value;
1840
1841 const qreal zero = 0.0;
1842 const qreal one = 1.0;
1843
1844 m_cdg000.setZero();
1845 m_cdg100.setZero();
1846 m_cdg010.setZero();
1847 m_cdg001.setZero();
1848 for (qint64 p = 0; p < m_nprim; ++p) {
1849 qreal xx0 = xyz(0) - m_X0(p);
1850 qreal yy0 = xyz(1) - m_Y0(p);
1851 qreal zz0 = xyz(2) - m_Z0(p);
1852
1853 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
1854
1855 if (b0arg > m_cutoff) {
1856 qint64 aax0 = 1;
1857 qint64 aay0 = 1;
1858 qint64 aaz0 = 1;
1859 qint64 aax1 = m_xamom(p);
1860 qint64 aay1 = m_yamom(p);
1861 qint64 aaz1 = m_zamom(p);
1862
1863 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
1864 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
1865 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
1866
1867 qreal ax1;
1868 qreal ay1;
1869 qreal az1;
1870 if (m_xamom(p) < 1) {
1871 ax1 = zero;
1872 } else if (m_xamom(p) == 1) {
1873 ax1 = one;
1874 } else {
1875 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
1876 }
1877
1878 if (m_yamom(p) < 1) {
1879 ay1 = zero;
1880 } else if (m_yamom(p) == 1) {
1881 ay1 = one;
1882 } else {
1883 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
1884 }
1885
1886 if (m_zamom(p) < 1) {
1887 az1 = zero;
1888 } else if (m_zamom(p) == 1) {
1889 az1 = one;
1890 } else {
1891 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
1892 }
1893
1894 qreal b0 = exp(b0arg);
1895
1896 qreal bx1 = -2 * m_alpha(p) * xx0;
1897 qreal by1 = -2 * m_alpha(p) * yy0;
1898 qreal bz1 = -2 * m_alpha(p) * zz0;
1899
1900 qreal dg000 = ax0 * ay0 * az0 * b0;
1901 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
1902 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
1903 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
1904
1905 for (qint64 m = 0; m < m_nmo; ++m) {
1906 m_cdg000(m) += m_coef(m, p) * dg000;
1907 m_cdg100(m) += m_coef(m, p) * dg100;
1908 m_cdg010(m) += m_coef(m, p) * dg010;
1909 m_cdg001(m) += m_coef(m, p) * dg001;
1910 }
1911 }
1912 }
1913
1914 value = zero;
1915 for (qint64 m = 0; m < m_nmo; ++m) {
1916 value +=
1917 (0.5) * (m_occno(m) * (ipow(m_cdg100(m), 2) + ipow(m_cdg010(m), 2) +
1918 ipow(m_cdg001(m), 2)));
1919 }
1920
1921 return value;
1922 }
1923
kineticEnergyDensityK(const Matrix<qreal,3,1> xyz)1924 qreal QTAIMWavefunctionEvaluator::kineticEnergyDensityK(
1925 const Matrix<qreal, 3, 1> xyz)
1926 {
1927
1928 qreal value;
1929
1930 const qreal zero = 0.0;
1931 const qreal one = 1.0;
1932
1933 m_cdg000.setZero();
1934 m_cdg200.setZero();
1935 m_cdg020.setZero();
1936 m_cdg002.setZero();
1937 for (qint64 p = 0; p < m_nprim; ++p) {
1938 qreal xx0 = xyz(0) - m_X0(p);
1939 qreal yy0 = xyz(1) - m_Y0(p);
1940 qreal zz0 = xyz(2) - m_Z0(p);
1941
1942 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
1943
1944 if (b0arg > m_cutoff) {
1945 qint64 aax0 = 1;
1946 qint64 aay0 = 1;
1947 qint64 aaz0 = 1;
1948 qint64 aax1 = m_xamom(p);
1949 qint64 aay1 = m_yamom(p);
1950 qint64 aaz1 = m_zamom(p);
1951 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
1952 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
1953 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
1954
1955 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
1956 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
1957 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
1958
1959 qreal ax1;
1960 qreal ay1;
1961 qreal az1;
1962 if (m_xamom(p) < 1) {
1963 ax1 = zero;
1964 } else if (m_xamom(p) == 1) {
1965 ax1 = one;
1966 } else {
1967 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
1968 }
1969
1970 if (m_yamom(p) < 1) {
1971 ay1 = zero;
1972 } else if (m_yamom(p) == 1) {
1973 ay1 = one;
1974 } else {
1975 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
1976 }
1977
1978 if (m_zamom(p) < 1) {
1979 az1 = zero;
1980 } else if (m_zamom(p) == 1) {
1981 az1 = one;
1982 } else {
1983 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
1984 }
1985
1986 qreal ax2;
1987 qreal ay2;
1988 qreal az2;
1989 if (m_xamom(p) < 2) {
1990 ax2 = zero;
1991 } else if (m_xamom(p) == 2) {
1992 ax2 = one;
1993 } else {
1994 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
1995 }
1996
1997 if (m_yamom(p) < 2) {
1998 ay2 = zero;
1999 } else if (m_yamom(p) == 2) {
2000 ay2 = one;
2001 } else {
2002 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
2003 }
2004
2005 if (m_zamom(p) < 2) {
2006 az2 = zero;
2007 } else if (m_zamom(p) == 2) {
2008 az2 = one;
2009 } else {
2010 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
2011 }
2012
2013 qreal b0 = exp(b0arg);
2014
2015 qreal bx1 = -2 * m_alpha(p) * xx0;
2016 qreal by1 = -2 * m_alpha(p) * yy0;
2017 qreal bz1 = -2 * m_alpha(p) * zz0;
2018 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
2019 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
2020 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
2021
2022 qreal dg000 = ax0 * ay0 * az0 * b0;
2023 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
2024 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
2025 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
2026
2027 for (qint64 m = 0; m < m_nmo; ++m) {
2028 m_cdg000(m) += m_coef(m, p) * dg000;
2029 m_cdg200(m) += m_coef(m, p) * dg200;
2030 m_cdg020(m) += m_coef(m, p) * dg020;
2031 m_cdg002(m) += m_coef(m, p) * dg002;
2032 }
2033 }
2034 }
2035
2036 value = 0.0;
2037 for (qint64 m = 0; m < m_nmo; ++m) {
2038 value +=
2039 (0.25) * (m_occno(m) *
2040 (2 * m_cdg000(m) * (m_cdg200(m) + m_cdg020(m) + m_cdg002(m))));
2041 }
2042
2043 return value;
2044 }
2045
quantumStressTensor(const Matrix<qreal,3,1> xyz)2046 const Matrix<qreal, 3, 3> QTAIMWavefunctionEvaluator::quantumStressTensor(
2047 const Matrix<qreal, 3, 1> xyz)
2048 {
2049
2050 Matrix<qreal, 3, 3> value;
2051
2052 const qreal zero = 0.0;
2053 const qreal one = 1.0;
2054
2055 m_cdg000.setZero();
2056 m_cdg100.setZero();
2057 m_cdg010.setZero();
2058 m_cdg001.setZero();
2059 m_cdg200.setZero();
2060 m_cdg020.setZero();
2061 m_cdg002.setZero();
2062 m_cdg110.setZero();
2063 m_cdg101.setZero();
2064 m_cdg011.setZero();
2065 for (qint64 p = 0; p < m_nprim; ++p) {
2066 qreal xx0 = xyz(0) - m_X0(p);
2067 qreal yy0 = xyz(1) - m_Y0(p);
2068 qreal zz0 = xyz(2) - m_Z0(p);
2069
2070 qreal b0arg = -m_alpha(p) * (xx0 * xx0 + yy0 * yy0 + zz0 * zz0);
2071
2072 if (b0arg > m_cutoff) {
2073 qint64 aax0 = 1;
2074 qint64 aay0 = 1;
2075 qint64 aaz0 = 1;
2076 qint64 aax1 = m_xamom(p);
2077 qint64 aay1 = m_yamom(p);
2078 qint64 aaz1 = m_zamom(p);
2079 qint64 aax2 = m_xamom(p) * (m_xamom(p) - 1);
2080 qint64 aay2 = m_yamom(p) * (m_yamom(p) - 1);
2081 qint64 aaz2 = m_zamom(p) * (m_zamom(p) - 1);
2082
2083 qreal ax0 = aax0 * ipow(xx0, m_xamom(p));
2084 qreal ay0 = aay0 * ipow(yy0, m_yamom(p));
2085 qreal az0 = aaz0 * ipow(zz0, m_zamom(p));
2086
2087 qreal ax1;
2088 qreal ay1;
2089 qreal az1;
2090 if (m_xamom(p) < 1) {
2091 ax1 = zero;
2092 } else if (m_xamom(p) == 1) {
2093 ax1 = one;
2094 } else {
2095 ax1 = aax1 * ipow(xx0, m_xamom(p) - 1);
2096 }
2097
2098 if (m_yamom(p) < 1) {
2099 ay1 = zero;
2100 } else if (m_yamom(p) == 1) {
2101 ay1 = one;
2102 } else {
2103 ay1 = aay1 * ipow(yy0, m_yamom(p) - 1);
2104 }
2105
2106 if (m_zamom(p) < 1) {
2107 az1 = zero;
2108 } else if (m_zamom(p) == 1) {
2109 az1 = one;
2110 } else {
2111 az1 = aaz1 * ipow(zz0, m_zamom(p) - 1);
2112 }
2113
2114 qreal ax2;
2115 qreal ay2;
2116 qreal az2;
2117 if (m_xamom(p) < 2) {
2118 ax2 = zero;
2119 } else if (m_xamom(p) == 2) {
2120 ax2 = one;
2121 } else {
2122 ax2 = aax2 * ipow(xx0, m_xamom(p) - 2);
2123 }
2124
2125 if (m_yamom(p) < 2) {
2126 ay2 = zero;
2127 } else if (m_yamom(p) == 2) {
2128 ay2 = one;
2129 } else {
2130 ay2 = aay2 * ipow(yy0, m_yamom(p) - 2);
2131 }
2132
2133 if (m_zamom(p) < 2) {
2134 az2 = zero;
2135 } else if (m_zamom(p) == 2) {
2136 az2 = one;
2137 } else {
2138 az2 = aaz2 * ipow(zz0, m_zamom(p) - 2);
2139 }
2140
2141 qreal b0 = exp(b0arg);
2142
2143 qreal bx1 = -2 * m_alpha(p) * xx0;
2144 qreal by1 = -2 * m_alpha(p) * yy0;
2145 qreal bz1 = -2 * m_alpha(p) * zz0;
2146 qreal bx2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(xx0, 2));
2147 qreal by2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(yy0, 2));
2148 qreal bz2 = -2 * m_alpha(p) + 4 * (ipow(m_alpha(p), 2) * ipow(zz0, 2));
2149
2150 qreal dg000 = ax0 * ay0 * az0 * b0;
2151 qreal dg100 = ay0 * az0 * b0 * (ax1 + ax0 * bx1);
2152 qreal dg010 = ax0 * az0 * b0 * (ay1 + ay0 * by1);
2153 qreal dg001 = ax0 * ay0 * b0 * (az1 + az0 * bz1);
2154 qreal dg200 = ay0 * az0 * b0 * (ax2 + 2 * ax1 * bx1 + ax0 * bx2);
2155 qreal dg020 = ax0 * az0 * b0 * (ay2 + 2 * ay1 * by1 + ay0 * by2);
2156 qreal dg002 = ax0 * ay0 * b0 * (az2 + 2 * az1 * bz1 + az0 * bz2);
2157 qreal dg110 = az0 * b0 * (ax1 + ax0 * bx1) * (ay1 + ay0 * by1);
2158 qreal dg101 = ay0 * b0 * (ax1 + ax0 * bx1) * (az1 + az0 * bz1);
2159 qreal dg011 = ax0 * b0 * (ay1 + ay0 * by1) * (az1 + az0 * bz1);
2160
2161 for (qint64 m = 0; m < m_nmo; ++m) {
2162 m_cdg000(m) += m_coef(m, p) * dg000;
2163 m_cdg100(m) += m_coef(m, p) * dg100;
2164 m_cdg010(m) += m_coef(m, p) * dg010;
2165 m_cdg001(m) += m_coef(m, p) * dg001;
2166 m_cdg200(m) += m_coef(m, p) * dg200;
2167 m_cdg020(m) += m_coef(m, p) * dg020;
2168 m_cdg002(m) += m_coef(m, p) * dg002;
2169 m_cdg110(m) += m_coef(m, p) * dg110;
2170 m_cdg101(m) += m_coef(m, p) * dg101;
2171 m_cdg011(m) += m_coef(m, p) * dg011;
2172 }
2173 }
2174 }
2175
2176 value.setZero();
2177 for (qint64 m = 0; m < m_nmo; ++m) {
2178 value(0, 0) +=
2179 (m_occno(m) * (2 * m_cdg000(m) * m_cdg200(m) - 2 * ipow(m_cdg100(m), 2)));
2180 value(0, 1) += (m_occno(m) * (2 * m_cdg000(m) * m_cdg110(m) -
2181 2 * m_cdg100(m) * m_cdg010(m)));
2182 value(0, 2) += (m_occno(m) * (2 * m_cdg000(m) * m_cdg101(m) -
2183 2 * m_cdg100(m) * m_cdg001(m)));
2184 value(1, 1) +=
2185 (m_occno(m) * (2 * m_cdg000(m) * m_cdg020(m) - 2 * ipow(m_cdg010(m), 2)));
2186 value(1, 2) += (m_occno(m) * (2 * m_cdg000(m) * m_cdg011(m) -
2187 2 * m_cdg010(m) * m_cdg001(m)));
2188 value(2, 2) +=
2189 (m_occno(m) * (2 * m_cdg000(m) * m_cdg002(m) - 2 * ipow(m_cdg001(m), 2)));
2190 }
2191 value(1, 0) = value(0, 1);
2192 value(2, 0) = value(0, 2);
2193 value(2, 1) = value(1, 2);
2194
2195 return 0.25 * value;
2196 }
2197
2198 } // namespace QtPlugins
2199 } // namespace Avogadro
2200