So I left only half of that DPA ligand structure on PDB file and changed the occupancy of N1 and C4 to 1/2 while keeping the occupancy of all other atoms of the half ligand structure at 1.00.

Keep the occupancy of N1 and C4 as 1.00, then phenix.refine should automatically detect those as special positions and keep them exactly on the two-fold axis.

The bonds and other restraints to the atoms on special positions will not work automatically. Currently, the only option is to define custom bonds (see phenix.refine docs) with symmetry_operation=x,y,z. You can use Coot to obtain the symmetry operation (turn on drawing of symmetry copies; when clicking atom Coot will show the sym op). Emulate angle and dihedral restraints with pseudo-bonds. Inspect the .geo file to verify that all restraints are correct.

You can also keep all DPA atoms, with all occupancies set to 0.5. Then the restraints will work as usual; probably you don't have to do anything special. However, the N1 and C4 are allowed to float around the special position. It may be OK though. I think I'd try this option first.

Ralf