SHARE
TWEET

Untitled

a guest Oct 20th, 2018 109 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. void LeapfrogStep(int part)
  2. {
  3.     int i;
  4.     if (part == 1) {
  5.         for (i = 0; i < nMol; i++) {
  6.             mol[i].rv.x += 0.5 * deltaT * mol[i].ra.x;
  7.             mol[i].rv.y += 0.5 * deltaT * mol[i].ra.y;
  8.             mol[i].r.x += deltaT * mol[i].rv.x;
  9.             mol[i].r.y += deltaT * mol[i].rv.y;
  10.         }
  11.     }
  12.     else {
  13.         for (i = 0; i < nMol; i++) {
  14.             mol[i].rv.x += 0.5 * deltaT * mol[i].ra.x;
  15.             mol[i].rv.y += 0.5 * deltaT * mol[i].ra.y;
  16.         }
  17.     }
  18. }
  19.  
  20. void ComputeForces()
  21. {
  22.     int i, j;
  23.     VecR dist, forces;
  24.     real r, force;
  25.     for (i = 0; i < nMol - 1; i++)
  26.         for (j = i + 1; j < nMol; j++) {
  27.             dist.x = mol[j].r.x - mol[i].r.x;
  28.             dist.y = mol[j].r.y - mol[i].r.y;
  29.             r = sqrt(dist.x * dist.x + dist.y * dist.y); //distance formula
  30.             force = 24 * (r*r*r*r*r*r - 2) / (r*r*r*r*r*r*r*r*r*r*r*r*r); //derivative of the Lennard-Jones potential
  31.             forces.x = dist.x / r * force;
  32.             forces.y = dist.y / r * force;
  33.             mol[i].rv.x += forces.x; //the program works correctly when I apply the force to rv (velocity) but not ra (acceleration)
  34.             mol[i].rv.y += forces.y;
  35.             mol[j].rv.x -= forces.x;
  36.             mol[j].rv.y -= forces.y;
  37.         }
  38. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top