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.

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 ;;]];

Can you publish the data used to create the film?

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"}];

Attachments:

sim-100-ELECTRON...txt

template.txt

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)