Molecular Orbitals in VR

A while ago I started working on a VR Molecular Visualizer. A feature I wanted to add to it was the ability to overlay a plot of a Molecular Orbital on top of the XYZ structure. I recently got my very own copy of Turbomole, so I thought this would be a great test case! I fired up an example of benzene, and had it output a generic xyz format for using an arbitrary program to plot 3D grids.

Hrm. Well then. That won't do at all! These programs output a 3D grid, where every point has an iso value associated with it. I usually plot MOs with an iso value of +- 0.01 in gOpenMol, so I whipped up a cleaver little awk script to weed out anything where (f > 0.01 || f < -0.01)!

Well, shit. That gave me the balloon I wanted, along with all the air in it... Looks like I'm going to have to make a proper program to do this right. Luckily, I write more Scala than Fortran now-a-days :-)

I know that I want a case class Orbital(fileLines: Seq[String]) that I can feed the whole file, and parse the lines. I also know the first 11 lines of this output have the info I need and that I want to limit the isosurface value to a magnitude of 0.01.

private val nHeaderLines: Int = 11
private val isoLimit: Double = 0.01

In these 11 lines, I will get the grid information about the x, y, and z axes.

val grid1: Grid = Grid(fileLines(4))
val grid2: Grid = Grid(fileLines(5))
val grid3: Grid = Grid(fileLines(6))

where each grid line looks like


case class Grid(line: String) {
      private val data = line.replaceAll(" +", " ").split(" ")
      val start: Double = data(2).toDouble
      val delta: Double = data(4).toDouble
      val nPoints: Int = data(6).toInt
}

We will then have x y*z blocks of data to parse through

private val nBlocks = grid2.nPoints * grid3.nPoints  
val data: Seq[GridData] = {  
    val allData = for {
      index <- 0 until nBlocks
    } yield {
      val skip = nHeaderLines + (grid1.nPoints * index) + index
      fileLines.slice(skip, skip + grid1.nPoints).map(GridData)
    }
    allData.flatten.filter(p=> Math.abs(p.f) >= isoLimit)
  }

with each line of grid data looking like

case class GridData(line: String) {  
  private val data = line.replaceAll(" +", " ").split(" ").filter(!_.equals(""))
  val x: Double = data(0).toDouble
  val y: Double = data(1).toDouble
  val z: Double = data(2).toDouble
  val f: Double = data(3).toDouble
}

If we look at the positive iso values, all we need to do is group by (x,y), which will select a plane, and then sort by the z value. Then we can simply take the first and last elements to get the outter top/bottom of our data

  val posDensity = {

    val xyPlane = data.filter(_.f > 0)
      .groupBy(d => (d.x, d.y))

    val outerPoints = xyPlane.flatMap{ d =>
      val sorted = d._2.sortBy(_.z)
      val maxZ = sorted.head
      val minZ = sorted.last
      Set(maxZ, minZ)
    }

    outerPoints.toSeq

  }

We do the same thing with the negative iso values, and then combine the data and get the following

Bingo!

Now our whole process is as simple as

object Main extends App {

  val fileName = "dnt/1e2u2.xyz"
  val fileLines = Source.fromFile(fileName).getLines().toSeq
  val orbital: Orbital = Orbital(fileLines)

  val posPw = new PrintWriter(new File("dnt/pos.xyz"))
  val negPw = new PrintWriter(new File("dnt/neg.xyz"))

  orbital.posDensity.foreach{ d =>
    posPw.println(s"${d.x} ${d.y} ${d.z} ${d.f}")
  }
  posPw.flush()
  posPw.close()

  orbital.negDensity.foreach{ d =>
    negPw.println(s"${d.x} ${d.y} ${d.z} ${d.f}")
  }
  negPw.flush()
  negPw.close()
}

Let's load these up in Unity and see what they look like in VR!

Pretty sweet!

If you want to play around with the code yourself, you can snag it from github