1 /******************************************************************************
2
3 This source file is part of the Avogadro project.
4
5 Copyright 2010 Eric C. Brown
6 Copyright 2013 Kitware, Inc.
7
8 This source code is released under the New BSD License, (the "License").
9
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
15
16 ******************************************************************************/
17
18 #include "qtaimextension.h"
19
20 #include <avogadro/qtgui/molecule.h>
21
22 #include <QAction>
23
24 #include <QDebug>
25 #include <QDir>
26 #include <QFileDialog>
27 #include <QList>
28 #include <QPair>
29 #include <QString>
30 #include <QVector3D>
31
32 #include <QThread>
33
34 #include "qtaimcriticalpointlocator.h"
35 #include "qtaimcubature.h"
36 #include "qtaimwavefunction.h"
37 #include "qtaimwavefunctionevaluator.h"
38
39 #include <QTime>
40
41 using namespace std;
42 using namespace Eigen;
43
44 namespace Avogadro {
45 namespace QtPlugins {
46
47 enum QTAIMExtensionIndex
48 {
49 FirstAction = 0,
50 SecondAction,
51 ThirdAction
52 };
53
QTAIMExtension(QObject * aParent)54 QTAIMExtension::QTAIMExtension(QObject* aParent)
55 : QtGui::ExtensionPlugin(aParent)
56 {
57 // create an action for our first action
58 QAction* action = new QAction(this);
59 action->setText(tr("Molecular Graph..."));
60 m_actions.append(action);
61 action->setData(FirstAction);
62 connect(action, SIGNAL(triggered()), SLOT(triggered()));
63
64 // create an action for our second action
65 action = new QAction(this);
66 action->setText(tr("Molecular Graph with Lone Pairs..."));
67 m_actions.append(action);
68 action->setData(SecondAction);
69 connect(action, SIGNAL(triggered()), SLOT(triggered()));
70
71 // create an action for our third action
72 action = new QAction(this);
73 action->setText(tr("Atomic Charge..."));
74 m_actions.append(action);
75 action->setData(ThirdAction);
76 connect(action, SIGNAL(triggered()), SLOT(triggered()));
77 }
78
~QTAIMExtension()79 QTAIMExtension::~QTAIMExtension()
80 {
81 }
82
actions() const83 QList<QAction*> QTAIMExtension::actions() const
84 {
85 return m_actions;
86 }
87
menuPath(QAction *) const88 QStringList QTAIMExtension::menuPath(QAction*) const
89 {
90 return QStringList() << tr("&Analysis") << tr("QTAIM");
91 }
92
setMolecule(QtGui::Molecule * molecule)93 void QTAIMExtension::setMolecule(QtGui::Molecule* molecule)
94 {
95 m_molecule = molecule;
96 }
97
triggered()98 void QTAIMExtension::triggered()
99 {
100 QAction* action = qobject_cast<QAction*>(sender());
101 if (!action)
102 return;
103
104 bool wavefunctionAlreadyLoaded;
105
106 if (m_molecule->property("QTAIMComment").isValid()) {
107 wavefunctionAlreadyLoaded = true;
108 } else {
109 wavefunctionAlreadyLoaded = false;
110 }
111
112 int i = action->data().toInt();
113
114 QTime timer;
115 timer.start();
116
117 QString fileName;
118 if (wavefunctionAlreadyLoaded) {
119 // do nothing
120 } else {
121 fileName = QFileDialog::getOpenFileName(
122 new QWidget, tr("Open WFN File"), QDir::homePath(),
123 tr("WFN files (*.wfn);;All files (*.*)"));
124
125 if (fileName.isNull()) {
126 qDebug() << "No such file.";
127 return;
128 }
129 }
130
131 // Instantiate a Wavefunction
132 bool success;
133 QTAIMWavefunction wfn;
134 if (wavefunctionAlreadyLoaded) {
135 success = wfn.initializeWithMoleculeProperties(m_molecule);
136 } else {
137 success = wfn.initializeWithWFNFile(fileName);
138 }
139
140 if (!success) {
141 if (wavefunctionAlreadyLoaded) {
142 qDebug() << "Error initializing wavefunction.";
143 } else {
144 qDebug() << "Error reading WFN file.";
145 }
146 return;
147 }
148
149 QtGui::Molecule::MoleculeChanges changes;
150 if (m_molecule->atomCount() > 0)
151 changes |= QtGui::Molecule::Atoms | QtGui::Molecule::Removed;
152 if (m_molecule->bondCount() > 0)
153 changes |= QtGui::Molecule::Bonds | QtGui::Molecule::Removed;
154
155 m_molecule->clearAtoms();
156 m_molecule->emitChanged(static_cast<unsigned int>(changes));
157
158 // Instantiate an Evaluator
159 QTAIMWavefunctionEvaluator eval(wfn);
160
161 switch (i) {
162 case FirstAction: // Molecular Graph
163 {
164 // Instantiate a Critical Point Locator
165 QTAIMCriticalPointLocator cpl(wfn);
166
167 // Locate the Nuclear Critical Points
168 cpl.locateNuclearCriticalPoints();
169
170 // QLists of results
171 QList<qint64> nucChargeList = wfn.nuclearChargesList();
172 QList<QVector3D> ncpList = cpl.nuclearCriticalPoints();
173
174 QVariantList xNCPsVariantList;
175 QVariantList yNCPsVariantList;
176 QVariantList zNCPsVariantList;
177 QVariantList nuclearChargesVariantList;
178
179 const qreal convertBohrToAngstrom = 0.529177249;
180
181 // Nuclear Critical Points
182 for (qint64 n = 0; n < ncpList.length(); ++n) {
183 QVector3D thisNuclearCriticalPoint = ncpList.at(n);
184
185 qreal x = thisNuclearCriticalPoint.x() * convertBohrToAngstrom;
186 qreal y = thisNuclearCriticalPoint.y() * convertBohrToAngstrom;
187 qreal z = thisNuclearCriticalPoint.z() * convertBohrToAngstrom;
188
189 xNCPsVariantList.append(x);
190 yNCPsVariantList.append(y);
191 zNCPsVariantList.append(z);
192 nuclearChargesVariantList.append(wfn.nuclearCharge(n));
193 }
194
195 m_molecule->setProperty("QTAIMXNuclearCriticalPoints", xNCPsVariantList);
196 m_molecule->setProperty("QTAIMYNuclearCriticalPoints", yNCPsVariantList);
197 m_molecule->setProperty("QTAIMZNuclearCriticalPoints", zNCPsVariantList);
198 m_molecule->setProperty("QTAIMNuclearCharges", nuclearChargesVariantList);
199
200 // Nuclei stored as Atoms
201 for (qint64 n = 0; n < wfn.numberOfNuclei(); ++n) {
202 qreal x = wfn.xNuclearCoordinate(n) * convertBohrToAngstrom;
203 qreal y = wfn.yNuclearCoordinate(n) * convertBohrToAngstrom;
204 qreal z = wfn.zNuclearCoordinate(n) * convertBohrToAngstrom;
205
206 int Z = (int)wfn.nuclearCharge(n);
207
208 m_molecule->addAtom(static_cast<unsigned char>(Z))
209 .setPosition3d(Vector3(static_cast<Real>(x), static_cast<Real>(y),
210 static_cast<Real>(z)));
211 }
212
213 if (m_molecule->atomCount() > 0) {
214 m_molecule->emitChanged(QtGui::Molecule::Atoms |
215 QtGui::Molecule::Added);
216 }
217
218 // Locate the Bond Critical Points and Trace Bond Paths
219 cpl.locateBondCriticalPoints();
220
221 // BCP and Bond Path Results
222 QList<QVector3D> bcpList = cpl.bondCriticalPoints();
223 QList<QList<QVector3D>> bondPathList = cpl.bondPaths();
224 QList<QPair<qint64, qint64>> bondedAtomsList = cpl.bondedAtoms();
225 QList<qreal> laplacianAtBondCriticalPoints =
226 cpl.laplacianAtBondCriticalPoints();
227 QList<qreal> ellipticityAtBondCriticalPoints =
228 cpl.ellipticityAtBondCriticalPoints();
229
230 QVariantList xBCPsVariantList;
231 QVariantList yBCPsVariantList;
232 QVariantList zBCPsVariantList;
233 QVariantList firstNCPIndexVariantList;
234 QVariantList secondNCPIndexVariantList;
235 QVariantList laplacianAtBondCriticalPointsVariantList;
236 QVariantList ellipticityAtBondCriticalPointsVariantList;
237
238 QVariantList bondPathSegmentStartIndexVariantList;
239 QVariantList bondPathSegmentEndIndexVariantList;
240 QVariantList xBondPathsVariantList;
241 QVariantList yBondPathsVariantList;
242 QVariantList zBondPathsVariantList;
243
244 // Connectivity stored as Bonds
245
246 qint64 bpCtr = 0;
247 qint64 numAtoms = static_cast<qint64>(m_molecule->atomCount());
248
249 for (qint64 atom0 = 0; atom0 < numAtoms - 1; ++atom0) {
250 for (qint64 atom1 = atom0 + 1; atom1 < numAtoms; ++atom1) {
251
252 bool areBonded = false;
253
254 for (qint64 bondPair = 0; bondPair < bondedAtomsList.length();
255 ++bondPair) {
256 if (atom0 == bondedAtomsList.at(bondPair).first &&
257 atom1 == bondedAtomsList.at(bondPair).second) {
258 areBonded = true;
259
260 if (areBonded) {
261
262 if ((wfn.nuclearCharge(atom0) == 1 ||
263 wfn.nuclearCharge(atom1) == 1) &&
264 laplacianAtBondCriticalPoints.at(bondPair) > 0.0) {
265 // do not draw Bond because it looks like hydrogen bond
266 } else {
267 m_molecule->addBond(
268 m_molecule->atom(static_cast<Index>(atom0)),
269 m_molecule->atom(static_cast<Index>(atom1)));
270 // bond->setAromaticity(isAromatic);
271 // bond->setOrder( (int) order);
272 }
273
274 qreal x = bcpList.at(bondPair).x() * convertBohrToAngstrom;
275 qreal y = bcpList.at(bondPair).y() * convertBohrToAngstrom;
276 qreal z = bcpList.at(bondPair).z() * convertBohrToAngstrom;
277
278 xBCPsVariantList.append(x);
279 yBCPsVariantList.append(y);
280 zBCPsVariantList.append(z);
281
282 firstNCPIndexVariantList.append(atom0);
283 secondNCPIndexVariantList.append(atom1);
284
285 laplacianAtBondCriticalPointsVariantList.append(
286 laplacianAtBondCriticalPoints.at(bondPair));
287 ellipticityAtBondCriticalPointsVariantList.append(
288 ellipticityAtBondCriticalPoints.at(bondPair));
289
290 bondPathSegmentStartIndexVariantList.append(bpCtr);
291 for (qint64 j = 0; j < bondPathList.at(bondPair).length();
292 ++j) {
293 x =
294 bondPathList.at(bondPair).at(j).x() * convertBohrToAngstrom;
295 y =
296 bondPathList.at(bondPair).at(j).y() * convertBohrToAngstrom;
297 z =
298 bondPathList.at(bondPair).at(j).z() * convertBohrToAngstrom;
299
300 xBondPathsVariantList.append(x);
301 yBondPathsVariantList.append(y);
302 zBondPathsVariantList.append(z);
303
304 bpCtr++;
305 }
306 bondPathSegmentEndIndexVariantList.append(bpCtr);
307 }
308 }
309 } // bond pairs
310 } // atom1
311 } // atom 0
312
313 m_molecule->setProperty("QTAIMXBondCriticalPoints", xBCPsVariantList);
314 m_molecule->setProperty("QTAIMYBondCriticalPoints", yBCPsVariantList);
315 m_molecule->setProperty("QTAIMZBondCriticalPoints", zBCPsVariantList);
316 m_molecule->setProperty("QTAIMFirstNCPIndexVariantList",
317 firstNCPIndexVariantList);
318 m_molecule->setProperty("QTAIMSecondNCPIndexVariantList",
319 secondNCPIndexVariantList);
320 m_molecule->setProperty("QTAIMLaplacianAtBondCriticalPoints",
321 laplacianAtBondCriticalPointsVariantList);
322 m_molecule->setProperty("QTAIMEllipticityAtBondCriticalPoints",
323 ellipticityAtBondCriticalPointsVariantList);
324
325 m_molecule->setProperty("QTAIMBondPathSegmentStartIndex",
326 bondPathSegmentStartIndexVariantList);
327 m_molecule->setProperty("QTAIMBondPathSegmentEndIndex",
328 bondPathSegmentEndIndexVariantList);
329 m_molecule->setProperty("QTAIMXBondPaths", xBondPathsVariantList);
330 m_molecule->setProperty("QTAIMYBondPaths", yBondPathsVariantList);
331 m_molecule->setProperty("QTAIMZBondPaths", zBondPathsVariantList);
332
333 if (m_molecule->bondCount()) {
334 m_molecule->emitChanged(QtGui::Molecule::Bonds |
335 QtGui::Molecule::Added);
336 }
337 } break;
338 case SecondAction: // Molecular Graph with Lone Pairs
339 {
340 // Instantiate a Critical Point Locator
341 QTAIMCriticalPointLocator cpl(wfn);
342
343 // Locate the Nuclear Critical Points
344 cpl.locateNuclearCriticalPoints();
345
346 // QLists of results
347 QList<qint64> nucChargeList = wfn.nuclearChargesList();
348 QList<QVector3D> ncpList = cpl.nuclearCriticalPoints();
349
350 QVariantList xNCPsVariantList;
351 QVariantList yNCPsVariantList;
352 QVariantList zNCPsVariantList;
353 QVariantList nuclearChargesVariantList;
354
355 const qreal convertBohrToAngstrom = 0.529177249;
356
357 // Nuclear Critical Points
358 for (qint64 n = 0; n < ncpList.length(); ++n) {
359 QVector3D thisNuclearCriticalPoint = ncpList.at(n);
360
361 qreal x = thisNuclearCriticalPoint.x() * convertBohrToAngstrom;
362 qreal y = thisNuclearCriticalPoint.y() * convertBohrToAngstrom;
363 qreal z = thisNuclearCriticalPoint.z() * convertBohrToAngstrom;
364
365 xNCPsVariantList.append(x);
366 yNCPsVariantList.append(y);
367 zNCPsVariantList.append(z);
368 nuclearChargesVariantList.append(wfn.nuclearCharge(n));
369 }
370
371 m_molecule->setProperty("QTAIMXNuclearCriticalPoints", xNCPsVariantList);
372 m_molecule->setProperty("QTAIMYNuclearCriticalPoints", yNCPsVariantList);
373 m_molecule->setProperty("QTAIMZNuclearCriticalPoints", zNCPsVariantList);
374 m_molecule->setProperty("QTAIMNuclearCharges", nuclearChargesVariantList);
375
376 // Nuclei stored as Atoms
377 for (qint64 n = 0; n < wfn.numberOfNuclei(); ++n) {
378 qreal x = wfn.xNuclearCoordinate(n) * convertBohrToAngstrom;
379 qreal y = wfn.yNuclearCoordinate(n) * convertBohrToAngstrom;
380 qreal z = wfn.zNuclearCoordinate(n) * convertBohrToAngstrom;
381
382 int Z = (int)wfn.nuclearCharge(n);
383
384 m_molecule->addAtom(static_cast<unsigned char>(Z))
385 .setPosition3d(Vector3(static_cast<Real>(x), static_cast<Real>(y),
386 static_cast<Real>(z)));
387 }
388
389 if (m_molecule->atomCount() > 0) {
390 m_molecule->emitChanged(QtGui::Molecule::Atoms |
391 QtGui::Molecule::Added);
392 }
393
394 // Locate the Bond Critical Points and Trace Bond Paths
395 cpl.locateBondCriticalPoints();
396
397 // BCP and Bond Path Results
398 QList<QVector3D> bcpList = cpl.bondCriticalPoints();
399 QList<QList<QVector3D>> bondPathList = cpl.bondPaths();
400 QList<QPair<qint64, qint64>> bondedAtomsList = cpl.bondedAtoms();
401 QList<qreal> laplacianAtBondCriticalPoints =
402 cpl.laplacianAtBondCriticalPoints();
403 QList<qreal> ellipticityAtBondCriticalPoints =
404 cpl.ellipticityAtBondCriticalPoints();
405
406 QVariantList xBCPsVariantList;
407 QVariantList yBCPsVariantList;
408 QVariantList zBCPsVariantList;
409 QVariantList firstNCPIndexVariantList;
410 QVariantList secondNCPIndexVariantList;
411 QVariantList laplacianAtBondCriticalPointsVariantList;
412 QVariantList ellipticityAtBondCriticalPointsVariantList;
413
414 QVariantList bondPathSegmentStartIndexVariantList;
415 QVariantList bondPathSegmentEndIndexVariantList;
416 QVariantList xBondPathsVariantList;
417 QVariantList yBondPathsVariantList;
418 QVariantList zBondPathsVariantList;
419
420 // Connectivity stored as Bonds
421
422 qint64 bpCtr = 0;
423 qint64 numAtoms = static_cast<qint64>(m_molecule->atomCount());
424
425 for (qint64 atom0 = 0; atom0 < numAtoms - 1; ++atom0) {
426 for (qint64 atom1 = atom0 + 1; atom1 < numAtoms; ++atom1) {
427
428 bool areBonded = false;
429
430 for (qint64 bondPair = 0; bondPair < bondedAtomsList.length();
431 ++bondPair) {
432 if (atom0 == bondedAtomsList.at(bondPair).first &&
433 atom1 == bondedAtomsList.at(bondPair).second) {
434 areBonded = true;
435
436 if (areBonded) {
437
438 if ((wfn.nuclearCharge(atom0) == 1 ||
439 wfn.nuclearCharge(atom1) == 1) &&
440 laplacianAtBondCriticalPoints.at(bondPair) > 0.0) {
441 // do not draw Bond because it looks like hydrogen bond
442 } else {
443 m_molecule->addBond(
444 m_molecule->atom(static_cast<Index>(atom0)),
445 m_molecule->atom(static_cast<Index>(atom1)));
446 // bond->setAromaticity(isAromatic);
447 // bond->setOrder( (int) order);
448 }
449
450 qreal x = bcpList.at(bondPair).x() * convertBohrToAngstrom;
451 qreal y = bcpList.at(bondPair).y() * convertBohrToAngstrom;
452 qreal z = bcpList.at(bondPair).z() * convertBohrToAngstrom;
453
454 xBCPsVariantList.append(x);
455 yBCPsVariantList.append(y);
456 zBCPsVariantList.append(z);
457
458 firstNCPIndexVariantList.append(atom0);
459 secondNCPIndexVariantList.append(atom1);
460
461 laplacianAtBondCriticalPointsVariantList.append(
462 laplacianAtBondCriticalPoints.at(bondPair));
463 ellipticityAtBondCriticalPointsVariantList.append(
464 ellipticityAtBondCriticalPoints.at(bondPair));
465
466 bondPathSegmentStartIndexVariantList.append(bpCtr);
467 for (qint64 j = 0; j < bondPathList.at(bondPair).length();
468 ++j) {
469 x =
470 bondPathList.at(bondPair).at(j).x() * convertBohrToAngstrom;
471 y =
472 bondPathList.at(bondPair).at(j).y() * convertBohrToAngstrom;
473 z =
474 bondPathList.at(bondPair).at(j).z() * convertBohrToAngstrom;
475
476 xBondPathsVariantList.append(x);
477 yBondPathsVariantList.append(y);
478 zBondPathsVariantList.append(z);
479
480 bpCtr++;
481 }
482 bondPathSegmentEndIndexVariantList.append(bpCtr);
483 }
484 }
485 } // bond pairs
486 } // atom1
487 } // atom 0
488
489 m_molecule->setProperty("QTAIMXBondCriticalPoints", xBCPsVariantList);
490 m_molecule->setProperty("QTAIMYBondCriticalPoints", yBCPsVariantList);
491 m_molecule->setProperty("QTAIMZBondCriticalPoints", zBCPsVariantList);
492 m_molecule->setProperty("QTAIMFirstNCPIndexVariantList",
493 firstNCPIndexVariantList);
494 m_molecule->setProperty("QTAIMSecondNCPIndexVariantList",
495 secondNCPIndexVariantList);
496 m_molecule->setProperty("QTAIMLaplacianAtBondCriticalPoints",
497 laplacianAtBondCriticalPointsVariantList);
498 m_molecule->setProperty("QTAIMEllipticityAtBondCriticalPoints",
499 ellipticityAtBondCriticalPointsVariantList);
500
501 m_molecule->setProperty("QTAIMBondPathSegmentStartIndex",
502 bondPathSegmentStartIndexVariantList);
503 m_molecule->setProperty("QTAIMBondPathSegmentEndIndex",
504 bondPathSegmentEndIndexVariantList);
505 m_molecule->setProperty("QTAIMXBondPaths", xBondPathsVariantList);
506 m_molecule->setProperty("QTAIMYBondPaths", yBondPathsVariantList);
507 m_molecule->setProperty("QTAIMZBondPaths", zBondPathsVariantList);
508
509 if (m_molecule->bondCount()) {
510 m_molecule->emitChanged(QtGui::Molecule::Bonds |
511 QtGui::Molecule::Added);
512 }
513
514 // Locate Electron Density Sources / Lone Pairs
515
516 cpl.locateElectronDensitySources();
517 QList<QVector3D> electronDensitySourcesList =
518 cpl.electronDensitySources();
519
520 QVariantList xElectronDensitySourcesVariantList;
521 QVariantList yElectronDensitySourcesVariantList;
522 QVariantList zElectronDensitySourcesVariantList;
523
524 for (qint64 n = 0; n < electronDensitySourcesList.length(); ++n) {
525 QVector3D thisCriticalPoint = electronDensitySourcesList.at(n);
526
527 qreal x = thisCriticalPoint.x() * convertBohrToAngstrom;
528 qreal y = thisCriticalPoint.y() * convertBohrToAngstrom;
529 qreal z = thisCriticalPoint.z() * convertBohrToAngstrom;
530
531 xElectronDensitySourcesVariantList.append(x);
532 yElectronDensitySourcesVariantList.append(y);
533 zElectronDensitySourcesVariantList.append(z);
534 }
535
536 m_molecule->setProperty("QTAIMXElectronDensitySources",
537 xElectronDensitySourcesVariantList);
538 m_molecule->setProperty("QTAIMYElectronDensitySources",
539 yElectronDensitySourcesVariantList);
540 m_molecule->setProperty("QTAIMZElectronDensitySources",
541 zElectronDensitySourcesVariantList);
542
543 // TODO need some way to indicate that the properties have changed:
544 // m_molecule->update();
545
546 } break;
547 case ThirdAction:
548 // perform third action
549 {
550 // Instantiate a Critical Point Locator
551 QTAIMCriticalPointLocator cpl(wfn);
552
553 // Locate the Nuclear Critical Points
554 cpl.locateNuclearCriticalPoints();
555
556 // QLists of results
557 QList<qint64> nucChargeList = wfn.nuclearChargesList();
558 QList<QVector3D> ncpList = cpl.nuclearCriticalPoints();
559
560 QVariantList xNCPsVariantList;
561 QVariantList yNCPsVariantList;
562 QVariantList zNCPsVariantList;
563 QVariantList nuclearChargesVariantList;
564
565 const qreal convertBohrToAngstrom = 0.529177249;
566
567 // Nuclear Critical Points
568 for (qint64 n = 0; n < ncpList.length(); ++n) {
569 QVector3D thisNuclearCriticalPoint = ncpList.at(n);
570
571 qreal x = thisNuclearCriticalPoint.x() * convertBohrToAngstrom;
572 qreal y = thisNuclearCriticalPoint.y() * convertBohrToAngstrom;
573 qreal z = thisNuclearCriticalPoint.z() * convertBohrToAngstrom;
574
575 xNCPsVariantList.append(x);
576 yNCPsVariantList.append(y);
577 zNCPsVariantList.append(z);
578 nuclearChargesVariantList.append(wfn.nuclearCharge(n));
579 }
580
581 m_molecule->setProperty("QTAIMXNuclearCriticalPoints",
582 xNCPsVariantList);
583 m_molecule->setProperty("QTAIMYNuclearCriticalPoints",
584 yNCPsVariantList);
585 m_molecule->setProperty("QTAIMZNuclearCriticalPoints",
586 zNCPsVariantList);
587 m_molecule->setProperty("QTAIMNuclearCharges",
588 nuclearChargesVariantList);
589
590 // Nuclei stored as Atoms
591 for (qint64 n = 0; n < wfn.numberOfNuclei(); ++n) {
592 qreal x = wfn.xNuclearCoordinate(n) * convertBohrToAngstrom;
593 qreal y = wfn.yNuclearCoordinate(n) * convertBohrToAngstrom;
594 qreal z = wfn.zNuclearCoordinate(n) * convertBohrToAngstrom;
595
596 int Z = (int)wfn.nuclearCharge(n);
597
598 m_molecule->addAtom(static_cast<unsigned char>(Z))
599 .setPosition3d(Vector3(static_cast<Real>(x), static_cast<Real>(y),
600 static_cast<Real>(z)));
601 }
602
603 if (m_molecule->atomCount() > 0) {
604 m_molecule->emitChanged(QtGui::Molecule::Atoms |
605 QtGui::Molecule::Added);
606 }
607
608 // Locate the Bond Critical Points and Trace Bond Paths
609 cpl.locateBondCriticalPoints();
610
611 // BCP and Bond Path Results
612 QList<QVector3D> bcpList = cpl.bondCriticalPoints();
613 QList<QList<QVector3D>> bondPathList = cpl.bondPaths();
614 QList<QPair<qint64, qint64>> bondedAtomsList = cpl.bondedAtoms();
615 QList<qreal> laplacianAtBondCriticalPoints =
616 cpl.laplacianAtBondCriticalPoints();
617 QList<qreal> ellipticityAtBondCriticalPoints =
618 cpl.ellipticityAtBondCriticalPoints();
619
620 QVariantList xBCPsVariantList;
621 QVariantList yBCPsVariantList;
622 QVariantList zBCPsVariantList;
623 QVariantList firstNCPIndexVariantList;
624 QVariantList secondNCPIndexVariantList;
625 QVariantList laplacianAtBondCriticalPointsVariantList;
626 QVariantList ellipticityAtBondCriticalPointsVariantList;
627
628 QVariantList bondPathSegmentStartIndexVariantList;
629 QVariantList bondPathSegmentEndIndexVariantList;
630 QVariantList xBondPathsVariantList;
631 QVariantList yBondPathsVariantList;
632 QVariantList zBondPathsVariantList;
633
634 // Connectivity stored as Bonds
635
636 qint64 bpCtr = 0;
637 qint64 numAtoms = static_cast<qint64>(m_molecule->atomCount());
638
639 for (qint64 atom0 = 0; atom0 < numAtoms - 1; ++atom0) {
640 for (qint64 atom1 = atom0 + 1; atom1 < numAtoms; ++atom1) {
641
642 bool areBonded = false;
643
644 for (qint64 bondPair = 0; bondPair < bondedAtomsList.length();
645 ++bondPair) {
646 if (atom0 == bondedAtomsList.at(bondPair).first &&
647 atom1 == bondedAtomsList.at(bondPair).second) {
648 areBonded = true;
649
650 if (areBonded) {
651
652 if ((wfn.nuclearCharge(atom0) == 1 ||
653 wfn.nuclearCharge(atom1) == 1) &&
654 laplacianAtBondCriticalPoints.at(bondPair) > 0.0) {
655 // do not draw Bond because it looks like hydrogen bond
656 } else {
657 m_molecule->addBond(
658 m_molecule->atom(static_cast<Index>(atom0)),
659 m_molecule->atom(static_cast<Index>(atom1)));
660 // bond->setAromaticity(isAromatic);
661 // bond->setOrder( (int) order);
662 }
663
664 qreal x = bcpList.at(bondPair).x() * convertBohrToAngstrom;
665 qreal y = bcpList.at(bondPair).y() * convertBohrToAngstrom;
666 qreal z = bcpList.at(bondPair).z() * convertBohrToAngstrom;
667
668 xBCPsVariantList.append(x);
669 yBCPsVariantList.append(y);
670 zBCPsVariantList.append(z);
671
672 firstNCPIndexVariantList.append(atom0);
673 secondNCPIndexVariantList.append(atom1);
674
675 laplacianAtBondCriticalPointsVariantList.append(
676 laplacianAtBondCriticalPoints.at(bondPair));
677 ellipticityAtBondCriticalPointsVariantList.append(
678 ellipticityAtBondCriticalPoints.at(bondPair));
679
680 bondPathSegmentStartIndexVariantList.append(bpCtr);
681 for (qint64 j = 0; j < bondPathList.at(bondPair).length();
682 ++j) {
683 x = bondPathList.at(bondPair).at(j).x() *
684 convertBohrToAngstrom;
685 y = bondPathList.at(bondPair).at(j).y() *
686 convertBohrToAngstrom;
687 z = bondPathList.at(bondPair).at(j).z() *
688 convertBohrToAngstrom;
689
690 xBondPathsVariantList.append(x);
691 yBondPathsVariantList.append(y);
692 zBondPathsVariantList.append(z);
693
694 bpCtr++;
695 }
696 bondPathSegmentEndIndexVariantList.append(bpCtr);
697 }
698 }
699 } // bond pairs
700 } // atom1
701 } // atom 0
702
703 m_molecule->setProperty("QTAIMXBondCriticalPoints", xBCPsVariantList);
704 m_molecule->setProperty("QTAIMYBondCriticalPoints", yBCPsVariantList);
705 m_molecule->setProperty("QTAIMZBondCriticalPoints", zBCPsVariantList);
706 m_molecule->setProperty("QTAIMFirstNCPIndexVariantList",
707 firstNCPIndexVariantList);
708 m_molecule->setProperty("QTAIMSecondNCPIndexVariantList",
709 secondNCPIndexVariantList);
710 m_molecule->setProperty("QTAIMLaplacianAtBondCriticalPoints",
711 laplacianAtBondCriticalPointsVariantList);
712 m_molecule->setProperty("QTAIMEllipticityAtBondCriticalPoints",
713 ellipticityAtBondCriticalPointsVariantList);
714
715 m_molecule->setProperty("QTAIMBondPathSegmentStartIndex",
716 bondPathSegmentStartIndexVariantList);
717 m_molecule->setProperty("QTAIMBondPathSegmentEndIndex",
718 bondPathSegmentEndIndexVariantList);
719 m_molecule->setProperty("QTAIMXBondPaths", xBondPathsVariantList);
720 m_molecule->setProperty("QTAIMYBondPaths", yBondPathsVariantList);
721 m_molecule->setProperty("QTAIMZBondPaths", zBondPathsVariantList);
722
723 if (m_molecule->bondCount() > 0) {
724 m_molecule->emitChanged(QtGui::Molecule::Bonds |
725 QtGui::Molecule::Added);
726 }
727
728 // Electron Density
729 qint64 mode = 0;
730
731 // All Atomic Basins
732 QList<qint64> basins;
733 for (qint64 j = 0; j < wfn.numberOfNuclei(); ++j) {
734 basins.append(j);
735 }
736
737 QTAIMCubature cub(wfn);
738
739 // QTime time;
740 // time.start();
741 QList<QPair<qreal, qreal>> results = cub.integrate(mode, basins);
742 // qDebug() << "Time Elapsed:" << time.elapsed();
743
744 for (qint64 j = 0; j < results.length(); ++j) {
745 qDebug() << "basin" << j << results.at(j).first
746 << results.at(j).second;
747 }
748
749 // TODO: Set the properties of the atoms.
750 // I don't know why this bombs.
751 for (qint64 j = 0; static_cast<qint64>(m_molecule->atomCount()); ++j) {
752 // Atom *atom=m_molecule->atoms().at(i);
753 // const qreal charge=results.at(i).first;
754 // atom->setPartialCharge( charge );
755 }
756 }
757 break;
758 }
759
760 emit requestActiveTool("Navigator");
761 emit requestActiveDisplayTypes(QStringList() << "QTAIMScenePlugin");
762
763 return;
764 }
765
766 } // end namespace QtPlugins
767 } // end namespace Avogadro
768