From c944183fe20a94bef6630c2cf62dcd09b16b2c6c Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 9 Sep 2023 11:39:21 -0400 Subject: [PATCH 1/2] Initial support for scripting surface generation Needs testing to make sure dialog isn't required Signed-off-by: Geoff Hutchison --- avogadro/core/cube.h | 1 + avogadro/qtplugins/surfaces/surfaces.cpp | 225 +++++++++++++++++++---- avogadro/qtplugins/surfaces/surfaces.h | 35 ++-- 3 files changed, 210 insertions(+), 51 deletions(-) diff --git a/avogadro/core/cube.h b/avogadro/core/cube.h index b33d63c50..357087157 100644 --- a/avogadro/core/cube.h +++ b/avogadro/core/cube.h @@ -41,6 +41,7 @@ class AVOGADROCORE_EXPORT Cube VdW, ESP, ElectronDensity, + SpinDensity, MO, FromFile, None diff --git a/avogadro/qtplugins/surfaces/surfaces.cpp b/avogadro/qtplugins/surfaces/surfaces.cpp index bff82f66f..7b6010b21 100644 --- a/avogadro/qtplugins/surfaces/surfaces.cpp +++ b/avogadro/qtplugins/surfaces/surfaces.cpp @@ -102,6 +102,115 @@ Surfaces::~Surfaces() delete m_cube; } +void Surfaces::registerCommands() +{ + // register some scripting commands + emit registerCommand("renderVDW", tr("Render the van der Waals surface.")); + emit registerCommand("renderVanDerWaals", + tr("Render the van der Waals molecular surface.")); + emit registerCommand("renderSolventAccessible", + tr("Render the solvent-accessible molecular surface.")); + emit registerCommand("renderSolventExcluded", + tr("Render the solvent-excluded molecular surface.")); + emit registerCommand("renderOrbital", tr("Render a molecular orbital.")); + emit registerCommand("renderElectronDensity", + tr("Render the electron density.")); + emit registerCommand("renderSpinDensity", tr("Render the spin density.")); + emit registerCommand("renderCube", + tr("Render a cube supplied with the file.")); +} + +bool Surfaces::handleCommand(const QString& command, const QVariantMap& options) +{ + if (m_molecule == nullptr) + return false; // No molecule to handle the command. + + // Set up some defaults for the options. + int index = -1; + int homo = -1; + float isoValue = 0.03; + float cubeResolution = resolution(); // use default resolution + + if (options.contains("resolution") && + options["resolution"].canConvert()) { + bool ok; + float res = options["resolution"].toFloat(&ok); + if (ok) + cubeResolution = res; + } + if (options.contains("isovalue") && + options["isolvalue"].canConvert()) { + bool ok; + float iso = options["isovalue"].toFloat(&ok); + if (ok) + isoValue = iso; + } + + if (m_basis != nullptr) + homo = m_basis->homo(); + if (options.contains("index")) { + // check if options contains "homo" or "lumo" + bool string = options["index"].canConvert(); + if (string) { + // should be something like "homo-1" or "lumo+2" + QString name = options["index"].toString(); + QString expression, modifier; + if (name.contains("homo", Qt::CaseInsensitive)) { + index = homo; // modified by the expression after the string + expression = name.remove("homo", Qt::CaseInsensitive); + } else if (name.contains("lumo", Qt::CaseInsensitive)) { + index = homo + 1; // again modified by the expression + expression = name.remove("homo", Qt::CaseInsensitive); + } + // modify HOMO / LUMO based on "+ number" or "- number" + if (expression.contains('-')) { + modifier = expression.remove('-'); + bool ok; + int n = modifier.toInt(&ok); + if (ok) + index = index - n; + } else if (expression.contains('+')) { + modifier = expression.remove('+'); + bool ok; + int n = modifier.toInt(&ok); + if (ok) + index = index + n; + } + + } else + index = options.value("index").toInt(); + } + bool beta = false; + if (options.contains("spin")) { + beta = options["spin"].toString().contains("beta"); + } + + if (command.compare("renderVanDerWaals", Qt::CaseInsensitive) == 0) { + calculateEDT(VanDerWaals, cubeResolution); + return true; + } else if (command.compare("renderSolventAccessible", Qt::CaseInsensitive) == + 0) { + calculateEDT(SolventAccessible, cubeResolution); + return true; + } else if (command.compare("renderSolventExcluded", Qt::CaseInsensitive) == + 0) { + calculateEDT(SolventExcluded, cubeResolution); + return true; + } else if (command.compare("renderOrbital", Qt::CaseInsensitive) == 0) { + calculateQM(MolecularOrbital, index, beta, isoValue, cubeResolution); + return true; + } else if (command.compare("renderElectronDensity", Qt::CaseInsensitive) == + 0) { + calculateQM(ElectronDensity, index, beta, isoValue, cubeResolution); + return true; + } else if (command.compare("renderSpinDensity", Qt::CaseInsensitive) == 0) { + calculateQM(SpinDensity, index, beta, isoValue, cubeResolution); + return true; + } + + return false; +} + void Surfaces::setMolecule(QtGui::Molecule* mol) { if (mol->basisSet()) { @@ -130,14 +239,6 @@ QStringList Surfaces::menuPath(QAction*) const void Surfaces::surfacesActivated() { - if (!m_dialog) { - m_dialog = new SurfaceDialog(qobject_cast(parent())); - connect(m_dialog, SIGNAL(calculateClickedSignal()), - SLOT(calculateSurface())); - connect(m_dialog, SIGNAL(recordClicked()), SLOT(recordMovie())); - connect(m_dialog, SIGNAL(stepChanged(int)), SLOT(stepChanged(int))); - } - if (m_basis) { // we have quantum data, set up the dialog accordingly auto gaussian = dynamic_cast(m_basis); @@ -150,10 +251,12 @@ void Surfaces::surfacesActivated() m_dialog->setupBasis(m_basis->electronCount(), m_basis->molecularOrbitalCount(), beta); } + // only show cubes from files so we don't duplicate orbitals if (m_cubes.size() > 0) { QStringList cubeNames; for (auto* cube : m_cubes) { - cubeNames << cube->name().c_str(); + if (cube->cubeType() == Core::Cube::Type::FromFile) + cubeNames << cube->name().c_str(); } m_dialog->setupCubes(cubeNames); } @@ -170,20 +273,25 @@ void Surfaces::surfacesActivated() m_dialog->show(); } -float Surfaces::resolution() +float Surfaces::resolution(float specified) { - if (!m_dialog->automaticResolution()) + if (specified != 0.0) + return specified; + + if (m_dialog != nullptr && !m_dialog->automaticResolution()) return m_dialog->resolution(); float r = 0.02 * powf(m_molecule->atomCount(), 1.0f / 3.0f); float minimum = 0.05; float maximum = 0.5; - switch (m_dialog->surfaceType()) { - case SolventExcluded: - minimum = 0.1; - break; - default:; + if (m_dialog != nullptr) { + switch (m_dialog->surfaceType()) { + case SolventExcluded: + minimum = 0.1; + break; + default:; + } } r = std::max(minimum, std::min(maximum, r)); @@ -226,11 +334,14 @@ float inline square(float x) return x * x; } -void Surfaces::calculateEDT() +void Surfaces::calculateEDT(Type type, float defaultResolution) { + if (type == Unknown && m_dialog != nullptr) + type = m_dialog->surfaceType(); + QFuture future = QtConcurrent::run([=]() { double probeRadius = 0.0; - switch (m_dialog->surfaceType()) { + switch (type) { case VanDerWaals: break; case SolventAccessible: @@ -248,7 +359,7 @@ void Surfaces::calculateEDT() QtGui::RWLayerManager layerManager; for (size_t i = 0; i < m_molecule->atomCount(); i++) { if (!layerManager.visible(m_molecule->layer(i))) - continue; + continue; // ignore invisible atoms auto radius = Core::Elements::radiusVDW(m_molecule->atomicNumber(i)) + probeRadius; atoms->emplace_back(atomPositions[i], radius); @@ -257,10 +368,10 @@ void Surfaces::calculateEDT() } double padding = max_radius + probeRadius; - m_cube->setLimits(*m_molecule, resolution(), padding); + m_cube->setLimits(*m_molecule, resolution(defaultResolution), padding); m_cube->fill(-1.0); - const float res = resolution(); + const float res = resolution(defaultResolution); const Vector3 min = m_cube->min(); // then, for each atom, set cubes around it up to a certain radius @@ -300,7 +411,7 @@ void Surfaces::calculateEDT() }); // SolventExcluded requires an extra pass - if (m_dialog->surfaceType() == SolventExcluded) { + if (type == SolventExcluded) { m_performEDTStepWatcher.setFuture(future); } else { m_displayMeshWatcher.setFuture(future); @@ -369,11 +480,18 @@ void Surfaces::performEDTStep() m_displayMeshWatcher.setFuture(future); } -void Surfaces::calculateQM() +void Surfaces::calculateQM(Type type, int index, bool beta, float isoValue, + float defaultResolution) { - if (!m_basis || !m_dialog) + if (!m_basis) return; // nothing to do + if (m_dialog != nullptr) { + beta = m_dialog->beta(); // we're using the GUI + } + + // TODO: check if we already calculated the requested cube + // Reset state a little more frequently, minimal cost, avoid bugs. m_molecule->clearCubes(); m_molecule->clearMeshes(); @@ -409,28 +527,45 @@ void Surfaces::calculateQM() if (!m_cube) m_cube = m_molecule->addCube(); - Type type = m_dialog->surfaceType(); - int index = m_dialog->surfaceIndex(); - m_isoValue = m_dialog->isosurfaceValue(); - m_cube->setLimits(*m_molecule, resolution(), 5.0); + if (type == Unknown) + type = m_dialog->surfaceType(); + + if (index == -1 && m_dialog != nullptr) + index = m_dialog->surfaceIndex(); + if (isoValue == 0.0 && m_dialog != nullptr) + m_isoValue = m_dialog->isosurfaceValue(); + else + m_isoValue = isoValue; + float cubeResolution = resolution(defaultResolution); + float padding = 5.0; + // TODO: allow extra padding for some molecules / properties + m_cube->setLimits(*m_molecule, cubeResolution, padding); QString progressText; if (type == ElectronDensity) { progressText = tr("Calculating electron density"); m_cube->setName("Electron Density"); + m_cube->setCubeType(Core::Cube::Type::ElectronDensity); if (dynamic_cast(m_basis)) { m_gaussianConcurrent->calculateElectronDensity(m_cube); } else { m_slaterConcurrent->calculateElectronDensity(m_cube); } - } - - else if (type == MolecularOrbital) { + } else if (type == ElectronDensity) { + progressText = tr("Calculating spin density"); + m_cube->setName("Spin Density"); + m_cube->setCubeType(Core::Cube::Type::SpinDensity); + if (dynamic_cast(m_basis)) { + m_gaussianConcurrent->calculateSpinDensity(m_cube); + } else { + m_slaterConcurrent->calculateSpinDensity(m_cube); + } + } else if (type == MolecularOrbital) { progressText = tr("Calculating molecular orbital %L1").arg(index); m_cube->setName("Molecular Orbital " + std::to_string(index + 1)); + m_cube->setCubeType(Core::Cube::Type::MO); if (dynamic_cast(m_basis)) { - m_gaussianConcurrent->calculateMolecularOrbital(m_cube, index, - m_dialog->beta()); + m_gaussianConcurrent->calculateMolecularOrbital(m_cube, index, beta); } else { m_slaterConcurrent->calculateMolecularOrbital(m_cube, index); } @@ -471,14 +606,22 @@ void Surfaces::calculateQM() } } -void Surfaces::calculateCube() +void Surfaces::calculateCube(int index, float isoValue) { - if (!m_dialog || m_cubes.size() == 0) + if (m_cubes.size() == 0) + return; + + if (index == -1 && m_dialog != nullptr) + index = m_dialog->surfaceIndex(); + if (index < 0 || index >= static_cast(m_cubes.size())) return; // check bounds - m_cube = m_cubes[m_dialog->surfaceIndex()]; - m_isoValue = m_dialog->isosurfaceValue(); + m_cube = m_cubes[index]; + if (isoValue == 0.0 && m_dialog != nullptr) + m_isoValue = m_dialog->isosurfaceValue(); + else + m_isoValue = isoValue; displayMesh(); } @@ -507,7 +650,10 @@ void Surfaces::displayMesh() // qDebug() << " running displayMesh"; - m_smoothingPasses = m_dialog->smoothingPassesValue(); + if (m_dialog != nullptr) + m_smoothingPasses = m_dialog->smoothingPassesValue(); + else + m_smoothingPasses = 0; if (!m_mesh1) m_mesh1 = m_molecule->addMesh(); @@ -520,6 +666,9 @@ void Surfaces::displayMesh() // TODO - only do this if we're generating an orbital // and we need two meshes // How do we know? - likely ask the cube if it's an MO? + qDebug() << "Cube " << m_cube->name().c_str() << " type " + << m_cube->cubeType(); + if (!m_mesh2) m_mesh2 = m_molecule->addMesh(); if (!m_meshGenerator2) { diff --git a/avogadro/qtplugins/surfaces/surfaces.h b/avogadro/qtplugins/surfaces/surfaces.h index c40472da8..d25df78c1 100644 --- a/avogadro/qtplugins/surfaces/surfaces.h +++ b/avogadro/qtplugins/surfaces/surfaces.h @@ -28,7 +28,7 @@ namespace Core { class BasisSet; class Cube; class Mesh; -} +} // namespace Core namespace QtPlugins { @@ -61,7 +61,7 @@ class Surfaces : public QtGui::ExtensionPlugin FromFile, Unknown }; - + enum ColorProperty { None, @@ -69,7 +69,10 @@ class Surfaces : public QtGui::ExtensionPlugin }; QString name() const override { return tr("Surfaces"); } - QString description() const override { return tr("Read and render surfaces."); } + QString description() const override + { + return tr("Read and render surfaces."); + } QList actions() const override; @@ -77,19 +80,26 @@ class Surfaces : public QtGui::ExtensionPlugin void setMolecule(QtGui::Molecule* mol) override; + void registerCommands() override; + +public slots: + bool handleCommand(const QString& command, + const QVariantMap& options) override; + private slots: void surfacesActivated(); void calculateSurface(); - void calculateEDT(); + void calculateEDT(Type type = Unknown, float defaultResolution = 0.0); void performEDTStep(); // EDT step for SolventExcluded - void calculateQM(); - void calculateCube(); + void calculateQM(Type type = Unknown, int index = -1, bool betaSpin = false, + float isoValue = 0.0, float defaultResolution = 0.0); + void calculateCube(int index = -1, float isoValue = 0.0); void stepChanged(int); void displayMesh(); void meshFinished(); - + void colorMesh(); void colorMeshByPotential(); @@ -97,10 +107,9 @@ private slots: void movieFrame(); private: - float resolution(); - Core::Color3f chargeGradient( - double value, double clamp, tinycolormap::ColormapType colormap - ) const; + float resolution(float specified = 0.0); + Core::Color3f chargeGradient(double value, double clamp, + tinycolormap::ColormapType colormap) const; tinycolormap::ColormapType getColormapFromString(const QString& name) const; QList m_actions; @@ -139,7 +148,7 @@ private slots: class PIMPL; PIMPL* d = nullptr; }; -} -} +} // namespace QtPlugins +} // namespace Avogadro #endif // AVOGADRO_QTPLUGINS_QUANTUMOUTPUT_H From c5e074110e53a18bc6484998d6d9e37eb3c7f95d Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Wed, 13 Sep 2023 22:53:35 -0400 Subject: [PATCH 2/2] Fix some crashes found with testing A few remaining spots required the dialog to exist Signed-off-by: Geoff Hutchison --- avogadro/core/cube.h | 2 ++ avogadro/qtplugins/surfaces/surfaces.cpp | 40 +++++++++++++++++++----- 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/avogadro/core/cube.h b/avogadro/core/cube.h index 357087157..b62709fe9 100644 --- a/avogadro/core/cube.h +++ b/avogadro/core/cube.h @@ -39,6 +39,8 @@ class AVOGADROCORE_EXPORT Cube enum Type { VdW, + SolventAccessible, + SolventExcluded, ESP, ElectronDensity, SpinDensity, diff --git a/avogadro/qtplugins/surfaces/surfaces.cpp b/avogadro/qtplugins/surfaces/surfaces.cpp index 7b6010b21..ad503955a 100644 --- a/avogadro/qtplugins/surfaces/surfaces.cpp +++ b/avogadro/qtplugins/surfaces/surfaces.cpp @@ -99,7 +99,7 @@ Surfaces::Surfaces(QObject* p) : ExtensionPlugin(p), d(new PIMPL()) Surfaces::~Surfaces() { delete d; - delete m_cube; + // delete m_cube; // should be freed by the molecule } void Surfaces::registerCommands() @@ -113,6 +113,7 @@ void Surfaces::registerCommands() emit registerCommand("renderSolventExcluded", tr("Render the solvent-excluded molecular surface.")); emit registerCommand("renderOrbital", tr("Render a molecular orbital.")); + emit registerCommand("renderMO", tr("Render a molecular orbital.")); emit registerCommand("renderElectronDensity", tr("Render the electron density.")); emit registerCommand("renderSpinDensity", tr("Render the spin density.")); @@ -125,6 +126,8 @@ bool Surfaces::handleCommand(const QString& command, const QVariantMap& options) if (m_molecule == nullptr) return false; // No molecule to handle the command. + qDebug() << "handle surface cmd:" << command << options; + // Set up some defaults for the options. int index = -1; int homo = -1; @@ -148,12 +151,12 @@ bool Surfaces::handleCommand(const QString& command, const QVariantMap& options) if (m_basis != nullptr) homo = m_basis->homo(); - if (options.contains("index")) { + if (options.contains("orbital")) { // check if options contains "homo" or "lumo" - bool string = options["index"].canConvert(); + bool string = options["orbital"].canConvert(); if (string) { // should be something like "homo-1" or "lumo+2" - QString name = options["index"].toString(); + QString name = options["orbital"].toString(); QString expression, modifier; if (name.contains("homo", Qt::CaseInsensitive)) { index = homo; // modified by the expression after the string @@ -185,7 +188,8 @@ bool Surfaces::handleCommand(const QString& command, const QVariantMap& options) beta = options["spin"].toString().contains("beta"); } - if (command.compare("renderVanDerWaals", Qt::CaseInsensitive) == 0) { + if ((command.compare("renderVanDerWaals", Qt::CaseInsensitive) == 0) || + (command.compare("renderVDW", Qt::CaseInsensitive) == 0)) { calculateEDT(VanDerWaals, cubeResolution); return true; } else if (command.compare("renderSolventAccessible", Qt::CaseInsensitive) == @@ -196,7 +200,8 @@ bool Surfaces::handleCommand(const QString& command, const QVariantMap& options) 0) { calculateEDT(SolventExcluded, cubeResolution); return true; - } else if (command.compare("renderOrbital", Qt::CaseInsensitive) == 0) { + } else if ((command.compare("renderOrbital", Qt::CaseInsensitive) == 0) || + (command.compare("renderMO", Qt::CaseInsensitive) == 0)) { calculateQM(MolecularOrbital, index, beta, isoValue, cubeResolution); return true; } else if (command.compare("renderElectronDensity", Qt::CaseInsensitive) == @@ -239,6 +244,14 @@ QStringList Surfaces::menuPath(QAction*) const void Surfaces::surfacesActivated() { + if (!m_dialog) { + m_dialog = new SurfaceDialog(qobject_cast(parent())); + connect(m_dialog, SIGNAL(calculateClickedSignal()), + SLOT(calculateSurface())); + connect(m_dialog, SIGNAL(recordClicked()), SLOT(recordMovie())); + connect(m_dialog, SIGNAL(stepChanged(int)), SLOT(stepChanged(int))); + } + if (m_basis) { // we have quantum data, set up the dialog accordingly auto gaussian = dynamic_cast(m_basis); @@ -339,14 +352,20 @@ void Surfaces::calculateEDT(Type type, float defaultResolution) if (type == Unknown && m_dialog != nullptr) type = m_dialog->surfaceType(); + if (!m_cube) + m_cube = m_molecule->addCube(); + QFuture future = QtConcurrent::run([=]() { double probeRadius = 0.0; switch (type) { case VanDerWaals: + m_cube->setCubeType(Core::Cube::Type::VdW); break; case SolventAccessible: + m_cube->setCubeType(Core::Cube::Type::SolventAccessible); case SolventExcluded: probeRadius = 1.4; + m_cube->setCubeType(Core::Cube::Type::SolventExcluded); break; default: break; @@ -618,6 +637,9 @@ void Surfaces::calculateCube(int index, float isoValue) // check bounds m_cube = m_cubes[index]; + if (m_cube == nullptr) + return; + if (isoValue == 0.0 && m_dialog != nullptr) m_isoValue = m_dialog->isosurfaceValue(); else @@ -764,6 +786,9 @@ void Surfaces::colorMeshByPotential() void Surfaces::colorMesh() { + if (m_dialog == nullptr) + return; + switch (m_dialog->colorProperty()) { case None: break; @@ -784,7 +809,8 @@ void Surfaces::meshFinished() m_molecule->emitChanged(QtGui::Molecule::Added); movieFrame(); } else { - m_dialog->reenableCalculateButton(); + if (m_dialog != nullptr) + m_dialog->reenableCalculateButton(); qDebug() << " mesh finished";