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

Plot UV and CD spectra #1872

Closed
wants to merge 2 commits into from
Closed

Plot UV and CD spectra #1872

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
fdd4e50
Make sure to store electronic and CD spectra
ghutchis Dec 19, 2024
8e4ea5a
Initial update for electronic and CD spectra plots
ghutchis Dec 19, 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
130 changes: 88 additions & 42 deletions avogadro/qtplugins/spectra/spectradialog.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,8 @@ SpectraDialog::SpectraDialog(QWidget* parent)
m_ui->dataTable->hide();
m_ui->push_exportData->hide();

readSettings();

// connections for options
connect(m_ui->push_options, SIGNAL(clicked()), this, SLOT(toggleOptions()));
connect(m_ui->push_colorBackground, SIGNAL(clicked()), this,
Expand All @@ -77,6 +79,17 @@ SpectraDialog::SpectraDialog(QWidget* parent)
SLOT(changeFontSize()));
connect(m_ui->lineWidthSpinBox, SIGNAL(valueChanged(double)), this,
SLOT(changeLineWidth()));
connectOptions();
}

SpectraDialog::~SpectraDialog()
{
writeSettings();
}

void SpectraDialog::connectOptions()
{
// connect (or reconnect) anything that calls change or update plot
connect(m_ui->combo_spectra, SIGNAL(currentIndexChanged(int)), this,
SLOT(changeSpectra()));
connect(m_ui->xAxisMinimum, SIGNAL(valueChanged(double)), this,
Expand All @@ -89,20 +102,32 @@ SpectraDialog::SpectraDialog(QWidget* parent)
SLOT(updatePlot()));
connect(m_ui->peakWidth, SIGNAL(valueChanged(double)), this,
SLOT(updatePlot()));

readSettings();
}

SpectraDialog::~SpectraDialog()
void SpectraDialog::disconnectOptions()
{
writeSettings();
// disconnect anything that calls change or update plot
disconnect(m_ui->combo_spectra, SIGNAL(currentIndexChanged(int)), this,
SLOT(changeSpectra()));
disconnect(m_ui->xAxisMinimum, SIGNAL(valueChanged(double)), this,
SLOT(updatePlot()));
disconnect(m_ui->xAxisMaximum, SIGNAL(valueChanged(double)), this,
SLOT(updatePlot()));
disconnect(m_ui->yAxisMinimum, SIGNAL(valueChanged(double)), this,
SLOT(updatePlot()));
disconnect(m_ui->yAxisMaximum, SIGNAL(valueChanged(double)), this,
SLOT(updatePlot()));
disconnect(m_ui->peakWidth, SIGNAL(valueChanged(double)), this,
SLOT(updatePlot()));
}

void SpectraDialog::changeSpectra()
{
// TODO: change the scale and offset based on defaults and settings
QSettings settings;

disconnectOptions();

// what type of spectra are we plotting?
SpectraType type =
static_cast<SpectraType>(m_ui->combo_spectra->currentData().toInt());
Expand Down Expand Up @@ -154,48 +179,22 @@ void SpectraDialog::changeSpectra()
break;
}

MatrixX& spectra =
m_spectra[m_ui->combo_spectra->currentText().toStdString()];
float maxIntensity = 1.0;
// update the data table
m_ui->dataTable->setRowCount(spectra.rows());
m_ui->dataTable->setColumnCount(spectra.cols());
for (auto i = 0; i < spectra.rows(); ++i) {
for (auto j = 0; j < spectra.cols(); ++j) {
QTableWidgetItem* item =
new QTableWidgetItem(QString::number(spectra(i, j), 'f', 4));
m_ui->dataTable->setItem(i, j, item);
}
}
// if there's a second column, check for intensities
if (spectra.cols() > 1) {
for (auto i = 0; i < spectra.rows(); ++i) {
if (spectra(i, 1) > maxIntensity)
maxIntensity = spectra(i, 1);
}
maxIntensity = maxIntensity * 1.25;
}
// if transmission for IR, set the max intensity to 100
if (type == SpectraType::Infrared)
maxIntensity = 100.0;

if (maxIntensity < 1.0)
maxIntensity = 1.0;

// update the spin box
m_ui->yAxisMaximum->setValue(maxIntensity);

updatePlot();
connectOptions();
}

void SpectraDialog::setSpectra(const std::map<std::string, MatrixX>& spectra)
{
m_spectra = spectra;

// update the combo box
disconnect(m_ui->combo_spectra, SIGNAL(currentIndexChanged(int)), this,
SLOT(changeSpectra()));

m_ui->combo_spectra->clear();
for (auto& spectra : m_spectra) {
QString name = QString::fromStdString(spectra.first);

if (name == "IR") {
name = tr("Infrared");
m_ui->combo_spectra->addItem(name,
Expand All @@ -222,7 +221,9 @@ void SpectraDialog::setSpectra(const std::map<std::string, MatrixX>& spectra)
}

changeSpectra();
updatePlot();
// connect again
connect(m_ui->combo_spectra, SIGNAL(currentIndexChanged(int)), this,
SLOT(changeSpectra()));
}

void SpectraDialog::writeSettings() const
Expand Down Expand Up @@ -319,6 +320,9 @@ void SpectraDialog::changeLineWidth()

void SpectraDialog::updatePlot()
{
// disconnect the options while we update the plot
disconnectOptions();

// the raw data
std::vector<double> transitions, intensities;
// for the plot
Expand Down Expand Up @@ -392,7 +396,11 @@ void SpectraDialog::updatePlot()
break;
case SpectraType::CircularDichroism:
transitions = fromMatrix(m_spectra["Electronic"].col(0));
intensities = fromMatrix(m_spectra["Electronic"].col(2));
// check if electronic has a third column
if (m_spectra["Electronic"].cols() > 2)
intensities = fromMatrix(m_spectra["Electronic"].col(2));
else // grab it from the CD data
intensities = fromMatrix(m_spectra["CircularDichroism"].col(1));
windowName = tr("Circular Dichroism Spectra");
xTitle = tr("eV)");
yTitle = tr("Intensity");
Expand All @@ -419,12 +427,41 @@ void SpectraDialog::updatePlot()
}
setWindowTitle(windowName);

// update the data table
m_ui->dataTable->setRowCount(transitions.size());
m_ui->dataTable->setColumnCount(2);
for (auto i = 0; i < transitions.size(); ++i) {
// frequency or energy
QTableWidgetItem* item =
new QTableWidgetItem(QString::number(transitions[i], 'f', 4));
m_ui->dataTable->setItem(i, 0, item);
// intensities
item = new QTableWidgetItem(QString::number(intensities[i], 'f', 4));
m_ui->dataTable->setItem(i, 1, item);
}

double maxIntensity = 0.0f;
for (auto intensity : intensities) {
if (intensity > maxIntensity)
maxIntensity = intensity;
}

// if transmission for IR, set the max intensity to 100
if (type == SpectraType::Infrared)
maxIntensity = 100.0;

if (maxIntensity < 1.0)
maxIntensity = 1.0;

// update the spin boxes
m_ui->yAxisMaximum->setValue(maxIntensity);
m_ui->yAxisMinimum->setMinimum(0.0);
// if CD, set the minimum too
if (type == SpectraType::CircularDichroism) {
m_ui->yAxisMinimum->setMinimum(-maxIntensity * 2.0);
m_ui->yAxisMinimum->setValue(-maxIntensity);
}

// now compose the plot data
float scale = m_ui->scaleSpinBox->value();
float offset = m_ui->offsetSpinBox->value();
Expand All @@ -433,11 +470,17 @@ void SpectraDialog::updatePlot()
float xMin = m_ui->xAxisMinimum->value();
float xMax = m_ui->xAxisMaximum->value();

int start = std::min(static_cast<int>(xMin), static_cast<int>(xMax));
int end = std::max(static_cast<int>(xMin), static_cast<int>(xMax));
float start = std::min(xMin, xMax);
float end = std::max(xMin, xMax);
// for some spectra, we need to take small steps, so we scale the x axis
float xScale = 1.0;
if (type == SpectraType::Electronic || type == SpectraType::CircularDichroism)
xScale = 1.0f / 0.05f;

float stickWidth = fwhm / (xScale * 30.0);

for (unsigned int x = start; x < end; ++x) {
float xValue = static_cast<float>(x);
for (unsigned int x = round(start * xScale); x < round(end * xScale); ++x) {
float xValue = static_cast<float>(x) / xScale;
xData.push_back(xValue);
yData.push_back(0.0f);
yStick.push_back(0.0f);
Expand All @@ -448,7 +491,7 @@ void SpectraDialog::updatePlot()
float peak = intensities[index];

float intensity = scaleAndBlur(xValue, freq, peak, scale, offset, fwhm);
float stick = scaleAndBlur(xValue, freq, peak, scale, offset, 1.0);
float stick = scaleAndBlur(xValue, freq, peak, scale, offset, stickWidth);

yData.back() += intensity;
yStick.back() += stick;
Expand Down Expand Up @@ -492,6 +535,9 @@ void SpectraDialog::updatePlot()

chart->setXAxisLimits(xAxisMin, xAxisMax);
chart->setYAxisLimits(yAxisMin, yAxisMax);

// re-enable the options
connectOptions();
}

VTK::ChartWidget* SpectraDialog::chartWidget()
Expand Down
3 changes: 3 additions & 0 deletions avogadro/qtplugins/spectra/spectradialog.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,9 @@ class SpectraDialog : public QDialog

VTK::ChartWidget* chartWidget();

void disconnectOptions();
void connectOptions();

private slots:
void changeBackgroundColor();
void changeForegroundColor();
Expand Down
3 changes: 3 additions & 0 deletions avogadro/qtplugins/spectra/spectradialog.ui
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,9 @@
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
<property name="minimum">
<double>-4000.000000000000000</double>
</property>
<property name="maximum">
<double>4000.000000000000000</double>
</property>
Expand Down
39 changes: 20 additions & 19 deletions avogadro/quantumio/orca.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,23 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
return false;
}

// this should be the final coordinate set (e.g. the optimized geometry)
molecule.setCoordinate3d(molecule.atomPositions3d(), 0);
if (m_coordSets.size() > 1) {
for (unsigned int i = 0; i < m_coordSets.size(); i++) {
Array<Vector3> positions;
positions.reserve(molecule.atomCount());
for (size_t j = 0; j < molecule.atomCount(); ++j) {
positions.push_back(m_coordSets[i][j] * BOHR_TO_ANGSTROM);
}
molecule.setCoordinate3d(positions, i + 1);
}
}

// guess bonds and bond orders
molecule.perceiveBondsSimple();
molecule.perceiveBondOrders();

if (m_frequencies.size() > 0 &&
m_frequencies.size() == m_vibDisplacements.size() &&
m_frequencies.size() == m_IRintensities.size()) {
Expand Down Expand Up @@ -98,23 +115,6 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
molecule.setSpectra("NMR", nmrData);
}

// this should be the final coordinate set (e.g. the optimized geometry)
molecule.setCoordinate3d(molecule.atomPositions3d(), 0);
if (m_coordSets.size() > 1) {
for (unsigned int i = 0; i < m_coordSets.size(); i++) {
Array<Vector3> positions;
positions.reserve(molecule.atomCount());
for (size_t j = 0; j < molecule.atomCount(); ++j) {
positions.push_back(m_coordSets[i][j] * BOHR_TO_ANGSTROM);
}
molecule.setCoordinate3d(positions, i + 1);
}
}

// guess bonds and bond orders
molecule.perceiveBondsSimple();
molecule.perceiveBondOrders();

// check bonds from calculated bond orders
if (m_bondOrders.size() > 0) {
for (unsigned int i = 0; i < m_bondOrders.size(); i++) {
Expand Down Expand Up @@ -251,7 +251,8 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)
getline(in, key); // skip header
}
// starts at the next line
} else if (Core::contains(key, "CD SPECTRUM")) {
} else if (Core::contains(key, "CD SPECTRUM") &&
!Core::contains(key, "TRANSITION VELOCITY DIPOLE")) {
m_currentMode = ECD;
for (int i = 0; i < 4; ++i) {
getline(in, key); // skip header
Expand Down Expand Up @@ -697,7 +698,7 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)

wavenumbers = Core::lexicalCast<double>(list[1]);
// convert to eV
m_electronicTransitions.push_back(wavenumbers / 8065.544);
// m_electronicTransitions.push_back(wavenumbers / 8065.544);
m_electronicRotations.push_back(Core::lexicalCast<double>(list[3]));

getline(in, key);
Expand Down
Loading