function [trj] = simple_md(nparticles,density,T,dt,nsteps,nrecord)

% nparticles = number of particles
% density = density of the unit cell
% T = temperature in Kelvin
% dt = step time
% nsteps = total number of steps
% nrecord = coordinates are stored every nrecord
        
    % Ensemble parameters
    mass = 1;            % Particles mass
    ndims = 3;           % Dimensions 
    
    % Simulation parameters
    dt2 = dt*dt;         
    T = 300.0;           % Temperature
    
    % Temperature control parameters
    nu = 0.1;                 % Frequency of collisions with the heat bath
    vel_std = sqrt(T/mass);   % standard deviation of the velocity
    
    % Unit cell
    [coords L] = unit_cell(nparticles,density);
    
    % Initial normally distributed velocities 
    vels = zeros(ndims,nparticles);
    for part = 1:nparticles
        vels(:,part) = randn(ndims,1) * vel_std;        
    end
        
    % Initial momentum correcting for any drift of the center of mass 
    V = sum(vels,2)/nparticles;     % Center of mass velocity
    for dim = 1:ndims
        vels(dim,:) = vels(dim,:) - V(dim);    
    end
    
    % Initial kinetic energy
    V2 = sum(dot(vels,vels))/nparticles;    % Mean-squared velocity
    vel_scale = sqrt(ndims*T/V2);           % Velocity scaling factor
    for dim = 1:ndims
        vels(dim,:) = vels(dim,:)*vel_scale;    
    end
    
    % Initial forces
    forces = lennard_jones(coords,L,ndims,nparticles);
    
    % MD. Coordinates are stored every nrecord steps
    trj = zeros(3,nparticles,floor(nsteps/nrecord)+1);
    trj(:,:,1) = coords;
    
    time = 0; 
    frame_ind = 1;
    for step = 1:nsteps
                
        % Update coordinates 
        coords = coords + dt*vels + 0.5*dt2*forces; 
        
        % Apply periodic boundary conditions
        for part = 1:nparticles
            coords(:,part) = periodic_boundaries(coords(:,part),L,ndims);
        end
        
        % Update trajectory 
        if mod(step,nrecord) == 0
            display(['frame no: ' num2str(step)]);
            frame_ind = frame_ind+1;
            trj(:,:,frame_ind) = coords;
        end
        
        % Update velocities
        vels = vels + 0.5*dt*forces;
         
        % Update forces
        forces = lennard_jones(coords,L,ndims,nparticles);
                
        % Update velocities 
        vels = vels + 0.5*dt*forces;
        
        % Temperature control. Test for collisions with the heat bath. If
        % the particle collided with the bath, draw a new velocity out of a
        % normal distribution.
        for part = 1:nparticles
            if (rand < nu*dt)
                vels(:,part) = randn(ndims,1) * vel_std;        
            end
        end
               
        % Update time
        time = time + dt;
                
    end
end    

    function [coords,L] = unit_cell(nparticles,density)
    
        coords = zeros(3,nparticles);
    
        % Cubic unit cell size
        L = (nparticles/density)^(1.0/3);
    
        % Lowest cube greater than or equal to the number of particles
        cube = 2;
        
        while (cube^3 < nparticles)
            cube = cube + 1;
        end
                
        % Coordinates
        ind = [0,0,0]';
        for i = 1:nparticles
            coords(:,i) = (ind+[0.5,0.5,0.5]')*(L/cube);
            ind(1) = ind(1) + 1;
            if (ind(1) == cube) 
                ind(1) = 0;
                ind(2) = ind(2) + 1;
                if (ind(2) == cube)
                    ind(2) = 0;
                    ind(3) = ind(3) + 1;
                end
            end
        end   
    end
        
    function coords = periodic_boundaries(coords,L,ndims)
    
        % Half box size in each direction
        hL = L/2;

        % Distance vector should be in the range -hLx -> hLx and -hLy ->
        % hLy
        for i = 1:ndims
            if (coords(i) > hL)
                coords(i) = coords(i)-L;
            elseif (coords(i) < -hL)
                coords(i) = coords(i)+L;
            end
        end   
    end
    
    function forces = lennard_jones(coords,L,ndims,nparticles)
    
        forces = zeros(size(coords));
        
        % Loop over all particle pairs
        for i = 1:nparticles-1
            for j = (i+1):nparticles
                
                % Particle-particle distance accounting for periodic
                % boundary conditions
                r = coords(:,i) - coords(:,j);
                r = periodic_boundaries(r,L,ndims);
                inv_r2 = 1/(r'*r);
   
                % Lennard-Jones potential:
                % U(r) = 4*epsilon* [(sigma/r)^12 - (sigma/r)^6]
                %
                % If sigma = 1, epsilon = 1 (reduced distance and energy
                % units) we derive:
                %
                % U(r) = 4 * [(1/r)^12 - (1/r)^6]
                
                U = 4*(inv_r2^6 - inv_r2^3);
                    
                % Action and reaction (Newton's third law)
                forces(:,i) = forces(:,i) + r*U;
                forces(:,j) = forces(:,j) - r*U;
               
            end
        end
    end

    
    
    