%two models of molecular dynamics simulated with a 2D potential energy
%on the domain [0,1]^2 with periodic boundary conditions
%blue curve = Langevin dynamics, red curve = Brownian dynamics

V = @(x) cos(4*pi*x(1))+cos(4*pi*x(2))+0.2*sin(2*pi*x(1))+0.2*sin(2*pi*x(2));    %potential energy (must be periodic)
grad_V = @(x) [-4*pi*sin(4*pi*x(1))+0.2*2*pi*cos(2*pi*x(1));-4*pi*sin(4*pi*x(2))+0.2*2*pi*cos(2*pi*x(2))];
dt = 0.01;        %discrete step size
dt2 = 0.0001;
beta = 0.5;
beta2 = 3; 

steps = 10^9;      %number of simulation steps  

close all
figure;
axis([0 1 0 1])
axis square
a = linspace(0,1,100);
b = linspace(0,1,100);
c = zeros(100,100);
for i=1:100
    for j = 1:100
        x = [a(i),b(j)];
        c(j,i) = V(x);   %order needs to be reversed for contourf for some reason
    end
end
contourf(a,b,c);
colorbar
colormap summer
hold on;

%initial points

q_old = [.4;.4];  
q_new = q_old;
p_old = [-.1;.2];
p_new = p_old;
x_old = [.6;.6];
x_new = x_old;

%parameters

gam = .1;   %damping
m_inv = 3; %inverse mass
beta = .2;  %inverse temperature
sig = sqrt(2*gam/beta);  %noise constant

for step=1:steps
    q_old = q_new;      %first Langevin dynamics
    p_old = p_new;
    u = rand;
    v = rand;
    mu1 = sqrt(-2*log(u))*cos(2*pi*v);
    nu1 = sqrt(-2*log(u))*sin(2*pi*v);
    u = rand;
    v = rand;
    mu2 = sqrt(-2*log(u))*cos(2*pi*v);
    nu2 = sqrt(-2*log(u))*sin(2*pi*v);
    p2 = p_old - 0.5*dt*grad_V(q_old) - 0.5*gam*m_inv*p_old + sqrt(dt/2)*sig*[mu1;nu1];
    q_new = mod(q_old + dt*m_inv*p2,1);
    p_new = (p2 - 0.5*dt*grad_V(q_new) + sqrt(dt/2)*sig*[mu2;nu2])/(1+0.5*dt*gam*m_inv);
    if step > 1
       children = get(gca,'children');
       delete(children(1:2));
       delete(children(400:length(children)-1));
    end
    x_old = x_new;       %now Brownian dynamics
    u = rand;
    v = rand;
    mu3 = sqrt(-2*log(u))*cos(2*pi*v);
    nu3 = sqrt(-2*log(u))*sin(2*pi*v);
    x_new = mod(x_old - dt2*grad_V(x_old) + sqrt(2*dt2/beta2)*[mu3;nu3],1);
    line([q_old(1),q_new(1)],[q_old(2),q_new(2)],'color','b');
    line([x_old(1),x_new(1)],[x_old(2),x_new(2)],'color','r');
    plot(q_new(1),q_new(2),'.','MarkerSize',14,'color','b');
    plot(x_new(1),x_new(2),'.','MarkerSize',14,'color','r');
    axis([0 1 0 1])
    title('blue = langevin, red = brownian dynamics','interpreter','latex')
    axis square
    M(step) = getframe(gca);
end

close(gcf);
