I am doing simulations on a gold sample (20x20x10 lattice constants) at 0K. I want to calculate surface stress on the surface normal to z axis. Surface stress is calculated as (surface energy) + (first derivative of the surface energy with respect to strain @ strain=0) . I can get the surface energy (after minimization) as the energy difference between the sample with pps boundary conditions and the sample with ppp boundary conditions.

For calculating the derivative of the surface energy w.r.t strain, I need to plot surface energy vs strain; for which i should do tensile test on samples with pps and ppp boundary conditions, then calculate the surface energy at desired strain.

Doing uniaxial tensile test on ppp sample is quite easy with box_relax command. For pps sample, i can specify a non-zero (say 10,000 bars) pressure along x and zero pressure along y directions to simulate the uniaxial tensile test. In this case, the pressure along the z-direction does not become to zero. The pressure values after minimization are given below.
pxx pyy pzz (in bars)
9999.0856 -0.83783101 3555.6296

1. You stated that you want derivatives w.r.t. strain, but your simulations
increment the stress and measure strain. I think you will have an easier
time incrementing the strain, using fix deform, and measuring stress.

2. Some of the qunatities you are measuring are not well-defined for pps
boundary conditions. Pzz is not well-defined, because the dimension Lz is
somewhat arbitrary for pps. Also Pxx and Pzz are not well-defined because
they vary with distance z from the surface. In your simulation, you are
computing an average of Pxx(z) and Pyy(z) over the thickness of the sample.

3. I am also somewhat surprised that Pzz is not close to zero in your
simulation. Clearly, given that Pzz is positive, the energy could be lowered
by allowing the atoms to move further apart in the z direction.
I have not performed minimization with free surfaces, but I think the answer
is that the system is in a local energy minimum that is not the global
minumum, or the minimizer may not have fully converged. I suggest you
examine the results of your simulations more closely, and see how results
depend on initial condtions and simulation settings.