molecularDynamics

classdef molecularDynamics < handle
    properties
        N;
        L;
        t = 0;
        dt;
        otime = 1;
        rc = 2.5;
        pos = [];
        vel;
        m;
        phand = [];
        timeser;
        Eser;
        KEser;
        PEser;
        KEfix = false;
        KEset = 0.2;
    end
    
    methods
        function obj = molecularDynamics(N, L, dt, KEset)
            switch nargin
                case 4
                    obj.N = N;
                    obj.L = L;
                    obj.dt = dt;
                    obj.KEset = KEset;
                case 0
                    obj.N = 36;
                    obj.L = 6;
                    obj.dt = 0.01;
                    obj.KEset = 0.2;
            end
            
            % Generate initial positions
            [X,Y] = meshgrid([0.5:1:obj.L, 0.5:1:obj.L]);
            for i = 1:obj.L
                for j = 1:obj.L
                    obj.pos = [obj.pos; X(i,j) Y(i,j)];
                end
            end
            obj.pos = obj.pos(1:obj.N, 1:2); % pos is only first N elements
            
            % Generate distribution of velocities with avgKE = KEset
            randVel = randn(obj.N, 1);
            randVel = (randVel - mean(randVel))/sqrt(sum(0.5.*(randVel - mean(randVel)).^2))*sqrt(obj.KEset);
            randDir = 2*pi*rand(obj.N, 1);
                        subplot(1,2,1);

            % Generate velocities and handles for circles
            for i = 1:obj.N
                obj.vel = [obj.vel; randVel(i)*cos(randDir(i)) randVel(i)*sin(randDir(i))];
                obj.phand(i) = rectangle('Position', [obj.pos(i,1)-0.5 obj.pos(i,2)-0.5 1 1],'Curvature',[1 1]); 
            end
        end
        
        function [] = draw(obj)
            % draw a frame every otime
            subplot(1,2,1);
            pbaspect([1 1 1]);
            axis([-1 7 -1 7]);
            for i = 1:obj.N
                set(obj.phand(i),'Position', [obj.pos(i,1)-0.5 obj.pos(i,2)-0.5 1 1]);
            end

            subplot(1,2,2);
            hold on
            plot(obj.timeser, obj.KEser, 'b');
            plot(obj.timeser, obj.Eser, 'g');
            plot(obj.timeser, obj.PEser, 'r');
            legend('Kinetic Energy', 'Total Energy', 'Potential Energy')
            hold off
        end
        
        function [] = step(obj)
            obj.t = obj.t + 1;
            obj.pos = obj.pos + obj.dt.*obj.vel;
            for i = 1:obj.N
                for dim = 1:2
                    if (obj.pos(i,dim) > 6)
                        obj.pos(i,dim) = obj.pos(i,dim) - 6;
                    elseif (obj.pos(i,dim) < 0)
                        obj.pos(i,dim) = obj.pos(i,dim) + 6;
                    end
                end
            end
            
            % Calculate forces between atoms
            for i = 1:obj.N
                for j = 1:obj.N
                    if (i ~= j)
                        x = obj.pos(i, 1) - obj.pos(j, 1);
                        y = obj.pos(i, 2) - obj.pos(j, 2);
                        
                        if (x > obj.L/2)
                            x = x - obj.L;
                        elseif (x < -obj.L/2)
                            x = x + obj.L;
                        end
                        if (y > obj.L/2)
                            y = y - obj.L;
                        elseif (y < -obj.L/2)
                            y = y + obj.L;
                        end
                        
                        r = sqrt(x^2 + y^2);
                        if (r < 2.5)
                            f = [24 * (2 * r^-14 - r^-8) * x, 24 * (2 * r^-14 - r^-8) * y];
                            obj.vel(i,:) = obj.vel(i,:) + obj.dt .* f;
                        end
                    end
                end
            end
            
            % Calculate data every otime seconds
            if (mod(obj.t, obj.otime) == 0)
                PE = 0;
                for i = 1:obj.N
                    for j = 1:obj.N
                        if (i ~= j)
                            x = obj.pos(i, 1) - obj.pos(j, 1);
                            y = obj.pos(i, 2) - obj.pos(j, 2);

                            if (x > obj.L/2)
                                x = x - obj.L;
                            elseif (x < -obj.L/2)
                                x = x + obj.L;
                            end
                            if (y > obj.L/2)
                                y = y - obj.L;
                            elseif (y < -obj.L/2)
                                y = y + obj.L;
                            end

                            r = sqrt(x^2 + y^2);
                            if (r < 2.5)
                                PE = PE + 4 * (r^-12 - r^-6);
                            end
                        end
                    end
                end
                obj.timeser = [obj.timeser; obj.t];
                obj.PEser = [obj.PEser; PE];
                KE = 0.5*sum(obj.vel(:,1).^2+obj.vel(:,2).^2);
                obj.KEser = [obj.KEser; KE];
                obj.Eser = [obj.Eser; KE+PE];
            end
        end
    end
end