Why is the energy decreasing in my atom simulation?

[ebrattr]

I'm trying to simulate a long chain of atoms where I restrict the bond (length between two consecutive atoms) and angle (between three consecutive atoms). Therefore, if I want to simulate N particles, then I have N+3 degrees of freedom. The only interactions between atoms is due to Lennard-Jones Potential.
I am using the classical Runge-Kutta fourth order method to integrate the equations. But, I got that the energy is always decreasing and then it keeps constant, like in the picture shown below.

http://img194.imageshack.us/img194/3757/examplexy.png

What can I do to keep constant the energy ?

What is happening ? I'm new with the topic. My english is not that good, so if you do not understand something, just let me know it.
Thanks !

 

[A.T.]

What velocity dependent forces do you have exactelly? If you have Coriolis forces then you are obviously in a rotating frame. To get energy conservation in a rotating frame you have to include the potential energy due to the inertial centrifugal force. If the angular velocity of the frame is not constant, you can't expect energy to be consereved at all.

 

[ebrattr]

What velocity dependent forces do you have exactelly? If you have Coriolis forces then you are obviously in a rotating frame. To get energy conservation in a rotating frame you have to include the potential energy due to the inertial centrifugal force. If the angular velocity of the frame is not constant, you can't expect energy to be consereved at all.

No, I don't have velocity dependent forces, because my only forces are due to a Lennard-Jones Potential. I was a little confused. What I do really have is a rotating frame. Got it. I was not including the potential energy due to inertial centrifugal force. I was just including, kinetick energies (\dfrac{1}{2}mv^2) and the corresponding LJ potential. However, these velocities are the global velocities of the particles.

Thanks !

 

[A.T.]

No, I don't have velocity dependent forces, because my only forces are due to a Lennard-Jones Potential. I was a little confused. What I do really have is a rotating frame.

In a rotating frame you have inertial forces too (Centrifugal and Coriolis). And Coriolis is velocity dependent, but it is not doing any work, so it is irrelevant for Energy conservation.

However, these velocities are the global velocities of the particles.

To make an energy analysis you have to transform every thing into one frame.

 

[DrClaude]

Can you post your program?

 

[ebrattr]

Yeah, I compute all the velocities to one frame, my global frame.

 

[A.T.]

Yeah, I compute all the velocities to one frame, my global frame.

And that global frame is rotating?

 

[ebrattr]

No. It is not rotating.

 

[ebrattr]

Can you post your program?

Yeah, off corse. But what part ? There many lines of code.

 

[DrClaude]

How do you implement the constraints?

 

[ebrattr]

How do you implement the constraints?

By setting every positiion in terms of those coordinates I mentioned before? I don't understand the question.

 

[DrClaude]

I was asking about this

I'm trying to simulate a long chain of atoms where I restrict the bond (length between two consecutive atoms) and angle (between three consecutive atoms).

I took "restrict" to mean that you had constraints of the bond lengths and angles.

Do I understand correctly that you don't have an angular potential? Do you calculate the Lennard-Jones potential over all particle pairs?

 

[ebrattr]

I was asking about this

I took "restrict" to mean that you had constraints of the bond lengths and angles.

Do I understand correctly that you don't have an angular potential? Do you calculate the Lennard-Jones potential over all particle pairs?

Yes. I do not have an angular potential. And I calculate Lennard Jones potential over all particle pairs, except those that have an angle or a bond.

In fact, in the most simple case, with M chains and 2 particles for every chain, I get the same thing.

In fact, can we only discuss this particular case ? To see what is going on here.

 

[DrClaude]

Before we get into the particulars, I have a couple more questions.

Do you use a cutoff for the LJ potential?

Have you tired changing the time step?

 

[ebrattr]

Before we get into the particulars, I have a couple more questions.

Do you use a cutoff for the LJ potential?

Have you tired changing the time step?

Yeah I am tired changin the time step. And I'm not using a cutoff for the LJ potential.

 

[DrClaude]

What is the smallest possible system you can simulate and what are the results then?

 

[DrClaude]

I'll try to have a look at all that tomorrow.

 

[A.T.]

Forces:
<br /> \mathbf{F}_{1}^{1}=f_{LJ}(\mathbf{r}_{1}^{2},\mathbf{r}_{2}^{2},\mathbf{r}_{1}^{1})/m<br />

What are F, f_LJ and m here?

 

[ebrattr]

What are F, f_LJ and m here?

F_i^j is the resultant force over the particle i in the j-chain. f_{LJ} is the corresponding magnitude.

 

[A.T.]

F_i^j is the resultant force over the particle i in the j-chain. f_{LJ} is the corresponding magnitude.

The magnitude of what? And what is m?

 

[ebrattr]

The magnitude of what? And what is m?

Mass. Magnitude of the force due to Lennard Jones potential.

 

[A.T.]

Mass. Magnitude of the force due to Lennard Jones potential.

So you divide force by mass and get force again?

 

[ebrattr]

So you divide force by mass and get force again?

Suppose I have mass of 1 kg. This not real important.

 

[ebrattr]

Thanks !
I have found the problem !, when I compute the Total energy, in the internal cycle, it is not k \neq j, it is instead k &gt; j.