Molecular Bond Visualization with Gizmo in Igor Pro 7

1.png

This blog post is meant to show how Igor Pro 7 can be a useful tool for 3D visualization of some molecular bond structures having too many bonds to be input by manual entry. An earlier WaveMetrics example, Molecule.pxp, made with Igor Pro 6, had few molecules, and could not use Igor Pro 7's new Gizmo Path tube-visualization capability. Another useful topic described below is the use of Path waves which contain data for disparate bond segments, separated by NaN rows. I wish to thank AG (WaveMetrics) for suggesting the NaN path-gap idea, and visually enhancing the final Gizmo display.

The post is not meant to be a Gizmo tutorial; some basic familiarity with Gizmo usage is assumed. WaveMetrics provides extensive Help files and Examples on all of Gizmo's aspects if you need to come up to speed.

The simplest example that interested me is the fullerene molecule C60 (shown above; also called buckminsterfullerene, for the famous architect who came up with the "bucky-ball" structure). It is simple because there is only one type of atomic constituent, carbon. The structure is highly symmetric and has only two types of C-C bonds and bond lengths, which are either on pentagonal rings, or shared hexagon edges. This molecule is receiving new attention for its potential usage in very small atomic clocks:
https://spectrum.ieee.org/semiconductors/materials/to-build-the-worlds-smallest-atomic-clock-trap-a-nitrogen-atom-in-a-carbon-cage

A necessary starting point for the bonds' description is a knowledge of their x, y, and z coordinates.
If you want to exercise the code below and need the Cartesian coordinates for the C60 structure you can use the following link:
http://www.ccl.net/chemistry/resources/messages/2000/12/07.005-dir/
Look for the section section headed by:
  Center     Atomic     Atomic              Coordinates (Angstroms)
  Number     Number      Type              X           Y           Z

Copy and paste the 60 numerical columns below the section header into a plain-text file, which will have 60 row entries. Use the Igor Pro 7 Data→Load Waves→Load General Text... menu item. Select your text file, skip the first three columns, and load the last three columns into waves named 'wx', 'wy', and 'wz'. It will also be necessary to have a triplet N=(60,3) wave to construct the Gizmo Path object. this is easily accomplished by using a WaveMetrics utility procedure in your Procedure file:

#include <XYZtoTripletToXYZ>

The included ipf adds a Macro in the main Menu "XYZ waves to XYZ Triplet". Combine the three waves 'wx', 'wy', and 'wz' into a new triplet wave 'wxyz'. Next add to your Procedure window a user utility function that calculates the Cartesian distance between any two C atoms. 

function dist(m,n)
    variable m, n
    
    WAVE wx, wy, wz
    if(m==n)
        return 0
    else
        return sqrt((wx[m]-wx[n])^2 + (wy[m]-wy[n])^2 + (wz[m]-wz[n])^2)
    endif
end

Now we have to figure out which pairs of C atoms are nearest neighbors. There are two possible bond lengths, 1.45331 or 1.36705 Angstrom, for the data file given above. (If you use a different data set, you should recalculate possibly different bond lengths.) A useful summary of the nearest-neighbor locations is an Adjacency Matrix. In the present example the full matrix has N=(60,60), with values either 1 or 0, depending on whether the atoms at the matrix[p][q] indices are nearest neighbors or not.
https://www.researchgate.net/publication/287198431_Characterization_of_Fullerene_Using_Adjacency_Matrix

Only the upper-right matrix corner is required, but for coding simplicity the full Adjacency Matrix is constructed in the following function. The test for being a nearest-neighbor pair is that the bond distance be >0 and less than the maximum bond length (plus a small tolerance factor). The final returned value is a check on the correct number of bonds (or 3D polygon edges).

function fadj() //  make and check full adjacency matrix

    WAVE wx, wy, wz
    make/O/N=(60,60) Adj=0
    Adj = (dist(p,q)>0 && dist(p,q)<1.4534) ? 1 : 0
    MatrixOP/O msum = sum(Adj)
    return msum[0][0]/2         //  check the number of edges 3N/2
end

Finally, the Path object for visualizing the nearest-neighbor bonds is calculated. The Pathwave, 'wpath', is a series of two triplets connecting nearest-neighbors, with each pair separated by an additional NaN triplet. This expains the wave dimensions    N = (3*(Number of bonds), 3 coordinates/atom). In this function, only the upper-right corner of the Adjacency Matrix is used.

function fpath()    //  create a path wave for NaN-separated bond paths

    WAVE wxyz   //  triplet coordinate wave
    WAVE Adj        //  Adjacency matrix
    
    make/O/N=((90*3), 3) wpath  //  2 rows + 1 Nan per bond
    variable m, n, ktr=0
    for(m=0;m<60;m+=1)
        for(n=m+1; n<60; n+=1)
            if(Adj[m][n]==1)
                wpath[ktr][]          = wxyz[m][q]
                wpath[ktr+1][]  = wxyz[n][q]
                wpath[ktr+2][]  =   NaN
                ktr+=3
            endif
        endfor
    endfor
end

Then the Gizmo display is created, based on the fundamental Scatter0 object using the triplet coordinate wave 'wxyz' and the Path0 object using the wave 'wpath'. The path uses Igor Pro 7's new Tube display mode. All other Gizmo additions are to enhance the display appearance. Although this example is simple in its symmetry and composition, it may suggest ways to apply similar Igor Pro 7 methods to visualizing more complicated molecules. For example, straightforward modification of the above functions enables creation of two separate adjacency matrices and paths for long (blue) and short (yellow) bonds. There are 60 long and 30 short bonds. Note that pentagonal rings contain only long bonds, while short bonds occur only at edges between shared hexagons. 

2.png

The details for the display shown at the top of the post are given in the following (simpler) recreation macro.

Window Gizmo0() : GizmoPlot
    PauseUpdate; Silent 1       // building window...
    // Building Gizmo 7 window...
    NewGizmo/W=(62.25,49.25,652.5,590)
    ModifyGizmo startRecMacro=700
    ModifyGizmo scalingMode=8
    ModifyGizmo setOuterBox={-4.5,4.5,-4.5,4.5,-4.5,4.5}
    ModifyGizmo scalingOption=0
    ModifyGizmo keepPlotSquare=1
    AppendToGizmo Scatter=root:wxyz,name=scatter0
    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ scatterColorType,0}
    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ markerType,0}
    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ sizeType,0}
    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ rotationType,0}
    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ Shape,2}
    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ size,1}
    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ color,1,0,0,1}
    ModifyGizmo modifyObject=scatter0,objectType=scatter,property={calcNormals,1}
    AppendToGizmo freeAxesCue={0,0,0,0.4},name=freeAxesCue0
    ModifyGizmo modifyObject=freeAxesCue0,objectType=freeAxesCue,property={calcNormals,1}
    AppendToGizmo Path=root:wpath,name=path0
    ModifyGizmo ModifyObject=path0,objectType=path,property={ pathColorType,1}
    ModifyGizmo ModifyObject=path0,objectType=path,property={ lineWidthType,0}
    ModifyGizmo ModifyObject=path0,objectType=path,property={ pathColor,0.466667,0.466667,0.466667,1}
    ModifyGizmo ModifyObject=path0,objectType=path,property={ drawTube,1}
    ModifyGizmo ModifyObject=path0,objectType=path,property={ fixedRadius,0.1000}
    ModifyGizmo ModifyObject=path0,objectType=path,property={ calcNormals,1}
    ModifyGizmo modifyObject=path0,objectType=Path,property={calcNormals,1}
        ModifyGizmo setObjectAttribute={path0,blendFunc0}
        ModifyGizmo setObjectAttribute={path0,specular0}
    AppendToGizmo light=Directional,name=light0
    ModifyGizmo modifyObject=light0,objectType=light,property={ position,-0.241800,-0.664500,0.707100,0.000000}
    ModifyGizmo modifyObject=light0,objectType=light,property={ direction,-0.241800,-0.664500,0.707100}
    ModifyGizmo modifyObject=light0,objectType=light,property={ ambient,0.133000,0.133000,0.133000,1.000000}
    ModifyGizmo modifyObject=light0,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
    AppendToGizmo attribute blendFunc={770,771},name=blendFunc0
    AppendToGizmo attribute specular={1,1,1,1,1032},name=specular0
    AppendToGizmo attribute specular={1,1,0,1,1032},name=specular1
    AppendToGizmo attribute shininess={5,20},name=shininess0
    AppendToGizmo attribute specular={1,1,0,1,1032},name=specular2
    AppendToGizmo attribute shininess={5,20},name=shininess1
    ModifyGizmo setDisplayList=0, object=freeAxesCue0
    ModifyGizmo setDisplayList=1, object=light0
    ModifyGizmo setDisplayList=2, attribute=shininess1
    ModifyGizmo setDisplayList=3, attribute=specular2
    ModifyGizmo setDisplayList=4, object=path0
    ModifyGizmo setDisplayList=5, object=scatter0
    ModifyGizmo setDisplayList=6, opName=clearColor, operation=clearColor, data={0.8,0.8,0.8,1}
    ModifyGizmo currentGroupObject=""
    ModifyGizmo showInfo
    ModifyGizmo infoWindow={911,4,1724,245}
    ModifyGizmo endRecMacro
    ModifyGizmo SETQUATERNION={-0.479222,-0.350733,-0.486510,-0.640841}
EndMacro

Stephen Chinn, Ph.D. (EE)

How to compute and display scalp potentials from EEG recordings

EEG researchers commonly create spatial maps of the electrical activity acros the scalp. Data are interpolated or triangulated from a discrete number of electrode placements distributed across the scalp (generally greater than 32 sites) often in fixed/standardized positions (see, e.g., http://www.ant-neuro.com/products/waveguard_original).  To do this, electrode coordinates, either as {x, y, z} triplets or in spherical coordinates are needed and can be obtained from the manufacturer of electrode caps (e.g., https://www.egi.com/knowledge-center/item/71-sensor-position-files-for-visualizing-data-in-third-party-software).  

Theoretical analysis typically involves modeling the head as a perfect sphere.  The mapping of the electrode coordinates to the surface of a sphere can be accomplished with the following function:

Function scaleToSpherical(inWave)
    Wave inWave
    
    // offset the coordinates to the  origin:
    MatrixOP/o/free wave0c=subtractMean(inWave,1)
    
    // compute the distance of electrodes from the origin:
    MatrixOP/O/free rawRad=sqrt(magsqr(col(wave0c,0))+magsqr(col(wave0c,1))+magsqr(col(wave0c,2)))
    
    // rescale the distance to a radius of 1:
    MatrixOP/O centeredNormalized=scalerows(wave0c,rec(rawRad))
End
Figure 1: Initial electrode positions relative to the surface of a sphere.

Figure 1: Initial electrode positions relative to the surface of a sphere.

Figure 2: Normalized electrode coordinates on the surface of a sphere.

Figure 2: Normalized electrode coordinates on the surface of a sphere.

Using the normalized set of coordinates we can triangulate the data on the surface of a sphere using the SphericalTriangulate operation. 

SphericalTriangulate centeredNormalized // this creates M_TriangulationData

Following the triangulation we can use the wave M_TriangulationData combined with the EEG data to interpolate the potentials at any point on the sphere.  In order to display the resulting interpolation in Gizmo we start by creating a parametric surface representing the sphere. 

Function makeParametricSphere(pointsx,pointsy)
    Variable pointsx,pointsy
    
    Variable i,j,rad
    Make/O/N=(pointsx,pointsy,3) parametricW
    Variable anglePhi,angleTheta
    Variable dPhi,dTheta
    
    
    dPhi=2*pi/(pointsx-1)
    dTheta=pi/(pointsy-1)
    Variable xx,yy,zz
    Variable sig
    
    for(j=0;j<pointsy;j+=1)
        angleTheta=j*dTheta
        zz=sin(angleTheta)
        if(angleTheta>pi/2)
            sig=-1
        else
            sig=1
        endif
        for(i=0;i<pointsx;i+=1)
            anglePhi=i*dPhi
            xx=zz*cos(anglePhi)
            yy=zz*sin(anglePhi)
            parametricW[i][j][0]=xx
            parametricW[i][j][1]=yy
            parametricW[i][j][2]=sig*sqrt(1-xx*xx-yy*yy)
        endfor
    endfor
End

To complete the interpolation we need to construct one more wave: dataPointsWave.  This wave has 4-columns consisting of the electrode coordinates in the first three columns and the potentials measured by each electrode at some fixed time t0 in the fourth column:

    Variable nRows=DimSize(centeredNormalized,0)
    Make/O/D/N=(nRows,4) dataPointsWave
    dataPointsWave[][0,2]=centeredNormalized
    dataPointsWave[][3]=electrodePotentials[p]

The interpolation is computed by:

    SphericalInterpolate M_TriangulationData,dataPointsWave,sphereData

The resulting interpolation is stored in the wave W_SphericalInterpolation where each row contains the scalar potential interpolated for the corresponding XYZ location in the wave sphereData.  To display these values on the parametric surface we need to construct a color wave

Function createParametricColorWave(pWave,interpWave,reverseCTAB)
    Wave pWave,interpWave
    Variable reverseCTAB            // set to 1 to reverse; 0 otherwise
    
    // the color wave must match the parametric surface:
    Variable rows=dimsize(pwave,0)
    Variable cols=dimSize(pWave,1)
    
    // Create the parametric color wave:
    Make/O/N=(rows,cols,4) pWaveColor=1     
    
    // find the range of values for scaling:
    Variable mmax=WaveMax(interpWave)
    Variable mmin=WaveMin(interpWave)
    
    // optionally replace the colortable here:
    ColorTab2Wave rainbow256
    Wave M_Colors
    // not all colortables have the same number of colors:
    Variable nTableCols=DimSize(M_Colors,0)-1

    if(reverseCTAB)
        MatrixOP/o/free aa=col(M_Colors,0)
        WaveTransform/o flip aa
        MatrixOP/O M_Colors=setCol(M_Colors,0,aa)
        MatrixOP/o/free aa=col(M_Colors,1)
        WaveTransform/o flip aa
        MatrixOP/O M_Colors=setCol(M_Colors,1,aa)
        MatrixOP/o/free aa=col(M_Colors,2)
        WaveTransform/o flip aa
        MatrixOP/O M_Colors=setCol(M_Colors,2,aa)
    endif
    
    MatrixOP/O M_Colors=fp32(M_Colors/65535)
    Variable nor=nTableCols/(mmax-mmin)
    Variable i, np=rows*cols,index,rr,cc
    for(i=0;i<np;i+=1)
        index=trunc((interpWave[i]-mmin)*nor)
        rr=mod(i,rows)
        cc=trunc(i/rows)
        pWaveColor[rr][cc][0]=M_Colors[index][0]
        pWaveColor[rr][cc][1]=M_Colors[index][1]
        pWaveColor[rr][cc][2]=M_Colors[index][2]
    endfor
End

An example of this approach is illustrated in a cross-lingustic speech discrimination study at the following webpage:

http://www.gc.cuny.edu/Page-Elements/Academics-Research-Centers-Initiatives/Doctoral-Programs/Speech-Language-Hearing-Sciences/Research-Labs/Developmental-Neurolinguistics-Laboratory/Research-Highlights

Valerie L. Shafer, Ph.D.
vshafer@gc.cuny.edu
Professor and Deputy Executive Officer
Ph.D. Program in Speech-Language-Hearing Sciences
The Graduate Center, CUNY
365 Fifth Avenue, NY, NY 10016

Macintosh 64-bit XOPs and Igor Pro 8

There is an issue that will affect Macintosh 64-bit XOPs running with Igor Pro 8. We don't know when we will start beta testing Igor8, but it will be no earlier than the fall of 2017.

The short story is that, because of a change in Igor8, existing Macintosh 64-bit XOPs need to be modified and recompiled with XOP Toolkit 7.01 or later to run with the 64-bit version of Igor Pro 8 on Macintosh. Macintosh 32-bit XOPs as well as Windows 32-bit and 64-bit XOPs will run with Igor8 without modification.

The 64-bit Macintosh version of Igor7 uses Mac OS handles both internally and to exchange data with XOPs. Mac OS handles are limited to roughly 2^31 bytes, even in 64-bit applications. To overcome this limit, the Macintosh 64-bit version of Igor8 will use WaveMetrics handles - that is, handles created and manipulated using WaveMetrics code, as all Windows versions have always done.

XOP Toolkit 7.01 provides the necessary support to allow Macintosh 64-bit XOPs to use WaveMetrics handles. It adds WaveMetrics memory management XOPSupport routines that allow an XOP to use compatible handles no matter what version of Igor it is running with.

Although this change is critical only for Macintosh 64-bit XOPs, we recommend that all actively-developed XOPs be changed, when it is convenient, to use the new WaveMetrics memory management XOPSupport routines so that the source code will be compatible with any version of Igor on any platform. The source code changes required are minimal and amount mostly to changing NewHandle to WMNewHandle, DisposeHandle to WMDisposeHandle, and so on.

When you are ready to update your XOP, here is what you need to do:

  1. Download the latest XOP Toolkit. If you have an XOP Toolkit 7 license, use the download instructions that you previously received from WaveMetrics. If you have an XOP Toolkit 6 license, contact WaveMetrics sales for a free upgrade to XOP Toolkit 7.
  2. If you have not yet updated your XOP Toolkit 6 XOP for XOP Toolkit 7, follow the instructions in Appendix A of the XOP Toolkit 7 manual.
  3. See Appendix C of the XOP Toolkit 7 manual for an overview of changes since the release of XOP Toolkit 7.00.
  4. Follow the link in Appendix C to the section "WM Memory XOPSupport Routines" in Chapter 12. This explains the issue in greater detail.
  5. Follow the instructions in the section "Updating Old Code to Use WM Memory XOPSupport Routines" in Chapter 12.

The XOP Toolkit 7 manual is included in the XOP Toolkit 7. It is also available from http://www.wavemetrics.net/doc/XOPMan7.pdf.

For an overview of the XOP Toolkit, see https://www.wavemetrics.com/products/xoptoolkit/xoptoolkit.htm.

Get to Know a Feature: Color Table Wave Creation

Get to Know a Feature: Color Table Wave Creation

In "Get to Know a Feature: Color Table Wave Basics", we showed how Igor 7's color table waves are more flexible than Igor 6's color index waves because multiple uses of the same color table wave can cover differing Z value ranges. This post discusses several ways that you can create your own color table waves such as including importing RGB values from a CSV file and recreating a color table from a screen capture,