headmyshoulder / odeint-v2

odeint - solving ordinary differential equations in c++ v2
http://headmyshoulder.github.com/odeint-v2/
Other
344 stars 101 forks source link

feature request: Velocity Verlet stepper #89

Open kirilligum opened 11 years ago

kirilligum commented 11 years ago

Velocity Verlet is the most used integration method and used in all/most molecular dynamics packages, physics simulations and games.

http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet

headmyshoulder commented 11 years ago

Should be quite easy. I already implemented once a velocity verlet method.

On 17.06.2013 13:44, kirilligum wrote:

Velocity Verlet is the most used integration method and used in all/most molecular dynamics packages, physics simulations and games.

http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet

kirilligum commented 11 years ago

I agree. I meant it as a part of a library. When I would have time, I would love to submit it but I''m not optimistic about finding time in near future. I just wanted to put it out there in case someone has already a working version. not a top priority since there are already a couple of simplistic algorithms.

headmyshoulder commented 11 years ago

Hi Kirill, sorry for the delay. I have an implementation of the velociry verlet now but I do not have any test ODE. Can you look in the velocity_verlet branch, velocity_verlet.hpp and libs/numeric/odeint/test/velocity_verlet.cpp if the code is working as expected? Step size control and dense output are not provided.

kirilligum commented 11 years ago

will do. step size is not possible with verlet. thank you.

headmyshoulder commented 11 years ago

On 04.08.2013 19:33, kirilligum wrote:

will do. step size is not possible with verlet. thank you.

just one remark: I added the time as a parameter to system function such that you can solve non-autonomous systems.

headmyshoulder commented 11 years ago

@kirilligum :I finished the stepper, even with docs but you have to build them manually right now, since the stepper is not in the master branch. Could you provide an example for the tutorial, or even write an entry in the tutorial?

kirilligum commented 11 years ago

looking into it. i'll reply in the next few days.

On Tue, Aug 27, 2013 at 1:04 AM, headmyshoulder notifications@github.comwrote:

@kirilligum https://github.com/kirilligum :I finished the stepper, even with docs but you have to build them manually right now, since the stepper is not in the master branch. Could you provide an example for the tutorial, or even write an entry in the tutorial?

— Reply to this email directly or view it on GitHubhttps://github.com/headmyshoulder/odeint-v2/issues/89#issuecomment-23319179 .

kirilligum commented 11 years ago

I have trouble compiling your code. I'm not sure how odeint-v2-velocity_verlet fits into the boost installation. I changed some paths in files but i get an error finding unit-testing framework. I'm gonna try to figure it out but it would also help if you can give a hint about your set up of boost 1.54 and odeint-velocity-verlet.

the simplest system is with potential energy V(x)=x_x_0.5, so the hamiltonian (total) energy is h(x,v)=v_v_0,5+x_x_0.5. the only thing you need for verlet is a(x)=-dV(x)/dt=-x

here is a simple code:

include

include

include

include

using namespace std;

struct harmonic_system { double energy(double p, double q){ return 0.5_p_p+0.5_q_q; } double a(double q){ return -q; } };

int main(int argc, char const _argv[]) { harmonic_system sys; double empty , t=0.0 , p=0 , q=1.0 , adt=sys.a(q) , at=adt , dt=1e-1 ; for (size_t i = 0; i < 100; ++i) { t+=dt; at=adt; q+=p_dt+0.5_at_dt_dt; adt=sys.a(q); p+=(at+adt)_0.5*dt; cout << t << " " << q << " " << p << " " << adt << " " << sys.energy(p,q) << endl; } return 0; }

I'll try to look into how to compile odeint-v2-velocity_verlet and see if i can write a sample program.

On Wed, Aug 28, 2013 at 12:24 AM, Kirill Igumenshchev bribeme@gmail.comwrote:

looking into it. i'll reply in the next few days.

On Tue, Aug 27, 2013 at 1:04 AM, headmyshoulder notifications@github.comwrote:

@kirilligum https://github.com/kirilligum :I finished the stepper, even with docs but you have to build them manually right now, since the stepper is not in the master branch. Could you provide an example for the tutorial, or even write an entry in the tutorial?

— Reply to this email directly or view it on GitHubhttps://github.com/headmyshoulder/odeint-v2/issues/89#issuecomment-23319179 .

headmyshoulder commented 11 years ago

Hi @kirilligum ? Are boundary condition important for? As fas as I know the Velocity Verlet is often used for molecular dyanmics simulations and here boundary conditions are crucial. I think of introducing an adapter for the BC.

kirilligum commented 11 years ago

boundary condition, like a periodic box?

that is usually done outside of the integrator. also other things can be involved to complicate a general integrator. for example interacting with other particles across the boundary. to implement the periodic box, I would just subtract the size of the box when the particle crosses the boundary (x=x-box*(int)(x/box)). also, in eom, the other particles that involved into calculating the energy derivative should be mirrored from the original box to the neighboring boxes.

but I guess, you can add integrate_periodic function for usability.

Understanding Molecular Simulation, Second Edition: From Algorithms to Applications (Computational Science) Daan Frenkel (Author), Berend Smit (Author) has an example on page 86 -- algorithm 7

On Wed, Sep 18, 2013 at 12:40 PM, headmyshoulder notifications@github.comwrote:

Hi @kirilligum https://github.com/kirilligum ? Are boundary condition important for? As fas as I know the Velocity Verlet is often used for molecular dyanmics simulations and here boundary conditions are crucial. I think of introducing an adapter for the BC.

— Reply to this email directly or view it on GitHubhttps://github.com/headmyshoulder/odeint-v2/issues/89#issuecomment-24693140 .

headmyshoulder commented 11 years ago

Ok, I think you are right. One should do the boundary conditions outside of the integrator. It has the consequence that one needs to take them into account in the system function. Here, points can possibly be outside of the domain, since a "streaming step" is done before the system function is called. But it is not very difficult.

On 09/18/2013 10:30 PM, kirilligum wrote:

boundary condition, like a periodic box?

that is usually done outside of the integrator. also other things can be involved to complicate a general integrator. for example interacting with other particles across the boundary. to implement the periodic box, I would just subtract the size of the box when the particle crosses the boundary (x=x-box*(int)(x/box)). also, in eom, the other particles that involved into calculating the energy derivative should be mirrored from the original box to the neighboring boxes.

but I guess, you can add integrate_periodic function for usability.

Understanding Molecular Simulation, Second Edition: From Algorithms to Applications (Computational Science) Daan Frenkel (Author), Berend Smit (Author) has an example on page 86 -- algorithm 7

On Wed, Sep 18, 2013 at 12:40 PM, headmyshoulder notifications@github.comwrote:

Hi @kirilligum https://github.com/kirilligum ? Are boundary condition important for? As fas as I know the Velocity Verlet is often used for molecular dyanmics simulations and here boundary conditions are crucial. I think of introducing an adapter for the BC.

— Reply to this email directly or view it on GitHubhttps://github.com/headmyshoulder/odeint-v2/issues/89#issuecomment-24693140

.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/89#issuecomment-24696700.

headmyshoulder commented 11 years ago

Hi @kirilligum I have added an example for molecular dynamics now. It looks reasonable for me. Maybe you have a look. If your work on Linux you can pipe the output to gnuplot which should show you a nice animation of dynamics.

kirilligum commented 11 years ago

@headmyshoulder awesome. I'll check it in the next few days