I am currently trying to simulate a polymer/MOF interface by means of classical MD. I have a single data file already containing information concerning position of the atoms belonging to both phases. The initial goal is to equilibrate the polymer in the presence of the MOF. This is to be achieved by doing a 21 NVT/NPzT steps cycle which cherishes high temperature and pressure conditions which I do not want to apply to the MOF phase simultaneously (its structure is not supposed to be stable on these conditions). To do so, I am setting up a set of commands line that resembles what is briefly shown below (the NVT/NVT/NPzT conditions are simply roughly repeated again for 7 more times with different pressure conditions):
(…) compute myTemp polymer temp compute myPress all pressure myTemp fix 2 mof setforce 0.0 0.0 0.0 fix 1 all nvt temp 600 600 100 fix_modify 1 temp myTemp run 50000 unfix 1 fix 1 all nvt temp 300 300 100 fix_modify 1 temp myTemp run 50000 unfix 1 fix 1 all npt temp 300 300 100 z 986.923 986.923 500 fix_modify 1 temp myTemp press myPress run 50000 unfix 1 (…)
From my understanding, with this set up, I am freezing the atoms of the MOF phase (group-ID: MOF) whilst at the same time considering their non-bonded interaction with the polymer’s atoms to exist throughout the dynamics since I run NVT and NPT in group “all”. I am also addressing the “partial temperature” issue by redefining a “compute temperature” that consider only atoms of the polymer phase (group-ID: polymer). My issue is because for some reason even though during the cycle I apply significantly high pressures, I end up getting a result after equilibration where the two phases are not close to one another, in the sense that there is a huge void in the z direction as seen in the data file output by the end of the simulation. I am attaching here the input script, the initial microstate and the final microstate. Is the setforce somehow preventing the existance of non-bonded potential between atoms of the polymer and the MOF ?