1D ising model - Helix-coil transistion

[Matt atkinson]

Homework Statement

A simple model of a polymer undergoing a helix-coil transition is to describe the polymer
in terms of N equal length segments, each of which can be in either a coil or
a helix state. A more realistic model also takes into account the energy cost associated
with a boundary between a helix section and a coil section which leads to
cooperativity. This is analogous to the Ising model.
Write a simple Monte Carlo simulation to model a polymer helix-coil transition. Each
of the N segments can be in a coil or helix state and can switch state with a probability
depending on the energy cost of switching state. You may assume that the change in
free energy of a segment switching from coil to helix in units of $k_BT$ is
$$\frac{∆Fh}
{k_BT}
= a(T − T_c)$$
where a is a constant, $T$ is the temperature and $T_c$ is the critical temperature for the
transition. Your simulation should compare generated pseudo random numbers with
the relevant probability to decide whether a segment should switch state or not. You
should start your simulation with an initial state for your polymer and then run over
many steps of trial segment state switches. You will need enough steps so that the
polymer has relaxed to equilibrium. Repeat your simulation experiment many times to
get good averages. Run your simulation at different temperatures to find the fraction
of helices as a function of temperature.
It is recommended to do the case with no cooperativity first. To introduce cooperativity,
i.e. a non-zero energy cost of junctions between segments of different states,
you will need to modify your code to check the state of a segment’s neighbours to
calculate the probability of switching state.

Homework Equations

The Attempt at a Solution

Just confused as to if I am approaching the problem correctly;
In c++ I set up at array of size Nx1 where N is the number of monomers, or elements.
I then have relevant code to flip the particle form a "helix" to "coil" state if the probability given by $e^{\Delta E}$ where for the non-cooperativity would be just $e^{-a(T-T_c)}$ right? this part I understand. Also if $\Delta E<0$ or the probability given above is less than a generated random number then the state change will occur because it minizes entropy.
In my simulation I have used the value T as $T-T_c$ because no specific polymer was given and i range the value of T from -10 to 10. The constant "a" is apparently a scale factor to fit with experimental data.
Here's my problem, i have no idea how to add cooperativity maybe some local energy of +1 or -1 if the adjacent elements are in the same or opposite state?
To get the fraction of helix's I looped a huge number of times over the flip probability and then add up the number of elements with the value I set for the for the "helix" state.

I can attach my code if that would help?

 

[Päällikkö]

First of all, the model you describe with no cooperativity is like the Ising model in an external field with no interactions between the sites. Adding cooperativity, then, is what adding the nearest-neighbor interactions is to the Ising model.

So basically you need to modify your energy calculation: the interaction energy between coil and helix is some E1, the energy between coil and coil or helix and helix some other number, say 0 (remember Flory-Huggins and how we get the chi parameter, so interaction energies between things of the same kind can effectively be put to 0?). Anyway, as you flip a an entity on the array, you'll need to check if this creates new helix-coil bonds (increasing the energy) or destroys them (decreasing the energy).

 

[hibphys]

hey how do i post a question ?

 

[Matt atkinson]

First of all, the model you describe with no cooperativity is like the Ising model in an external field with no interactions between the sites. Adding cooperativity, then, is what adding the nearest-neighbor interactions is to the Ising model.

So basically you need to modify your energy calculation: the interaction energy between coil and helix is some E1, the energy between coil and coil or helix and helix some other number, say 0 (remember Flory-Huggins and how we get the chi parameter, so interaction energies between things of the same kind can effectively be put to 0?). Anyway, as you flip a an entity on the array, you'll need to check if this creates new helix-coil bonds (increasing the energy) or destroys them (decreasing the energy).

Ah thank you that cleared a lot up, i think I've pretty much got the code working now. Just a question though so basically if before the flip of an entity on the array the state is coil - coil - coil i would be increasing the energy when it is flipped to coil - helix - coil, but if the opposite occurs it would be a decrease in energy? so for the first transition would be a +1 to the local energy but the second would be a -1?

and i could have some code like
if(particles[right_column].spinvalue()==1 && particles[current_column].spinvalue()==1) L_energy+=0;
if(particles[right_column].spinvalue()==1 && particles[current_column].spinvalue()==-1) L_energy+=1;
if(particles[right_column].spinvalue()==-1 && particles[current_column].spinvalue()==1) L_energy-=1;
if(particles[right_column].spinvalue()==-1 && particles[current_column].spinvalue()==-1) L_energy+=0;

for the local energy of the interaction with the current entity and the one right of it?
And finally was it correct to use $T-T_c$ as a separate variable say $B=T-T_c$ and vary B as I don't have a specific polymer to study?

I haven't really had much experience with biology/biological physics so sorry if these are simple questions.

 

[Päällikkö]

Note that with the flip you describe, you are creating two H-C bonds, so the energy increase should probably be +2 (and yes, -2 when you flip it back). An easier way to check whether the spins align is to do what the Ising Hamiltonian does: play with signs, i.e. s*s[i+1] == 1 if the spins are the same (or == -1 if the spins are different).

As for the new variable, I think it was fine for the first part, but now that you have cooperativity, you need to think about how these scale with temperature (in that the Boltzmann term has E/kT, where T is the absolute temperature, not a relative one). You could just include the 1/kT in the constant (as in, the energy of the bond between H-C is b/kT), like is done in the Flory-Huggins chi parameter, but you then need to remember the meaning of the parameter (that it is in fact highly temperature-dependent).

As for your implementation in code, and this is a matter of taste, but I would probably like a method in the Particle class that returns an array of pointers to its neighbours. You could then very easily just sum over all neighbours rather than worrying about right_column, left_column (up and down columns if you decide to generalize the code into multiple dimensions) in the computational part.

 

[Matt atkinson]

Okay thanks got it to work and almost give the expected output, just got two more questions if that's fine:
-you said +2 and -2 for flips but if i have my code like below, that takes into account the +2 and -2 right?

if(particles[right_column].spinvalue()==1 && particles[current_column].spinvalue()==1) L_energy+=0;
if(particles[right_column].spinvalue()==1 && particles[current_column].spinvalue()==-1) L_energy+=1;
if(particles[right_column].spinvalue()==-1 && particles[current_column].spinvalue()==1) L_energy-=1;
if(particles[right_column].spinvalue()==-1 && particles[current_column].spinvalue()==-1) L_energy+=0;

if(particles[left_column].spinvalue()==1 && particles[current_column].spinvalue()==1) L_energy+=0;
if(particles[left_column].spinvalue()==1 && particles[current_column].spinvalue()==-1) L_energy+=1;
if(particles[left_column].spinvalue()==-1 && particles[current_column].spinvalue()==1) L_energy-=1;
if(particles[left_column].spinvalue()==-1 && particles[current_column].spinvalue()==-1) L_energy+=0;

-also I am pretty new to montecarlo simulations and for the error calculation is it expected to be tiny for example, we were given a temperature of $T_c=427 ^o K$ for a polymer of $N=10$ so if i run the program say 10000 times for each temperature and only take into account the last 5000 (after the system got into a equilibrium state) and i was getting errors of $\approx 10^{-9}$. using this formula:
$$ \sigma=\sqrt{\frac{<a^2>-<a>^2}{5000}} $$

 

[Päällikkö]

I don't quite understand your code example. Has the current_column already been flipped or not? I'm concerned by these two lines:
if(particles[right_column].spinvalue()==1 && particles[current_column].spinvalue()==-1) L_energy+=1;
if(particles[right_column].spinvalue()==-1 && particles[current_column].spinvalue()==1) L_energy-=1;

If the right_column has the same spin as the current_column, should the energy not always decrease (or, increase, if you have not yet flipped current_column)?

Suppose that you in your code pick a current_column.
Now you flip the spin.
Then you compute the energy of this new system and compare it to the previous state (with no flip). This means that if particles[right_column].spinvalue() == particles[current_column].spinvalue(), then you have broken a H-C bond and formed a H-H or C-C bond instead. This is to say that the energy should decrease. If those two spins are not equal, your spin flip has generated a H-C bond where there used to be a C-C or H-H. This is to say that the energy of the conformation should increase. Same goes for the comparison with the left_column.
Now you do the Monte Carlo comparison (energy vs. a random number). If it fails, you flip the current_column back to its previous state, if it succeeds, you keep the new conformation.
Finally, you start again from the beginning.

 

[Matt atkinson]

Oh okay, my thinking was if the spins are different then there was a energy in between then when i flipped it it would decrease the energy, maybe i have my codes backwards, where the one where the spin=spin should be plus or minus and the different spins should be 0.
As for the calculation I use this code in my main body:

particle.choose();
Eb = particle.localEnergy();
particle.perturb();
Ef = particle.localEnergy();
particle.perturb();
if(a==1 && b==-1)DeltaE2=(T-Tc);
if(a==-1 && b==1) DeltaE2=(Tc-T);
DeltaE=DeltaE2+(Ef-Eb)/T;
r=(double)rand()/(double)RAND_MAX;
e=exp(-DeltaE);
if(r < e || DeltaE<0)
{
particle.perturb();
}

Which would calculate the local energy before and after a flip and them compare with random number as you said.
the code which I pasted the local energy calculation is part of the particle class and i run it before and after.
I appear to get results which agree with theory, but that means nothing if i have my physics wrong in calculating the local energy. I wonder if it should be like this instead;

if(particles[right_column].spinvalue()==1 && particles[current_column].spinvalue()==1) L_energy+=+1;
if(particles[right_column].spinvalue()==1 && particles[current_column].spinvalue()==-1) L_energy+=-;
if(particles[right_column].spinvalue()==-1 && particles[current_column].spinvalue()==1) L_energy-=0;
if(particles[right_column].spinvalue()==-1 && particles[current_column].spinvalue()==-1) L_energy-=1;

if(particles[left_column].spinvalue()==1 && particles[current_column].spinvalue()==1) L_energy+=1;
if(particles[left_column].spinvalue()==1 && particles[current_column].spinvalue()==-1) L_energy+=0;
if(particles[left_column].spinvalue()==-1 && particles[current_column].spinvalue()==1) L_energy-=0;
if(particles[left_column].spinvalue()==-1 && particles[current_column].spinvalue()==-1) L_energy-=1;

I attached a file of the results plotted with the theory that my group member calculated.

 

[Päällikkö]

What I don't quite understand is what the term DeltaE2 in your code is. I also see a typo in your equation ( L_energy+=-; ). Are you sure this is the part that you meant to paste here?

So I'm supposing the figure you derived comes from the system with no interactions, so that the Hamiltonian may be written \mathcal{H} = h\sum_{i=1}^N s_i. It is then trivial to show that the partition function \mathcal{Z} = (2\cosh \beta h)^N, so the portion of helices is \frac{1}{2}(1-\tanh\beta h) (the more general case with interactions between sites is a bit more complicated, but using the transfer matrix method, nothing that couldn't be handled).

Your theoretical result may be written in the form \frac{1}{2}(1-\tanh\frac{T-T_c}{2}), which is to say that I think you are off by a factor of \frac{1}{2\beta} somewhere in your calculations. Well, not quite: you probably did mean to divide by beta, but I still am unsure what's up with the 1/2.

If you want something to compare with, I wrote a quick and dirty version that implements your interface:

Class def'n

class Particle {
public:
        Particle(float T);
  void  choose();
  void  perturb();
  float localEnergy() const;
private:
  std::vector<int> neighbors() const;
private:
  std::uniform_int_distribution<int>    mDistI{0, N-1};
  std::mt19937                          mGen;
  int                                   mChosen {0};
  float                                 mT;
  std::array<int, N>                    mArray;
};

The implementation of localEnergy():

float Particle::localEnergy() const {
  auto ener = 0.f;
  ener += (mT-Tc)*mArray[mChosen];
  auto neighs = neighbors();
  for(const auto n : neighs)
    ener += -mArray[mChosen]*mArray[n];
  return ener*mT;
}

And the indices to neighboring sites (note that I assume periodic boundary conditions; the return type is std::vector rather than std::array in case you want to change this, or have larger neighborhoods or whatever, in case you were wondering)

std::vector<int> Particle::neighbors() const {
  return {(mChosen-1 < 0) ? N-1 : mChosen-1,  
          (mChosen+1 > N-1) ? 0 : mChosen+1};
}

Note that I divide by beta (multiply by the temperature) at the end of localEnergy(). This is to make the constant depend on temperature, as I'm guessing you've done in your code. My main loop then is virtually the same as yours, except the DeltaE2 term (and I'm using the STL functions from <random> rather than the C implementation, so you need to include those and compile with the -std=c++11 or similar flag if you wish to just copy paste what I wrote above).

 

[Matt atkinson]

Ah yes sorry it should be L_energy+=0; the "-" was a typo.
As for the DeltaE2 I just took it from the question for the free energy in units of $k_B T$ and my supervisor said that the flip back would give a negative energy, that is why one is $T-T_c$ and the other is $T_c-T$.
Im pretty sure that my group derived the equation with interactions between sites, but they did struggle with it so I'm not entirely sure.
Your code looks much better, I'm pretty new to c++ (not had much c experience either) so ill have to get my head around it and try it.

It appears that my code is wrong, but I must have gotten lucky to get a similar shape to the solution the theory gave.

Really appreciate the help, and I understand I'm new to coding (as you can see from how messy it is) and also biological physics (no experience what so ever).

 

[Päällikkö]

You already compute the flip energy in localEnergy(), and it should be negative/positive depending on whether the spin aligns with the external field. Suppose the external field points to (T-Tc), meaning it is positive when T>Tc. This means that -(T-Tc)*Si minimizes energy when the spin is to the same direction, i.e. positive (meaning that when T>Tc the positive spin, i.e. coil, is favored, and the other way around). In my last post I wrote the Hamiltonian without the minus sign in front, which would mean that the spins align opposite to the external field, I hope you won't find this terribly confusing (the only effect of the "wrong" sign is a different definition of spins: +1 = helix, -1 = coil, nothing else changes; This is why I had 1-tanh instead of 1+tanh in the final equation).

That wasn't very clear. What I wanted to say is that DeltaE2 should not enter the picture: The energy is computed by the localEnergy() method.

I am sure that the equation your group members derived does not have interactions in it (there is no term accounting for them). If you want to know what the correct equation looks like, here it is (not promising there are no typos):
\frac{1}{2}\left(1-\frac{\sinh(\beta h)}{\cosh(\beta h) + \sqrt{\sinh^2(\beta h)+\exp(-4\beta J)}}\left(1+\frac{\cosh(\beta h)}{\sqrt{\sinh^2(\beta h)+\exp(-4\beta J)}}\right) \right)

J here is the interaction between the sites, and h the external field, as before (so with the "wrong" sign); I've again assumed periodic boundaries, so especially for shorter polymers you will get results that slightly vary.

Anyway, if you have more questions, feel free to ask. The C++11 features I was using in my code include "auto", "braced initialization, nd "ranged for loops" (and from the new STL, <array> and <random>).

 

[Matt atkinson]

Okay so I think I understand, in my code I am basically accounting for the energy of the flips twice? once with DeltaE2 and the and again in the local energy?

Also if i wanted to try your code also just for a comparison point of view would I just have to write a main function similar to my current one?

 

[Päällikkö]

Yes, the role of DeltaE2 could be regarded as counting the energy a second time. Here's my main:

int main() {
  std::uniform_real_distribution<float> dist{0.f, 1.f};
  std::mt19937                          gen;

  std::vector<float> Ts;
  for(float T = 426.f; T < 428.f; T+=.025f)
    Ts.push_back(T);
  std::vector<float> res;

  for(auto T : Ts) {
    Particle particle(T);
    std::vector<float> medres;
    medres.reserve(10000);
    for(int i = 0; i < 100000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy(); 
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%10) && i > 1000)
        medres.push_back(particle.mean());
    }
    float xx = 0.f;
    for(const auto x : medres)
      xx += x;
    res.push_back(xx/medres.size());
  }
  for(unsigned int i = 0; i < Ts.size(); ++i)
    std::cout << Ts[i] << " " << (1+res[i])*.5 << std::endl;
  std::cout << std::endl;
}

So this is pretty much the same as yours (I used your code as a model, so I could easily replace parts of my code with yours or the other way around, if need be).
The vectors medres and res are there just to compute some statistics (and note that the Particle class has a mean() method that I did not list earlier).

 

[Matt atkinson]

Okay so I essential put these into separate files in a project, do i need any preprocessor directives specially?, and the mean method what would that do? just find the mean fraction of helix's?
Again sorry if these are really basic questions I've had to look up a lot of the code you send because i haven't gone really far into c++ as of yet.

 

[Päällikkö]

The method mean() would return the average spin (not quite the same thing as the proportion of helices, but similar). You can make it work by returning the proportion of helices, but need to modify one line at the std::cout line where stuff are printed out.
I have everything in the same file (I don't use headers for small projects like this; well under 100 lines), and just compile with g++ -std=c++11 -O3 file.cpp
The -std flag is necessary (as is a compiler that supports the C++11 standard). The -O3 makes the compiler optimize the code for speed (and it gives a very, very noticeable performance boost), but you can drop it if you so desire.

Here are my includes:

#include <cmath>
#include <array>
#include <vector>
#include <random>
#include <iostream>

I don't particularly like the code's structure, and especially not how I have been rather inconsistent in my typing (for example I sometimes use auto, sometimes don't), so don't use this as an example for your future projects.

 

[Matt atkinson]

I won't thank you, I'm going to take the time to learn C++ in more depth, and get to grips with it after this project is over.
As for a compiler or where you said "g++ -std=c++11 -O3 file.cpp" I use Dev c++ at the moment and have no idea if it had C++11 standard and I am not sure what to do with that line "g++ -std=c++11 -O3 file.cpp".

 

[Päällikkö]

Which Dev-C++ are you using? If it is the one available from http://orwelldevcpp.blogspot.com/, you have GCC 4.8.1 which should support most of the C++11 standard (use -std=c++0x if -std=c++11 is not recognized). As for how to add compiler flags (-std=c++11) to the particular program, I wouldn't know, I use Linux and prefer to compile everything from the command line. I'm sure that Dev-C++ does have a compiler flags option somewhere (or if it creates Makefiles, then you could just add it there, to CFLAGS). You would want to use -O3 or -O2 for all your projects, anyway (except when debugging).

When using an IDE, I prefer C::B (not that there is anything wrong with Dev-C++, I remember testing it some 7 years ago and it was fine), available at http://www.codeblocks.org/downloads/26. If you do download this, take the one with GCC 4.8.1 bundled with it (or alternatively, take the one without GCC and install MinGW/TDM-GCC separately).

 

[Matt atkinson]

Okay yes, i found compiler option that make's the compiler use C++11 but i get the errors that " [Error] 'vector' in namespace 'std' does not name a type"

 

[Päällikkö]

Did you forget to
#include <vector>
?
I believe you'll need to include this both in the header and in the file containing main (supposing that you do have separate files).

 

[Matt atkinson]

Ah yea sorry, completely missed all my includes...
May I ask what your mean function was?

 

[Päällikkö]

This is already implicitly in my main() function (computing to res from medres), but here:

float mean() const {
  float xx = 0.f;
  for(const auto x: mArray)
    xx += x;
  return xx/mArray.size();
}

Might as well include the remaining functions:

void Particle::choose() {
  mChosen = mDistI(mGen);
}

void Particle::perturb() {
  mArray[mChosen] *= -1;
}

 

[Matt atkinson]

I feel like there's some stuff missing, and also get quite a few errors like "[Error] 'mArray' was not declared in this scope"
Again I apologise for not being very good at coding.

#include <cmath>
#include <array>
#include <vector>
#include <random>
#include <iostream>

class Particle {
public:
        Particle(float T);
  void  choose();
  void  perturb();
  float localEnergy() const;
private:
  std::vector<int> neighbors() const;
private:
  std::uniform_int_distribution<int>    mDistI{0, N-1};
  std::mt19937                          mGen;
  int                                   mChosen {0};
  float                                 mT;
  std::array<int, N>                    mArray;
};

int main() {
  std::uniform_real_distribution<float> dist{0.f, 1.f};
  std::mt19937                          gen;

  std::vector<float> Ts;
  for(float T = 426.f; T < 428.f; T+=.025f)
    Ts.push_back(T);
  std::vector<float> res;

  for(auto T : Ts) {
    Particle particle(T);
    std::vector<float> medres;
    medres.reserve(10000);
    for(int i = 0; i < 100000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy();
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%10) && i > 1000)
        medres.push_back(particle.mean());
    }
    float xx = 0.f;
    for(const auto x : medres)
      xx += x;
    res.push_back(xx/medres.size());
  }
  for(unsigned int i = 0; i < Ts.size(); ++i)
    std::cout << Ts[i] << " " << (1+res[i])*.5 << std::endl;
  std::cout << std::endl;
}

float Particle::localEnergy() const {
  auto ener = 0.f;
  ener += (mT-Tc)*mArray[mChosen];
  auto neighs = neighbors();
  for(const auto n : neighs)
    ener += -mArray[mChosen]*mArray[n];
  return ener*mT;
}
float mean() const {
  float xx = 0.f;
  for(const auto x: mArray)
    xx += x;
  return xx/mArray.size();
}

std::vector<int> Particle::neighbors() const {
  return {(mChosen-1 < 0) ? N-1 : mChosen-1, 
          (mChosen+1 > N-1) ? 0 : mChosen+1};
}
void Particle::choose() {
  mChosen = mDistI(mGen);
}

void Particle::perturb() {
  mArray[mChosen] *= -1;
}

 

[Päällikkö]

Please try to debug the programs a bit before asking. You have not defined N or Tc anywhere in the file (I use global constexpr). The former is needed for the array to be allocated. The mean() method ought to be a part of the Particle class (and as it isn't, it doesn't know what mArray is, creating errors).

 

[Matt atkinson]

Sorry about that I've fixed the bugs' defining Tc and N, but now get an error;
"C:\Users\Matt\AppData\Local\Temp\ccI6g56X.o Helixfraction.cpp:(.text+0x130): undefined reference to `Particle::Particle(float)' "

 

[Päällikkö]

You are missing the definition of the constructor (i.e. Particle::particle(float T)). I picked +1, -1 randomly, but you can just assign +1 spin to everything in mArray (this works as there should be no danger of metastable states in the transition you are looking at). Remember to assign mT from the T input.

Here's my constructor (I decided to show it here as it does highlight a couple of C++ features not used elsewhere in the code: the initialization of mT and how to use ranged loops for assigning back to the looped variable)

Particle::Particle(float T) 
: mT{T} {
  std::uniform_real_distribution<float> distR{0.f, 1.f};
  for(auto& x: mArray) 
    x = (distR(mGen) > .5f) ? 1 : -1; 
}

 

[Matt atkinson]

Okay, I have to be honest I have no idea, I tried doing;

int spin;

to the class.

Particle::Particle(float T){
    spin=-1;
    }

but I know this is wrong cause it didnt work.

 

[Päällikkö]

I edited my previous post with my constructor already, but am writing again as there seems to be some confusion about terminology. When I talk about spins I mean that the helix might be represented by the value -1 and the coil by +1. I call these spins (because historically that's what the Ising model is about, and the terminology is very convenient even where the spins are not magnetic spins). Right, so mArray is an array that contains all the spins of all the particles, i.e. mArray[pos] gives the spin (either +1 or -1) at position pos.

So a better name for the class might be Polymer and instead of mArray, one might have mSpins or mPartsOfPolymer or something along those lines. (I start the names of all variables that are members of the class with an m; I've seen alternative conventions where people use this->variable when they refer to a variable that is defined in the class, Sometimes the variables are not distinguished at all, and sometimes prefixed with an "m_". I prefer the "m"-style. If you are asking why bother, it is easy to imagine writing a method that computes the temperature somehow. Now you might accidentally assign over the temperature in the class or shadow it: It is not clear from the code what you are trying to do. The "m"-style documents itself and you know when a developer is trying to change the member values and when they are just creating local data)

 

[Matt atkinson]

Oh, I understand now. It works, and i think its the c++ that's gotten me so confused I should probably take time to learn it thoroughly.

One last question how could I print to a file in this example, I tried how i would normally would but it doesn't work.

Thank you so much for your help, i know it must have been frustrating considering my lack of experience.

 

[Päällikkö]

Well you can run the program and output to a file "./prog.exe > file". If you want the program to do this automatically, see for example http://www.cplusplus.com/doc/tutorial/files/.

 

[Matt atkinson]

Thank you so much! i managed to get something like the attached file.

The errors would be calculated like this right?
$$\sigma=\sqrt{\frac{<a^2>-<a>^2}{N}}$$

 

[Päällikkö]

Great! If you comment out everything dealing with neighbors in the localEnergy() function (the loop), and multiply the variable ener by .5f, you should get a curve that goes on top of the green one.

Finally, there are a couple of issues in the code you should be aware of. Most critically the way the random generator is instantiated: Each time a particle is created, it constructs a new (Mersenne-twister) random generator, std::mt19937, with default arguments. This means that the seed is always the same, which is a bad thing (each temperature is going to get exactly the same random numbers). To avoid this, one could make the variable static: "static std::mt19937 mGen;". This requires an instantiation, "std::mt19937 Particle::mGen = std::mt19937{};" somewhere in the code (usually in the .cpp file where you define the class, i.e. does not have to be inside a function). What static does is "share" the variable among all variables that are of type Particle (so mGen will be destructed only at program exit rather than when each Particle is destructed).

There might also be a correlation between the random number generator inside the particle and the random number generator that is used to sample the Monte Carlo energy condition. One might avoid this by seeding one of the generators from system time or something else (like a Singleton class for a global random generator etc).

These are issues you don't need to worry about for simple coursework, but if you start doing serious research and/or start noticing curious non-random patterns, these might turn out to be the culprits.

This might have been a bit difficult to follow especially if you're not very comfortable with C++.

 

[Matt atkinson]

Once again, thank you so much you really helped me, helping me with my program and pointing out some of the pitfall's with it! :D
Now to calculate errors!

 

[Matt atkinson]

Just a question I am trying to work out errors with your code using a statistics class,

class statistics{
      private:
              int n;
              double sum;
              double sumsq;
      public:
             statistics();
             int getNumber() const;
             double getAverage() const;
             double getSqAverage() const;
             void add(double x);
};

int statistics::getNumber()             const{
                                              return n;
                                              }
double statistics::getAverage()         const {
                                              if(n==0) return 0.0;
                                              return sum/n;
                                              }
double statistics::getSqAverage()       const{
                                              if(n==0) return -1;
                                              return sumsq/n;
                                              }
void statistics::add(double x)                {
                                               n++;
                                               sum += x;
                                               sumsq += x*x;
       }

But I am not sure where to implement it in the main i tried doing something like;

res.add();
error=sqrt((res.getSqAverage()-res.getAverage()*res.getAverage())/(res.getNumber()));

Im pretty sure that res isn't the thing i should be using I am just strugglign to figure out what i should.

 

[Päällikkö]

Well, you're not really specifying what the errors might be, nor have you supplied a minimal example (i.e. please do not include Particle in the example: These objects are supposed to work independently) of using the Statistics class (I prefer user defined type names to be written with the first letter capitalized).

You did not write out the constructor: I hope you are initializing n, sum, and sumsq to zero. I also hope that res.add() is not an actual use case, but rather you are adding something (you did declare Statistics res; and didn't just try to use the std::vector<float> res; I had used, right? This part of your text was a bit ambiguous.).

If you are using the STL's sqrt from <cmath>, it might be better to be explicit about it: std::sqrt. I don't think the standard mandates the STL implementation to declare the function in the global namespace (but allows it to), whereas it must be present in std:: (don't quote me on this, though).

 

[Matt atkinson]

Oh okay, and my interpretation of the errors was that each value was an average of many measurements so each of the points i plot would have a spread of values that could be represented as an error bar?

Also no I didn't declare statistics res, I realize that was a mistake now.

Also sorry i forgot to paste the constructor,

statistics::statistics()                     {
                                              n=0;
                                              sum=0.0;
                                              sumsq=0.0;
                                              }

 

[Päällikkö]

So medres in my code collects several data points at each temperature, say 10000. res, on the other hand, stores an average of these 10000 values to represent the result at each particular temperature. What you want is not only the mean of medres (which is res), but also the standard deviation, sigma (which is an error bar estimate, although quite often people use 2*sigma as their error estimate, or at LHC/CERN, 5*sigma when they announced they had found the Higgs boson). From what I understand, you are trying to achieve a universal way of going about this with your statistics class, which is a good idea and ought to clarify the code.

If you paste your main, or better yet, a minimal example of how you want to use your statistics class (that is, one without interaction with the Particle class etc.), things would probably be easier to debug (if there sitll are bug). This might also clarify to yourself what the limits of your clss design might be (for example, it is not immediately clear to me how you suppose to use the class for multiple temperatures, maybe it is just a dropin for medres, rather than res?).

Finally, I'd recommend using

statistics::statistics() : n(0), sum(0.0), sumsq(0.0) {}

Here it doesn't really matter, but suppose n were const int rather than an int. Only the command I listed would then work (you can initialize with a value, but cannot assign one; same is true if you tried const int x; x = 1; versus const int x = 1; The latter is cheaper too when not talking about consts).

 

[Matt atkinson]

well for example i want to calculate $\sigma$ for each mean value of medres, for each temperature. So inside the loop for temperature, i could put something like this ;

A=medres.getAverage();
B=medres.getSqAverage();
Cnum=medres.getNumber();
Sigma=std::sqrt((B-A*A)/(Cnum));

Not sure if I am understanding it right.

 

[Päällikkö]

That'd be correct, yes, if you define "statistics medres;" and rather than "medres.push_back(particle.mean());" you use "medres.add(particle.mean());"
Sigma should be an array like res (res stores the means, Sigma should store the sigmas at each temperature), or rather res could be replaced by A should A also be made into an array. Given that you compute A, B, C using the methods of statistics, you'd do well to make Sigma a method of statistics, too (say, getError(), or whatever you consider appropriate).

 

[Matt atkinson]

Oh okay thankyou! ill try that when i get back to my code.
Define sigma, in the statistics class then have something like;

sigma=std::sqrt((sumsq-sum*sum)/(n))

i believe this should work
also with

class statistics{
      private:
              int n;
              double sum;
              double sumsq;
              double sigma;
              double medres;
      public:
             statistics();
             int getNumber() const;
             double getAverage() const;
             double getSqAverage() const;
             void add(double x);
};

If i have understood correctly

 

[Matt atkinson]

thankyou for all your help, just checked and all the errors i get are the exact same (~0.00042) have i done something incorrectly?
I think its my definition of medres2.

int main() {
  std::uniform_real_distribution<float> dist{0.f, 1.f};
  std::mt19937                          gen;

  std::vector<float> Ts;
  std::ofstream ofile;
  std::array<float, N>   A,B,C;
  statistics medres2;
  float error;
  ofile.open ("helixfrac.txt");
  for(float T = 420.f; T < 434.f; T+=.025f)
    Ts.push_back(T);
  std::vector<float> res;
  for(auto T : Ts) {
    Particle particle(T);
    std::vector<float> medres;
    medres.reserve(10000);
    for(int i = 0; i < 100000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy();
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%10) && i > 1000){
        medres.push_back(particle.mean());
        medres2.add(particle.mean());
                            }
    }
    float xx = 0.f;
    for(const auto x : medres)
      xx += x;
    res.push_back(xx/medres.size());
  }
  for(unsigned int i = 0; i < Ts.size(); ++i){
    A[i]=medres2.getSqAverage();
    B[i]=medres2.getAverage();
    C[i]=1/(double) medres2.getNumber();
    std::cout << Ts[i] << "," << (1+res[i])*.5 << std::endl;
    ofile << Ts[i] << "," << (1+res[i])*.5 << "," << std::sqrt((A[i]*C[i]-B[i]*B[i]*C[i])) << std::endl;
    }
    ofile.close();
  std::cout << std::endl;
}

 

[Päällikkö]

It sounds more or less correct. There are caveats in the way that you compute the error (it gives an underestimate), but it is not really your fault, it's a simplification made in your course materials. The reason is that the measurements are highly correlated with each other (the length of the polymer is, perhaps, N = 512, while the !(i%10) means that you are storing an average after every 10 attempts at flipping: Obviously even if all the flips were successful, you'd only move the average computer by Particle::mean() a few percentages at best). If you are interested in better methods, look up block-averaging and other autocorrelation based calculations (the idea is that you wait for as long as there is no correlation in the mean data to take averages). Vary N and the number of flips per temperature (i < 10000) and see how they affect the error estimates.

Finally, there's no reason for you to keep the medres's in your code. I cleaned it up a bit (still not pretty, but now it doesn't have anything that's not needed). I also use a couple of different loop styles, so you can compare and see what you feel comfortable with. For example, stats is an array (std::vector) of Statistics objects, so that each temperature gets its own Statistics. Now I initialize stat as the iterator (if you're not familiar with iterators, just think of it as a pointer) to the first Statistics object in the array, and then increment the iterator at each temperature. When I print out, I use a more traditional loop, where everything is indexed (note, though, that many data types, such as std::list, do not permit direct indexing, i.e. random access, but rather you must traverse through them).

int main() {
  auto dist = std::uniform_real_distribution<float>{0.f, 1.f};
  auto gen  = std::mt19937{};

  auto Ts = std::vector<float>{};
  for(float T = 426.f; T < 428.f; T+=.05f)
    Ts.push_back(T);

  std::vector<Statistics> stats(Ts.size());

  auto stat = stats.begin();
  for(auto T : Ts) {
    auto particle = Particle(T);
    for(int i = 0; i < 10000; ++i) {
      particle.choose();
      auto Eb = particle.localEnergy(); 
      particle.perturb();
      auto Ef = particle.localEnergy();
      if(std::exp(-(Ef-Eb)/T) < dist(gen))
        particle.perturb();
      if(!(i%100) && i > 5000)
        stat->add(particle.mean());
    }
    ++stat;
  }
  for(unsigned int i = 0; i < Ts.size(); ++i) {
    auto mean  = stats[i].getAverage();
    auto stdev = stats[i].getStdDeviation();
    std::cout << Ts[i] << " " 
              << (1+mean)*.5 << " " 
              << (1+mean-2*stdev)*.5 << " " 
              << (1+mean+2*stdev)*.5 << std::endl;
  }
  std::cout << std::endl;
}

You'll want to add the file outstream (i.e. std:fstream) instead of std::cout to write directly to a file. The reason I'm using 2*stdev is that I consider this a better estimate of the error than stdev (and this is indeed more commonly used). You might want to change this.

Finally, as you may notice, I made some modifications to your statistics class:

class Statistics {
public:
  int    getNumber() const;
  float  getAverage() const;
  float  getStdDeviation() const;
  void   add(double x);
private:
  int    n     = 0;
  double sum   = 0.f;
  double sumsq = 0.f;
};

int Statistics::getNumber() const {
  return n;
}

float Statistics::getAverage() const {
  if(n==0) 
    return 0.0;
  return sum/n;
}

float Statistics::getStdDeviation() const {
  if(n==0) 
    return 0.0;
  return std::sqrt(sumsq - sum*sum/n)/n;
}

void Statistics::add(double x) {
   n++;
   sum += x;
   sumsq += x*x;
}

[/code]

 

[Matt atkinson]

Wow thank you, your help has been greatly appreciated with this. I learned quite a lot about coding and physics too, while also doing the problem.
A lot of the code you used, i didn't understand at first but your explanations were great and anything else I looked up.

Thanks again!