diff --git a/report/ChainMCReport.tex b/report/ChainMCReport.tex
new file mode 100644
index 0000000..656fe11
--- /dev/null
+++ b/report/ChainMCReport.tex
@@ -0,0 +1,260 @@
+\documentclass{article}
+\usepackage{amsmath}
+\usepackage{float}
+\usepackage{amssymb}
+\usepackage{graphicx}
+\usepackage[margin=1in]{geometry}
+\usepackage{multicol,caption,listings}
+\usepackage{cite}
+\usepackage[stable]{footmisc}
+
+\bibliographystyle{plain}
+
+\begin{document}
+
+\title{Polymer Chain Monte Carlo}
+\author{A. Lovell}
+\maketitle
+
+\noindent \textbf{Abstract}  Polymer chains are of great interest in biophysics, and while one might think to use a molecular dynamics calculation to simulate their configurations, for practical reasons, Monte Carlo techniques are more widely used.  In this project, we will use the Rosenbluth algorithm within a Monte Carlo framework to simulate polymer chains of varying lengths around $10^4$ atoms at several temperatures to investigate the correlation between the size of the polymer structure and the number of atoms within the chain.  This will be compared to theoretical models to investigate the systems we can create using a fixed distance between sequential atoms and the Lennard-Jones potential between all other pairs.
+
+\begin{multicols}{2}
+
+\section{Introduction}
+
+Polymers are chains of atoms or molecules of the same type that interact directly between sequential atoms through a spring force and through all pairs via some other long- or short-range potential (given, for example, by hydrogen bonds, van der Waals interactions, or the Lennard-Jones interaction).  The typical number of atoms or molecules in a polymer is between $10^3$ and $10^5$ which can all be of the same type or different types.  \\
+
+One would often think that a molecular dynamics calculation would be used to simulate a polymer chain.  In this case, the movement of the chain could be followed by calculating the forces between each pair of atoms for each time step and then used to update the chain's next position.  However, because of the different types of motion involved in the problem - ranging from fast motion of the individual atomic units to slow motion of the chain as a whole - this process can be inefficient, or even impossible, for long chains.  \cite{PhilNotes} Instead, Monte Carlo simulations are used to discover properties of the polymer.  \\
+
+This report is broken up into the following sections.  In Section \ref{theory}, we will discuss the theory behind the problem, including the interaction between atoms and the algorithm that was used to run the Monte Carlo simulation.  Section \ref{IC} gives a brief discussion of the initial conditions of the simulation and a note on the units used in the calculation.  The results of the simulation are discussed in Section \ref{discuss}.  Finally, in Section \ref{concl}, we provide a summary of our work, along with further progress that could be made on this project in the future.  The appendix then provides some detail about our error calculations.  \\
+
+\section{Theory}
+\label{theory}
+
+In the following section, we will describe some of the concepts that will be needed to simulate a polymer chain as well as briefly explain the path of the simulation code.
+
+\subsection{Interaction}
+
+The interactions between the atoms of polymer chains can be modeled in several ways.  One of the most common ways is to describe the interaction between sequential pairs of atoms as a stiff spring and put in a Lennard-Jones potential between every other pair of atoms.  \cite{PhilNotes}  However, in this project we will take a somewhat simplified approach.  Instead of stiff springs, we will fix the length between each successive pair of atoms (infinitely stiff springs), while keeping the interaction between all other pairs a Lennard-Jones interaction.  There are, of course, more complicated potentials that could be implemented to fold the polymer into more complex shapes (such as proteins or even origami ducks), but these will not be explored here.\\
+
+The Lennard-Jones potential is 
+
+\begin{equation}
+\label{VLJ}
+V_{LJ}(r) = 4\epsilon \left [ \left ( \frac{\sigma}{r} \right ) ^{12} - \left ( \frac{\sigma}{r}) \right ) ^6 \right ]
+\end{equation}
+
+\noindent where $\epsilon$ is the depth of the interaction, $\sigma$ is related to the particle size, and $r$ is the distance between the interacting pairs.  When the two atoms are far enough apart, they can be viewed as non-interacting.  For this reason, we would be able to introduce a cut-off radius beyond which the potential between two particles is taken to be zero.  \\
+
+The form of (\ref{VLJ}) takes into account the short-range repulsive force between the two atoms due to Coulomb repulsion and the longer-range attractive forces from dipole-dipole and dipole-induced dipole forces.  \cite{verlet}  An example of the Lennard-Jones potential can be found in Figure \ref{VLJfig}.
+
+\begin{figure}[H]
+\begin{center}
+\includegraphics[width=\linewidth]{Figures/VLJ.pdf}
+\caption{The shape of the Lennard-Jones potential.  Here, $\sigma=1$ fm (the distance where the potential crosses the r-axis) and $\epsilon =1$ MeV (the maximum depth of the interaction).}
+\label{VLJfig}
+\end{center}
+\end{figure}
+
+\subsection{Algorithm}
+\label{RA}
+
+In this Monte Carlo simulation, we build up the polymer chain by adding one atom at a time.  The algorithm that we use to pick the next atom's position is the Rosenbluth algorithm, proposed in 1956 by Marshall and Arianna Rosenbluth.  As compared to simple sampling - which generates all configurations with equal probability - the Rosenbluth algorithm generates configurations with different probabilities.  These probabilities cause the chain to be self-avoiding and to terminate when it is trapped in a configuration where every next possible addition would be on top of an already-placed atom.  \cite{Rosenbluth}  (This is a more straightforward concept for polymers created on a grid, which have a fixed number of nearest neighbors.)  \\
+
+On a square (or cubic) lattice, this means that each nearest neighbor is tested as to whether or not it is already occupied, and the the probability that the next atom is placed at that site is based on the energy of the potential new configuration using the statistical factor
+
+\begin{equation}
+\label{prob}
+P_i = \frac{1}{Z} e^{-V_i/T}
+\end{equation}
+
+\noindent where
+
+\begin{equation}
+\label{part}
+Z = \sum \limits _{i=1}^{N_a} e^{-V_i/T}
+\end{equation}
+
+\noindent Here, $V_i$ is the potential energy between the $i^{th}$ candidate configuration and the rest of the chain, $N_a$ is the number of unoccupied nearest neighbors, and T is the temperature of the system.  These probabilities are mapped onto the range $[0,1]$, and a random number is picked within that same range.  If this random number falls within the range $P_i$, then the next atom is added at the position from which $V_i$ was calculated.  \\
+
+However, we can also grow our polymer freely, without a grid.  In this case we have to sample from a circle in two-dimensions or a sphere in three-dimensions with the radius of the distance between each successive pair of atoms.  To do this, we sample from $\mathrm{cos}\theta \in [-1,1]$ and $\phi \in [0,2\pi]$ (taking $\theta = \pi/2$ for the two-dimensional case, to fix the third dimension).  This samples over the circumference of a unit circle in two-dimensions and the surface of a sphere in three-dimensions.  Some number, $N_a$, of points are randomly chosen, and new positions are calculated from
+
+\begin{equation}
+\begin{split}
+x^{j}_i & = x^{j-1} + \mathrm{sin}\theta _i \mathrm{cos}\phi _i \\
+y^{j}_i & = y^{j-1} + \mathrm{sin} \theta _i \mathrm{sin} \phi _i \\
+z^{j}_i & = z^{j-1} + \mathrm{cos} \theta _i
+\end{split}
+\end{equation}
+
+\noindent where we are placing the $j$th atom in the $i$th candidate configuration.  The probabilities for each $\textbf{r}_i$ to be chosen is the same as described above with equations (\ref{prob}) and (\ref{part}).
+
+\subsection{Configuration}
+
+The main quantity that we want to explore is the relationship between the size of the system and the number of atoms in the polymer, which should be related to the dimensionality of the simulation (e.g. two-, three-, or more dimensional space).  To do this, we relate the gyration radius, $R_g$, to the number of atoms in the chain, $N$, where 
+
+\begin{equation}
+\label{rg}
+R_g^2 = \frac{1}{N} \sum \limits _{i=1}^{N} \langle (r_i - R_{cm})^2 \rangle
+\end{equation} 
+
+\noindent and 
+
+\begin{equation}
+\label{RNcomparison}
+R \sim N^{\nu}
+\end{equation}
+
+The temperature, $T_{\Theta}$, defines the transition between a chain-like polymer ($T<T_{\Theta}$) and a dense, ball-like polymer ($T>T_{\Theta}$).  For temperatures below $T_{\Theta}$, $\nu = 1/ 2$, at $T_{\Theta}$, $\nu = 1/3$, and above $T_{\Theta}$, $\nu$ scales with the dimensionality of the simulation - namely, $\nu = 3/(d+2)$, where $d$ is the dimensionality of the simulation.  \cite{PhilNotes}\\
+
+$R_g$ is the main observable we will be calculating for this Monte Carlo simulation.  From this, we can also define - or at least put some constraint on - the $\Theta $-temperature, $T_{\Theta}$.
+
+\subsection{Code}
+
+In this subsection, we will briefly describe the structure of the code and where the above ideas are implemented. \\
+
+Initially, two atoms are placed at a distance of one unit length apart from each other.  The potential energy of the system is calculated.  To prepare to place the next atom, some number of points on the unit sphere are randomly chosen from a distribution of $cos\theta d\theta d\phi$.  The potential energy of each of these configurations is calculated, and the Rosenbluth algorithm, described in \ref{RA}, is implemented to choose the position of the next atom.  Once the next position has been chosen, the potential energy of the new system is calculated and the method repeats:  choose points on the unit sphere, implement the Rosenbluth algorithm, calculate the potential energy of the system.  \\
+
+Once the entire chain (a given number of atoms, $N$) is constructed, the radius of gyration, from (\ref{rg}), is calculated for this polymer.  This chain-building process is repeated again and again to build up enough statistics to calculate a mean and and standard deviation (for the error).  
+
+\section{Initial Conditions}
+\label{IC}
+
+The first two atoms of the chain were positioned at (0,0,0) and (1,0,0).  Essentially, these positions are arbitrary (only constrained to be 1 unit length apart), but it gives a good starting point for both the polymer on a grid and in free space.   Because this calculation is not performed on a grid, for each new position, 90 points on the unit sphere are tested.  This more easily prevents the chain from becoming stuck by only having options to place the next atom too close to previously placed atoms.  \\
+
+In this calculation, we also chose $\epsilon =1$ and $\sigma=1$ for the Lennard-Jones potential parameters in (\ref{VLJ}).  In order to compare to data from the lab, we would just need to change our potentials to include realistic numbers for these two values.  Along with this, we also are using units where $k_B=1$.  This gives us temperature in eV, with $1$ eV $\approx 11600$ K.  Thus, even seemingly small temperatures quickly become large - and potentially unphysical.\\
+
+\section{Discussion}
+\label{discuss}
+
+In this section, we present the main result of our calculation.  We were able to vary the temperature of the polymer, from 1.0 eV to 5.0 eV, as well as the number of atoms, $N=5000, 10000, 15000$, to calculate the radius of gyration and see how well our simulation compares to (\ref{RNcomparison}) for known values of $\nu$.  These results are summarized in Table 1 and the method of error, $\sigma_{R_g}$, calculation is discussed in the appendix.  \\
+
+\begin{table*}
+\begin{center}
+\begin{tabular}{| c | c | c | c | c | c | c |}
+\hline & $T=1.0$ eV &  & $T=2.0$ eV &  & $T=5.0$ eV & \\ \hline
+ & \textbf{$log(R_g)$} & \textbf{$\sigma _{R_g}$} & \textbf{$log(R_g)$} & \textbf{$\sigma _{R_g}$} & \textbf{$log(R_g)$} & \textbf{$\sigma _{R_g}$} \\ \hline
+$N=5000$ & 2.5119 & 0.343484 & 3.109664 & 0.312228 & 3.308621 & 0.330394 \\ \hline
+$N=10000$ & 2.783211 & 0.295998 & 3.656107 & 0.335429 & 3.642807 & 0.321192 \\ \hline
+$N=15000$ & 2.994968 & 0.376915 & 3.722428 & 0.343053 & 3.872835 & 0.318437 \\ \hline
+\end{tabular}
+\label{results}
+\caption{Results for a three-dimension polymer chain in free space (not on a grid) for three different temperatures at three different polymer lengths.  Both the mean value ($log(R_g)$) and the standard deviation ($\sigma_{R_g}$) for each set of points is given.  Here, and throughout the report "log" means log base e - the natural log.}
+\end{center}
+\end{table*}
+
+\begin{figure}[H]
+\begin{center}
+\includegraphics[width=\linewidth]{Figures/T1plot.pdf}
+\caption{Results for $T=1.0$ eV.  Black dots give the Monte Carlo values calculated for this project, and the dotted-dashed line gives the best fit line through the simulation points.  The solid, dotted, and dashed lines show possible, interesting values of $\nu$ from (\ref{RNcomparison}).  Each of these three lines has the same y-intercept as the best fit line; thus, if one of these lines was the correct description of the simulation, it would go through the black dots and be aligned with the best fit line.}
+\label{T1fig}
+\end{center}
+\end{figure}
+
+In Figures \ref{T1fig}, \ref{T2fig}, and \ref{T5fig}, we show the results of the simulation plotted along with the best fit line for each set of data.  These fits are given in Table 2.  \\
+
+\begin{figure}[H]
+\begin{center}
+\includegraphics[width=\linewidth]{Figures/T2plot.pdf}
+\caption{Same as Figure \ref{T1fig} but for $T=2.0$ eV.  We can see that $\nu = 3/5$ best describes the simulation.}
+\label{T2fig}
+\end{center}
+\end{figure}
+
+From these fits, we find that $\nu = 0.4345 \pm 0.03440$ for $T=1.0$ eV, $\nu = 0.5826 \pm 0.1642$ for $T = 2.0$ eV, and $\nu = 0.5102 \pm 0.02239$ for $T = 5.0$ eV.\footnote{These errors are discussed in more detail in the Appendix.}  The first two cases match up well with a temperature just below $T_{\Theta}$ for $T=1.0$ eV and for a temperature above $T_{\Theta}$ that follows the pattern $\nu = 3/(d+2)$ for $T = 2.0$ eV.  (Recall, for this simulation $d=3$ so $\nu = 0.6$.)  However, we see that for $T = 5.0$ eV, $\nu \approx 0.5$ which would seem to indicate that this is the $\Theta$-temperature. This cannot be the case.  To postulate an explanation for this, we note that $5.0$ eV $\approx 58000$ K, which is an order of magnitude hotter than the coolest part of the sun.  \cite{sun}  We assume that at these temperatures, our simple model is no longer valid, and we get nonsensical results from our calculations.  \\
+
+From Table 2, we can also see that the $\chi^2$ value for each of these fits is extremely low.  This is a product of the fact that only three data points are being fit with by linear regression, causing a misleadingly small $\chi^2$ value.  Ideally, we want a $\chi^2$ value of 1 for every degree of freedom in the problem (e.g. $\chi ^2/dof =1$) which would indicate a perfect fit.  A very low $\chi ^2$, such as we have here, indicates that we have under-constrained our model by allowing too many free parameters for the number of data points. Thus, having run the simulation for more $N$ and $T$ values would serve to better constrain our fit, as well as give us more confidence that the parameters we calculated adequately described the results of our simulation.\\
+
+\begin{figure}[H]
+\begin{center}
+\includegraphics[width=\linewidth]{Figures/T5plot.pdf}
+\caption{Same as Figure \ref{T1fig} but for $T=5.0$ eV.  We can see that $\nu = 1/2$ best describes the simulation.  This unexpected result is discussed more within the text.}
+\label{T5fig}
+\end{center}
+\end{figure}
+
+\begin{table*}
+\begin{center}
+\begin{tabular}{| c | c | c | c | c |}
+\hline \textbf{T (eV)} & \textbf{m} & \textbf{b} & \textbf{s$_m$} & \textbf{$\chi ^2$} \\ \hline
+1.0 & 0.4345 & -1.1968 & 0.03440 & 0.007327 \\ \hline
+2.0 & 0.5826 & -1.8139 & 0.1642 & 0.1483 \\ \hline
+5.0 & 0.5102 & -1.042 & 0.02239 & 0.003000 \\ \hline
+\end{tabular}
+\caption{Slopes ($m$) and y-intercepts ($b$) for the best fit lines for the simulation data in Table 1.  Best fit lines are given as $log(R_g) = m*log(N) + b$.  Here $log$ indicates the natural log.  In this case, the slope $m$ gives the value for $\nu$ from (\ref{RNcomparison}).  $s_m$ gives the error on the fit - only describing the difference between the values calculated from the fit and the simulation values, and $\chi^2$ takes into account, somewhat, the error from the simulation.}
+\end{center}
+\end{table*}
+
+We can see from Figures \ref{T1fig}, \ref{T2fig}, and \ref{T5fig} that the error bars are significant for the range of data displayed.  These large errors did not decrease with the number of statistics, a trend which can be seen in Table 3.  From these figures, we can see that these large errors contribute to some uncertainty around the slope of the best fit lines, as the calculated value for $\nu$ does not take into account the error from the original simulation.  This uncertainty makes it unclear as to whether or not we can confidently describe the behavior of the chain.  Even though the best fit lines have a very low $\chi ^2$ value, this value is misleading, as described previously.    \\
+  
+\begin{table*}
+\begin{center}
+\begin{tabular}{| c | c | c | c |}
+\hline & $T=1.0$ eV & $T = 2.0$ eV & $T = 5.0$ eV \\ \hline
+$N=5000$ & 64 (0.343484) & 500 (0.312228) & 1000 (0.330394) \\ \hline
+$N=10000$ & 43 (0.295998)  & 500 (0.335429) & 810 (0.321192) \\ \hline
+$N=15000$ & 45 (0.376915) & 290 (0.343053) & 1147 (0.318437) \\ \hline
+\end{tabular}
+\caption{Number of chains produced for each combination of number of atoms in the chain (first column) and temperature (top row).  The number in parentheses gives the error calculated for each of these points, also given in Table 1.  It is clear to see that the errors do not decrease as the number of calculations increases, even though, in some cases, the number of simulations increases by orders of magnitude.}
+\end{center}
+\end{table*}
+
+Still, in the region where our simulation is valid, we have calculated reasonable results for our three-dimensional simulation of a polymer chain, even if the error analysis leaves somewhat of a question as to the behavior of the chain at our various temperatures of interest.   \\
+
+\section{Conclusion}
+\label{concl}
+
+In this work, we were able to construct a Monte Carlo simulation of a three-dimensional polymer chain of approximately $10^4$ atoms at temperatures of 1.0 eV, 2.0 eV, and 5.0 eV.  The polymer was simulated without being confined to a grid.  By comparing the natural logarithm of the gyration radius and number of atoms in the chain, we were able to discover properties of the chain itself and find the exponent, $\nu$.  We found that $\nu = 0.4345 \pm 0.03440$ for $T = 1.0 $ eV, $\nu = 0.5826 \pm 0.1642$ for $T =2.0$ eV, and $\nu = 0.5102 \pm 0.02239$ for $T=5.0$ eV.  Through this, we were able to compare to known behavior above and below the $\Theta$-temperature.  \\
+
+However, there is still future work that can be accomplished.  While we calculated a polymer chain in three-dimensions, it would also be interesting to calculate higher (or lower) dimensional chains and examine the scaling of the size of the chain with the number of atoms added to the chain, $\nu$.  The same scaling ($\nu = 1/3$ below $T_\Theta$, $\nu = 1/2$ at $T_\Theta$, and $\nu = 3/(d+2)$ above $T_\Theta$) should hold, but it would be an interesting project nonetheless to generalize the code.  It would also be interesting to examine the differences between this type of self-avoiding walk that is not confined to a grid and a self-avoid walk confined to a cubic grid (or other shape grid).  From this, we could see how the dimensionality affects various grids on which the chain is created.  For all of this, we could also run the calculation in smaller temperature steps to find the $\Theta$-temperature.  \\
+
+There are also several thermodynamical quantities that could be calculated from this system, including specific heat and thermal energy.  We could also change the interaction between the atoms within each chain.  Although we used a fixed bond length between each successive atom and a Lennard-Jones interaction between all pairs of atoms, there are other interactions that could describe a more realistic polymer, including using a stiff spring to model the interaction between successive atoms instead of a fixed bond length.  We also could have changed the strength (and type) of the interaction between various pairs to see what kinds of shapes we could form with our polymer chain.  \\
+
+There have also been developments within the last two decades of the Pruned-enriched Rosenbluth method which combines the Rosenbluth method with recursive enrichment.  This algorithm allows for simulations of up to $N=10^6$ atoms in a single chain with high statistics, applicable to calculations both on and off of a lattice.  \cite{PERM}  In future work, this could be an interesting algorithm to implement, especially if it is more efficient and can give higher statistics in a shorter amount of time than the original Rosenbluth algorithm.  \\
+
+Still, without all of this, we were able to model a polymer chain of $10^4$ atoms in three-dimensions without a grid.  \\
+
+\appendix
+
+\section{Error Analysis}
+
+Just as experimentally measured data points should quote errors - whether due to systematics, timing, or unknown quantities - theoretical simulations and calculations should also give error bars.  To do this, we run our simulation several times to find a mean value and a standard deviation.  The mean is calculated by
+
+\begin{equation}
+\bar x = \frac{1}{N} \sum \limits _{i=1}^N x_i
+\end{equation}
+
+\noindent where $\bar x$ is the mean value, $x_i$ are the values of a given quantity from each calculation (for instance the length of the chain), and $N$ is the number of simulations that were performed.\\
+
+To find the error on each data point, we average the standard deviation for each run.
+
+\begin{equation}
+\label{stddev}
+\sigma _x^2 = \frac{1}{N} \sum \limits _{i=1} ^N (x_i - \bar x)^2
+\end{equation}
+
+\noindent Here, $\sigma _x$ is the error for the variable $x$ from each calculation.  We can see that the errors decrease as $1/\sqrt{N}$.  \\
+
+Because we are calculating a fit to data, we should also report a $\chi ^2$ value for that fit, to assess its "goodness" - or the quality of the fit when compared to data.  \cite{Nunes}  In this case, we would calculate the $\chi^2$ of each of the linear fits that are made to the simulation data.
+
+\begin{equation}
+\chi ^2 = \sum \limits _{i=1} ^N \frac{(x^{sim} - x^{fit})^2}{\sigma _x ^2}
+\end{equation}
+
+\noindent In this case, $x^{sim}$ is the parameter value from the simulation, $x^{fit}$ is the parameter value from the fit through the simulation data, and $\sigma _x ^2$ is the square of the standard deviation (the variance) from (\ref{stddev}). \\
+
+Although this does not give us an error on the fits that are being made, it does give some indication as to the quality of the fit.\\
+
+Using similar ideas, errors for the fits can be calculated (although these errors still do not take into account the error from the simulations).  For these errors (such as those in Table 1), we take
+
+\begin{equation}
+s^2 = \frac{1}{n} \sum \limits _{i=1} ^{n} (y^{data}_i - y^{fit}_i)^2
+\end{equation}
+
+\noindent where $s$ is the error on the fit, $n$ is the number of data points fit, $y^{data}$ is the result of the simulation, and $y^{fit}$ is the result from the fit (both calculated at the same point).  \cite{staterr}
+
+\end{multicols}
+
+\bibliography{chainbib}
+
+\end{document}
\ No newline at end of file
diff --git a/report/chainbib.bib b/report/chainbib.bib
new file mode 100644
index 0000000..4131872
--- /dev/null
+++ b/report/chainbib.bib
@@ -0,0 +1,78 @@
+@article{PhilNotes,
+      title = {Monte Carlo Simulation of Polymers:  Coarse-Grained Models},
+      author = {J\"org Baschnagel, Joachim P. Wittmer, Hendrik Meyer},
+      publisher = {},
+      year = 2004,
+      journal = {John von Neumann Institute for Computing},
+      volume = 23,
+      pages = {83-140},
+      url = {http://www.pa.msu.edu/~duxbury/courses/phy480/BaschnagelReview2004.pdf} 
+}
+@article{Rosenbluth,
+      title = {Monte Carlo Calculation of the Average Extension of Molecular Chains},
+      author = {Marshall N. Rosenbluth and Arianna W. Rosenbluth},
+      publisher = {},
+      year = 1955,
+      journal = {The Journal of Chemical Physics},
+      volume = 23,
+      pages = 356
+}
+@article{sun,
+      title = {The Photosphere},
+      author = {University of St. Andrews},
+      publisher = {AWH/JOC},
+      year = 1995,
+      journal = {Lecture notes at University of St. Andrews},
+      url = {www-solar.mcs.st-andrews.ac.uk/~alan/sun_course/Introduction/Photosphere.html} 
+}
+@article{stuff,
+      title = {From Rosenbluth Sampling to PERM - rare event sampling with stochastic growth algorithms},
+      author = {Thomas Prellberg},
+      publisher = {Oldenburg, Germany},
+      year = 2012,
+      journal = {Lecture given at the International Summer School, Modern Computational Science 2012 - Optimization},
+      url = {http://www.maths.qmul.ac.uk/~tp/papers/pub084pre.pdf} 
+}
+@book{Nunes,
+	title = {Nuclear Reactions for Astrophysics},
+	author = {Ian J. Thompson and Filomena M. Nunes},
+	publisher = {Cambridge University Press},
+	year = 2009,
+	address = {New York}
+}
+@article{verlet,
+  title = {Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules},
+  author = {Verlet, Loup},
+  journal = {Phys. Rev.},
+  volume = {159},
+  issue = {1},
+  pages = {98--103},
+  numpages = {0},
+  year = {1967},
+  month = {Jul},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRev.159.98},
+  url = {http://link.aps.org/doi/10.1103/PhysRev.159.98}
+}
+@article{PERM,
+  title = {Pruned-enriched Rosenbluth method: Simulations of $\theta$ polymers of chain length up to 1 000 000},
+  author = {Grassberger, Peter},
+  journal = {Phys. Rev. E},
+  volume = {56},
+  issue = {3},
+  pages = {3682--3693},
+  numpages = {0},
+  year = {1997},
+  month = {Sep},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevE.56.3682},
+  url = {http://link.aps.org/doi/10.1103/PhysRevE.56.3682}
+}
+@book{staterr,
+	title = {Numerical Recipes in FORTRAN:  The Art of Scientific Computing},
+	author = {W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling},
+	publisher = {Cambridge University Press},
+	year = 1992,
+	edition = 2nd
+	address = {Cambridge}
+}
diff --git a/report/figures/T1plot.pdf b/report/figures/T1plot.pdf
new file mode 100644
index 0000000..70e9498
Binary files /dev/null and b/report/figures/T1plot.pdf differ
diff --git a/report/figures/T2plot.pdf b/report/figures/T2plot.pdf
new file mode 100644
index 0000000..a2f60a4
Binary files /dev/null and b/report/figures/T2plot.pdf differ
diff --git a/report/figures/T5plot.pdf b/report/figures/T5plot.pdf
new file mode 100644
index 0000000..91ab4eb
Binary files /dev/null and b/report/figures/T5plot.pdf differ
diff --git a/report/figures/VLJ.pdf b/report/figures/VLJ.pdf
new file mode 100644
index 0000000..2562799
Binary files /dev/null and b/report/figures/VLJ.pdf differ
diff --git a/report/figures/fig_1.pdf b/report/figures/fig_1.pdf
deleted file mode 100644
index 629164b..0000000
Binary files a/report/figures/fig_1.pdf and /dev/null differ
diff --git a/report/report.tex b/report/report.tex
deleted file mode 100644
index 72d1253..0000000
--- a/report/report.tex
+++ /dev/null
@@ -1,156 +0,0 @@
-\documentclass[]{article}
-
-\usepackage{amsmath}  % AMS math package
-\usepackage{amssymb}  % AMS symbol package
-\usepackage{bm}       % bold math
-\usepackage{graphicx} % Include figure files
-\usepackage{dcolumn}  % Align table columns on decimal point
-\usepackage{multirow} % Multirow/column tables
-\usepackage{hyperref} % Hyperlinks
-
-\begin{document}
-
-\title{Manuscript Title:\\with Forced Linebreak}% Force line breaks with \\
-\author{Ann Author}
-\date{\today}% It is always \today, today, but you can specify other dates manually 
-\maketitle
-
-\begin{abstract}
-An article usually includes an abstract, a concise summary of the work
-covered at length in the main body of the article. 
-\end{abstract}
-
-
-\section{First-level heading} %Title for the section
-\label{sec:level1} %Label for the section, to be used for referencing in other parts of the document
-
-This sample document demonstrates some common uses of \LaTeX\ in documents you'll write for ICCP. Further information can be found online at the \href{http://en.wikibooks.org/wiki/LaTeX}{\LaTeX\ wikibook} or in the distribution documentation. 
-
-When we refer to commands in this example file, they always have their required formal arguments in normal \TeX{} format. In this format, \verb+#1+, \verb+#2+, etc. stand for required author-supplied arguments to commands. For example, in \verb+\section{#1}+ the \verb+#1+ stands for the title text of the author's section heading, and in \verb+\title{#1}+ the \verb+#1+ stands for the title text of the paper.
-
-Line breaks in section headings at all levels can be introduced using \textbackslash\textbackslash. A blank input line tells \TeX\ that the paragraph has ended. 
-
-\subsection{\label{sec:level2} Second-level heading: Formatting}
-
-This file may be formatted in either the \texttt{preprint} or \texttt{reprint} style. \texttt{reprint} format mimics final journal output. The paper size may be specified with the option \texttt{letter}. These options are inserted in the square brackets inside the \texttt{\textbackslash documentclass[]\{article\}} command.
-
-\section{Math and Equations}
-Inline math may be typeset using the \verb+$+ delimiters. Bold math symbols may be achieved using the \verb+bm+ package and the \verb+\bm{#1}+ command it supplies. For instance, a bold $\alpha$ can be typeset as \verb+$\bm{\alpha}$+ giving $\bm{\alpha}$. Fraktur and Blackboard (or open face or double struck) characters should be typeset using the \verb+\mathfrak{#1}+ and \verb+\mathbb{#1}+ commands respectively. Both are supplied by the \texttt{amssymb} package. For
-example, \verb+$\mathbb{R}$+ gives $\mathbb{R}$ and \verb+$\mathfrak{G}$+ gives $\mathfrak{G}$
-
-In \LaTeX\ there are many different ways to display equations, and a few preferred ways are noted below. Displayed math will center by default. Use the class option \verb+fleqn+ to flush equations left.
-
-Below we have numbered single-line equations; this is generally the most common type of equation in \LaTeX: 
-\begin{equation} 
-\label{eq:one} %Label for the equation, to be used for referencing in other parts of the document
-  i \hbar \frac{\partial \Psi}{\partial t} = -\frac{\hbar^2}{2m}\nabla^2 \psi + U(\mathbf{x}) \psi
-\end{equation}
-
-When the \verb+\label{#1}+ command is used [cf. input for Eq.~(\ref{eq:one})], the equation can be referred to in text without knowing the equation number that \TeX\ will assign to it. Just use \verb+\ref{#1}+, where \verb+#1+ is the same name that was used in the \verb+\label{#1}+ command.
-
-Unnumbered single-line equations can be typeset using the \verb+equation*+ environment:
-\begin{equation*}
-  \hat{f}(\xi) = \int_{-\infty}^\infty f(x) e^{-2\pi i x \xi} \, \mathrm{d}x
-\end{equation*}
-
-\subsection{Multiline equations}
-
-You'll generally want to use the \verb+align+ environment for multiline equations. Use the \verb+\nonumber+ command at the end of each line (again denoted with \textbackslash\textbackslash) to avoid assigning a number, or \verb+align*+ to disable numbering entirely:
-\begin{align}
-  \sin(2\theta) &= 2 \sin(\theta) \cos(\theta) \nonumber \\
-                &= \frac{2\tan(\theta)}{1 + \tan^2(\theta)}
-\end{align}
-Note how alignment occurs at the (unprinted) \verb+&+ character.
-
-Do not use \verb+\label{#1}+ on a line of a multiline equation if \verb+\nonumber+ is also used on that line. Incorrect cross-referencing will result. Notice the use \verb+\text{#1}+ for using a Roman font within a math environment.
-
-\section{Cross-referencing}
-\LaTeX\ will automatically number such things as sections, footnotes, equations, figure captions, and table captions for you. In order to reference them in text, use the \verb+\label{#1}+ and \verb+\ref{#1}+ commands\footnote{Check out the \texttt{cleveref} (one r) package, too!}. To reference a particular page, use the \verb+\pageref{#1}+ command.
-
-The \verb+\label{#1}+ should appear within the section heading (or immediately after), within the footnote text, within the equation, or within the table or figure caption. The \verb+\ref{#1}+ command is used in text at the point where the reference is to be displayed.  Some examples: Section~\ref{sec:level1} on page~\pageref{sec:level1}, Table~\ref{tab:table1}, and Fig.~\ref{fig:epsart}.
-\begin{figure}[b]
-  \centering
-  \includegraphics{figures/fig_1}% Imports a figure - does not automatically scale
-  \caption{\label{fig:epsart} A figure caption. The figure captions are automatically numbered.}
-\end{figure}
-
-\section{Floats: Figures, Tables, etc.}
-Figures (images) and tables are usually allowed to ``float'', which means that their placement is determined by \LaTeX, while the document is being typeset. 
-
-Use the \texttt{figure} environment for a figure, the \texttt{table} environment for a table. In each case, use the \verb+\caption+ command within to give the text of the figure or table caption along with the \verb+\label+ command to provide a key for referring to this figure or table. Insert an image using either the \texttt{graphics} or \texttt{graphix} packages, which define the \verb+\includegraphics{#1}+ command. (The two packages differ in respect of the optional arguments used to specify the orientation, scaling, and translation of the image.) To create an alignment, use the \texttt{tabular} environment. 
-
-\begin{table}
-  \centering
-  \begin{tabular}{ |l|l|l| }
-  \hline
-  \multicolumn{3}{ |c| }{Team sheet} \\
-  \hline
-  Goalkeeper & GK & Paul Robinson \\ \hline
-  \multirow{4}{*}{Defenders} & LB & Lucus Radebe \\
-   & DC & Michael Duburry \\
-   & DC & Dominic Matteo \\
-   & RB & Didier Domi \\ \hline
-  \multirow{3}{*}{Midfielders} & MC & David Batty \\
-   & MC & Eirik Bakke \\
-   & MC & Jody Morris \\ \hline
-  Forward & FW & Jamie McMaster \\ \hline
-  \multirow{2}{*}{Strikers} & ST & Alan Smith \\
-   & ST & Mark Viduka \\
-  \hline
-  \end{tabular}
-  \caption{\label{tab:table1}A table, demonstrating the use of the \texttt{multirow} package for spanning rows and columns.}
-\end{table}
-
-The best place for the \texttt{figure} or \texttt{table} environment is immediately following its first reference in text; this sample document illustrates this practice for Fig.~\ref{fig:epsart}, which shows a figure that is small enough to fit in a single column. In exceptional cases, you will need to move the float earlier in the document: \LaTeX's float placement algorithms need to know about a full-page-width float sooner.
-
-The content of a table is typically a \texttt{tabular} environment, giving rows of type in aligned columns. Column entries separated by \verb+&+'s, and 
-each row ends with \textbackslash\textbackslash. The required argument for the \texttt{tabular} environment specifies how data are aligned in the columns. 
-For instance, entries may be centered, left-justified, right-justified, aligned on a decimal point. 
-
-Extra column-spacing may be be specified as well, although \LaTeX sets this spacing so that the columns fill the width of the table. Horizontal rules are typeset using the \verb+\hline+ command. Rows with that columns span multiple columns can be typeset using the \verb+\multicolumn{#1}{#2}{#3}+ command (for example, see the first row of Table~\ref{tab:table1}).
-
-\appendix
-
-\section{Appendixes}
-
-To start the appendixes, use the \verb+\appendix+ command. This signals that all following section commands refer to appendixes instead of regular sections. Therefore, the \verb+\appendix+ command should be used only once---to setup the section commands to act as appendixes. Thereafter normal section commands are used. The heading for a section can be left empty. For example, 
-\begin{verbatim}
-\appendix
-\section{}
-\end{verbatim}
-will produce an appendix heading that says ``APPENDIX A'' and
-\begin{verbatim}
-\appendix
-\section{Background}
-\end{verbatim}
-will produce an appendix heading that says ``APPENDIX A: BACKGROUND'' (note that the colon is set automatically).
-
-If there is only one appendix, then the letter ``A'' should not appear. This is suppressed by using the star version of the appendix command (\verb+\appendix*+ in the place of \verb+\appendix+).
-
-\section{A little more on appendixes}
-
-Observe that this appendix was started by using
-\begin{verbatim}
-\section{A little more on appendixes}
-\end{verbatim}
-
-Note the equation number in an appendix:
-\begin{equation}
-  E^2=p^2c^2 + m^2c^4.
-\end{equation}
-
-\subsection{\label{app:subsec}A subsection in an appendix}
-
-You can use a subsection or subsubsection in an appendix. Note the numbering: we are now in Appendix~\ref{app:subsec}.
-
-Note the equation numbers in this appendix, produced with the subequations environment:
-\begin{subequations}
-\begin{align}
-  E^2 &= m^2c^4 \quad \text{(rest)}, \label{appa} \\
-  E^2 &= m^2c^4 + p^2 c^2 \quad \text{(relativistic)}, \label{appb} \\
-  E^2 &\approx p^2 c^2 \quad \text{(ultrarelativistic)} \label{appc}
-\end{align}
-\end{subequations}
-They turn out to be Eqs.~(\ref{appa}), (\ref{appb}), and (\ref{appc}).
-
-\end{document}
diff --git a/src/Makefile b/src/Makefile
deleted file mode 100644
index eb24188..0000000
--- a/src/Makefile
+++ /dev/null
@@ -1,22 +0,0 @@
-FC      = gfortran
-FFLAGS  = -Wall -Wextra -march=native -O3 -fimplicit-none 
-#FFLAGS += -pedantic -fbounds-check -fmax-errors=1 -g
-FFLAGS += $(shell pkg-config --cflags plplotd-f95)
-LDFLAGS = $(shell pkg-config --libs plplotd-f95)
-LIBS    =
-
-COMPILE = $(FC) $(FFLAGS)
-LINK = $(FC) $(LDFLAGS)
-
-TARGET = ising.exe			# Name of final executable to produce
-OBJS = plot.o ising.o	# List of object dependencies
-
-$(TARGET): $(OBJS)
-	$(LINK) -o $@ $^ $(LIBS)
-
-%.o:%.f95
-	$(COMPILE) -c $<
-
-.PHONY:clean
-clean:
-	rm *.exe *.o *.mod
diff --git a/src/R_code/mag_critical.R b/src/R_code/mag_critical.R
deleted file mode 100644
index 4a4878b..0000000
--- a/src/R_code/mag_critical.R
+++ /dev/null
@@ -1,51 +0,0 @@
-wolff <- read.table("temperature_sweep.dat",header=T)
-attach(wolff)
-
-steps = 14
-
-mag.crit <- which.max(Susceptibility)
-
-x1 = log(abs(Temperature - Temperature[mag.crit]))
-y1 = log(Magnetization - Magnetization[mag.crit])
-
-plot(x1,y1,type="n",
-    xlim=c(-3.5,-1),
-    ylim=c(-1,0),
-    xlab = "Log(Temperature - Tc)",
-    ylab ="Log(Magnetization)",
-    main = expression(paste("Magnetic Exponent ",beta,sep=" ")))
-
-grid()
-
-bound = c(120,120)
-tc.guesses = Temperature[seq(mag.crit - bound[1], mag.crit + bound[2],
-    length.out = steps)]
-
-mag.guesses= Magnetization[seq(mag.crit - bound[1], mag.crit + bound[2],
-    length.out = steps)]
-
-colors = rainbow(steps)
-colors[7] = "black"
-
-for(i in 1:steps ){
-    d1 = log(-Temperature + tc.guesses[i])
-    d2 = log(Magnetization)
-
-    lines(d1,d2,col = colors[i])
-
-    fitrange = which(-3.5 < d1 & d1 < -2)
-
-    fit = lm(d2[fitrange]~d1[fitrange])
-    print(fit)
-
-    if (i == 7) abline(fit, col = colors[i],lty = 4)
-    print(paste(i,tc.guesses[i]))
-}
-
-
-
-lines(log(-Temperature + tc.guesses[7]),d2)
-fit <- lm(d2[fitrange] ~ log(-Temperature + tc.guesses[7])[fitrange])
-abline(fit,lty = 4)
-
-detach(wolff)
diff --git a/src/R_code/make_plots.R b/src/R_code/make_plots.R
deleted file mode 100644
index 4780b30..0000000
--- a/src/R_code/make_plots.R
+++ /dev/null
@@ -1,39 +0,0 @@
-wolff <- read.table("temperature_sweep.dat",header=T)
-attach(wolff) #Attaches wolff to R search path; wolff$Value becomes Value
-
-
-sus.max = which.max(Susceptibility)
-cv.max = which.max(Cv)
-print(paste("Susceptibility max occurs at:",Temperature[sus.max]) ,sep = " ")
-print(paste("Cv max occurs at:",Temperature[cv.max]) ,sep = " ")
-
-par(mfrow = c(2,2), oma = c(0,0,2,0), mar = c(5,5,2,2) - .3)
-#opens graphics device with 2x2 matrix for plots, 2 line outer margin
-
-plot(Temperature,Magnetization,type="l",
-    xlab="Temperature",ylab="<Magnetization>")
-abline(v = Temperature[sus.max],col = "red",lty = 3)
-abline(v = Temperature[cv.max],col = "green",lty = 3)
-
-plot(Temperature,Susceptibility,type="l",
-    xlab="Temperature",ylab="Susceptibility")
-abline(v = Temperature[sus.max],col = "red",lty = 3)
-abline(v = Temperature[cv.max],col = "green",lty = 3)
-
-plot(Temperature,Energy,type="l",
-    xlab="Temperature",ylab="<Energy>")
-abline(v = Temperature[sus.max],col = "red",lty = 3)
-abline(v = Temperature[cv.max],col = "green",lty = 3)
-
-plot(Temperature,Cv,type="l",
-    xlab="Temperature",ylab="<Heat Capacity>")
-abline(v = Temperature[sus.max],col = "red",lty = 3)
-abline(v = Temperature[cv.max],col = "green",lty = 3)
-
-
-mtext("Wolff Algorithm", line = .5, outer = T)
-
-
-detach(wolff)
-
-#dev.copy2eps(file="system_plots.eps")
diff --git a/src/R_code/sus_critical.R b/src/R_code/sus_critical.R
deleted file mode 100644
index bccbbee..0000000
--- a/src/R_code/sus_critical.R
+++ /dev/null
@@ -1,33 +0,0 @@
-wolff <- read.table("temperature_sweep.dat",header=T)
-attach(wolff)
-
-steps = 14
-
-sus.crit <- which.max(Susceptibility)
-
-x1 = log(Temperature - Temperature[sus.crit])
-y1 = log(Susceptibility)
-
-plot(x1,y1,type="n",
-    ylim = c(-20,-12),
-    xlab = "Log(Temperature - Tc)",
-    ylab = "Log(Susceptibility)",
-    main = expression(paste("Susceptibility Exponent ",gamma," (tail)")))
-
-grid()
-
-colors = rainbow(steps)
-tmps = seq(2, 3, length.out = steps)
-for( i in steps:1 ){
-    lines(log(Temperature - tmps[i]),log(Susceptibility),col=colors[i])
-    }
-#plots warm colors first, warm colors represent "guesses" < t_c
-lines(x1,y1)
-
-pow.rng <- which(-3 < x1 & x1 < -1)
-
-sus.fit <- lm(y1[pow.rng] ~ x1[pow.rng])
-
-abline(sus.fit,lty = 4)
-
-detach(wolff)
diff --git a/src/R_code/visualize_lattice.R b/src/R_code/visualize_lattice.R
deleted file mode 100644
index 53c2f73..0000000
--- a/src/R_code/visualize_lattice.R
+++ /dev/null
@@ -1,15 +0,0 @@
-spin.file <- "array.dat"
-#change to path of spin output file
-
-spin.colors <- c("purple","red","green") #downspin, upspin
-spin.matrix <- as.matrix(read.table(spin.file))
-
-par(mar=c(1,1,3,1)+.1)
-
-image(spin.matrix,col=spin.colors,axes=FALSE)
-title(main="Wolff Cluster Algorithm")
-
-#axis(1,at=seq(0,1,by=(1/(dim(spin.matrix)[1]-1))*1),labels=F)
-#axis(2,at=seq(0,1,by=(1/(dim(spin.matrix)[2]-1))*1),labels=F)
-box()
-dev.copy2pdf(file="vis_lattice.pdf") #uncomment to generate output file
diff --git a/src/calcerror.f90 b/src/calcerror.f90
new file mode 100644
index 0000000..1301d0c
--- /dev/null
+++ b/src/calcerror.f90
@@ -0,0 +1,43 @@
+   ! calculate the mean and error for each temperature/chain length set
+   
+   program error
+   implicit none
+   integer i,nruns,io
+   real*8 logn,total,mean,variance,err
+   real*8, allocatable :: logrg(:)
+   
+   ! read in file until eof to get number of data points
+   nruns = 0
+   do
+      read(5,*,iostat=io) 
+      if (io/=0) exit
+      nruns = nruns + 1
+   enddo 
+   
+   rewind(5)
+   
+   allocate(logrg(nruns))
+   
+   ! read file to get logR and logN
+   total=0
+   do i=1,nruns
+      read(5,*) logrg(i),logn
+      total = total + logrg(i)
+   enddo  
+   
+   ! calculate mean
+   mean = total/real(nruns)
+   
+   ! calculate variance
+   total=0
+   do i=1,nruns
+      variance = (mean - logrg(i))**2
+      total = total + variance
+   enddo
+   
+   ! calculate the error
+   err = sqrt(total/real(nruns))
+   
+   print *, mean, err
+   
+   end program error
diff --git a/src/chainmc.f90 b/src/chainmc.f90
new file mode 100644
index 0000000..7f87e76
--- /dev/null
+++ b/src/chainmc.f90
@@ -0,0 +1,73 @@
+! main program for advanced monte carlo
+! creates a chain polymer
+! in 3 dimensions
+
+program chainmc
+   implicit none
+   integer N,i,nangles,j,k
+   real*8 T,pos,Vtot
+   real*8, allocatable, dimension(:,:) :: r
+   real*8, allocatable, dimension(:,:) :: V
+   real*8, allocatable, dimension(:,:) :: temppos
+   
+   ! seed random number generator
+   !call random_seed()
+   call init_random_seed()
+   
+   ! number atoms in the chain
+   N=10000
+   
+   ! Temperature
+   T = 2.d0
+   
+   ! allocate memory for the positions in the chain
+   allocate (r(N,3))
+   ! allocate memory for the potentials between pairs
+   allocate (V(N,N))
+   ! initilize to zero
+   r=0.d0
+   V=0.d0
+   
+   ! set positions of the first two atoms (0,0,0) and (1,0,0)
+   r(1,1)=0
+   r(1,2)=0
+   r(1,3)=0
+   r(2,1)=1
+   r(2,2)=0
+   r(2,3)=0
+   
+   ! add the rest of the atoms
+   ! define a number of sample points
+   nangles = 90
+   do i=3,N
+      ! calculate the energy of the previous configuration
+      call calcpot(i,V,Vtot,N,r)
+      
+      ! pick new direction, this will fill positions (r)
+      ! first pick possible directions
+      allocate (temppos(nangles,3))
+      call pickdir(i,r,N,temppos,nangles)
+      
+      ! then calculate the energy of each configuration
+      ! use metropolis algorithm to decide "best" configuration
+      call pickconfig(i,r,N,T,temppos,Vtot,nangles)
+      deallocate(temppos)
+   enddo
+   
+   ! check the lengths of all adjacent segements
+   do i=1,N-1
+      pos = 0
+      pos = (r(i,1)-r(i+1,1))**2 + (r(i,2)-r(i+1,2))**2 + (r(i,3)-r(i+1,3))**2
+      pos = sqrt(pos)
+   enddo 
+   
+   ! calculate the size of the system
+   call calcsize(r,N)
+   
+   ! prepare to run program many times to get many chains
+   ! for some sort of error
+   deallocate(r)
+   deallocate(V)
+
+
+end program chainmc
diff --git a/src/direction.f90 b/src/direction.f90
new file mode 100644
index 0000000..489e74b
--- /dev/null
+++ b/src/direction.f90
@@ -0,0 +1,32 @@
+! choose a direction for the next atom
+! randomly choose theta and phi
+
+subroutine pickdir(i,r,N,pos,nang)
+   implicit none
+   integer, intent(in) :: i,N,nang
+   real*8, intent(inout) :: r(N,3),pos(nang,3)
+   real*8 theta,phi,pi
+   integer j
+   
+   ! define pi
+   pi = 4.d0 * datan(1.d0)
+   
+   !random number
+   do j=1,nang
+      call random_number(theta)
+      call random_number(phi)
+      !print *, theta,phi
+   
+      ! sample theta from cos(theta)
+      ! theta between 0 and pi
+      ! phi between 0 and 2pi
+      theta = 2*theta - 1 ! between 1 and -1
+      theta = acos(theta) ! becomes an angle
+      phi = phi*2.d0*pi
+      ! define array of new possible positions
+      pos(j,1) = r(i-1,1)+sin(theta)*cos(phi)
+      pos(j,2) = r(i-1,2)+sin(theta)*sin(phi)
+      pos(j,3) = r(i-1,3)+cos(theta)
+    enddo
+   
+end subroutine pickdir
diff --git a/src/energycalc.f90 b/src/energycalc.f90
new file mode 100644
index 0000000..239f085
--- /dev/null
+++ b/src/energycalc.f90
@@ -0,0 +1,47 @@
+   subroutine pickconfig(i,r,N,T,pos,Vtot,nang)
+   implicit none
+   integer, intent(in) :: i,N,nang
+   real*8, intent(in) :: T,pos(nang,3),Vtot
+   real*8, intent(inout) :: r(N,3)
+   real*8 Z,Vtemp(20),dtemp,rand,prob
+   integer j,k,count,pickk,l
+   
+   ! calculate potential for each possible configuration
+   ! for the ith atom
+   do k=1,nang
+      Vtemp(k) = 0
+      do j=1,i-1
+         dtemp = 0
+         dtemp = (r(j,1)-pos(k,1))**2 + (r(j,2)-pos(k,2))**2 + (r(j,3)-pos(k,3))**2
+         Vtemp(k) = Vtemp(k) + 4*((1.d0/dtemp)**6 - (1.d0/dtemp)**3)
+      enddo 
+      Vtemp(k) = Vtemp(k) !+ Vtot
+   enddo 
+   
+   ! calculate partition function
+   Z = 0
+   do l=1,nang
+      Z = Z + exp(-Vtemp(l)/T)
+   enddo 
+   
+   ! generate random number 
+   call random_number(rand)
+   ! pick new position based on random number
+   count = 1
+   prob = 0
+   do while (count < nang+1)
+      prob = prob + exp(-Vtemp(count)/T)/Z
+      if (rand <= prob) then
+         pickk = count
+	 count = 1000
+      else
+         count = count + 1
+      end if 
+   enddo 
+   
+   ! define the new position
+   r(i,1) = pos(pickk,1)
+   r(i,2) = pos(pickk,2)
+   r(i,3) = pos(pickk,3)   
+   
+   end subroutine pickconfig
diff --git a/src/makefile b/src/makefile
new file mode 100644
index 0000000..9363476
--- /dev/null
+++ b/src/makefile
@@ -0,0 +1,19 @@
+FC      = ifort
+FFLAGS = -cm -w -02
+#FFLAGS  = -Wall -Wextra -O3 -fimplicit-none -march=native 
+#FFLAGS += -pedantic -fbounds-check -fmax-errors=1 -g
+#FFLAGS += $(shell pkg-config --cflags plplotd-f95)
+#LDFLAGS = $(shell pkg-config --libs plplotd-f95)
+LIBS    =
+
+COMPILE = $(FC) $(FFLAGS)
+LINK = $(FC) $(LDFLAGS)
+
+TARGET = createchain       # Name of final executable to produce
+OBJS = direction.o sizecalc.o energycalc.o potcalc.o chainmc.o seedrn.o # List of object dependencies
+
+$(TARGET): $(OBJS)
+	$(LINK) -o $@ $^ $(LIBS)
+
+%.o:%.f90
+	$(COMPILE) -c $<
diff --git a/src/plot.f95 b/src/plot.f95
deleted file mode 100644
index 1fa8c67..0000000
--- a/src/plot.f95
+++ /dev/null
@@ -1,64 +0,0 @@
-module plot
-  use plplot
-  implicit none
-  private
-
-  public plot_init, plot_close, plot_lattice
-
-contains
-
-  subroutine plot_init(latticeSize)
-    integer, intent(in) :: latticeSize
-    call plsdev("xcairo")
-    call plinit()
-
-    !You can find default colors at
-    !http://plplot.sourceforge.net/docbook-manual/plplot-html-5.9.9/plcol0.html
-
-    !call plscol0(0, 255, 255, 255)  ! white
-    !call plscol0(1, 255, 0, 0)      ! red
-    !call plscol0(2, 255, 77, 0)     ! orange
-    !call plscol0(3, 255, 255, 0)    ! yellow
-    !call plscol0(4, 0, 255, 0)      ! green
-    !call plscol0(5, 0, 0, 255)      ! blue
-    !call plscol0(6, 0, 255, 255)    ! cyan
-    !call plscol0(7, 255, 0, 255)    ! magenta
-    !call plscol0(8, 128, 128, 128)  ! gray
-    !call plscol0(9, 0, 0, 0)        ! black
-
-    call plenv(0d0, latticeSize + 1d0, 0d0, latticeSize + 1d0, 0, 0)
-  end subroutine plot_init
-
-  subroutine plot_close()
-    call plspause(.false.)
-    call plend()
-  end subroutine plot_close
-
-  subroutine plot_lattice(lattice)
-    ! Draw a rectangular grid of up- and down-arrows corresponding to Ising
-    ! states. Because this redraws the entire screen, it is *very* slow for
-    ! Metropolis models - instead, consider a routine that only redraws the
-    ! sites that change.
-    integer, intent(in) :: lattice(:,:)
-    integer :: i, j
-    real(8) :: x, y
-    !Assumes 1 corresponds to up-spin, -1 to down-spin.
-
-    call plclear()
-    do i = 1, size(lattice, 1)
-      do j = 1, size(lattice, 2)
-        x = i; y = j
-        if (lattice(i, j) .eq. 1) then
-          call plcol0(1)            !default red
-          call plpoin([x], [y], 30) !30 denotes up-arrow glyph
-        else ! Comment out the else clauses to speed drawing
-          call plcol0(11)           !default cyan
-          call plpoin([x], [y], 31) !31 denotes down-arrow glyph
-        end if
-      end do
-    end do
-
-    call plflush()
-  end subroutine plot_lattice
-
-end module plot
diff --git a/src/potcalc.f90 b/src/potcalc.f90
new file mode 100644
index 0000000..8cc8769
--- /dev/null
+++ b/src/potcalc.f90
@@ -0,0 +1,37 @@
+   subroutine calcpot(i,V,Vtot,N,r)
+   implicit none
+   integer, intent(in) :: i,N
+   real*8, intent(inout) :: Vtot
+   real*8, intent(inout) :: V(N,N)
+   real*8, intent(in) :: r(N,3)
+   integer j,k
+   real*8 distsq
+   
+   ! calculate pair-wise potentials
+   if (i==3) then
+      ! calculate the potential between atoms 1 and 2
+      distsq = (r(1,1)-r(2,1))**2 + (r(1,2)-r(2,2))**2 + (r(1,3)-r(2,3))**2
+      V(1,2) = 4*((1.d0/distsq)**6 - (1.d0/distsq)**3)
+      V(2,1) = V(1,2)
+      !print *, distsq,V(1,2),"1,2"
+      ! should be zero unless we change sigma
+   else 
+      ! calculate the potential between the newest atom and the previous ones
+      do j=1,i-2
+         distsq = (r(j,1)-r(i-1,1))**2 + (r(j,2)-r(i-1,2))**2 + (r(j,3)-r(i-1,3))**2
+	 V(i-1,j) = 4*((1.d0/distsq)**6 - (1.d0/distsq)**3)
+	 V(j,i-1) = V(i-1,j)
+      enddo 
+   end if 
+   
+   ! calculate the total potential energy of the system
+   ! this is for the old configuration - before the ith atom is added
+   if (i==3) then
+      Vtot = V(1,2)
+   else
+      do j=1,i-2
+         Vtot = Vtot + V(i-1,j)
+      enddo 
+   end if 
+   
+   end subroutine calcpot
diff --git a/src/seedrn.f90 b/src/seedrn.f90
new file mode 100644
index 0000000..53018d1
--- /dev/null
+++ b/src/seedrn.f90
@@ -0,0 +1,18 @@
+   ! seed the random number generator
+   ! from:  fortranwiki.org/fortran/show/random_seed
+   
+   subroutine init_random_seed()
+      integer i,n,clock
+      integer, dimension(:), allocatable :: seed
+      
+      call random_seed(size = n)
+      allocate(seed(n))
+      
+      call system_clock(count=clock)
+      
+      seed = clock + 37 * (/ (i-1, i=1, n) /)
+      call random_seed(put=seed)
+      
+      deallocate(seed)
+      
+   end subroutine
diff --git a/src/sizecalc.f90 b/src/sizecalc.f90
new file mode 100644
index 0000000..7e90032
--- /dev/null
+++ b/src/sizecalc.f90
@@ -0,0 +1,34 @@
+   subroutine calcsize(r,N)
+   implicit none
+   integer, intent(in) :: N
+   real*8, intent(in) :: r(N,3)
+   real*8 xcm,ycm,zcm,rcm,rg,ri
+   integer j
+   
+   ! calculate the radius of gyration
+   ! R_g^2 = (1/N) \sum <(r_i - R_cm)^2>
+   
+   ! calculate the center of the distribution, R_cm
+   xcm = 0
+   ycm = 0
+   zcm = 0
+   do j=1,N
+      xcm = xcm + r(j,1)
+      ycm = ycm + r(j,2)
+      zcm = zcm + r(j,3)
+   enddo
+   
+   rcm = sqrt(xcm**2 + ycm**2 + zcm**2)/N
+   
+   ! R_g calculation
+   rg = 0
+   do j=1,N
+      ri = sqrt(r(j,1)**2 + r(j,2)**2 + r(j,3)**2)
+      rg = rg + (ri - rcm)**2
+   enddo
+   rg = sqrt(rg/N)
+   open(unit=20,file="RvsNT2010000.txt",access="append")
+   write(20,*) log(rg), log(real(N))
+   close(20)
+      
+   end subroutine calcsize