2
$\begingroup$

UPDATED

I've coded a multiparticle MD simulation in 3D. It is based on Langevin Dynamics, with random impulse and dissipation. I think the program works correctly now? I have attached the plots of kinetic, potential, and total energy. I believe that in the case of no Langevin Dynamics, total energy should be conserved as dissipation is turned off. But when dissipation is turned on, one shouldn't expect total energy conservation anymore. At this point, I would like to know if my results seem reasonable. I have included some of the more important code snippets (please ask for more if needed)

Physical parameters:

//V0 is the potential depth, r0 is the effective radius for the potential
double const m = 1., V0 = 1., r0 = 1., boxLength = 15., kT = 0.5;
int const n_atoms = 12;
double gamma = 0.5;   //damping coeff.
double var = 2.*gamma*m*kT*dt;  //variance of the gaussian distributed I
double c = (2.-gamma*dt)/(2.+gamma*dt);    //due to damping
double dt = 0.005;   //time step
int N = (int)(1000 / dt);   //simulation time

Initializing the system

void init(struct Atom system[n_atoms]){
    double rc = 0.9;
    double xx = 0., yy = boxLength, zz = boxLength;
    int atoms_in_row = (int)(boxLength / rc);
    int n = 0;

    for(int k = 0; k < (int)boxLength; k++) {
        for(int i = 0; i <= atoms_in_row; i++) {
            for(int j = 0; j <= atoms_in_row; j++) {
                if(n >= n_atoms) break;
                system[n].x = xx;
                system[n].y = yy;
                system[n].z = zz;
                system[n].vx = system[n].vy = system[n].vz = 0.;
                system[n].ax = system[n].ay = system[n].az = 0.;
                xx += rc;
                n++;
            }
            yy -= rc;
            xx = 0.;
        }
        yy = 0.;
        zz -= rc;
    }
}

Here are the energy plots with dissipation/random impulse turned off and on, respectively. The first one shows the plot of energy for two atoms only, in a non LD case. The next two are for 12 atoms. In particular, I think the result for non LD cases is fine as total energy is more or less conserved. And maybe the results are good for LD case as well? I'm not quite sure on this one.

2 atoms no LD No LD LD

Improve this question
$\endgroup$
  • $\begingroup$ Yes, the total energy is the sum of contributions of each pair. Regarding your other question, what happens when you decrease your timestep? $\endgroup$ – nicoguaro♦ Dec 24 '17 at 2:55
  • $\begingroup$ So then is the way I have implemented it correct? It's goes like $PE_{12}+PE_{13}+ ... +PE_{1N}+PE_{23}+PE_{24}+ ... +PE_{2N}+...+PE_{1-NN}$. And right now I'm using a time step of 0.02. My kinetic energy rises to about 0.3 and oscillates with small amplitudes about that point which is promising, but I'm not sure if it's right. In a sense, I'm not sure if my entire code is right :/ $\endgroup$ – Ptheguy Dec 24 '17 at 5:38
  • $\begingroup$ The number of the timestep is not that relevant right now. The question is about the behavior when you decrease that value, e.g., halve it every time $\endgroup$ – nicoguaro♦ Dec 24 '17 at 13:37
  • $\begingroup$ Potential Energy can change, total energy cannot. Also, ensure you are setting the net drift momentum to constant initially. Make sure you are using your if statement for the computation for potential where particle 1 and particle 2 are calculated in the same loop. $\endgroup$ – Ayush Agrawal Dec 26 '17 at 5:49
  • 1
    $\begingroup$ Please avoid making updates that change the original question. That would render the answer obsolete. You should ask new questions instead. $\endgroup$ – nicoguaro♦ Jan 5 '18 at 16:18
1
$\begingroup$

From your comments I think you should spend more time debugging your code. You have said that sometimes the program get trapped into an infinite loop: this is the first thing you need to fix.

I would suggest you to read three books that are very clear and well written in which concerns Molecular Dynamics method:

1 D. Frenkel & B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd Edition.

[2] Allen & Tildesley, Computer Simulation of Liquids.

[3] D.C. Rapaport, The Art of Molecular Dynamics Simulation

They bring routines implemented where you can compare to your code and verify if you are translating well the method to the programming. Follow the first one as the main reference for Molecular Dynamics and take a look into the others if you may have some doubt not clarified in the Frenkel.

Lennard-Jones Molecular Dynamics parametrisation [1]:

$U(r) = 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]$

tstep = 0.01; // time step
temperature = 0.728; 
epsilon = 1.0; // potential depth
sigma = 1.0; // effective radius
cut_off = 2.5*sigma; // LJ cut-off
min_distance = 0.87; // minimum initial distance between atoms
Improve this answer
$\endgroup$
  • $\begingroup$ The first suggestion was very helpful. It allowed me to spot some significant mistakes and correct them. But, now I think my code meets the requirements mentioned in the book, and still I'm facing issues. I've updated the post with my recent progress, if you are willing to take another look. $\endgroup$ – Ptheguy Jan 3 '18 at 14:22
  • $\begingroup$ Are you testing already for 100 atoms? Why don't you start testing for 2 atoms only? If it works for a small group of atoms it should work for a high number too. I don't think you are supposed to use a cut off to avoid divergence due to small distances, the hardcore component of the potential must be enough to repeal the atoms then avoiding the divergence. $\endgroup$ – The Doctor Jan 3 '18 at 15:31
  • $\begingroup$ I started with a 100 because it was easier to write the init method. I admit it's not written robustly at all, but it does the job for this configuration. I wrote a different MD simulation for 2 atoms and that was all good. I the expanded on that and the result is the above code that doesn't work properly. Exactly, I agree that the potential should repel the atoms, preventing them from being too close (or overlapping). I don't understand why its not doing that. That's exactly what's confusing me. $\endgroup$ – Ptheguy Jan 3 '18 at 16:04
  • $\begingroup$ Have you considered that the initial configuration you are imposing is already a divergent condition? Try to increase the distance between the atoms and verify if you are using a valid set of parameters (I'm concerned that the effective radius you are using is too small, be sure that this not the case). Re-write the init function to a more robust form, it will facilitate a lot the debugging. Debugging with 100 atoms is madness. $\endgroup$ – The Doctor Jan 3 '18 at 16:59
  • $\begingroup$ One of my problems is I do not know what constitutes a valid set of parameters. I've never done MD before so I've been using values which I obtained from a 1D MD simulation I had done before. How would I check for their validity? Trial and error? Is my initial condition divergent because there's not enough space between the atoms? Some other reason? I will work on the init and try much fewer atoms. I will update post accordingly. Thank you for following and helping. $\endgroup$ – Ptheguy Jan 3 '18 at 20:33
0
$\begingroup$

Based on a cursory look at your code, I think you might want to change the following lines in calculateForce:

        system[i].ax += f * dx / d;
        system[i].ay += f * dy / d;
        system[i].az += f * dz / d;
        system[j].ax += -system[i].ax;
        system[j].ay += -system[i].ay;
        system[j].az += -system[i].az;

and replace them with something to the effect of:

        system[i].ax += f * dx / d;
        system[i].ay += f * dy / d;
        system[i].az += f * dz / d;
        system[j].ax -= f * dx / d;
        system[j].ay -= f * dy / d;
        system[j].az -= f * dz / d;

Irrespective of your current problem, a good debugging strategy is to turn off the fluctuations and the dissipation introduced by Langevin dynamics (i.e., simply use Verlet) and see if you conserve energy. If the answer is negative, it means there is a bug that you need to fix. If you don't do this, it is possible that the stochastic dynamics will make it harder to spot a problem in your code.

Improve this answer
$\endgroup$
  • $\begingroup$ I can't believe I made such obvious mistake. I fixed that part of the code. I turned off dissipation and random force, and still the results are wrong. In particular, I'm thinking of the following: (1) are my physical parameters such as boxLength set to physically meaningful values? (2) is my PBC1D method corred? Because sometimes my program goes into an infinite loop, and the only place it could do this is in this method. $\endgroup$ – Ptheguy Jan 2 '18 at 11:43
  • $\begingroup$ Concerning your second point, you are using PBC1D to compute nearest neighbors distance. It's not correct. To compute the distance relative to x, you should compute dx=atm1.x-atm2.x then withdraw the correct number of boxlength, for example with a floor function, dx = dx-boxlength*floor(dx/boxlength). Plus Juan is right, you have to be in NVE ensemble to see if energy is conserved. Langevin thermostat sets you in NVT ensemble where temperature of the system is the controled quantity, not the energy. Not knowing how you initialize the system, we can be of no help concerning its convergence. $\endgroup$ – G.Clavier Jan 2 '18 at 22:55
  • $\begingroup$ Slight mistake in my above comment, I was more thinking about the round() function, not floor(). $\endgroup$ – G.Clavier Jan 2 '18 at 23:48
  • $\begingroup$ Isn't 'dx=dx-boxLength*round(dx/boxLength)' the same thing as my routine though? If dx/boxLength > 0.5 which means dx > boxLength/2, then dx=dx-boxLength. This is what I have, isn't it? $\endgroup$ – Ptheguy Jan 3 '18 at 5:40
  • $\begingroup$ Also, I'm starting with this configuration: all atoms on a 10*10*10 lattice separated by one unit, with velocities and accelerations of zero. $\endgroup$ – Ptheguy Jan 3 '18 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.