Setforce and lack of non-bonded interaction (?)

Hey there,

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 ?

Thanks a lot in advance!
final_pos_21.dat (6.6 MB)
in.input (10.1 KB)
merged_pos_polymer_oncee.dat (4.1 MB)

There are multiple issues here to keep in mind:

  • using fix setforce 0.0 0.0 0.0 is not sufficient to immobilize atoms, you also need to set their velocities to zero. in general, it is much easier to simply not include those atoms in the time integration group (i.e. don’t use group “all” for fix nvt)
  • there are some special problems with applying fix npt and immobilized atoms. generally, when you change the box dimensions all atom positions are scaled accordingly, thus your immobilized atoms are not fully immobile. Also their temperature (or lack thereof) and pressure contributions are included in the pressure value that controls the volume change. in general, it is therefor usually better to not consider using fix npt, but rather use fix nvt and manually adjust the box dimension occasionally and also only compute force, energy, pressure contributions from mobile atoms and exclude the intra-group interactions of the immobilized atoms. If you prefer to use fix npt regardless, then it is strongly advised to not immobilize your MOF phase, but make it rigid with fix npt rigid (and allow only translation in z-direction as well as exclude the intra-group interactions to avoid bogus contribution to pressure) and then apply fix nvt to the remaining atoms.

No, but the pressure contributions from the intra-MOF forces can be bogus.

2 Likes

Thanks Axel, I managed to fix it with the possible issues you highlighted in your second comment :))