0
$\begingroup$

Now there are two particles $j^{th}$ and $i^{th}$, Their coordinates are $\mathbf{r}_j (1.0, 0.0)$, $ \mathbf{r}_i (3.0, 1.0)$ respectively. Their speed $\mathbf{v}_j (-0.3, 0.0)$, $ \mathbf{v}_i (-0.5, -0.3)$ respectively. Their acceleration $\mathbf{a}_j (-0.1, 0.0)$, $ \mathbf{a}_i (-0.2, -0.01)$ respectively. How long does it take for two particles to be separated by $rcut = 1.30 $? My approach is to first find the components of the relative velocity and relative acceleration in the direction of $\mathbf{r}_{ij}$. $$ \begin{align} \mathbf{r}_{ij} =& \mathbf{r}_j - \mathbf{r}_i \\ \mathbf{v}_{ji} =& \mathbf{v}_i - \mathbf{v}_j \\ \mathbf{r}_{ji} =& \mathbf{a}_i - \mathbf{a}_j \\ \mathrm{dr} = &| \mathbf{r}_{ij} | \\ \mathrm{dv} = &\dfrac{ \mathbf{v}_{ji} \cdot \mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} \\ \mathrm{da} = &\dfrac{ \mathbf{a}_{ji} \cdot \mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} \end{align} $$ Then use Newton's equations of motion to find dt: $$ rcut = dr + dv*dt + \dfrac{1}{2} da dt^2 $$ My problem is that when I write a program to find dt, when $j^{th}$ and $i^{th}$ move through dt time, the distance between them is not rcut, and there is a big error. How to eliminate this error, at least to $10^{-10}$ level of accuracy.
This is my code：

program main
  implicit none

  real(8), dimension(1:2) :: rj, ri, rij
  real(8), dimension(1:2) :: vj, vi, vji
  real(8), dimension(1:2) :: aj, ai, aji
  real(8) :: rcut, dr, dv, da, dt

  rcut = 1.30d0

  rj = (/1.0d0, 0.0d0 /)    ;  ri = (/3.0d0, 1.0d0 /)
  vj = (/-0.30d0, 0.0d0 /)  ;  vi = (/-0.50d0, -0.30d0 /)
  aj = (/-0.10d0, 0.0d0 /)  ;  ai = (/-0.20d0, -0.0010d0 /)

  rij = rj - ri
  vji = vi - vj
  aji = ai - aj

  dr = sqrt( rij(1)*rij(1) + rij(2)*rij(2) )
  dv = dot_product(vji, rij)/dr
  da = dot_product(aji, rij)/dr

  dt = ( -dv + sqrt( dv*dv-2.0d0*da*(rcut-dr) ) )/da

  write(*,*) 'dr=', dr
  write(*,*) 'dv=', dv
  write(*,*) 'da=', da
  write(*,*) 'dt=', dt

  rj = rj + vj*dt + 0.50d0*aj*dt*dt
  ri = ri + vi*dt + 0.50d0*ai*dt*dt

  rij = rj - ri

  write(*,*) 'new dr=', sqrt( rij(1)*rij(1) + rij(2)*rij(2) )

end

The result of running is as follows:

 dr=   2.2360679774997898
 dv=  0.31304951684997051
 da=   8.9889932695491545E-002
 dt=   2.2580918192732384
 new dr=   1.3324349346471371

It can be seen that the error between new dr and rcut is very large. How to deal with it?

$\endgroup$
8
  • $\begingroup$ You're missing some terms in $da$, the derivative of $dv$: $$da=\frac{\mathbf a\cdot\mathbf r}{\lVert\mathbf r\rVert}+\frac{\mathbf v\cdot\mathbf v}{\lVert\mathbf r\rVert}-\frac{(\mathbf v\cdot\mathbf r)^2}{\lVert\mathbf r\rVert^3}$$ And $da$ is not constant. The equation of motion you're using is only valid for constant acceleration. $\endgroup$
    – mr_e_man
    Commented Jun 29, 2022 at 12:12
  • $\begingroup$ @mr_e_man $da$ is not the derivative of $dv$, but the acceleration difference between $j^{th}$ and $i^{th}$. $\endgroup$
    – Zhao Dazhuang
    Commented Jun 29, 2022 at 12:19
  • $\begingroup$ The vector acceleration $\mathbf a=\mathbf a_i-\mathbf a_j$ is constant (I assume). So $$\mathbf r(t)=\mathbf r(t_0)+\mathbf v(t_0)\,(t-t_0)+\frac12\mathbf a(t_0)\,(t-t_0)^2$$ which describes a parabola in the plane. And $\lVert\mathbf r\rVert=r_{cut}$ describes a circle with radius $r_{cut}$. So you want to find the intersection of these two curves. $\endgroup$
    – mr_e_man
    Commented Jun 29, 2022 at 12:29
  • 1
    $\begingroup$ Taking the dot product of my vector equation with $\mathbf r(t_0)/\lVert\mathbf r(t_0)\rVert$, I get $$\mathbf r(t)\cdot\mathbf r(t_0)/\lVert\mathbf r(t_0)\rVert=dr+dv\,(t-t_0)+\frac12da\,(t-t_0)^2$$ Why are you thinking that $\mathbf r(t)\cdot\mathbf r(t_0)/\lVert\mathbf r(t_0)\rVert=r_{cut}$ ? $\endgroup$
    – mr_e_man
    Commented Jun 29, 2022 at 12:33
  • $\begingroup$ @mr_e_man I got it, thank you very much for your answer $\endgroup$
    – Zhao Dazhuang
    Commented Jun 29, 2022 at 12:35

0

Reset to default

You must log in to answer this question.