Thanks for your help.
1, I fixed VGG’s parameters and trained the last layer. It outputs about 35%. Not as expected ( 60%-70%)
2, I use VGG16 without batch-normalization, which has the same structure with its original MatConvnet model.
3, I follow the official tw-step training: first to extract deep features, and then train the last layer. It  helps and boost the accuracy.( But I don’t know why. In my mind, freezing VGG layers is equal to extracting features first and training the last layer later.) The accuracy of first stage is 59% (80.1% reported in the paper), the accuracy of the second stage is 78%(84.1% reported).
4, I use the SGD optimizer with momentum (same as the paper).
5, Interestingly, the Tensorflow version 9 works well. (69% for the 1st stage and 84% for 2rd stage). Is the pretrained VGG16 model’s fault in different DL frameworks? I love Pytorch, but it does not work well this time…