Don't use ISIF=3 but use  ISIF=4  instead! Then the cell shape only (b/a, c/a) will be relaxed at  *fixed volume* and the volume itself can be given in POSCAR on the second line as a negative value (i.e., minus volume). From that and the initial three lattice vectors below VASP calculates internally the correct a, b, c. After relaxation you will find changes in these three lattice vectors (lenghts) but such that volume will remain conserved. At the end, on CONTCAR you can then find as usual a "scaling parameter" in the second line and the three altered lattice vectors from which you can (after multiplication with the "scaling parameter") determine the (relaxed) lattice constants a, b, and c for the given volume. And of course you will also know the total energy then but even hydrostatic pressure printed in OUTCAR ("pressure = ..."). So you could not only determine an energy-volume but even a pressure-volume curve!

Only if you have no idea at all about the equilibrium lattice constants / volume you might consider a "pre-relaxation" run with ISIF=3 first before you begin to sample different volumes around the (approximate) equilibrium volume with ISIF=4 ...

Just take care that your plane-wave cutoff is high enough (in order to avoid too large Pullay stresses -- I usually recommend to approximately double the cutoff with respect to what is set by VASP for PREC=high/accurate) and that you perform enough relaxation steps to get all stress tensor components converged to at least better than order of 0.1 kBar. And also use a rather small EDIFF (use maybe values of about 1-9 to 1-6, in order to obtain well-converged wave functions which are important for numerically accurate forces and stresses). And overall try to be as accurate as possible (e.g. PREC=high or PREC=accurate is recommended, maybe also ADDGRID=.TRUE.).

This strategy works by the way for all types of (non-cubic) lattices, also monoclinic or triclinic ones (getting not only the three lattice constants but also the one or three angles). And even for tetragonal or hexagonal ones it can still be useful.

Don't use ISIF=3 but use  ISIF=4  instead! Then the cell shape only (b/a, c/a) will be relaxed at  *fixed volume* and the volume itself can be given in POSCAR on the second line as a negative value (i.e., minus volume). From that and the initial three lattice vectors below VASP calculates internally the correct a, b, c. After relaxation you will find changes in these three lattice vectors (lenghts) but such that volume will remain conserved. At the end, on CONTCAR you can then find as usual a "scaling parameter" in the second line and the three altered lattice vectors from which you can (after multiplication with the "scaling parameter") determine the (relaxed) lattice constants a, b, and c for the given volume. And of course you will also know the total energy then but even hydrostatic pressure printed in OUTCAR ("pressure = ..."). So you could not only determine an energy-volume but even a pressure-volume curve!

Only if you have no idea at all about the equilibrium lattice constants / volume you might consider a "pre-relaxation" run with ISIF=3 first before you begin to sample different volumes around the (approximate) equilibrium volume with ISIF=4 ...

Just take care that your plane-wave cutoff is high enough (in order to avoid too large Pullay stresses -- I usually recommend to approximately double the cutoff with respect to what is set by VASP for PREC=high/accurate) and that you perform enough relaxation steps to get all stress tensor components converged to at least better than order of 0.1 kBar. And also use a rather small EDIFF (use maybe values of about 1-9 to 1-6, in order to obtain well-converged wave functions which are important for numerically accurate forces and stresses). And overall try to be as accurate as possible (e.g. PREC=high or PREC=accurate is recommended, maybe also ADDGRID=.TRUE.).

This strategy works by the way for all types of (non-cubic) lattices, also monoclinic or triclinic ones (getting not only the three lattice constants but also the one or three angles). And even for tetragonal or hexagonal ones it can still be useful.

Many thanks Jürgen for the nice reply. I would come back to you after trying your method. Thanks once again for sharing your ideas.

I was trying something different. E.g., for a 'fixed volume', keeping b/a constant I was calculating E-V curve for a set of c/a values and same for keeping c/a constant using ISIF=0 i.e., instead of varying volume automatically by VASP,  I was changing the volume by-hand and then plot a 2-D contour plot for b/a vs c/a and find the optimal b/a and c/a for that particular volume and thus a,b,c. This is bit exhaustive calculations.

Well, that's a kind of brute-force alternative which also works, but as you said it's very exhaustive (if one uses just 7 "grid points" for each "degree of freedom" it would be 49 calculations for 7 volumes ...). Of course, 7 relaxation runs taking something like 100 steps (or more) will also need some time. But latest if you have internal structural parameters which require an  ionic (re-)relaxation in each calculation with a few ionic relaxation steps, I think the all-in-one relaxation jobs for ionic coordinates plus cell shape are at the end (much) more efficient (therefore I prefer this). But let's see; at least comparing the two results you may also get an impression about potential errors due to Pullay stresses (even for high cutoffs small errors of the order few 10^-4 Angstroem will occur in the lattice constants -- on the other hand fit procedures for finding the minimum of your total-energy surface may result in similar errors).

I am trying something : first of all,  i am writing my POSCAR file from experimental lattice parameter and  I allow all relaxation with ISIF=3, and then I get CONTCAR file. I copy my CONTCAR file as POSCAR file, after then I apply some specific strain while ISIF=4. That's how I can get volume-energy or strain-energy curve. Do I follow wrong way? 

Otherwise, If we you ISIF =4 how we can get a, b, c lattice parameters? 

If you share your comment it will be so helpful for me. 

Many Thanks.

Well as already Jürgen said in his previous comments that using ISIF 4 you get a volume-energy curve and also you get b/a and c/a ratios for the different volumes. After that you fit your data with some equation of state and get the optimum volume. Once you get the optimum volume you can then use your already volume vs b/a and volume vs c/a plot to get the actual b/a and c/a ratios for your optimal volume.

Please note that you have optimized b/a and c/a for each and every volume. You then plot those as a function of volume and then use those fitting equation for your optimal volume to get your optimal b/a and c/a. Once you know b/a and c/a and volume you do simple mathematics to calculate a.b and c.

I have no experience with monoclinic structures but if you look at Jürgen's first comment (more precisely last paragraph of first comment), there he mentioned that this ISIF=4 tag holds for any non-cubic e.g., including monoclinic and triclinic structures.

Well, in a sense, yes you've to provide data for lattice parameters. Let's say you know (very roughly) your a,b and c and thus V0 (you can also use ISIF 3 if you don't have any idea). Then take 5/7 volumes around V0(±) and use ISIF 4 to get the energy for each volume (yes, here you should provide a,b and c to get V, but you don't need the correct ratios. VASP will do that).