Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Initial support for scripting surface generation #1344

Merged
merged 2 commits into from
Sep 14, 2023
Merged

Initial support for scripting surface generation #1344

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c944183
Initial support for scripting surface generation
ghutchis Sep 9, 2023
c5e0741
Fix some crashes found with testing
ghutchis Sep 14, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions avogadro/core/cube.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,11 @@ class AVOGADROCORE_EXPORT Cube
enum Type
{
VdW,
SolventAccessible,
SolventExcluded,
ESP,
ElectronDensity,
SpinDensity,
MO,
FromFile,
None
Expand Down
239 changes: 207 additions & 32 deletions avogadro/qtplugins/surfaces/surfaces.cpp
View file
Open in desktop
ghutchis marked this conversation as resolved.
Show resolved Hide resolved
ghutchis marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,121 @@
Surfaces::~Surfaces()
{
delete d;
delete m_cube;
// delete m_cube; // should be freed by the molecule
}

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("renderMO", 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.

qDebug() << "handle surface cmd:" << command << options;

// 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<float>()) {
bool ok;
float res = options["resolution"].toFloat(&ok);
if (ok)
cubeResolution = res;
}
if (options.contains("isovalue") &&
options["isolvalue"].canConvert<float>()) {
bool ok;
float iso = options["isovalue"].toFloat(&ok);
if (ok)
isoValue = iso;
}

if (m_basis != nullptr)
homo = m_basis->homo();
if (options.contains("orbital")) {
// check if options contains "homo" or "lumo"
bool string = options["orbital"].canConvert<QString>();
if (string) {
// should be something like "homo-1" or "lumo+2"
QString name = options["orbital"].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) ||
(command.compare("renderVDW", 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) ||
(command.compare("renderMO", 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)
Expand Down Expand Up @@ -150,10 +264,12 @@
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);
}
Expand All @@ -170,20 +286,25 @@
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:;
}
Comment on lines +302 to +307

Check notice

Code scanning / CodeQL

No trivial switch statements Note

This switch statement should either handle more cases, or be rewritten as an if statement.
}

r = std::max(minimum, std::min(maximum, r));
Expand Down Expand Up @@ -226,16 +347,25 @@
return x * x;
}

void Surfaces::calculateEDT()
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 (m_dialog->surfaceType()) {
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;
Expand All @@ -248,7 +378,7 @@
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);
Expand All @@ -257,10 +387,10 @@
}

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
Expand Down Expand Up @@ -300,7 +430,7 @@
});

// SolventExcluded requires an extra pass
if (m_dialog->surfaceType() == SolventExcluded) {
if (type == SolventExcluded) {
m_performEDTStepWatcher.setFuture(future);
} else {
m_displayMeshWatcher.setFuture(future);
Expand Down Expand Up @@ -369,11 +499,18 @@
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();
Expand Down Expand Up @@ -409,28 +546,45 @@
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<GaussianSet*>(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<GaussianSet*>(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<GaussianSet*>(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);
}
Expand Down Expand Up @@ -471,14 +625,25 @@
}
}

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<int>(m_cubes.size()))
return;

// check bounds
m_cube = m_cubes[m_dialog->surfaceIndex()];
m_isoValue = m_dialog->isosurfaceValue();
m_cube = m_cubes[index];
if (m_cube == nullptr)
return;

if (isoValue == 0.0 && m_dialog != nullptr)
m_isoValue = m_dialog->isosurfaceValue();
else
m_isoValue = isoValue;
displayMesh();
}

Expand Down Expand Up @@ -507,7 +672,10 @@

// 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();
Expand All @@ -520,6 +688,9 @@
// 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) {
Expand Down Expand Up @@ -615,6 +786,9 @@

void Surfaces::colorMesh()
{
if (m_dialog == nullptr)
return;

switch (m_dialog->colorProperty()) {
case None:
break;
Expand All @@ -635,7 +809,8 @@
m_molecule->emitChanged(QtGui::Molecule::Added);
movieFrame();
} else {
m_dialog->reenableCalculateButton();
if (m_dialog != nullptr)
m_dialog->reenableCalculateButton();

qDebug() << " mesh finished";

Expand Down
Loading
Loading