# Homework Help: 1D ising model - Helix-coil transistion

1. Dec 3, 2014

### Matt atkinson

1. The problem statement, all variables and given/known data
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.

2. Relevant equations

3. The attempt at a solution
Just confused as to if im 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?

Last edited: Dec 3, 2014
2. Dec 4, 2014

### 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).

3. Dec 4, 2014

### hibphys

hey how do i post a question ?

4. Dec 5, 2014

### Matt atkinson

Ah thank you that cleared a lot up, i think ive 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.

Last edited: Dec 5, 2014
5. Dec 5, 2014

### 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.

6. Dec 6, 2014

### 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 im 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}}$$

7. Dec 6, 2014

### 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.

8. Dec 6, 2014

### 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.

File size:
4 KB
Views:
140
9. Dec 6, 2014

### 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
Code (Text):

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():
Code (Text):

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)
Code (Text):

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).

10. Dec 7, 2014

### 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).

11. Dec 7, 2014

### 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>).

Last edited: Dec 7, 2014
12. Dec 7, 2014

### Matt atkinson

Okay so I think I understand, in my code im 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 jsut for a comparison point of view would I just have to write a main function similar to my current one?

13. Dec 7, 2014

### Päällikkö

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

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).

14. Dec 7, 2014

### 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 alot of the code you send because i havent gone really far into c++ as of yet.

15. Dec 7, 2014

### 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:
Code (Text):

#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.

16. Dec 7, 2014

### 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 im not sure what to do with that line "g++ -std=c++11 -O3 file.cpp".

17. Dec 7, 2014

### 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).

18. Dec 7, 2014

### 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"

19. Dec 7, 2014

### 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).

20. Dec 7, 2014

### Matt atkinson

Ah yea sorry, completely missed all my includes....