00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include <globalsearch/structure.h>
00017 #include <globalsearch/macros.h>
00018
00019 #include <avogadro/primitive.h>
00020 #include <avogadro/molecule.h>
00021 #include <avogadro/atom.h>
00022
00023 #include <openbabel/generic.h>
00024 #include <openbabel/forcefield.h>
00025
00026 #include <QtCore/QFile>
00027 #include <QtCore/QDebug>
00028 #include <QtCore/QTimer>
00029 #include <QtCore/QRegExp>
00030 #include <QtCore/QStringList>
00031 #include <QtCore/QtConcurrentMap>
00032
00033 #define KCAL_PER_MOL_TO_EV 0.043364122
00034
00035 using namespace Avogadro;
00036 using namespace OpenBabel;
00037 using namespace Eigen;
00038 using namespace std;
00039
00040 namespace GlobalSearch {
00041
00042 Structure::Structure(QObject *parent) :
00043 Molecule(parent),
00044 m_hasEnthalpy(false),
00045 m_updatedSinceDupChecked(true),
00046 m_histogramGenerationPending(false),
00047 m_hasBestOffspring(false),
00048 m_generation(0),
00049 m_id(0),
00050 m_rank(0),
00051 m_jobID(0),
00052 m_PV(0),
00053 m_optStart(QDateTime()),
00054 m_optEnd(QDateTime()),
00055 m_index(-1)
00056 {
00057 m_currentOptStep = 1;
00058 setStatus(Empty);
00059 resetFailCount();
00060 }
00061
00062 Structure::Structure(const Structure &other) :
00063 Molecule(other),
00064 m_histogramGenerationPending(false),
00065 m_updatedSinceDupChecked(true),
00066 m_hasBestOffspring(false),
00067 m_generation(0),
00068 m_id(0),
00069 m_rank(0),
00070 m_jobID(0),
00071 m_PV(0),
00072 m_optStart(QDateTime()),
00073 m_optEnd(QDateTime()),
00074 m_index(-1)
00075 {
00076 *this = other;
00077 }
00078
00079 Structure::Structure(const Avogadro::Molecule &other) :
00080 Molecule(other),
00081 m_histogramGenerationPending(false),
00082 m_updatedSinceDupChecked(true),
00083 m_hasBestOffspring(false),
00084 m_generation(0),
00085 m_id(0),
00086 m_rank(0),
00087 m_jobID(0),
00088 m_PV(0),
00089 m_optStart(QDateTime()),
00090 m_optEnd(QDateTime()),
00091 m_index(-1)
00092 {
00093 *this = other;
00094 }
00095
00096 void Structure::setupConnections()
00097 {
00098
00099 QList<const char*> slotlist;
00100 slotlist.append(SLOT(structureChanged()));
00101 for (int i = 0; i < slotlist.size(); i++) {
00102 connect(this, SIGNAL(updated()),
00103 this, slotlist.at(i),
00104 Qt::QueuedConnection);
00105 connect(this, SIGNAL(atomAdded(Atom*)),
00106 this, slotlist.at(i),
00107 Qt::QueuedConnection);
00108 connect(this, SIGNAL(atomUpdated(Atom*)),
00109 this, slotlist.at(i),
00110 Qt::QueuedConnection);
00111 connect(this, SIGNAL(atomRemoved(Atom*)),
00112 this, slotlist.at(i),
00113 Qt::QueuedConnection);
00114 connect(this, SIGNAL(bondAdded(Bond*)),
00115 this, slotlist.at(i),
00116 Qt::QueuedConnection);
00117 connect(this, SIGNAL(bondUpdated(Bond*)),
00118 this, slotlist.at(i),
00119 Qt::QueuedConnection);
00120 connect(this, SIGNAL(bondRemoved(Bond*)),
00121 this, slotlist.at(i),
00122 Qt::QueuedConnection);
00123 connect(this, SIGNAL(primitiveAdded(Primitive*)),
00124 this, slotlist.at(i),
00125 Qt::QueuedConnection);
00126 connect(this, SIGNAL(primitiveUpdated(Primitive*)),
00127 this, slotlist.at(i),
00128 Qt::QueuedConnection);
00129 connect(this, SIGNAL(primitiveRemoved(Primitive*)),
00130 this, slotlist.at(i),
00131 Qt::QueuedConnection);
00132 }
00133 }
00134
00135 void Structure::enableAutoHistogramGeneration(bool b)
00136 {
00137 if (b) {
00138 connect(this, SIGNAL(updated()),
00139 this, SLOT(requestHistogramGeneration()),
00140 Qt::QueuedConnection);
00141 connect(this, SIGNAL(atomAdded(Atom*)),
00142 this, SLOT(requestHistogramGeneration()),
00143 Qt::QueuedConnection);
00144 connect(this, SIGNAL(atomUpdated(Atom*)),
00145 this, SLOT(requestHistogramGeneration()),
00146 Qt::QueuedConnection);
00147 connect(this, SIGNAL(atomRemoved(Atom*)),
00148 this, SLOT(requestHistogramGeneration()),
00149 Qt::QueuedConnection);
00150 connect(this, SIGNAL(bondAdded(Bond*)),
00151 this, SLOT(requestHistogramGeneration()),
00152 Qt::QueuedConnection);
00153 connect(this, SIGNAL(bondUpdated(Bond*)),
00154 this, SLOT(requestHistogramGeneration()),
00155 Qt::QueuedConnection);
00156 connect(this, SIGNAL(bondRemoved(Bond*)),
00157 this, SLOT(requestHistogramGeneration()),
00158 Qt::QueuedConnection);
00159 connect(this, SIGNAL(primitiveAdded(Primitive*)),
00160 this, SLOT(requestHistogramGeneration()),
00161 Qt::QueuedConnection);
00162 connect(this, SIGNAL(primitiveUpdated(Primitive*)),
00163 this, SLOT(requestHistogramGeneration()),
00164 Qt::QueuedConnection);
00165 connect(this, SIGNAL(primitiveRemoved(Primitive*)),
00166 this, SLOT(requestHistogramGeneration()),
00167 Qt::QueuedConnection);
00168 } else {
00169 disconnect(this, SIGNAL(updated()),
00170 this, SLOT(requestHistogramGeneration()));
00171 disconnect(this, SIGNAL(atomAdded(Atom*)),
00172 this, SLOT(requestHistogramGeneration()));
00173 disconnect(this, SIGNAL(atomUpdated(Atom*)),
00174 this, SLOT(requestHistogramGeneration()));
00175 disconnect(this, SIGNAL(atomRemoved(Atom*)),
00176 this, SLOT(requestHistogramGeneration()));
00177 disconnect(this, SIGNAL(bondAdded(Bond*)),
00178 this, SLOT(requestHistogramGeneration()));
00179 disconnect(this, SIGNAL(bondUpdated(Bond*)),
00180 this, SLOT(requestHistogramGeneration()));
00181 disconnect(this, SIGNAL(bondRemoved(Bond*)),
00182 this, SLOT(requestHistogramGeneration()));
00183 disconnect(this, SIGNAL(primitiveAdded(Primitive*)),
00184 this, SLOT(requestHistogramGeneration()));
00185 disconnect(this, SIGNAL(primitiveUpdated(Primitive*)),
00186 this, SLOT(requestHistogramGeneration()));
00187 disconnect(this, SIGNAL(primitiveRemoved(Primitive*)),
00188 this, SLOT(requestHistogramGeneration()));
00189 }
00190 }
00191
00192 Structure::~Structure()
00193 {
00194 }
00195
00196 Structure& Structure::operator=(const Structure& other)
00197 {
00198 this->copyStructure(other);
00199
00200
00201 m_histogramGenerationPending = other.m_histogramGenerationPending;
00202 m_hasEnthalpy = other.m_hasEnthalpy;
00203 m_generation = other.m_generation;
00204 m_id = other.m_id;
00205 m_rank = other.m_rank;
00206 m_jobID = other.m_jobID;
00207 m_currentOptStep = other.m_currentOptStep;
00208 m_failCount = other.m_failCount;
00209 m_parents = other.m_parents;
00210 m_dupString = other.m_dupString;
00211 m_rempath = other.m_rempath;
00212 m_enthalpy = other.m_enthalpy;
00213 m_PV = other.m_PV;
00214 m_status = other.m_status;
00215 m_optStart = other.m_optStart;
00216 m_optEnd = other.m_optEnd;
00217 m_index = other.m_index;
00218
00219 return *this;
00220 }
00221
00222 Structure& Structure::operator=(const Avogadro::Molecule &mol)
00223 {
00224 this->copyStructure(mol);
00225 return *this;
00226 }
00227
00228 Structure& Structure::copyStructure(const Structure &other)
00229 {
00230 Molecule::operator=(other);
00231
00232 if (other.OBUnitCell()) {
00233 OpenBabel::OBUnitCell *cell = new OpenBabel::OBUnitCell(*other.OBUnitCell());
00234 setOBUnitCell(cell);
00235 }
00236
00237 return *this;
00238 }
00239
00240 void Structure::writeStructureSettings(const QString &filename)
00241 {
00242 SETTINGS(filename);
00243 const int VERSION = 2;
00244 settings->beginGroup("structure");
00245 settings->setValue("version", VERSION);
00246 settings->setValue("generation", getGeneration());
00247 settings->setValue("id", getIDNumber());
00248 settings->setValue("index", getIndex());
00249 settings->setValue("rank", getRank());
00250 settings->setValue("jobID", getJobID());
00251 settings->setValue("currentOptStep", getCurrentOptStep());
00252 settings->setValue("parents", getParents());
00253 settings->setValue("rempath", getRempath());
00254 settings->setValue("status", int(getStatus()));
00255 settings->setValue("failCount", getFailCount());
00256 settings->setValue("startTime", getOptTimerStart().toString());
00257 settings->setValue("endTime", getOptTimerEnd().toString());
00258 settings->setValue("needsPreoptimization", needsPreoptimization());
00259 settings->setValue("hasBestOffspring", hasBestOffspring());
00260
00261
00262 settings->beginGroup("history");
00263
00264 settings->beginWriteArray("atomicNums");
00265 for (int i = 0; i < m_histAtomicNums.size(); i++) {
00266 settings->setArrayIndex(i);
00267 const QList<unsigned int> *ptr = &(m_histAtomicNums.at(i));
00268 settings->beginWriteArray(QString("atomicNums-%1").arg(i));
00269 for (int j = 0; j < ptr->size(); j++) {
00270 settings->setArrayIndex(j);
00271 settings->setValue("value", ptr->at(j));
00272 }
00273 settings->endArray();
00274 }
00275 settings->endArray();
00276
00277
00278 settings->beginWriteArray("coords");
00279 for (int i = 0; i < m_histCoords.size(); i++) {
00280 settings->setArrayIndex(i);
00281 const QList<Eigen::Vector3d> *ptr = &(m_histCoords.at(i));
00282 settings->beginWriteArray(QString("coords-%1").arg(i));
00283 for (int j = 0; j < ptr->size(); j++) {
00284 settings->setArrayIndex(j);
00285 settings->setValue("x", ptr->at(j).x());
00286 settings->setValue("y", ptr->at(j).y());
00287 settings->setValue("z", ptr->at(j).z());
00288 }
00289 settings->endArray();
00290 }
00291 settings->endArray();
00292
00293
00294 settings->beginWriteArray("energies");
00295 for (int i = 0; i < m_histEnergies.size(); i++) {
00296 settings->setArrayIndex(i);
00297 settings->setValue("value", m_histEnergies.at(i));
00298 }
00299 settings->endArray();
00300
00301
00302 settings->beginWriteArray("enthalpies");
00303 for (int i = 0; i < m_histEnthalpies.size(); i++) {
00304 settings->setArrayIndex(i);
00305 settings->setValue("value", m_histEnthalpies.at(i));
00306 }
00307 settings->endArray();
00308
00309
00310 settings->beginWriteArray("cells");
00311 for (int i = 0; i < m_histCells.size(); i++) {
00312 settings->setArrayIndex(i);
00313 const Eigen::Matrix3d *ptr = &m_histCells.at(i);
00314 settings->setValue("00", (*ptr)(0, 0));
00315 settings->setValue("01", (*ptr)(0, 1));
00316 settings->setValue("02", (*ptr)(0, 2));
00317 settings->setValue("10", (*ptr)(1, 0));
00318 settings->setValue("11", (*ptr)(1, 1));
00319 settings->setValue("12", (*ptr)(1, 2));
00320 settings->setValue("20", (*ptr)(2, 0));
00321 settings->setValue("21", (*ptr)(2, 1));
00322 settings->setValue("22", (*ptr)(2, 2));
00323 }
00324 settings->endArray();
00325
00326 settings->endGroup();
00327
00328
00329 settings->beginWriteArray("optimizerLookup");
00330 QList<int> optKeys = m_optimizerLookup.keys();
00331 for (int i = 0; i < optKeys.size(); ++i) {
00332 settings->setArrayIndex(i);
00333 settings->setValue("key", optKeys.at(i));
00334 settings->setValue("value", m_optimizerLookup.value(optKeys.at(i)));
00335 }
00336 settings->endArray();
00337
00338 settings->endGroup();
00339 DESTROY_SETTINGS(filename);
00340 }
00341
00342 void Structure::readStructureSettings(const QString &filename)
00343 {
00344 SETTINGS(filename);
00345 settings->beginGroup("structure");
00346 int loadedVersion = settings->value("version", 0).toInt();
00347 if (loadedVersion >= 1) {
00348 setGeneration( settings->value("generation", 0).toInt());
00349 setIDNumber( settings->value("id", 0).toInt());
00350 setIndex( settings->value("index", 0).toInt());
00351 setRank( settings->value("rank", 0).toInt());
00352 setJobID( settings->value("jobID", 0).toInt());
00353 setCurrentOptStep( settings->value("currentOptStep", 0).toInt());
00354 setFailCount( settings->value("failCount", 0).toInt());
00355 setParents( settings->value("parents", "").toString());
00356 setRempath( settings->value("rempath", "").toString());
00357 setStatus( State(settings->value("status", -1).toInt()));
00358
00359 setOptTimerStart( QDateTime::fromString(settings->value("startTime", "").toString()));
00360 setOptTimerEnd( QDateTime::fromString(settings->value("endTime", "").toString()));
00361
00362 setNeedsPreoptimization(settings->value("needsPreoptimization", true).toBool());
00363 setHasBestOffspring(settings->value("hasBestOffspring", false).toBool());
00364
00365
00366 settings->beginGroup("history");
00367
00368 int size, size2;
00369 size = settings->beginReadArray("atomicNums");
00370 m_histAtomicNums.clear();
00371 for (int i = 0; i < size; i++) {
00372 settings->setArrayIndex(i);
00373 size2 = settings->beginReadArray(QString("atomicNums-%1").arg(i));
00374 QList<unsigned int> cur;
00375 for (int j = 0; j < size2; j++) {
00376 settings->setArrayIndex(j);
00377 cur.append(settings->value("value").toUInt());
00378 }
00379 settings->endArray();
00380 m_histAtomicNums.append(cur);
00381 }
00382 settings->endArray();
00383
00384
00385 size = settings->beginReadArray("coords");
00386 m_histCoords.clear();
00387 for (int i = 0; i < size; i++) {
00388 settings->setArrayIndex(i);
00389 size2 = settings->beginReadArray(QString("coords-%1").arg(i));
00390 QList<Eigen::Vector3d> cur;
00391 for (int j = 0; j < size2; j++) {
00392 settings->setArrayIndex(j);
00393 double x = settings->value("x").toDouble();
00394 double y = settings->value("y").toDouble();
00395 double z = settings->value("z").toDouble();
00396 cur.append(Eigen::Vector3d(x, y, z));
00397 }
00398 settings->endArray();
00399 m_histCoords.append(cur);
00400 }
00401 settings->endArray();
00402
00403
00404 size = settings->beginReadArray("energies");
00405 m_histEnergies.clear();
00406 for (int i = 0; i < size; i++) {
00407 settings->setArrayIndex(i);
00408 m_histEnergies.append(settings->value("value").toDouble());
00409 }
00410 settings->endArray();
00411
00412
00413 size = settings->beginReadArray("enthalpies");
00414 m_histEnthalpies.clear();
00415 for (int i = 0; i < size; i++) {
00416 settings->setArrayIndex(i);
00417 m_histEnthalpies.append(settings->value("value").toDouble());
00418 }
00419 settings->endArray();
00420
00421
00422 size = settings->beginReadArray("cells");
00423 m_histCells.clear();
00424 for (int i = 0; i < size; i++) {
00425 settings->setArrayIndex(i);
00426 Eigen::Matrix3d cur;
00427 cur(0, 0) = settings->value("00").toDouble();
00428 cur(0, 1) = settings->value("01").toDouble();
00429 cur(0, 2) = settings->value("02").toDouble();
00430 cur(1, 0) = settings->value("10").toDouble();
00431 cur(1, 1) = settings->value("11").toDouble();
00432 cur(1, 2) = settings->value("12").toDouble();
00433 cur(2, 0) = settings->value("20").toDouble();
00434 cur(2, 1) = settings->value("21").toDouble();
00435 cur(2, 2) = settings->value("22").toDouble();
00436 m_histCells.append(cur);
00437 }
00438 settings->endArray();
00439
00440 settings->endGroup();
00441
00442
00443 size = settings->beginReadArray("optimizerLookup");
00444 m_optimizerLookup.clear();
00445 for (int i = 0; i < size; ++i) {
00446 settings->setArrayIndex(i);
00447 int key = settings->value("key", -1).toInt();
00448 int value = settings->value("value", -1).toInt();
00449 if (key < 0 || value < 0) {
00450 qDebug() << "Missing key (" << key << ") value (" << value << ")"
00451 "pair for structure" << this->getIDString();
00452 continue;
00453 }
00454 m_optimizerLookup.insert(key, value);
00455 }
00456 settings->endArray();
00457
00458 }
00459 settings->endGroup();
00460
00461
00462 switch (loadedVersion) {
00463 case 0: {
00464
00465 qDebug() << "Updating "
00466 << filename
00467 << " from Version 0 -> 1";
00468 QFile file (filename);
00469 file.open(QIODevice::ReadOnly);
00470 QTextStream stream (&file);
00471 load(stream);
00472 }
00473 case 1:
00474
00475 case 2:
00476 default:
00477 break;
00478 }
00479 }
00480
00481 void Structure::structureChanged()
00482 {
00483 m_updatedSinceDupChecked = true;
00484 }
00485
00486 void Structure::updateAndSkipHistory(const QList<unsigned int> &atomicNums,
00487 const QList<Eigen::Vector3d> &coords,
00488 const double energy,
00489 const double enthalpy,
00490 const Eigen::Matrix3d &cell)
00491 {
00492 Q_ASSERT_X(atomicNums.size() == coords.size() && coords.size() == numAtoms(),
00493 Q_FUNC_INFO,
00494 "Lengths of atomicNums and coords must match numAtoms().");
00495
00496
00497 Atom *atm;
00498 bool useLUT = !m_optimizerLookup.isEmpty();
00499 for (int optIndex = 0; optIndex < numAtoms(); ++optIndex) {
00500 int atomIndex = optIndex;
00501 if (useLUT) {
00502 atomIndex = m_optimizerLookup.value(optIndex);
00503 }
00504 atm = atom(atomIndex);
00505 atm->setAtomicNumber(atomicNums.at(optIndex));
00506 atm->setPos(coords.at(optIndex));
00507 }
00508
00509
00510 if (fabs(enthalpy) < 1e-6) {
00511 m_hasEnthalpy = false;
00512 m_PV = 0.0;
00513 }
00514 else {
00515 m_hasEnthalpy = true;
00516 m_PV = enthalpy - energy;
00517 }
00518 m_enthalpy = enthalpy;
00519 setEnergy(energy * EV_TO_KJ_PER_MOL);
00520
00521
00522 if (!cell.isZero()) {
00523 OpenBabel::matrix3x3 obmat;
00524 for (int i = 0; i < 3; i++) {
00525 for (int j = 0; j < 3; j++) {
00526 obmat.Set(i,j,cell(i,j));
00527 }
00528 }
00529 OpenBabel::OBUnitCell *obcell = OBUnitCell();
00530 obcell->SetData(obmat);
00531 setOBUnitCell(obcell);
00532 }
00533
00534 emit moleculeChanged();
00535 update();
00536 }
00537
00538 void Structure::updateAndAddToHistory(const QList<unsigned int> &atomicNums,
00539 const QList<Eigen::Vector3d> &coords,
00540 const double energy,
00541 const double enthalpy,
00542 const Eigen::Matrix3d &cell)
00543 {
00544 Q_ASSERT_X(atomicNums.size() == coords.size() && coords.size() == numAtoms(),
00545 Q_FUNC_INFO,
00546 "Lengths of atomicNums and coords must match numAtoms().");
00547
00548
00549 m_histAtomicNums.append(atomicNums);
00550 m_histCoords.append(coords);
00551 m_histEnergies.append(energy);
00552 m_histEnthalpies.append(enthalpy);
00553 m_histCells.append(cell);
00554
00555
00556 Atom *atm;
00557 bool useLUT = !m_optimizerLookup.isEmpty();
00558 for (int optIndex = 0; optIndex < numAtoms(); ++optIndex) {
00559 int atomIndex = optIndex;
00560 if (useLUT) {
00561 atomIndex = m_optimizerLookup.value(optIndex);
00562 }
00563 atm = atom(atomIndex);
00564 atm->setAtomicNumber(atomicNums.at(optIndex));
00565 atm->setPos(coords.at(optIndex));
00566 }
00567
00568
00569 if (fabs(enthalpy) < 1e-6) {
00570 m_hasEnthalpy = false;
00571 m_PV = 0.0;
00572 }
00573 else {
00574 m_hasEnthalpy = true;
00575 m_PV = enthalpy - energy;
00576 }
00577 m_enthalpy = enthalpy;
00578 setEnergy(energy * EV_TO_KJ_PER_MOL);
00579
00580
00581 if (!cell.isZero()) {
00582 OpenBabel::matrix3x3 obmat;
00583 for (int i = 0; i < 3; i++) {
00584 for (int j = 0; j < 3; j++) {
00585 obmat.Set(i,j,cell(i,j));
00586 }
00587 }
00588 OpenBabel::OBUnitCell *obcell = OBUnitCell();
00589 obcell->SetData(obmat);
00590 setOBUnitCell(obcell);
00591 }
00592
00593 emit moleculeChanged();
00594 update();
00595 }
00596
00597 void Structure::deleteFromHistory(unsigned int index) {
00598 Q_ASSERT_X(index <= sizeOfHistory() - 1, Q_FUNC_INFO,
00599 "Requested history index greater than the number of available entries.");
00600
00601 m_histAtomicNums.removeAt(index);
00602 m_histEnthalpies.removeAt(index);
00603 m_histEnergies.removeAt(index);
00604 m_histCoords.removeAt(index);
00605 m_histCells.removeAt(index);
00606 }
00607
00608 void Structure::retrieveHistoryEntry(unsigned int index,
00609 QList<unsigned int> *atomicNums,
00610 QList<Eigen::Vector3d> *coords,
00611 double *energy,
00612 double *enthalpy,
00613 Eigen::Matrix3d *cell)
00614 {
00615 Q_ASSERT_X(index <= sizeOfHistory() - 1, Q_FUNC_INFO,
00616 "Requested history index greater than the number of available entries.");
00617
00618 if (atomicNums != NULL) {
00619 *atomicNums = m_histAtomicNums.at(index);
00620 }
00621 if (coords != NULL) {
00622 *coords = m_histCoords.at(index);
00623 }
00624 if (energy != NULL) {
00625 *energy = m_histEnergies.at(index);
00626 }
00627 if (enthalpy != NULL) {
00628 *enthalpy = m_histEnthalpies.at(index);
00629 }
00630 if (cell != NULL) {
00631 *cell = m_histCells.at(index);
00632 }
00633
00634 }
00635
00636 bool Structure::addAtomRandomly(uint atomicNumber, double minIAD, double maxIAD, int maxAttempts, Atom **atom)
00637 {
00638 INIT_RANDOM_GENERATOR();
00639 double IAD = -1;
00640 int i = 0;
00641 vector3 coords;
00642
00643
00644 if (numAtoms() == 0) {
00645 coords = vector3 (0,0,0);
00646 }
00647 else {
00648 do {
00649
00650 IAD = -1;
00651 double x = RANDDOUBLE() * radius() + maxIAD;
00652 double y = RANDDOUBLE() * radius() + maxIAD;
00653 double z = RANDDOUBLE() * radius() + maxIAD;
00654
00655 coords.Set(x,y,z);
00656 if (minIAD != -1) {
00657 getNearestNeighborDistance(x, y, z, IAD);
00658 }
00659 else { break;};
00660 i++;
00661 } while (i < maxAttempts && IAD <= minIAD);
00662
00663 if (i >= maxAttempts) return false;
00664 }
00665
00666 Atom *atm = addAtom();
00667 atom = &atm;
00668 Eigen::Vector3d pos (coords[0],coords[1],coords[2]);
00669 (*atom)->setPos(pos);
00670 (*atom)->setAtomicNumber(static_cast<int>(atomicNumber));
00671 return true;
00672 }
00673
00674 void Structure::setOBEnergy(const QString &ff) {
00675 OBForceField* pFF = OBForceField::FindForceField(ff.toStdString().c_str());
00676 if (!pFF) return;
00677 OpenBabel::OBMol obmol = OBMol();
00678 if (!pFF->Setup(obmol)) {
00679 qWarning() << "Structure::setOBEnergy: could not setup force field " << ff << ".";
00680 return;
00681 }
00682 std::vector<double> E;
00683 E.push_back(pFF->Energy());
00684 setEnergies(E);
00685 }
00686
00687 QString Structure::getResultsEntry() const
00688 {
00689 QString status;
00690 switch (getStatus()) {
00691 case Optimized:
00692 status = "Optimized";
00693 break;
00694 case Killed:
00695 case Removed:
00696 status = "Killed";
00697 break;
00698 case Duplicate:
00699 status = "Duplicate";
00700 break;
00701 case Error:
00702 status = "Error";
00703 break;
00704 case StepOptimized:
00705 case WaitingForOptimization:
00706 case InProcess:
00707 case Empty:
00708 case Updating:
00709 case Submitted:
00710 case Restart:
00711 status = "In progress";
00712 break;
00713 case Preoptimizing:
00714 status = "Preoptimizing";
00715 break;
00716 default:
00717 status = "Unknown";
00718 break;
00719 }
00720 return QString("%1 %2 %3 %4 %5")
00721 .arg(getRank(), 6)
00722 .arg(getGeneration(), 6)
00723 .arg(getIDNumber(), 6)
00724 .arg(getEnthalpy(), 10)
00725 .arg(status, 11);
00726 };
00727
00728 bool Structure::getNearestNeighborDistances(QList<double> * list) const
00729 {
00730 list->clear();
00731 QList<Atom *> atomList = atoms();
00732 if (atomList.size() < 2) return false;
00733
00734 #if QT_VERSION >= 0x040700
00735 list->reserve(atomList.size());
00736 #endif // QT_VERSION
00737
00738 double shortest;
00739 for (QList<Atom*>::const_iterator it = atomList.constBegin(),
00740 it_end = atomList.constEnd(); it != it_end; ++it) {
00741 getNearestNeighborDistance((*it), shortest);
00742 list->append(shortest);
00743 }
00744 return true;
00745 }
00746
00747 bool Structure::getShortestInteratomicDistance(double & shortest) const
00748 {
00749 QList<Atom*> atomList = atoms();
00750 if (atomList.size() <= 1) return false;
00751 QList<Vector3d> atomPositions;
00752 for (int i = 0; i < atomList.size(); i++)
00753 atomPositions.push_back(*(atomList.at(i)->pos()));
00754
00755 Vector3d v1= atomPositions.at(0);
00756 Vector3d v2= atomPositions.at(1);
00757 shortest = abs((v1-v2).norm());
00758 double distance;
00759
00760
00761 for (int i = 0; i < atomList.size(); i++) {
00762 v1 = atomPositions.at(i);
00763 for (int j = i+1; j < atomList.size(); j++) {
00764 v2 = atomPositions.at(j);
00765
00766 distance = abs((v1-v2).norm());
00767 if (distance < shortest) shortest = distance;
00768 }
00769 }
00770 return true;
00771 }
00772
00773 bool Structure::getNearestNeighborDistance(const double x,
00774 const double y,
00775 const double z,
00776 double & shortest) const
00777 {
00778 QList<Atom*> atomList = atoms();
00779 if (atomList.size() < 2) return false;
00780 QList<Vector3d> atomPositions;
00781 for (int i = 0; i < atomList.size(); i++)
00782 atomPositions.push_back(*(atomList.at(i)->pos()));
00783
00784 Vector3d v1 (x, y, z);
00785 Vector3d v2 = atomPositions.at(0);
00786 shortest = fabs((v1-v2).norm());
00787 double distance;
00788
00789
00790 for (int j = 0; j < atomList.size(); j++) {
00791 v2 = atomPositions.at(j);
00792
00793 distance = abs((v1-v2).norm());
00794 if (distance < shortest) shortest = distance;
00795 }
00796 return true;
00797 }
00798
00799 bool Structure::getNearestNeighborDistance(const Avogadro::Atom *atom,
00800 double & shortest) const
00801 {
00802 const Eigen::Vector3d *position = atom->pos();
00803 return getNearestNeighborDistance(position->x(),
00804 position->y(),
00805 position->z(), shortest);
00806 }
00807
00808 QList<Atom*> Structure::getNeighbors(const double x,
00809 const double y,
00810 const double z,
00811 const double cutoff,
00812 QList<double> *distances) const
00813 {
00814 QList<Atom*> neighbors;
00815 if (distances) {
00816 distances->clear();
00817 }
00818 Eigen::Vector3d vec (x,y,z);
00819 double cutoffSquared = cutoff*cutoff;
00820 for (QList<Atom*>::const_iterator it = m_atomList.constBegin(),
00821 it_end = m_atomList.constEnd(); it != it_end; ++it) {
00822 double distSq = ((*(*it)->pos()) - vec).squaredNorm();
00823 if (distSq <= cutoffSquared) {
00824 neighbors << *it;
00825 if (distances) {
00826 *distances << sqrt(distSq);
00827 }
00828 }
00829 }
00830 return neighbors;
00831 }
00832
00833 QList<Atom*> Structure::getNeighbors(const Avogadro::Atom *atom,
00834 const double cutoff,
00835 QList<double> *distances) const
00836 {
00837 const Eigen::Vector3d *position = atom->pos();
00838 return getNeighbors(position->x(),
00839 position->y(),
00840 position->z(), cutoff, distances);
00841 }
00842
00843 void Structure::requestHistogramGeneration()
00844 {
00845 if (!m_histogramGenerationPending) {
00846 m_histogramGenerationPending = true;
00847
00848 QTimer::singleShot(250, this, SLOT(generateDefaultHistogram()));
00849 }
00850 }
00851
00852 void Structure::generateDefaultHistogram()
00853 {
00854 generateIADHistogram(&m_histogramDist, &m_histogramFreq, 0, 10, 0.01);
00855 m_histogramGenerationPending = false;
00856 }
00857
00858 void Structure::getDefaultHistogram(QList<double> *dist, QList<double> *freq) const
00859 {
00860 dist->clear();
00861 freq->clear();
00862 for (int i = 0; i < m_histogramDist.size(); i++) {
00863 dist->append(m_histogramDist.at(i).toDouble());
00864 freq->append(m_histogramFreq.at(i).toDouble());
00865 }
00866 }
00867
00868 void Structure::getDefaultHistogram(QList<QVariant> *dist, QList<QVariant> *freq) const
00869 {
00870 (*dist) = m_histogramDist;
00871 (*freq) = m_histogramFreq;
00872 }
00873
00874 bool Structure::generateIADHistogram(QList<double> * dist,
00875 QList<double> * freq,
00876 double min,
00877 double max,
00878 double step,
00879 Avogadro::Atom *atom) const
00880 {
00881 QList<QVariant> distv, freqv;
00882 if (!generateIADHistogram(&distv, &freqv, min, max, step, atom)) {
00883 return false;
00884 }
00885 dist->clear();
00886 freq->clear();
00887 for (int i = 0; i < distv.size(); i++) {
00888 dist->append(distv.at(i).toDouble());
00889 freq->append(freqv.at(i).toDouble());
00890 }
00891 return true;
00892 }
00893
00894
00895
00897 struct NNHistMap {
00898 int i;
00899 double halfstep;
00900 QList<Vector3d> *atomPositions;
00901 QList<QVariant> *dist;
00902 };
00903
00904
00906
00907
00908 QList<int> calcNNHistChunk(const NNHistMap &m)
00909 {
00910 const Vector3d *v1 = &(m.atomPositions->at(m.i));
00911 const Vector3d *v2;
00912 QList<int> freq;
00913 double diff;
00914 for (int ind = 0; ind < m.dist->size(); ind++) {
00915 freq.append(0);
00916 }
00917 for (int j = m.i+1; j < m.atomPositions->size(); j++) {
00918 v2 = &(m.atomPositions->at(j));
00919 diff = fabs(((*v1)-(*v2)).norm());
00920 for (int k = 0; k < m.dist->size(); k++) {
00921 if (fabs(diff-(m.dist->at(k).toDouble())) < m.halfstep) {
00922 freq[k]++;
00923 }
00924 }
00925 }
00926 return freq;
00927 }
00928
00929 void reduceNNHistChunks(QList<QVariant> &final, const QList<int> &tmp)
00930 {
00931 if (final.size() != tmp.size()) {
00932 final.clear();
00933 for (int i = 0; i < tmp.size(); i++) {
00934 final.append(tmp.at(i));
00935 }
00936 }
00937 else {
00938 double d;
00939 for (int i = 0; i < final.size(); i++) {
00940 d = final.at(i).toDouble();
00941 d += tmp.at(i);
00942 final.replace(i, d);
00943 }
00944 }
00945 return;
00946 }
00947
00948 bool Structure::generateIADHistogram(QList<QVariant> * distance,
00949 QList<QVariant> * frequency,
00950 double min, double max, double step,
00951 Avogadro::Atom *atom) const
00952 {
00953 distance->clear();
00954 frequency->clear();
00955
00956 if (min > max && step > 0) {
00957 qWarning() << "Structure::getNearestNeighborHistogram: min cannot be greater than max!";
00958 return false;
00959 }
00960 if (step <= 0) {
00961 qWarning() << "Structure::getNearestNeighborHistogram: invalid step size:" << step;
00962 return false;
00963 }
00964
00965 double halfstep = step / 2.0;
00966
00967 double val = min;
00968 do {
00969 distance->append(val);
00970 frequency->append(0);
00971 val += step;
00972 } while (val < max);
00973
00974 QList<Atom*> atomList = atoms();
00975 QList<Vector3d> atomPositions;
00976 for (int i = 0; i < atomList.size(); i++)
00977 atomPositions.push_back(*(atomList.at(i)->pos()));
00978
00979 Vector3d v1= atomPositions.at(0);
00980 Vector3d v2= atomPositions.at(1);
00981 double diff;
00982
00983
00984
00985 if (atom == 0) {
00986 QList<NNHistMap> ml;
00987 for (int i = 0; i < atomList.size(); i++) {
00988 NNHistMap m;
00989 m.i = i; m.halfstep = halfstep; m.atomPositions = &atomPositions; m.dist = distance;
00990 ml.append(m);
00991 }
00992 (*frequency) = QtConcurrent::blockingMappedReduced(ml, calcNNHistChunk, reduceNNHistChunks);
00993 }
00994
00995 else {
00996 v1 = *atom->pos();
00997 for (int j = 0; j < atomList.size(); j++) {
00998 if (atomList.at(j) == atom) continue;
00999 v2 = atomPositions.at(j);
01000
01001 diff = fabs((v1-v2).norm());
01002 for (int k = 0; k < distance->size(); k++) {
01003 double radius = distance->at(k).toDouble();
01004 double d;
01005 if (diff != 0 && fabs(diff-radius) < halfstep) {
01006 d = frequency->at(k).toDouble();
01007 d++;
01008 frequency->replace(k, d);
01009 }
01010 }
01011 }
01012 }
01013
01014 return true;
01015 }
01016
01017 bool Structure::compareIADDistributions(const std::vector<double> &d,
01018 const std::vector<double> &f1,
01019 const std::vector<double> &f2,
01020 double decay,
01021 double smear,
01022 double *error)
01023 {
01024
01025 if (smear != 0 && d.size() <= 1) {
01026 qWarning() << "Cluster::compareNNDist: Cannot smear with 1 or fewer points.";
01027 return false;
01028 }
01029
01030 if (d.size() != f1.size() || f1.size() != f2.size()) {
01031 qWarning() << "Cluster::compareNNDist: Vectors are not the same size.";
01032 return false;
01033 }
01034
01035
01036
01037 double stepSize = fabs(d.at(1) - d.at(0));
01038 int boxSize = ceil(smear/stepSize);
01039 if (boxSize > d.size()) {
01040 qWarning() << "Cluster::compareNNDist: Smear length is greater then d vector range.";
01041 return false;
01042 }
01043
01044 vector<double> f1s, f2s, ds;
01045 if (smear != 0) {
01046 double f1t, f2t, dt;
01047 for (int i = 0; i < d.size() - boxSize; i++) {
01048 f1t = f2t = dt = 0;
01049 for (int j = 0; j < boxSize; j++) {
01050 f1t += f1.at(i+j);
01051 f2t += f2.at(i+j);
01052 }
01053 f1s.push_back(f1t / double(boxSize));
01054 f2s.push_back(f2t / double(boxSize));
01055 ds.push_back(dt / double(boxSize));
01056 }
01057 } else {
01058 for (int i = 0; i < d.size() - boxSize; i++) {
01059 f1s.push_back(f1.at(i));
01060 f2s.push_back(f2.at(i));
01061 ds.push_back(d.at(i));
01062 }
01063 }
01064
01065
01066 vector<double> diff;
01067 for (int i = 0; i < ds.size(); i++) {
01068 diff.push_back(fabs(f1s.at(i) - f2s.at(i)));
01069 }
01070
01071
01072
01073 double decayFactor = 0;
01074
01075 if (decay != 0) {
01076 decayFactor = 0.69314718055994530941723 / decay;
01077 }
01078
01079
01080 (*error) = 0;
01081 for (int i = 0; i < ds.size(); i++) {
01082 (*error) += exp(-decayFactor * ds.at(i)) * diff.at(i);
01083 }
01084
01085 return true;
01086 }
01087
01088 bool Structure::compareIADDistributions(const QList<double> &d,
01089 const QList<double> &f1,
01090 const QList<double> &f2,
01091 double decay,
01092 double smear,
01093 double *error)
01094 {
01095
01096 if (d.size() != f1.size() || f1.size() != f2.size()) {
01097 qWarning() << "Cluster::compareIADDist: Vectors are not the same size.";
01098 return false;
01099 }
01100 vector<double> dd, f1d, f2d;
01101 for (int i = 0; i < d.size(); i++) {
01102 dd.push_back(d.at(i));
01103 f1d.push_back(f1.at(i));
01104 f2d.push_back(f2.at(i));
01105 }
01106 return compareIADDistributions(dd, f1d, f2d, decay, smear, error);
01107 }
01108
01109 bool Structure::compareIADDistributions(const QList<QVariant> &d,
01110 const QList<QVariant> &f1,
01111 const QList<QVariant> &f2,
01112 double decay,
01113 double smear,
01114 double *error)
01115 {
01116
01117 if (d.size() != f1.size() || f1.size() != f2.size()) {
01118 qWarning() << "Cluster::compareIADDist: Vectors are not the same size.";
01119 return false;
01120 }
01121 vector<double> dd, f1d, f2d;
01122 for (int i = 0; i < d.size(); i++) {
01123 dd.push_back(d.at(i).toDouble());
01124 f1d.push_back(f1.at(i).toDouble());
01125 f2d.push_back(f2.at(i).toDouble());
01126 }
01127 return compareIADDistributions(dd, f1d, f2d, decay, smear, error);
01128 }
01129
01130 QList<QString> Structure::getSymbols() const {
01131 QList<QString> list;
01132 for (QList<Atom*>::const_iterator
01133 it = m_atomList.begin(),
01134 it_end = m_atomList.end();
01135 it != it_end;
01136 ++it) {
01137 QString symbol = QString(OpenBabel::etab.GetSymbol((*it)->atomicNumber()));
01138 if (!list.contains(symbol)) {
01139 list.append(symbol);
01140 }
01141 }
01142 qSort(list);
01143 return list;
01144 }
01145
01146 QList<uint> Structure::getNumberOfAtomsAlpha() const {
01147 QList<uint> list;
01148 QList<QString> symbols = getSymbols();
01149 for (int i = 0; i < symbols.size(); i++)
01150 list.append(0);
01151 int ind;
01152 for (QList<Atom*>::const_iterator
01153 it = m_atomList.begin(),
01154 it_end = m_atomList.end();
01155 it != it_end;
01156 ++it) {
01157 QString symbol = QString(OpenBabel::etab.GetSymbol((*it)->atomicNumber()));
01158 Q_ASSERT_X(symbols.contains(symbol), Q_FUNC_INFO,
01159 "getNumberOfAtomsAlpha found a symbol not in getSymbols.");
01160 ind = symbols.indexOf(symbol);
01161 ++list[ind];
01162 }
01163 return list;
01164 }
01165
01166 QString Structure::getOptElapsed() const {
01167 int secs;
01168 if (m_optStart.toString() == "") return "0:00:00";
01169 if (m_optEnd.toString() == "")
01170 secs = m_optStart.secsTo(QDateTime::currentDateTime());
01171 else
01172 secs = m_optStart.secsTo(m_optEnd);
01173 int hours = static_cast<int>(secs/3600);
01174 int mins = static_cast<int>( (secs - hours * 3600) / 60);
01175 secs = secs % 60;
01176 QString ret;
01177 ret.sprintf("%d:%02d:%02d", hours, mins, secs);
01178 return ret;
01179 }
01180
01181 void Structure::load(QTextStream &in) {
01182 QString line, str;
01183 QStringList strl;
01184 setStatus(InProcess);
01185 in.seek(0);
01186 while (!in.atEnd()) {
01187 line = in.readLine();
01188 strl = line.split(QRegExp(" "));
01189
01190
01191 if (line.contains("Generation:") && strl.size() > 1)
01192 setGeneration( (strl.at(1)).toUInt() );
01193 if (line.contains("ID#:") && strl.size() > 1)
01194 setIDNumber( (strl.at(1)).toUInt() );
01195 if (line.contains("Index:") && strl.size() > 1)
01196 setIndex( (strl.at(1)).toUInt() );
01197 if (line.contains("Enthalpy Rank:") && strl.size() > 2)
01198 setRank( (strl.at(2)).toUInt() );
01199 if (line.contains("Job ID:") && strl.size() > 2)
01200 setJobID( (strl.at(2)).toUInt() );
01201 if (line.contains("Current INCAR:") && strl.size() > 2)
01202 setCurrentOptStep( (strl.at(2)).toUInt() );
01203 if (line.contains("Current OptStep:") && strl.size() > 2)
01204 setCurrentOptStep( (strl.at(2)).toUInt() );
01205 if (line.contains("Ancestry:") && strl.size() > 1) {
01206 strl.removeFirst();
01207 setParents( strl.join(" ") );
01208 }
01209 if (line.contains("Remote Path:") && strl.size() > 2) {
01210 strl.removeFirst(); strl.removeFirst();
01211 setRempath( strl.join(" ") );
01212 }
01213 if (line.contains("Start Time:") && strl.size() > 2) {
01214 strl.removeFirst(); strl.removeFirst();
01215 str = strl.join(" ");
01216 setOptTimerStart(QDateTime::fromString(str));
01217 }
01218 if (line.contains("Status:") && strl.size() > 1) {
01219 setStatus( Structure::State((strl.at(1)).toInt() ));
01220 }
01221 if (line.contains("failCount:") && strl.size() > 1) {
01222 setFailCount((strl.at(1)).toUInt());
01223 }
01224 if (line.contains("End Time:") && strl.size() > 2) {
01225 strl.removeFirst(); strl.removeFirst();
01226 str = strl.join(" ");
01227 setOptTimerEnd(QDateTime::fromString(str));
01228 }
01229 }
01230 }
01231
01232 QHash<QString, QVariant> Structure::getFingerprint() const
01233 {
01234 QHash<QString, QVariant> fp;
01235 fp.insert("enthalpy", getEnthalpy());
01236 return fp;
01237 }
01238
01239 void Structure::sortByEnthalpy(QList<Structure*> *structures)
01240 {
01241 uint numStructs = structures->size();
01242
01243 if (numStructs <= 1) return;
01244
01245
01246 Structure *structure_i=0, *structure_j=0, *tmp=0;
01247 for (uint i = 0; i < numStructs-1; i++) {
01248 structure_i = structures->at(i);
01249 structure_i->lock()->lockForRead();
01250 for (uint j = i+1; j < numStructs; j++) {
01251 structure_j = structures->at(j);
01252 structure_j->lock()->lockForRead();
01253 if (structure_j->getEnthalpy() < structure_i->getEnthalpy()) {
01254 structures->swap(i,j);
01255 tmp = structure_i;
01256 structure_i = structure_j;
01257 structure_j = tmp;
01258 }
01259 structure_j->lock()->unlock();
01260 }
01261 structure_i->lock()->unlock();
01262 }
01263 }
01264
01265 void rankInPlace(const QList<Structure*> &structures)
01266 {
01267 if (structures.size() == 0) return;
01268 Structure *s;
01269 for (uint i = 0; i < structures.size(); i++) {
01270 s = structures.at(i);
01271 s->lock()->lockForWrite();
01272 s->setRank(i+1);
01273 s->lock()->unlock();
01274 }
01275 }
01276
01277 void Structure::rankByEnthalpy(const QList<Structure*> &structures)
01278 {
01279 uint numStructs = structures.size();
01280
01281 if (numStructs == 0) return;
01282
01283 QList<Structure*> rstructures;
01284
01285
01286 for (uint i = 0; i < numStructs; i++)
01287 rstructures.append(structures.at(i));
01288
01289
01290 Structure *structure_i=0, *structure_j=0, *tmp=0;
01291 for (uint i = 0; i < numStructs-1; i++) {
01292 structure_i = rstructures.at(i);
01293 structure_i->lock()->lockForRead();
01294 for (uint j = i+1; j < numStructs; j++) {
01295 structure_j = rstructures.at(j);
01296 structure_j->lock()->lockForRead();
01297 if (structure_j->getEnthalpy() < structure_i->getEnthalpy()) {
01298 rstructures.swap(i,j);
01299 tmp = structure_i;
01300 structure_i = structure_j;
01301 structure_j = tmp;
01302 }
01303 structure_j->lock()->unlock();
01304 }
01305 structure_i->lock()->unlock();
01306 }
01307
01308 rankInPlace(rstructures);
01309 }
01310
01311 void Structure::sortAndRankByEnthalpy(QList<Structure*> *structures)
01312 {
01313 sortByEnthalpy(structures);
01314 rankInPlace(*structures);
01315 }
01316
01317
01318 }
