Abaqus VUMAT doesn't given identical results for orthotropic linear elastic material

Dear iMechanica community,

Hi! I am trying to verify Abaqus VUMAT with simple examples so I wrote a VUMAT for orthotropic material and compare the result with Abaqus built-in material model (Elastic -> Engineering constants). The material constants are given in 1.png and 2.png. All other settings from part to mesh is identical. The stress distribution of built-in material model is in 3.png and the stress distribution of VUMAT is in 4.png. The discrepancy is not ignorable. The VUMAT code is also attached as a .txt file.

Is there anything that I miss in my code? I will appreciate your help. Thanks! If the problem were solved, I will write a tutorial for this and post it on iMechanica and Github.

Notes: The Fortran code is based on VUMAT from Rutger's repo: https://github.com/rutger-kok/composite_cdm_tan