How to Explain the Changes in the isometric heat capacity Cv ?

Hi everyone;

Water is driving me crazy, I'm trying to obtain isentropic heat capacity Cv of water using the mesocite module, :

I used Shinoda forcefield, all parameters left at their default values unless otherwise stated.

I used 1 bead = 3 molecule of water (mass=54amu and radius=2,78A), then I have ran a geometry optimisation (Ultra-fine), then the Dynamic task, temperature=298 K,  and finally I ran the simulation in the NVT ensemble with 300ps, the timestep is 10fs.

To calculate the isometric heat capacity, Cv, I launched the fluctuation properties analysis of mesocite, of course I used the NVT ensemble here too, and I filtered the frames when the structure is equilibrating.

I made the same calculations with three different configurations: 50x50x50A, then 70x70x70A and finally 80x80x80A, and here are the results;

Configuration

50x50x50A

70x70x70A

80x80x80A

Cv (kcal/molK)

8,37756210

20,87436761

260,28542699

As you can see when the configuration is 70x70x70A, the result (20,8744 kcal/molK) is comparable with the experimental value of Cv (17,86 kcal/molK), however, normally (unless I'm mistaken), if I increase the dimensions of the configuration, the results should become more accurate which is not the case here!, and I don't understand why...

Strangely, when in the three configurations I have obtained the correct value of solubility parameter of water (~47 (MPa)^(1/2)), which mean my simulation parameters are good.

I  tried  larger simulation time, I also tried to relax the configuration in other ensembles first (NVE and NPT), but nothing really changes.

here are my files (without the dynamics)

thank you.

JARRAY