You should determine the water box sizes based on the protein unfolded structure, using for example one of the unfolded configurations. Based on this unfolded structure, you can determine (x_max, x_min), (y_max, y_min) and (z_max,z_min) of protein. Then, construct a water layer around the protein with thickness at least the cutoff of non-bonded interactions size: let's say if the cutoff=11-12 Angstrom, then the thickness of the layer is 15 Angstrom (assuming that the water box may shrink a little if the water box is not initially equilibrated at normal conditions, i.e., density = 1 gr/cm^3). Then, the big water box sizes are defined as:

big_sizeX = x_max-x_min + 2*thickness

big_sizeY = y_max-y_min + 2*thickness

big_sizeZ = z_max-z_min + 2*thickness

Then, after building a water box with these dimensions (big_sizeX, big_sizeY, big_sizeZ), you can solvate the folded protein structure in this 'big' water box. With this, then, you are sure that even if the protein gets unfolded, the images will not see each other.

Thank you Hiqmet Kamberaj.

But where do I get this information from? How can I determine (x_max, x_min), (y_max, y_min) and (z_max,z_min) of protein from trajectory files?

You can use vmd, for example, with command 'measure minmax' in TkConsole.

Thank you very much.

I calculated minmax using VMD as you suggested and calculated big_sizeX, big_sizeY and big_sizeZ taking the thickness of 10. How do I correlate these values to size of the cubic box (-d option) I need to use in the following command:

gmx editconf -f protein.gro -o protein_newbox.gro -c -d 1.0 -bt cubic

Thickness 10 angstrom with a cutoff of at least 10 angstrom (as usually suggested for biomolecular systems) is a little small. It is recommended that the thickness should be greater than cutoff of non-bonded interactions. As I said in the previous post,  first, 

big_sizeX = x_max-x_min + 2*thickness

big_sizeY = y_max-y_min + 2*thickness

big_sizeZ = z_max-z_min + 2*thickness

then,

d = max(big_sizeX, big_sizeY, big_sizeZ)

where (x_min,x_max,y_min,y_max,z_min,z_max) are calculated using the most unfolded protein structure that you may have to see during the simulations.

Thank you very much, I have changed the thickness value to 15 and now my simulation works fine. Thank you very much.