Home   Manual

How to display vibration (nomal mode) on Avogadro 1


#1

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


#2

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.


#3

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


#4

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.


#5

Glad to know you got it to work. I guess I’m curious how this is different from the existing functionality?


#6

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.