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

Progressive refinement for surface #1692

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Progressive refinement for surface #1692

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
9215fb0
progressive-refinement-wip
perminder-17 Aug 23, 2024
ee85451
fixing things
perminder-17 Aug 26, 2024
13d844e
fixes
perminder-17 Aug 28, 2024
9279225
fixes
perminder-17 Aug 31, 2024
bf89942
Done
perminder-17 Aug 31, 2024
10a2dc0
minor fixes
perminder-17 Aug 31, 2024
d5961e8
minor fixes
perminder-17 Sep 1, 2024
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
304 changes: 155 additions & 149 deletions avogadro/qtplugins/surfaces/surfaces.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -332,6 +332,12 @@
return r;
}


float inline square(float x)
{
return x * x;
}

void Surfaces::calculateSurface()
{
if (!m_dialog)
Expand Down Expand Up @@ -363,99 +369,166 @@
}
}

float inline square(float x)
{
return x * x;
}

void Surfaces::calculateEDT(Type type, float defaultResolution)
{
if (type == Unknown && m_dialog != nullptr)
if (type == Unknown && m_dialog != nullptr){
type = m_dialog->surfaceType();
}
if (!m_cube){
m_cube = m_molecule->addCube();
}

auto* m_molecule2 = new QtGui::Molecule(*m_molecule);
m_cube2 = m_molecule2->addCube();

auto* m_temp = new Core::Cube(*m_cube);
m_cube = m_cube2;

// first pass low res.
QFuture futureLowRes = QtConcurrent::run([=]() {
calculateEDTpass(m_molecule2, m_cube, type, 0.5);
});
m_displayMeshWatcher.setFuture(futureLowRes);
m_molecule->emitChanged(QtGui::Molecule::Added);


m_cube = m_temp;
// // 2nd pass high res.
QFuture futureHighRes = QtConcurrent::run([=]() {
calculateEDTpass(m_molecule, m_cube, type, 0.1);
});

m_displayMeshWatcher.setFuture(futureHighRes);
m_molecule->emitChanged(QtGui::Molecule::Added);

}


void Surfaces::displayMesh()
{
if (!m_cube)
m_cube = m_molecule->addCube();
return;

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;
}
if (m_dialog != nullptr)
m_smoothingPasses = m_dialog->smoothingPassesValue();
else
m_smoothingPasses = 0;

if (!m_mesh1)
m_mesh1 = m_molecule->addMesh();
if (!m_meshGenerator1) {
m_meshGenerator1 = new QtGui::MeshGenerator;
connect(m_meshGenerator1, SIGNAL(finished()), SLOT(meshFinished()));
}
m_meshGenerator1->initialize(m_cube, m_mesh1, m_isoValue, m_smoothingPasses);

// first, make a list of all atom positions and radii
Array<Vector3> atomPositions = m_molecule->atomPositions3d();
auto* atoms = new std::vector<std::pair<Vector3, double>>();
double max_radius = probeRadius;
QtGui::RWLayerManager layerManager;
for (size_t i = 0; i < m_molecule->atomCount(); i++) {
if (!layerManager.visible(m_molecule->layer(i)))
continue; // ignore invisible atoms
auto radius =
Core::Elements::radiusVDW(m_molecule->atomicNumber(i)) + probeRadius;
atoms->emplace_back(atomPositions[i], radius);
if (radius > max_radius)
max_radius = radius;
// bool isMO = false;
// // if it's from a file we should "play it safe"
// if (m_cube->cubeType() == Cube::Type::MO ||
// m_cube->cubeType() == Cube::Type::FromFile) {
// isMO = true;
// }
Fixed Show fixed Hide fixed

// if (isMO) {
Fixed Show fixed Hide fixed
if (!m_mesh2)
m_mesh2 = m_molecule->addMesh();
if (!m_meshGenerator2) {
m_meshGenerator2 = new QtGui::MeshGenerator;
connect(m_meshGenerator2, SIGNAL(finished()), SLOT(meshFinished()));
}
m_meshGenerator2->initialize(m_cube2, m_mesh2, -m_isoValue,
m_smoothingPasses, true);
// }
Fixed Show fixed Hide fixed

double padding = max_radius + probeRadius;
m_cube->setLimits(*m_molecule, resolution(defaultResolution), padding);
m_cube->fill(-1.0);

const float res = resolution(defaultResolution);
const Vector3 min = m_cube->min();

// then, for each atom, set cubes around it up to a certain radius
QFuture innerFuture =
QtConcurrent::map(*atoms, [=](std::pair<Vector3, double>& in) {
double startPosX = in.first(0) - in.second;
double endPosX = in.first(0) + in.second;
int startIndexX = (startPosX - min(0)) / res;
int endIndexX = (endPosX - min(0)) / res + 1;
for (int indexX = startIndexX; indexX < endIndexX; indexX++) {
double posX = indexX * res + min(0);
double radiusXsq = square(in.second) - square(posX - in.first(0));
if (radiusXsq < 0.0)
continue;
double radiusX = sqrt(radiusXsq);
double startPosY = in.first(1) - radiusX;
double endPosY = in.first(1) + radiusX;
int startIndexY = (startPosY - min(1)) / res;
int endIndexY = (endPosY - min(1)) / res + 1;
for (int indexY = startIndexY; indexY < endIndexY; indexY++) {
double posY = indexY * res + min(1);
double lengthXYsq = square(radiusX) - square(posY - in.first(1));
if (lengthXYsq < 0.0)
continue;
double lengthXY = sqrt(lengthXYsq);
double startPosZ = in.first(2) - lengthXY;
double endPosZ = in.first(2) + lengthXY;
int startIndexZ = (startPosZ - min(2)) / res;
int endIndexZ = (endPosZ - min(2)) / res + 1;
m_cube->fillStripe(indexX, indexY, startIndexZ, endIndexZ - 1,
1.0f);
}
}
});
// Start the mesh generation - this needs an improved mutex with a read lock
// to function as expected. Write locks are exclusive, read locks can have
// many read locks but no write lock.
m_meshGenerator1->start();
// if (isMO)
m_meshGenerator2->start();

innerFuture.waitForFinished();
});
// Track how many meshes are left to show.
// if (isMO)
m_meshesLeft = 2;
// else
// m_meshesLeft = 1;
}

// SolventExcluded requires an extra pass
if (type == SolventExcluded) {
m_performEDTStepWatcher.setFuture(future);
} else {
m_displayMeshWatcher.setFuture(future);

void Surfaces::calculateEDTpass(QtGui::Molecule* mol, Core::Cube* cube, Type type, float defaultResolution)
{
double probeRadius = 0.0;
switch (type) {
case VanDerWaals:
cube->setCubeType(Core::Cube::Type::VdW);
break;
case SolventAccessible:
probeRadius = 1.4;
cube->setCubeType(Core::Cube::Type::SolventAccessible);
break;
case SolventExcluded:
probeRadius = 1.4;
cube->setCubeType(Core::Cube::Type::SolventExcluded);
break;
default:
break;
}

// first, make a list of all atom positions and radii
Array<Vector3> atomPositions = mol->atomPositions3d();
auto* atoms = new std::vector<std::pair<Vector3, double>>();
double max_radius = probeRadius;
QtGui::RWLayerManager layerManager;
for (size_t i = 0; i < mol->atomCount(); i++) {
if (!layerManager.visible(mol->layer(i)))
continue; // ignore invisible atoms
auto radius =
Core::Elements::radiusVDW(mol->atomicNumber(i)) + probeRadius;
atoms->emplace_back(atomPositions[i], radius);
if (radius > max_radius)
max_radius = radius;
}

double padding = max_radius + probeRadius;

const float res = resolution(defaultResolution);
cube->setLimits(*mol, res, padding);
cube->fill(-1.0);
const Vector3 min = cube->min();

// then, for each atom, set cubes around it up to a certain radius
QFuture innerFuture =
QtConcurrent::map(*atoms, [=](std::pair<Vector3, double>& in) {
double startPosX = in.first(0) - in.second;
double endPosX = in.first(0) + in.second;
int startIndexX = (startPosX - min(0)) / res;
int endIndexX = (endPosX - min(0)) / res + 1;
for (int indexX = startIndexX; indexX < endIndexX; indexX++) {
double posX = indexX * res + min(0);
double radiusXsq = square(in.second) - square(posX - in.first(0));
if (radiusXsq < 0.0)
continue;
double radiusX = sqrt(radiusXsq);
double startPosY = in.first(1) - radiusX;
double endPosY = in.first(1) + radiusX;
int startIndexY = (startPosY - min(1)) / res;
int endIndexY = (endPosY - min(1)) / res + 1;
for (int indexY = startIndexY; indexY < endIndexY; indexY++) {
double posY = indexY * res + min(1);
double lengthXYsq = square(radiusX) - square(posY - in.first(1));
if (lengthXYsq < 0.0)
continue;
double lengthXY = sqrt(lengthXYsq);
double startPosZ = in.first(2) - lengthXY;
double endPosZ = in.first(2) + lengthXY;
int startIndexZ = (startPosZ - min(2)) / res;
int endIndexZ = (endPosZ - min(2)) / res + 1;
cube->fillStripe(indexX, indexY, startIndexZ, endIndexZ - 1,
1.0f);
}
}
});
innerFuture.waitForFinished();
}

void Surfaces::performEDTStep()
Expand Down Expand Up @@ -686,56 +759,6 @@
}
}

void Surfaces::displayMesh()
{
if (!m_cube)
return;

if (m_dialog != nullptr)
m_smoothingPasses = m_dialog->smoothingPassesValue();
else
m_smoothingPasses = 0;

if (!m_mesh1)
m_mesh1 = m_molecule->addMesh();
if (!m_meshGenerator1) {
m_meshGenerator1 = new QtGui::MeshGenerator;
connect(m_meshGenerator1, SIGNAL(finished()), SLOT(meshFinished()));
}
m_meshGenerator1->initialize(m_cube, m_mesh1, m_isoValue, m_smoothingPasses);

bool isMO = false;
// if it's from a file we should "play it safe"
if (m_cube->cubeType() == Cube::Type::MO ||
m_cube->cubeType() == Cube::Type::FromFile) {
isMO = true;
}

if (isMO) {
if (!m_mesh2)
m_mesh2 = m_molecule->addMesh();
if (!m_meshGenerator2) {
m_meshGenerator2 = new QtGui::MeshGenerator;
connect(m_meshGenerator2, SIGNAL(finished()), SLOT(meshFinished()));
}
m_meshGenerator2->initialize(m_cube, m_mesh2, -m_isoValue,
m_smoothingPasses, true);
}

// Start the mesh generation - this needs an improved mutex with a read lock
// to function as expected. Write locks are exclusive, read locks can have
// many read locks but no write lock.
m_meshGenerator1->start();
if (isMO)
m_meshGenerator2->start();

// Track how many meshes are left to show.
if (isMO)
m_meshesLeft = 2;
else
m_meshesLeft = 1;
}

Core::Color3f Surfaces::chargeGradient(double value, double clamp,
ColormapType colormap) const
{
Expand Down Expand Up @@ -829,27 +852,10 @@

void Surfaces::meshFinished()
{
--m_meshesLeft;
if (m_meshesLeft == 0) {
colorMesh();

// finished, so request to enable the mesh display type
QStringList displayTypes;
displayTypes << tr("Meshes");
requestActiveDisplayTypes(displayTypes);

if (m_recordingMovie) {
// Move to the next frame.
qDebug() << "Let's get to the next frame…";
m_molecule->emitChanged(QtGui::Molecule::Added);
movieFrame();
} else {

if (m_dialog != nullptr)
m_dialog->reenableCalculateButton();

m_molecule->emitChanged(QtGui::Molecule::Added);
}
}
}

void Surfaces::recordMovie()
Expand Down
10 changes: 10 additions & 0 deletions avogadro/qtplugins/surfaces/surfaces.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,8 @@
void surfacesActivated();
void calculateSurface();
void calculateEDT(Type type = Unknown, float defaultResolution = 0.0);
// void meshDisplaying(Type type = Unknown);
Fixed Show fixed Hide fixed
void calculateEDTpass(QtGui::Molecule* mol, Core::Cube* cube, Type type = Unknown, float defaultResolution = 0.0);
void performEDTStep(); // EDT step for SolventExcluded
void calculateQM(Type type = Unknown, int index = -1, bool betaSpin = false,
float isoValue = 0.0, float defaultResolution = 0.0);
Expand All @@ -117,17 +119,25 @@
QProgressDialog* m_progressDialog = nullptr;

QtGui::Molecule* m_molecule = nullptr;
// QtGui::Molecule* m_molecule2 = nullptr;
Fixed Show fixed Hide fixed

Core::BasisSet* m_basis = nullptr;

GaussianSetConcurrent* m_gaussianConcurrent = nullptr;
SlaterSetConcurrent* m_slaterConcurrent = nullptr;

Core::Cube* m_cube = nullptr;
Core::Cube* m_cube2 = nullptr;
// Core::Cube* m_temp = nullptr;

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

std::vector<Core::Cube*> m_cubes;
/* One QFutureWatcher per asynchronous slot function, e.g.:*/
/* calculateEDT() -> [performEDTStep()] -> displayMesh() */
QFutureWatcher<void> m_performEDTStepWatcher;
QFutureWatcher<void> m_displayMeshWatcher;
QFutureWatcher<void> m_pataNahi;
QFutureWatcher<void> m_displayMeshWatcherLowRes;
QFutureWatcher<void> m_displayMeshWatcherHighRes;
Core::Mesh* m_mesh1 = nullptr;
Core::Mesh* m_mesh2 = nullptr;
/* displayMesh() -> meshFinished() */
Expand Down
Loading