Update file Checkpoint6.ipynb

88e15cc9 · cprutean · 38fe6fe4 · 88e15cc9
Commit 88e15cc9 authored 7 months ago by cprutean
--- a/Checkpoints/Checkpoint6.ipynb
+++ b/Checkpoints/Checkpoint6.ipynb
@@ -24,7 +24,7 @@
    "\n",
    "$$ A=4r^ 2 \\frac{N_i}{N} $$\n",
    "\n",
-    "Noting that A=$\\pi r^2$, we get an estimate for the value of π given by,\n",
+    "Noting that $A=\\pi r^2$, we get an estimate for the value of π given by,\n",
    "\n",
    "$$ π≈4 \\frac{N_i}{N} $$\n",
    "\n",

 %% Cell type:markdown id:51574bca tags:
 ## Checkpoint 6
 ### Aim
 To write a Python program to calculate the volume of molecules made up of overlapping spherical atoms; it is also an example in using OOP to simpify the task. Don't get frightened, for this checkpoint you will be given part of the code to help you progress (***this also simulates pretty well how multi-person software development works***).
 ## Monte Carlo estimate of π
 Consider a circle of radius r within a square of side $2r \times 2r$ as shown in the figure below.
 ![dynamicpng-width_300.png](dynamicpng-width_300.png)
 ### Calculation of π by estimation of area.
 If we randomly pick N points evenly distributed within the square and $N_{i}$ of them are inside the circle, then an estimate of the area of the circle A is given by the fraction falling inside the circle multiplied by the area of the square, so
 $$ A=4r^ 2 \frac{N_i}{N} $$
-Noting that A=$\pi r^2$, we get an estimate for the value of π given by,
+Noting that $A=\pi r^2$, we get an estimate for the value of π given by,
 $$ π≈4 \frac{N_i}{N} $$
 So provided the selection of points is evenly distributed over the whole square, than as N becomes large you will get a good estimate for $\pi$. This is the basis of Monte Carlo integration and can be used to estimate the internal area of any closed object.
 See an example Python code that uses the random.uniform() function to form random points:
 - Code example to calculate $\pi$ [here](../CodeExamples/PiExample.ipynb)
 Take a look at this and experiment, you will find that for N>200,000 you get a reasonable estimate for $\pi$, and that as N gets larger the estimate gets better.
 This is clearly a "very bad" way to calcuate π, but is very useful for finding the area of overlapping objects in two-dimensions as shown in the figure below; which is otherwise analytically tricky.
 ![dynamicpng-width_250.png](dynamicpng-width_250.png)
 This can also easily be extended to three-dimensions so for example overlapping spheres as shown below.
 ![dynamicpng-width_2502.png](dynamicpng-width_2502.png)
 Here the method is to form a bounding box of dimensions Vb=a×b×c that totally encloses the object. If N points are in the bounding box and N$_i$ are inside the object, then the volume of the object is simply given by
 $$ V_0=V_b \frac{N_i}{N} $$
 This technique is used in computer simulations of molecules to calculate the volume of molecules that are made up of overlapping spherical atoms.
 #### Code Structure
 This is rather tricky program but it can be vastly simplied by use of OOP; the suggested structure is
 1. Form an Atom class by extending the Vector3d class given by Vector3d Class
 2. The Atom class should add an extra internal variable radius and have a boolean method isInside(v) which takes argument v being a Vector3d and returns True if v is inside the atom. (Note the distance(v) method of Vector3d is useful here.)
 3. Form a Molecule class that has a list of objects of type Atom and a boolean method isInside(v), which takes argument v being a Vector3d and returns True if v is inside the molecule.
 - You also need to consider how to read in the location of the atoms; this is similar to Checkpoint 4.
 4. Form a BoundingBox class where the constructor takes a populated Molecule and finds and stores internally the minima / maxima of the box round the molecule (you need 6 floats). This class should have two methods:
 - volume() that gets the volume of the box.
 - getRandomPoint() which should return a random Vector3d that is randomly distributed inside the box.
    (This is the most “messy” bit of the code to write.)
 If you have this structure, then what your main program has to do is:
 - Ask for the file containing the Atom data and read this into a Molecule.
 - Construct a BoundingBox from the Molecule
 - Go round a while loop,
 - Getting a random Vector3d from the BoundingBox
 - Testing if it isInside() the Molecule and incrementing a counter if it is.
 - Finally calculating the estimate for the Molecule volume and printing it out.
 It is also useful to monitor the progress of the estimate as was done in the example code for estimation of π to check how the integration is going.
 #### Format of Data
 There are two data files suppied (in the current directory, so they can be read in and executed directly upon here), these being :
 - [oxygen.data](oxygen.data) : single oxygen atom
 - [water.data](water.data) : water molecule
 each consists of a comma delimited CSV file with one atom per line in the format,
 Name , r , x , y , z
 with comment lines delimited by # being the same format as coded in [Parsing Input](../CourseNotes/ParsingInput.ipynb) and used in [Checkpoint 4](Checkpoint4.ipynb).
 #### Computing Task
 Write a Python program based on the above structure to:
 1. Read in the data file defining the atom locations.
 2. Prompt the user for the number of points to run in the calculation.
 3. Calculate the volume of the molecule using Monte Carlo integration and output the final value.
 4. Optionally plot a graph of the convergence of the algorithm, being the estimate against number of points. The graph should have 100 data points at regular intervals, see Pi example for how to to this easily.
 5. Test your code with a single oxygen atom and with water. ***For a pass, only the atom is strictly necessary to run.***
 > #### Note
 > The volume of a single oxygen atom of radius 1.4 Å is 11.494 Å$^3$. The volume of water is $\approx$ 15.5 Å$^3$. Use these values to check your results before calling a demonstrator.
 6. Have another careful read of [final weekly task](../WeeklyTasks/Week6.ipynb).
 %% Cell type:code id:9376506d tags:
 ``` python
 """"
    A more developed Vector class for simple 3d vectors
 """
 import math
 import random
 class Vector3d(object):
    """
    Simple vector 3d class
    """
    def __init__(self, x, y, z):
        """
        Constructor to form a vector
        """
        self.x = float(x)
        self.y = float(y)
        self.z = float(z)
    def __str__(self):
        """
        Method to format a vector as a string (implments str())
        """
        return "[" + str(self.x) + ", " + str(self.y) + ", " + str(self.z) + "]"
    def copy(self):
        """
        Method to return a copy of current vector
        """
        return Vector3d(self.x,self.y,self.z)
    def magnitude(self):
        """
        Return the magnitude
        """
        return math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z)
    def distance(self,b):
        """
        Distance from self to vector b
        """
        dx = self.x - b.x
        dy = self.y - b.y
        dz = self.z - b.z
        return math.sqrt(dx*dx + dy*dy + dz*dz)
    def normalise(self):
        """
        Normalise current vector so its magnitude is 1.0 (unit vector)
        """
        mag = self.magnitude()
        if mag != 0.0 :         # if not already zero
            self.x /= mag       # dveive throgh by mag
            self.y /= mag
            self.z /= mag
    def add(self,b):
        """
        Add a vector b to self and return a new vector.
        """
        x = self.x + b.x
        y = self.y + b.y
        z = self.z + b.z
        return Vector3d(x,y,z)
    def addTo(self,b):
        """
        Add a vector b to self in place.
        """
        self.x += b.x
        self.y += b.y
        self.z += b.z
    def dot(self,b):
        """
        Form dot product between self and vector b
        """
        return self.x*b.x + self.y*b.y + self.z*b.z
    def cross(self,b):
        """
        Form cross product between self and vector b.
        """
        x = self.y*b.z - self.z*b.y
        y = self.z*b.x - self.x*b.z
        z = self.x*b.y - self.y*b.x
        return Vector3d(x,y,z)
 def main():
    """
    Simple test program to test the vectors
    """
    a = Vector3d(1,2,3)
    b = Vector3d(4,5,6)
    print("Magnitude of a is " + str(a.magnitude()))
    c = a.cross(b)           # Form cross product
    #
    #         the str() function call the __str__(self) method
    print("Cross product of a x b is : " + str(c))
 #main()
 ```
 %% Cell type:code id:e0aad9ca tags:
 ``` python
 ### Atom class
 class Atom(Vector3d):
    """
    The Atom class to extends the Vector3d that hold position
    """
    def __init__(self,x,y,z,r):
        """
        Constructor to form an atom with x,y,z,r parameters
        """
        # Your code here
    def __str__(self):
        """
        Implemnts str()
        """
        # Your code here
    def isInside(self,v):
        """
        isInside method, returns True of v in inside atom
        """
        # Your code here
 ```
 %% Cell type:code id:fccc8a4e tags:
 ``` python
 ### Molecule class
 class Molecule(list):
    """
    Molecule, being a list to hold Atoms
    """
    def __init__(self,*args):
        """
        Constructor to form a Molecule with optional arguments each one assumed to be an Atom
        """
        list.__init__(self)      # Init underlying list
        for a in args:           # Append any args to list
            self.append(a)
    def __str__(self):
        """      Format string of all atoms
        """
        s = ""
        for a in self:
            s += str(a) + "\n"
        return s
    def inInside(self, v):
        """
        Method to check if vector is inside moleule.
        """
        # Your code here
        return False           # if here then outside
    def fromFile(self,file):
        """
        Read Atoms from a file
        """
        # Your code here
 ```
 %% Cell type:code id:663fe126 tags:
 ``` python
 ### BoundingBox class
 class BoundingBox(object):
    """
    Class to define a  bounding box class round the molecule
    """
    def __init__(self, molecule):
        """
        The Constructor, which takes a molecule
        """
        self.xmin = float("inf")             # Initially set limits to +/- infty
        self.ymin = float("inf")
        self.zmin = float("inf")
        self.xmax = float("-inf")
        self.ymax = float("-inf")
        self.zmax = float("-inf")
        for a in molecule:                  # update limits for each atom in molecule
           # Your code here
    def volume(self):
        """
        Get the volume of the box
        """
        # Your code here
    def getRandomPoint(self):
        """
            Get a random point in the bounding box as a Vector3d.
        """
        # Your code here
 ```
 %% Cell type:code id:3e01a08e tags:
 ``` python
 """    Main program to do the simulation
 """
 file = open(str(input("File : ")),"r")   # Open file
 mol = Molecule().fromFile(file)              # Read in molecule
 print("Molecule is \n" + str(mol))           # Print out info
 box = BoundingBox(mol)                       # Form Bounding box
 maxpoint = float(input("Number of points : ")) # Max number of points
 plotInterval = maxpoint/100                    # Plot internal for monitoring
 # Your code here
 while p < maxpoint:                          # Loop counting internal points
    # Your code here
 ...
 print("Final estimate is : " + str(estimate))
 plt.plot(xData,yData)                       # Draw graph (with default plot)
 plt.ylim(0.9*estimate,1.1*estimate)
 plt.title("Estimate of volume for " + str(maxpoint) + " points.")
 plt.show()
 ```