Animation of Dupliverts with Shape Keys with Python api Description of Problem

I would like to create an animation of 100,000 atoms with 10,000 frames from my molecular dynamics simulations. Currently I am using dupliverts with a single uv sphere for each atom type. This has worked great in creating a single frame using the python api. However, I am unsure on the correct way to specify the change in vertices at each frame in Blender. I have data that will give me the position of each atom at each timestep/frame.



What I have Tried

What would be the ideal way to animate this simulation? From reading I have seen that using shape keys would allow me to specify the change in the vertices from initial positions (Basis). However, this would mean that I would have to create 10,000 shape keys and change the value from 0 to 1 for each step but would allow me to take advantage of interpolation. The best example I can find showing how this can be done using the blender python api simple example of shape key python api. However, when I used this I was only able to create one shape key not the 10,000 that I would need.

edit:
So I got the functionality that I need using this code. But it still does not answer are shape keys the correct way to handle 10,000 steps?

obj = bpy.context.object
shape_keys = ['step ' + str(i) for i in range(10000)]
for shape_key in shape_keys:
    obj.shape_key_add(name=shape_key, from_mix=False)

for shape_key in shape_keys:
    bm = bmesh.new() 
    bm.from_mesh(obj.data)
    shape = bm.verts..layers.shape[shape_key]
    for vert in bm.verts:
        vert[shape] = 'change vertex position here'
    bm.to_mesh(obj.data)


Probably very inefficient would be interested if there are better ways.
It is possible to keyframe each vertex. While by default Blender does not provide a GUI way to do this, it does include the AnimAll addon that provides a way to visually keyframe various components and values like vertices and edge creases. A python script can also easily keyframe the vertices.

import bpy
from random import random

obj = bpy.context.object
for f in range(bpy.context.scene.frame_end):
    for v in obj.data.vertices:
        m = random() * 0.1
        if f%2:
            v.co.z += m
        else:
            v.co.z -= m
        # index of 2 means only the z-axis
        # -1 (x,y, and z change), 0 x-axis, 1 y-axis
        # you will get weird behavior if the index and 
        # axis changed do not match
        v.keyframe_insert('co', index=2, frame=f)


Once keyframed they show in the graph editor and dopesheet like any other animation, meaning modifiers can be added and interpolation types changed. You will find the fcurves in obj.data.animation_data.action.fcurves.