For problems involving large material deformation rate, the material deformation time scale can be shorter than the material takes to reach a thermodynamical equilibrium. For such problems, it is difficult to obtain a constitutive relation. Furthermore, history dependency are usually important in these problems because of the thermodynamic non-equilibrium. A numerical method capable of tracking material deformation history is needed in a numerical simulation effort. A multi-scale numerical method considering non-equilibrium thermodynamics is developed based on dual domain material point method (DDMP). The DDMP method uses Lagrangian material points to track the history of the material whereas Eulerian grids are used to calculate the gradients in continuum level. Molecular dynamics (MD) calculations are performed in the material points and to calculate the closure quantities such as stress, under the condition of non-equilibrium thermodynamics, bypassing the need for a constitutive relation. Since the material points only communicate with mesh nodes and do not communicate among themselves, the MD calculations can be done in embarrassingly parallel. An efficient CUDA code is developed to accelerate MD calculations in GPU. Examples of shock wave propagation in Cerium and Copper single crystal are presented.