barostat simulation with/without rigid body

Hi all,

I’m doing a simulation in which I there is a wall, solvent and a rigid body (picture is attached). I have added the rigid body (fullerene) to the equilibrated wall and solvent system and I want to equilibrate this new system with keeping fullerene stationary. After equilibration I want to let the fullerene loose to do some free energy calculations on that. I’ve tried two methods for equilibrating the system:

  1. defining fullerene as a rigid body and do neigh_modify exclude to turn off interactions within that. I also turn off forces in rigid option to keep it stationary.

  2. I simply excluded fullerene from time integration without defining a rigid body (However, this is not what I want after equilibration because I want to let it loose and this option doesn’t work for that).

The problem that I have is that looking at the first few thousands time steps of the first scenario I get negative pressure around (-600 atm) where roughly +2400 is the kinetic part contribution and -3000 is the virial contribution. However, for the second scenario I get almost the same contribution for the kinetic part but the virial is increased to something around -1600 which causes a large positive pressure around 800. (my set pressure is 1atm in z direction).

Is this difference comes from turned off forces in rigid that don’t contribute to virial but they are counted in the second scenario?

I would appreciate it if you can help me with this.

Cheers,
Kasra.

nph.jpg

Hi all,

I'm doing a simulation in which I there is a wall, solvent and a rigid body
(picture is attached). I have added the rigid body (fullerene) to the
equilibrated wall and solvent system and I want to equilibrate this new
system with keeping fullerene stationary. After equilibration I want to let
the fullerene loose to do some free energy calculations on that. I've tried
two methods for equilibrating the system:

1. defining fullerene as a rigid body and do neigh_modify exclude to turn
off interactions within that. I also turn off forces in rigid option to keep
it stationary.

2. I simply excluded fullerene from time integration without defining a
rigid body (However, this is not what I want after equilibration because I
want to let it loose and this option doesn't work for that).

sure it does, you just have to add the atoms of the fullerene to the
group that you time integrate. but you also have to relax them.

The problem that I have is that looking at the first few thousands time
steps of the first scenario I get negative pressure around (-600 atm) where
roughly +2400 is the kinetic part contribution and -3000 is the virial
contribution. However, for the second scenario I get almost the same
contribution for the kinetic part but the virial is increased to something
around -1600 which causes a large positive pressure around 800. (my set
pressure is 1atm in z direction).

Is this difference comes from turned off forces in rigid that don't
contribute to virial but they are counted in the second scenario?

yes. you are comparing apples and oranges. you have forces (and thus
stress) from atoms in your fullerene, that you didn't relax properly
before starting the run.

axel.

Thank you Axel, that was so helpful.
Now if I want to treat fullerene as rigid (as some papers do): I do “neigh_exclude group c60 c60” and “delete_bonds c60 multi remove” which keeps the fullerene structure intact, right? if that’s correct then the fullerene is just an object having pairwise interaction with its surroundings. In this case: is that right if I claim that:

The only difference between a) using “fix rigid” on fullerene with turned off forces and torques with the case of b) not time integrate fullerene is that the first approach includes only the pairwise forces from fullerene on the surrounding atoms in the virial but the latter also includes pairwise forces from surrounding atoms on the fullerene?

if this statement is wrong, would you tell me what I’m missing here?

Cheers,
Kasra.

Thank you Axel, that was so helpful.
Now if I want to treat fullerene as rigid (as some papers do): I do
"neigh_exclude group c60 c60" and "delete_bonds c60 multi remove" which
keeps the fullerene structure intact, right? if that's correct then the
fullerene is just an object having pairwise interaction with its
surroundings. In this case: is that right if I claim that:

The only difference between a) using "fix rigid" on fullerene with turned
off forces and torques with the case of b) not time integrate fullerene is
that the first approach includes only the pairwise forces from fullerene on
the surrounding atoms in the virial but the latter also includes pairwise
forces from surrounding atoms on the fullerene?

you can easily verify this statement by either reading the source code
or printing out the data that LAMMPS is computing. why guess? or
depend on somebody's say-so?

if this statement is wrong, would you tell me what I'm missing here?

there are two things that *i* don't get:
- why don't you relax the fullerene structure (bonds and all) so that
it has negligible contributions? you will need to do that when you
re-add it to the system anyway, so why not get it over with right
away?

- what is the point of doing an NPT run on your system anyway? it will
equilibrate just as well with fixed volume.

axel.