We present results from our lattice QCD study of the contribution of the isovector quark chromoelectric dipole moment operator to the neutron EDM. The calculation was carried out on four 2+1+1-flavor highly improved staggerred quark ensembles (provided to us by MILC collaboration) using Wilson-clover quarks to construct correlation functions. We use the nonsinglet axial Ward identity including corrections up to O(a) to show how to control the power-divergent mixing with the lower dimensional pseudoscalar operator. Results are presented after conversion to the modified minimal-subtraction scheme at the leading-log order.