CODEWITHSUNDEEP
pythonsimulationphysicsmontecarlo
wiedzminYo
wiedzminYo

Reputation: 561

Metropolis monte carlo ideal gas simulation volume goes to zero

I'm trying to do simple simulation of ideal gas according to Clapeyron equation `pv=nkbT' using Metropolis Monte Carlo algorithm.This is very simple example,where I consider molecules in 2D with no interactions with each other and energy is equeal to E=pV wher V is area of circle containing all molecules.

My problem is that after very few monte carlo steps volume of my gas goes always to almost zero,no matter how many molecules or pressure I put.I can't figure out if I have bug in my code,or it is becouse I don't have any molecules interactions. All help will be much appriciated,here is my code

My results are shown in plot bellow,x-axis are monte carlo steps and y-axis is volume,i would expected as a result some none zero constant value of volume after greater number of steps.plot of volume vs monte carlo steps

import numpy as np
import random
import matplotlib.pyplot as plt


def centroid(arr):
    length = arr.shape[0]
    sum_x = sum([i.x for i in arr])
    sum_y = sum([i.y for i in arr])
    return sum_x/length, sum_y/length


class Molecule:
    def __init__(self, xpos, ypos):
        self.pos = (xpos, ypos)
        self.x = xpos
        self.y = ypos


class IdealGas:
    # CONSTANTS
    def __init__(self, n,full_radius, pressure, T, number_of_runs):
        gas = []
        for i in range(0, n):
            gas.append(Molecule(random.random() * full_radius,
                                random.random() * full_radius))
        self.gas = np.array(gas)
        self.center = centroid(self.gas)
        self.small_radius = full_radius/4.
        self.pressure = pressure
        self.kbT = 9.36E-18 * T

        self.n = n
        self.number_of_runs = number_of_runs

    def update_pos(self):
        random_molecule = np.random.choice(self.gas)
        old_state_x = random_molecule.x
        old_state_y = random_molecule.y
        old_radius = np.linalg.norm(np.array([old_state_x,old_state_y])-np.array([self.center[0],self.center[1]]))
        energy_old = np.pi * self.pressure * old_radius**2
        random_molecule.x = old_state_x + (random.random() * self.small_radius) * np.random.choice([-1, 1])
        random_molecule.y = old_state_y + (random.random() * self.small_radius) * np.random.choice([-1, 1])
        new_radius =  np.linalg.norm(np.array([random_molecule.x,random_molecule.y])-np.array([self.center[0],self.center[1]]))
        energy_new = np.pi * self.pressure * new_radius**2
        if energy_new - energy_old <= 0.0:
            return random_molecule
        elif np.exp((-1.0 * (energy_new - energy_old)) / self.kbT) > random.random():
            return random_molecule
        else:
            random_molecule.x = old_state_x
            random_molecule.y = old_state_y
            return random_molecule

    def monte_carlo_step(self):
        gas = []
        for molecule in range(0, self.n):
            gas.append(self.update_pos())
        self.gas = np.array(gas)
        #self.center = centroid(self.gas)
        return self.gas

    def simulation(self):
        volume = []
        for run in range(self.number_of_runs):
            step_gas = self.monte_carlo_step()
            step_centroid = centroid(step_gas)
            step_radius = max([np.linalg.norm(np.array([gas.x,gas.y])-np.array([step_centroid[0],step_centroid[1]]))
                               for gas in step_gas])
            step_volume = np.pi * step_radius**2
            volume.append(step_volume)
        return volume


Gas = IdealGas(500, 50, 1000, 300, 100)
vol = Gas.simulation()
plt.plot(vol)
plt.show()

Upvotes: 0

Views: 529

Answers (1)

gaFF
gaFF

Reputation: 757

You only allow your molecules to move if the new radius is inferior to the old radius:

if energy_new - energy_old <= 0.0:

is equivalent to:

np.pi * self.pressure * new_radius**2 <= np.pi * self.pressure * old_radius**2

that is:

abs(new_radius) <= abs(old_radius)

So all molecules goes to the centroïd.

To me your hypothesis are too strong: you fix temperature, pressure and number of molecules. According to the ideal gas equation, it means volumne v=nRT/p is constant too. If the volume can change, then pressure or temperature has to change. In your simulation, allowing pressure to change would mean that the product of pressure and volume is constant, so that energy is constant, so molecules can move freely in an arbitraty large volume.

By the way I think molecules should be initialized with:

(random.random() - 0.5) * full_radius

so that the occupy all the plane around zero.

Upvotes: 1

Related Questions