13.7 Concurrent coupling of length scales: FEM and MD 443
Obviously, we could take periodic boundary conditions in the transverse direction,
which implies a cylindrical description (this could be useful for describing a carbon
nanotube). However, this is not compatible with longitudinal waves, as the Poisson
ratio causes the system to expand where it is longitudinally compressed and vice
versa. A periodic boundary condition would not allow for this to happen and cause
unphysical, strong stresses to build up.
Once a FEM and a MD program are working, it is not so much work to couple
them along the lines decribed above. We use the velocity Verlet algorithm which
naturally splits into two steps, separated by a force calculation.
First step:
pi(t+h)=pi(t)+
h
2
Fi[r(t)]; (13.58a)
ri(t+h)=ri(t)+hpi(t); (13.58b)
CalculateFi[r(t)];
Second step:
pi(t+h)=pi(t+h)+
h
2
Fi[r(t+h)]
These steps must be kept in mind when setting up the algorithm for the full system.
This algorithm looks as follows:
Calculate MD forces;
Calculate FEM forces;
Copy locations of leftmost MD points to a shadow array
in the FEM region;
Copy locations of rightmost FEM points to a shadow array
in the MD region;
FOR TimeStep = 1, MaxStep DO
Set Initial values of boundary points;
Do first integration step (seeEq. (13.58a));
Copy locations of leftmost MD points to a shadow array
in the FEM region;
Copy locations of rightmost FEM points to a shadow array
in the MD region;
Calculate forces in FEM region, including those on the MD particles;
Calculate forces in MD region, including those on the FEM particles;
Add FEM forces acting on MD particles to MD forces;
Add MD forces acting on FEM particles to FEM forces;
Do second integration step [see Eq. (13.58b)];
END FOR.