If only mode_compute=6 can be used to find mode, since 9 or 4 has ‘‘hessian matrix at the “mode” is not positive definite’’ issue. Then after I start at the mode found by 6 and run MCMC, say , 1 million draws, and finally all
estimated variables converge. Then would this result reliable? Should I do other things to enhance the probability of finding a global mode?

Problems like the one your report often signify some deeper underlying problem. Have you compared the posterior density found by the respective mode-finders? Is the one by mode_compute=6 really the best?
In general, finding the global mode is hard and often requires a lot of trial and error. Looking at the mode_check plots after mode-finding often reveals whether something is wrong.