Possible way to counter displacement due to volume-rescalling in NPT

Hey everyone,

I am currently simulating a system that consists of (1) a porous material (PM) and (2) a polymer whose atoms are initially surrounding the PM. Despite aware that in order to get a “realistic final state” within the context of the force-fields I am using for my classical dynamics, let’s suppose that I want to forbid the polymer of getting inside the PM artificially. I managed to achieve this successfully with ‘region sphere’ + ‘fix wall/region’.

But let’s suppose now that I am going to start the dynamics and use nvt/rigid to move the PM, with intra-PM interactions turned off by neigh_modify + delete_bonds, and npt for the polymer. Naturally, the volume of my system will fluctuate and the position of all atoms, including the PM, will be rescalled accordingly. In my case, so will the center of my spherical region since I set to always be centered in the center of mass of the PM.

Within this rescaling, it is quite possible to, at some point, happen to have atoms of the polymer falling inside the sphere (which is a “forbidden” region since I am considering the outside of the sphere in my ‘region’ command). Does anyone have a suggestion of some command/strategy within LAMMPS that I could do to counter the problem that would allow me to directly use NPT? (I am aware of the possibility of using NVT and rescaling the box manually to ultimately attain a given desired pressure, but would there be a way to actually use NPT in this scenario?)


Ah, btw I forgot to share the part of the code in which I define the region and everything (in case it is useful):


compute compm pm com

thermo_style custom step temp etotal press vol c_compm[1] c_compm[2] c_compm[3]
thermo 10

fix 1 all nvt temp 300 300 $(100*dt)
run 0
unfix 1

variable xcom equal c_compm[1]
variable ycom equal c_compm[2]
variable zcom equal c_compm[3]

region ball sphere v_xcom v_ycom v_zcom 14.0 side out
fix reflection polymer wall/region ball harmonic 10.0 0.0 2.5


Hi @ceciliaalvares,

Without the specific options of your commands it is quite hard to answer you question.
When you say:

Do you mean that the NVT and NPT integrations are separated but that the dilation only applies to the polymer coordinates? Is there a reason for that?

Else you can use a variable to compute the radius of the required sphere (largest distance between PM atoms and COM of the PM for example) and it will be updated every timestep.

1 Like

Hey Germain,

Without the specific options of your commands it is quite hard to answer you question.

Yeah, sorry, I posted it a little late in a second comment above (at least part of the code). Indeed, the dynamics would be something as illustrated below (just for the sake of illustration).

fix 1 polymer npt temp 300 300 100 iso 3000 3000 1000
fix 1 PM nvt/rigid temp 300 300 100

Is there a reason for that?

Somewhat, yes (there is this published paper on how to computationally get rightful structures for polymers, which are amorphous in general, and it includes some NPT cycles which I would like not to apply for the PM (which is why I am using fix rigid to avoid some weird deformation on it)).
The fact that the center of mass of the PM would also be rescaled (and I understand it should) was messing up a bit my “made-up” strategy to counter the problem of the polymer getting inside the pore and I was out of ideas for an alternative.

In any case actually your idea is quite good, I had not thought of trying it :grimacing:
I will try it, thanks a lot! :slight_smile:

Else you can use a variable to compute the radius of the required sphere (largest distance between PM atoms and COM of the PM for example) and it will be updated every timestep.

Just a question on this, if I may. In my case, to put your idea in practice, I would need to, in each timestep, perform a calculation of distance between the “new center of mass” of my PM and every single existing atom composing the polymer in my system. And then I would need to find, within this big array of distances, what is the polymer atom sitting closest to the center of mass and set up the value of radius of my sphere to be that distance. This is what you mean, right? Cause this could work.

But would you mind referring me to a particular “section” of the manual which I can read to learn how to do this? Cause I dont know how to, in LAMMPS, calculate a distance between a point in space and atoms and also find which is the smallest distance within an array? Just the name of some command lines would help a lot for me to search the information about. Or would this be sort of python interface that I would need to learn how to do?

Couple of comments on your answers:

With these settings the scaling of coordinates happens to all the atoms of your box so it is unlikely a scaling step make you polymer atoms go inside your sphere if the radius also scales. I was curious to know if you used the dilate keyword to exclude atoms from the scaling operation. See the end section of fix rigid on hybrid simulations and fix npt pages to see why this is important.

This is something you can do using variables and computes in LAMMPS. If I am not mistaking, an atom style variable will compute a quantity for each atoms in your system if you do something like:

variable dist atom sqrt((x-c_compm[1])^2+(y-c_compm[2])^2+(z-c_compm[3])^2)

and then pass it to a compute reduce command that would find you the max value. The group used in the compute can restrict the search to only polymer atoms. You can the pass it to an equal style variable to be evaluated by your region command.

This is the simplest solution that comes to my mind. Python scripting is also a solution be might be an overkill in this case.

1 Like

I wonder if you would be better off just not integrating the porous material at all, using something like

fix integrator polymer npt ... dilate polymer

After all, from a physical point of view, the quantities of interest are the relative displacement between polymer and porous material, and you should generate the same dynamics simply keeping the porous material immobile (not just rigid) and letting the polymer come to it.

1 Like

To be very honest, what I had tried at first when setting up the input script was the setup below:

fix		 3 pm rigid single
fix		 4 polymer npt temp 300 300 $(100*dt) aniso 986.923 986.923 $(1000*dt)

Meaning that I had not assigned any equation of motion to my PM, but allowed it to move as a whole due to the rescalling of coordinates as the box shrinks/expands. Then, at some point during the dynamics, I happen to have an atom of the polymer inside the spherical region I had previously set up (which gives an error). From my judgment, this is coming from the rescalling the coordinates within the NPT (supported by the notion that if I try to run NVT I have no such error).
And in principle, in my understanding, I will continue to have this problem in my case, since it is intrinsic to the fact that I will have some rescalling of atomic positions, regardless if I am applying it only for the polymer or not. This is because, regardless, at some point I am always vulnerable to have atoms of the polymer lying within the region due to the rescalling of positions anyways. So this was why I was looking at fancier ways to counter the problem, particularly applied to my case. I havent tried the idea @Germain had yet but I think it will work!

@srtee I tried your idea, by changing my script to do be:

fix		 3 pm rigid single
fix		 4 polymer npt temp 300 300 $(100*dt) aniso 986.923 986.923 $(1000*dt) dilate polymer

but unfortunately I still seem to reach a point within the simulation where polymer atoms find themselves inside the sphere :frowning:
So indeed I think my problem is intrinsic to the rescalling. Thanks anyways for the suggestion!

I am now actually trying to assign the PM some dynamics by means of nvt/rigid as I had initially mentioned in my initial post and so far the simulation is still running (although in my mind, the suceptibility to the problem is still there and I believe the simulation will collapse at some point unless I am extremely lucky). If it also fails, then I will try the fancy idea Germain had.

PS: for what is worth it, in the cases for which I have no dynamics for my PM (in which case the velocity of all its atoms are 0) I also have a fix_modify that calls a compute command in which I define temperature calculations and etc to make my the “tempreature of the polymer phase” be indeed 300K (I am just not sharing every single detail of the script to save some time and energy of ppl reading).