Skip to content

Commit

Permalink
Merge pull request #1344 from ghutchis/add-surface-scripting
Browse files Browse the repository at this point in the history 
Initial support for scripting surface generation
  • Loading branch information
@ghutchis
ghutchis authored Sep 14, 2023
2 parents ab9e39a + c5e0741 commit 4f37ad4
Show file tree
Hide file tree
Showing 3 changed files with 232 additions and 45 deletions.
3 changes: 3 additions & 0 deletions avogadro/core/cube.h
View file
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
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,121 @@ 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()
{
// 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 @@ 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);
}
Expand All @@ -170,20 +286,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));
Expand Down Expand Up @@ -226,16 +347,25 @@ 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();

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 @@ 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);
Expand All @@ -257,10 +387,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
Expand Down Expand Up @@ -300,7 +430,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);
Expand Down Expand Up @@ -369,11 +499,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();
Expand Down Expand Up @@ -409,28 +546,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<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::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<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 @@ 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();
Expand All @@ -520,6 +688,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) {
Expand Down Expand Up @@ -615,6 +786,9 @@ void Surfaces::colorMeshByPotential()

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

switch (m_dialog->colorProperty()) {
case None:
break;
Expand All @@ -635,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";

Expand Down
Loading

0 comments on commit 4f37ad4

Please sign in to comment.