Hi
I’m trying to display CP2K vibration calculation result on Avogadro 1.
It is apparently inadequate only to do OBMol::SetData( OBVibrationData* ) and m_molecule->setOBMol(OBMol*) for this purpose.
This led only to display molecular geometry, but not vibration data.
Any help to display vibration would be appreciated.
Aoyama Iwao
The question is whether Open Babel supports reading CP2K vibrations. I don’t believe that’s the case yet.
Since Avogadro (both versions) use Babel to read data, if it’s not supported, Avogadro can’t visualize them.
Thank you so much for quick reply.
I’m sorry for the previous inadequate explanation.
Although my final goal is to visualize CP2K calculation result, I’ve only executed a preliminary test code as follows,
using namespace OpenBabel;
// Create a methane instance
OBMol* methane = new OBMol;
methane->BeginModify();
OBAtom a;
a.Clear();
a.SetIdx(1);
a.SetAtomicNum(6);
a.SetVector(-1.211371, 1.665435, 0.000000);
methane->AddAtom(a);
a.Clear();
a.SetIdx(2);
a.SetAtomicNum(1);
a.SetVector(0.810636, 1.665435, 0.000000);
methane->AddAtom(a);
a.Clear();
a.SetIdx(3);
a.SetAtomicNum(1);
a.SetVector(-1.885361, 3.223136, 1.099008);
methane->AddAtom(a);
a.Clear();
a.SetIdx(4);
a.SetAtomicNum(1);
a.SetVector(-1.885361, -0.065177, 0.799505);
methane->AddAtom(a);
a.Clear();
a.SetIdx(5);
a.SetAtomicNum(1);
a.SetVector(-1.885380, 1.838363, -1.898513);
methane->AddAtom(a);
methane->AddBond(1,2,1,0,-1);
methane->AddBond(1,3,1,0,-1);
methane->AddBond(1,4,1,0,-1);
methane->AddBond(1,5,1,0,-1);
// Set vibration data
std::vector< std::vector< vector3 > > vLx; // Normal modes.
uint nmode = 3; // Number of normal modes. Only 3 modes for simplicity.
uint natom = 5; // Number of atoms.
std::vector< double > vF; // Frequencies
std::vector< double > vI; // IR intensities
for (uint i=0;i<nmode; i++)
vLx.push_back(std::vector<vector3>(natom));
// 1st at() is iMode, 2nd at() is iAtom
vLx.at(0).at(0) = vector3(0.116389, -0.014994, -0.003072);
vLx.at(0).at(1) = vector3(0.148184, 0.076496, 0.015665);
vLx.at(0).at(2) = vector3(-0.483241, -0.178965, -0.164553);
vLx.at(0).at(3) = vector3(-0.536029, 0.224786, -0.052825);
vLx.at(0).at(4) = vector3(-0.514727, 0.056218, 0.238296);
vF.push_back( 1301.397928 );
vI.push_back( 0.000788 );
vLx.at(1).at(0) = vector3(0.005296, 0.017344, 0.115982);
vLx.at(1).at(1) = vector3(0.006741, -0.088488, -0.591676);
vLx.at(1).at(2) = vector3(-0.185609, 0.278061, -0.333112);
vLx.at(1).at(3) = vector3(-0.089159, -0.248082, -0.517980);
vLx.at(1).at(4) = vector3(0.204967, -0.147996, 0.061801);
vF.push_back( 1301.429171 );
vI.push_back( 0.000788 );
vLx.at(2).at(0) = vector3(-0.014361, -0.115130, 0.017872);
vLx.at(2).at(1) = vector3(-0.018284, 0.587355, -0.091164);
vLx.at(2).at(2) = vector3(0.230925, 0.223000, -0.348249);
vLx.at(2).at(3) = vector3(-0.161312, -0.014971, 0.187069);
vLx.at(2).at(4) = vector3(0.119660, 0.575441, 0.039548);
vF.push_back( 1301.437251 );
vI.push_back( 0.000788 );
OBVibrationData* vd = new OBVibrationData;
vd->SetData( vLx, vF, vI );
methane->SetData(vd);
methane->EndModify();
// Set methane to Avogadro::Molecule* m_molecule
m_molecule->setOBMol(methane);
// Show the vibration dock if it is hidden
foreach (QWidget *widget, QApplication::allWidgets())
{
if (widget->objectName().contains("vibrationDock",Qt::CaseInsensitive))
{
if (!widget->isVisible())
widget->show();
}
}
This code displays only methane geometry, but dose not vibartion data on Avogadro.
I’m afraid that this lacks some code (such as manipulating the vibration dock widget) for displaying vibration.
The used values were taken from a CP2K calculation result for methane, printed as Molden format.:
[Molden Format]
[FREQ]
1301.397928
1301.429171
1301.437251
1507.385663
1507.390295
2809.466776
2969.262023
2969.269193
2969.273286
[FR-COORD]
C -1.211371 1.665435 0.000000
H 0.810636 1.665435 0.000000
H -1.885361 3.223136 1.099008
H -1.885361 -0.065177 0.799505
H -1.885380 1.838363 -1.898513
[FR-NORM-COORD]
vibration 1
0.116389 -0.014994 -0.003072
0.148184 0.076496 0.015665
-0.483241 -0.178965 -0.164553
-0.536029 0.224786 -0.052825
-0.514727 0.056218 0.238296
vibration 2
0.005296 0.017344 0.115982
0.006741 -0.088488 -0.591676
-0.185609 0.278061 -0.333112
-0.089159 -0.248082 -0.517980
0.204967 -0.147996 0.061801
vibration 3
-0.014361 -0.115130 0.017872
-0.018284 0.587355 -0.091164
0.230925 0.223000 -0.348249
-0.161312 -0.014971 0.187069
0.119660 0.575441 0.039548
vibration 4
0.000002 0.000001 -0.000004
0.000003 -0.154951 0.475411
0.139013 0.315571 -0.362033
0.320580 -0.256618 -0.285237
-0.459616 0.095983 0.171911
vibration 5
-0.000000 -0.000001 -0.000000
-0.000001 -0.475379 -0.154944
-0.450437 -0.045120 -0.212295
0.345613 0.031676 0.359932
0.104828 0.488839 0.007310
vibration 6
-0.000001 0.000000 0.000002
-0.499993 -0.000000 -0.000000
0.166670 -0.385195 -0.271768
0.166670 0.427953 -0.197705
0.166665 -0.042760 0.469447
vibration 7
0.007508 -0.096350 -0.028492
-0.064175 0.018435 0.005452
-0.258279 0.612038 0.424257
0.189907 0.509749 -0.221522
0.043157 0.006997 0.131058
vibration 8
-0.070182 -0.025475 0.067654
0.599883 0.004874 -0.012947
0.126342 -0.256086 -0.197061
0.213822 0.519424 -0.250655
-0.104410 0.035110 -0.344874
vibration 9
0.071900 -0.014805 0.069012
-0.614551 0.002833 -0.013204
-0.007805 -0.010925 -0.022911
0.030784 0.117204 -0.066042
-0.264526 0.067172 -0.719553
[INT]
0.000788
0.000788
0.000788
0.000000
0.000000
0.000000
0.000688
0.000688
0.000688
Iwao
I examined the issue and concluded that it is necessary to manipulate the vibration widget for vibration display.
Pointer variables to components of the vibration widget can be retrieved as follows:
QWidget* vibrationDock;
QWidget* vibrationWidget;
QTableWidget* vibrationTable;
QPushButton* spectraButton;
QCheckBox* displayForcesCheckBox;
QPushButton* animationButton;
QPushButton* pauseButton;
foreach (QWidget *widget, QApplication::allWidgets())
{
if (widget->objectName().contains("vibrationDock",Qt::CaseInsensitive))
vibrationDock = widget;
vibrationWidget = vibrationDock->findChild<QWidget*>("VibrationWidget");
vibrationTable = vibrationDock->findChild<QTableWidget*>("vibrationTable");
spectraButton = vibrationDock->findChild<QPushButton*>("spectraButton");
displayForcesCheckBox = vibrationDock->findChild<QCheckBox*>("displayForcesCheckBox");
animationButton = vibrationDock->findChild<QPushButton*>("animationButton");
pauseButton = vibrationDock->findChild<QPushButton*>("pauseButton");
}
Some operations are necessary beforehand.
// Show the vibration dock if it is hidden
if (!vibrationDock->isVisible())
vibrationDock->show();
// Activate the vibration widget
vibrationWidget->setEnabled(true);
// Activate animation buttons
animationButton->setEnabled(true);
pauseButton->setEnabled(true);
// Disconnect SIGNAL/SLOTs to avoid exception when a component is clicked.
spectraButton->disconnect();
displayForcesCheckBox->disconnect();
animationButton->disconnect();
pauseButton->disconnect();
Set frequencis and intensities stored in vectors vF and vI on the vibration table.
(vF and vI were defined in the previous post.)
vibrationTable->horizontalHeader()->show();
vibrationTable->setColumnCount(2);
vibrationTable->setRowCount(vF.size());
for (unsigned int row = 0; row < vF.size(); ++row)
{
char buf[64];
sprintf(buf, "%lf", vF[row]);
QString Freq(buf);
sprintf(buf, "%lf", vI[row]);
QString Int(buf);
QTableWidgetItem* newFreq = new QTableWidgetItem;
newFreq->setText(Freq);
QTableWidgetItem* newInten = new QTableWidgetItem;
newInten->setText(Int);
vibrationTable->setItem(row, 0, newFreq);
vibrationTable->setItem(row, 1, newInten);
}
Create vibration animation for a selected mode.
Animation* m_animation = new Animation;
unsigned int m_framesPerStep = 8; // Frame number per a vibration motion.
qreal m_vibScale = 0.7; // Vibration amplitude scale factor
// Number of atoms consists in Avogadro::Molecule* m_molecule
unsigned int nAtoms = m_molecule->numAtoms();
// Atom list consists in Avogadro::Molecule* m_molecule
QList<Atom *> atoms = m_molecule->atoms();
m_molecule->setConformer(0);
std::vector<std::vector<Eigen::Vector3d> *> m_curFrames; // Conformers representing frames for a vibration
for(int i = 0; i< m_framesPerStep ; i++ )
m_curFrames.push_back(new std::vector<Eigen::Vector3d>(nAtoms));
vibrationTable->setCurrentCell(0,0); // Select a cell (a row, a column)
int curRow = vibrationTable->currentRow();
// This selects always 0th row.
// It needs use of SIGNAL/SLOT for dynamic row (or mode) selection.
for(int i = 0; i< m_framesPerStep ; i++ )
{
for(int j=0; j<nAtoms; j++)
{
m_curFrames.at(i)->at(j).x() = atoms.at(j)->pos()->x() + vLx.at(curRow).at(j).x() * m_vibScale * sin(2*M_PI*i/(m_framesPerStep-1));
m_curFrames.at(i)->at(j).y() = atoms.at(j)->pos()->y() + vLx.at(curRow).at(j).y() * m_vibScale * sin(2*M_PI*i/(m_framesPerStep-1));
m_curFrames.at(i)->at(j).z() = atoms.at(j)->pos()->z() + vLx.at(curRow).at(j).z() * m_vibScale * sin(2*M_PI*i/(m_framesPerStep-1));
}
}
m_animation->setFrame(1); // Set the current frame
m_animation->setFps(10);
m_animation->setLoopCount(0); // argument 0 means "repeat forever".
m_animation->setFrames(m_curFrames);
m_animation->setMolecule(m_molecule);
connect(animationButton, SIGNAL(clicked()), m_animation, SLOT(start()) );
connect(pauseButton, SIGNAL(clicked()), m_animation, SLOT(pause()) );
Plot IR spectrum with vectors vF and vI which defined as in the previous post.
SpectraPlotDialog *m_dialog = new SpectraPlotDialog(qobject_cast<QWidget*>(parent()));
if (m_molecule)
m_dialog->setIRSpectra(m_molecule, vF, vI);
connect( spectraButton, SIGNAL(clicked()), m_dialog, SLOT(show()) );
Oops, New class SpectraPlotDialog must be defined in advance.
I defined it as a same manner of the example, “libavogadro/examples/thirdPartyExtensions/05-ConformerPlot/”.
#include "ui_spectraplotdialog.h"
class SpectraPlotDialog : public QDialog
{
Q_OBJECT
public:
explicit SpectraPlotDialog(QWidget *parent = 0, Qt::WindowFlags f = 0 );
virtual ~SpectraPlotDialog(){}
public slots:
void setIRSpectra(Molecule *mol, std::vector< double > vF, std::vector< double > vI);
private:
Ui::SpectraPlotDialog ui;
Molecule *m_molecule;
};
SpectraPlotDialog::SpectraPlotDialog(QWidget *parent, Qt::WindowFlags f )
: QDialog( parent, f )
{ ui.setupUi(this); }
void SpectraPlotDialog::setIRSpectra(Molecule *mol, std::vector< double > vF, std::vector< double > vI)
{
ui.plot->resetPlot();
m_molecule = mol;
if (!m_molecule) return;
PlotObject *data = new PlotObject (Qt::red, // the color for the plotting object
PlotObject::Bars, // Plotting object. Points, Lines or Bars
2);
for(int i = 0; i< vF.size() ; i++)
{ data->addPoint( vF.at(i), vI.at(i), QString::number(vF.at(i)), 1.0); }
ui.plot->setDefaultLimits(500,
3500,
0.0,
0.01);
ui.plot->setAntialiasing(true);
ui.plot->setMouseTracking(true);
ui.plot->axis(PlotWidget::BottomAxis)->setLabel(tr("Wavenumber(cm-1)"));
ui.plot->axis(PlotWidget::LeftAxis)->setLabel(tr("Intensity (km/mol)"));
ui.plot->addPlotObject(data);
}
spectraplotdialog.ui is ( almost identical to
libavogadro/examples/thirdPartyExtensions/05-ConformerPlot/conformerplotdialog.ui):
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>SpectraPlotDialog</class>
<widget class="QDialog" name="SpectraPlotDialog">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>908</width>
<height>453</height>
</rect>
</property>
<property name="windowTitle">
<string>Spetrum Plot</string>
</property>
<layout class="QGridLayout" name="gridLayout">
<item row="0" column="0">
<widget class="Avogadro::PlotWidget" name="plot">
<property name="frameShape">
<enum>QFrame::StyledPanel</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
</widget>
</item>
</layout>
</widget>
<customwidgets>
<customwidget>
<class>Avogadro::PlotWidget</class>
<extends>QFrame</extends>
<header location="global">avogadro/plotwidget.h</header>
<container>1</container>
</customwidget>
</customwidgets>
<resources/>
<connections/>
</ui>
Although there are some problems such as no association with check boxes and the amplitude slider,IR spectrum lines have no band width and turn upward etc., basic functionality was impremented in this code.
Thank you.
Glad to know you got it to work. I guess I’m curious how this is different from the existing functionality?
I’d like to implement vibration visualization on my own Avogadro extension.
The vibration visualization functionality has existed on Avogadro in itself.
I only use it.
I added a function visualizing Molpro format vibration file on my extension
which can be found in https://github.com/brhr-iwao/libavogadro1cp2k.
CP2K can generate a Molpro vibration file during vibration calculation.
Visualizing Molpro vibration file has not been implemented in Avogadro 1 as far as I know.