I don't know much about molecular dynamics, but I am suspicious of my pressure being equal to kinetic energy. I am using
// sum the velocities squared:
double v2sum = 0.0;
for (int i = 0; i < N; i++) {
v2sum += vx[i]*vx[i] + vy[i]*vy[i];
}
KE = 0.5 * v2sum;
KE /= N; // kinetic energy per particle
Tp = KE; // reduced temperature (boltzmann's (kB) = 1)
for temperature and
//instantaneous pressure calcs
double w = 0.0;
w += 0.5 * (dx*fx + dy*fy);
Pr = N*Tp+w;
for pressure.
Hi Kaley Collins , the equation for the instantaneous pressure estimator is:
Pressure = (1/(3*Vol))*sum([(P(i)^2/(2*m(i)))+R(i).F(i)])
where I am using capital bold letters to notate vector quantities: R(i) is the position vector (x(i), y(i),z(i)) of atom i, F(i) is a vector containing the x, y, and z components of the forces acting on atom i, and R(i).F(i) is a dot product, m(i) is the mass of atom i, and P(i)^2 is essentially the dotted product of the momentum vector with itself i.e. P(i).P(i). Furthermore, the sum should run over all atoms (i) in your simulation from one to N, and "Vol" is the volume of your simulation cell. So when I look at your implementation, I am assuming you are working in two dimensions since you don't have a z-component anywhere. I also notice that you are not dividing by (3*Vol) as you should. I also see that your implementation of the second term is incorrect. First of all, you don't need to multiply by (0.5). Also, the equation indicates R(i) and not delta(R(i)) and it seems that you are using a displacement in your equation. So after every time step in your MD simulation you should calculate the total position vector for each atom and use that in your equation. It's a little hard to clearly notate mathematical equations in this textbox, but I hope that is helpful.
Adib Samin because I am running a 2D simulation what value do I use for volume?
- Badges
- Science topic
- Similar topics
- Chemistry
- Theoretical Chemistry
- Computational Chemistry
- Molecular Dynamics