Dear @Lulu.Liu,
I agree that the WAMIT .hst file to be used in HydroDyn should equal [C_ij]^Hydrostatic from your first equation. And I agree that this restoring matrix should not include the effects of body weight because body weight is intrinsically included in the structural modules of OpenFAST (ElastoDyn, SubDyn, BeamDyn), which would otherwise be double counted if included in [C_ij]^Hydrostatic. This restoring matrix will be different than the restoring matrix in your last equation (from Newman) because the restoring in your last equation must include the contribution from both buoyancy and body weight–i.e., the total restoring.
The undisplaced buoyancy term (rho * g * V_0) does not appear in your last equation (from Newman) because it should cancel with the weight of the floating system and vertical undisplaced mooring pretension.
Best regards,