Message Boards Message Boards

Quantum Chemistry Animations

Posted 6 years ago

A couple of weeks ago I experimented with automating inputs and visualizing results from the CP2K quantum chemistry software via Mathematica. My goal was to visualize the formation of a water molecule. I knew just enough about computational chemistry to stumble my way to making an animation that sort of looked like I wanted, although I'm sure it's horribly inaccurate. Are there any other people in the community who have worked on projects like this?

Video

enter image description here

POSTED BY: Michael Hale
13 Replies
Anonymous User
Anonymous User
Posted 6 years ago

https://sourceforge.net/projects/mathematica-spice/

this is a good example of a potential developer that started with foreign source and file format, relied on coding the format in another language and viewing in mathematica. in the end (or rather presently), there is no interface only graphics that is difficult to utilize

chemistry is a better situation. there are standard formats of "chemical information systems", mathematica offers live download of certain sources (data with pertinent information specific formatting, making it potable for mathematica to offer simulations and coding with), and mathematica has examples of minior simulation. if you compared it to solar information systems and data you might find the chemical simulations are ... not as well formated. (but compare mathematica's solar simulation with professional "space" solar simulators, you'll find there's also allot to be done comparably)

so these carefully laid out formats (containing data meant to be utilized in many kinds of simulations easily): they go so far but tend to stop before becoming as useful as the simulators they come from (i should rather say mathematica is like a swiss army knife tool that can do particular research the simulators cannot do because they have no foundation for doing things other than "mostly" pre-conceived actions)

POSTED BY: Anonymous User
Anonymous User
Anonymous User
Posted 6 years ago

I spent about two evenings on it that's how technology navigation is learned, the hard way

CP2K is one of the better free dynamic chemical simulators.

Always on my mind is: chemical simulators, electronics simulators, building design simulators (including many simulators, rich in pertinent information), raytrace simulation (physical light appearance simulators / 3D cad design even game design), city simulators. they are too large to just "release in mathematica" (would have to be just one installed add-on, for most)

they have huge interfaces (popular 3D game design interfaces are a good example) and mathematica is poor at that. huge code bases (mathematica is good at that but isn't connected for many reasons)

the simulator code is custom and the rich information locked in it. the math is poor in these (often so is ability to modify). there is little way to interact with Mathematica (a few file format converters - usually very limited in ability, cut&past - very limited). many ask about MathLink but few or none have ever produced any product that links up with any large simulator.

mathematica, despite not being "higher level than C", can compile and run code "as good or better". but it's ability to offer "large interfaces" is poor. the people willing to tackle "rewriting a simulator in mathematica" and provide a thorough MathLink interface: is about none. we don't have a mathematica download area with mathlinke'd simulators, not even one. CDF format is really only good at (trivial) problem exhibition.

this is a predicament that i frequently think about

a ton of custom simulators that can't do math and have no "central language"

a best mathematics, algorithms, code program which can't do "interfaces" (interfaces much due a cross platform lib issues, software wars, etc)

mapping is a good example of a technology that was simulated by custom apps that somehow gets "active" in mathematica

mapping companies made special large interfaces and multi-layered maps and mapping format files. these became orchestrated systems of standards of "how to transmit various mapping information". finally mathematica offers mapping data that has drivers capable of reading that information which is rich which (simulator specific information of how it is used)

the result: mathematica has powerful mapping simulation (though little interface). no circuit simulation.

but what about the rest? good question.

POSTED BY: Anonymous User

Can you publish the data used to create the film?

Posted 6 years ago

I overwrote the cube files from the first clip. I still have them from the second clip, but they are about 500 MB in total. I've attached an example of a random one (I had to change extension to .txt for the upload to work) and the template I used for the input file to CP2K. The initial conditions I used for the second run were just:

initPos = N@{{4, 4, 4}, {6, 6, 4}, {6, 6, 6}};

To generate the input file then you just use:

FileTemplateApply[File["template.txt"], 
 Append[Flatten@initPos, frame], "input.txt"];

Then run the calculation with:

DeleteFile["testoutput.txt"];
RunProcess[{"cp2k.sopt", "-i", "input.txt", "-o", "testoutput.txt"}];
POSTED BY: Michael Hale

Michael, thank you. In general, this is a very interesting approach to solving the problem of the dynamics of atoms and molecules.

@Michael Hale This is great! Does CP2K run the dynamics as well - so you only need to input the initial configuration? How does Blender come in to play? If you CP2K makes cube files, you can read those into Mathematica to plot them.

@Alexander Trounev It is definitely possible to code these simulations up in Wolfram Language alone - search "Hartree Fock mathematica" for some examples. But for anything but the smallest system, it would be terribly slow compared to the dedicated electronic structure packages available. A more viable solution is to make a linker package, one that would create input files for CP2K (or psi4, or Gaussian, molpro, etc), call the external program as an asynchronous task, and read in the results back into Mathematica for visualization or further computation.

POSTED BY: Jason Biggs
Posted 6 years ago

They do have options to calculate the dynamics in CP2K. And it seems like they have support for features like Langevin dynamics. But at a glance I only saw constant time stepping, and variable time stepping seems pretty important. Also monitoring the progress of CP2K calculations across multiple time steps is harder than just doing one at a time so you know when to look at new output. I'm not sure what the ideal middle ground is yet.

Here is essentially the code I used to read the cube files.

data = Import[
   cp2kDir <> "sim-" <> ToString@frame <> 
    "-ELECTRON_DENSITY-1_0.cube", "Table"];
output = Import[cp2kDir <> "testoutput.txt", "Table"];
atomCount = data[[3, 1]];
numX = data[[4, 1]];
basisX = data[[4, 2 ;;]];
numY = data[[5, 1]];
basisY = data[[5, 2 ;;]];
numZ = data[[6, 1]];
basisZ = data[[6, 2 ;;]];
volume = data[[7 + atomCount ;;]] // Flatten // Rescale // 
   Transpose[ArrayReshape[#, {numX, numY, numZ}], {3, 2, 1}] &;
atomTypes = data[[7 ;; 7 + atomCount - 1, 1]];
atomPos = 
  Rescale[#[[3 ;;]], {0, numX*basisX[[1]]}, {0, numX}] & /@ 
   data[[7 ;; 7 + atomCount - 1]];
atomForces = 
  Take[output, 
    Position[
       output, {"#", "Atom", "Kind", "Element", "X", "Y", "Z"}][[1, 
      1]] + {1, 3}][[All, 4 ;;]];
POSTED BY: Michael Hale
Posted 6 years ago

Oh, and I typically use Blender just to assemble large numbers of frames into a video.

POSTED BY: Michael Hale

Dear Jason Biggs, although the task is very difficult, but we see nothing complicated in the film itself.

I understand that you worked hard to get these pictures and create a video. But all this looks not very interesting. It is possible that not the best visualization tool is used. What is computed in this task and how? I looked at the code cp2k-6.1, but it's like a black box.

Posted 6 years ago

I spent about two evenings on it, so no, criticism won't bother me. But I would argue that it's already more interesting for a general audience than most quantum chemistry calculations. Most people have probably seen molecular dynamics simulations, but those are classical and focus on electrostatic interactions. They assume covalent bonds to be constant, so they can't show the formation of a water molecule, for example. It seems most quantum chemistry computations typically calculate electron configurations around stationary molecules, which they then use in classical molecular dynamics simulations. Or they are just calculating vibration spectra or something. I haven't seen many interesting quantum chemistry animations. I guess the workflow is still not simple enough, but it would be cool to see a lot of these animations used in a chemistry textbook. I know I have terrible intuition about chemistry, and maybe that would help.

For each frame in this animation, the inputs are the atom positions (two hydrogens and an oxygen), and the outputs are a 3D grid representing the probability of finding an electron at each point in the grid (the square of the wavefunction) as well as 3 vectors representing the force of the electron wavefunction and the other atoms acting on each atom. The reason you can separate the calculation of the electron wavefunction from the movement of the atoms is because it adjusts much more quickly than the atom positions. This is called the Born-Oppenheimer approximation. The method used by CP2K is called density functional theory. There are whole books written on it, but I think at a high level it works by taking a bunch of basis functions (think spherical harmonic functions from intro hydrogen orbital calculations) and adjusts their configurations and weights until an energy minimum is achieved. And importantly, this energy calculation takes into account effects like electron correlation and exchange, which is harder to do in Hartree-Fock methods.

POSTED BY: Michael Hale

The visualization is done somehow. You can do much better by using the graphical tools of Mathematica. Give the system of equations that must be solved to visualize the formation of the water molecule. Why do you need the CP2K package? At first glance, this task can be solved using only the code written in the Wolfram Language.

Posted 6 years ago

The visualization is done with ListDensityPlot3D in Mathematica. CP2K has an option to output "cube" files of the electron density after it finishes its density functional theory calculations, which are straightforward to parse.

The reason I used an external tool for the calculations is just because I was able to figure that out faster than figuring out how to implement it myself. Several years ago, I tried just "typing the equation in" for the time-dependent Schroedinger equation with a simple 3D spherical potential well, and I struggled quite a bit. So I would have even less success with a multi-atom potential and multiple electrons. I have a general understanding of the two main algorithms used in computational chemistry to approximate the time independent solutions (Hartree-Fock and density functional theory), but not enough to implement them myself and surely not enough to implement them better than a popular third party library. Although I'm sure there is room to improve the algorithms for various cases, and Mathematica would be a good tool to research that if I had the expertise. But I haven't even investigated most of the options available in CP2K yet. There are lots of methods, and options for those methods, and different sets of basis functions available.

To make those videos I generate the input file for CP2K, which contains the atom coordinates. I have it set to calculate the energy and forces, using a low number of iterations (for speed and because the atoms aren't near an equilibrium configuration, so it wouldn't converge anyway), and to output the cube file of the electron density. Then I import the forces from the general output file and update the atom positions and velocities. I use some simple adjustable time stepping logic, so when the atoms are moving quickly I update the velocities and positions using a smaller time step, so that I don't accidentally overshoot any important interactions. Then at the end I calculate the time at which I want each frame in the animation to occur, choose the cube file that occurred closest to that time, and parse and draw the result. Then I save the animation from the individual frames in Blender. So each one of those two clips took about an hour and I had a live view of the most recent output, which seems pretty fast and user friendly compared to most quantum chemistry software. Although like I said, I wasn't optimizing for accuracy.

POSTED BY: Michael Hale
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract