Rotation of molecules after using fix rigid

[RezaNamakian]

Hello LAMMPS community,

I’ve got a situation I could use your help with. My system involves organic-inorganic materials that have been equilibrated for 100 ns at 125 K, using ‘fix npt’ to control all pressure components in a triclinic simulation cell.

Now, I’m trying to tackle a specific issue: preventing the rotation of a group of methylammonium molecules (MA). To do this, I’ve been experimenting with the ‘fix rigid’ command to turn off torques. Here’s a snippet of the script I’ve been working on:

LAMMPS (28 Mar 2023 - Development)

read_restart 125.restart
special_bonds amber
pair_style hybrid buck/coul/long 17.0 17.0 lj/cut/coul/long 17.0 17.0
pair_modify shift yes
bond_style harmonic
angle_style harmonic
dihedral_style fourier
kspace_style pppm 1.0e-7

group MA type 3:10
group Thermalized subtract all MA

compute TempCOM Thermalized temp/com
compute Press all pressure TempCOM
compute Rotation MA temp/rotate

include force_field.data
timestep 0.1
run_style verlet
neigh_modify every 1 delay 0 check yes cluster yes one 10000 exclude molecule/intra MA
delete_bonds MA multi any

reset_timestep 0 time 0
velocity MA set 0 0 0 units box
balance 1.0 shift xyz 10 1.0 weight time 0.8 weight neigh 0.8
run 0 post no

fix Frozen MA rigid/nvt molecule temp 125 150 $(100.0*dt) torque * off off off tparam 10 100 5
fix_modify Frozen temp Rotation

velocity MA zero angular rigid Frozen

fix Mom Thermalized momentum 1000 linear 1 1 1

fix Equil Thermalized npt temp 125 150 (100.0*dt) x 0 0 (1000.0dt) y 0 0 (1000.0*dt) z 0 0 (1000.0dt) xy 0 0 (1000.0*dt) xz 0 0 (1000.0dt) yz 0 0 $(1000.0dt) nreset 1000 couple none tchain 5 tloop 100 pchain 5 ploop 100 nreset 1000
fix_modify Equil temp TempCOM press Press

run 5000000

However, I’ve encountered a couple of issues:

1. When I use ‘exclude molecule/intra MA’ and ‘delete_bonds MA multi any,’ I run into a problem with non-numeric pressure values.

2. Even when I don’t use the options mentioned above, I’m noticing that the molecules start rotating after 60000 timesteps. To give you a better idea, I’ve attached two snapshots illustrating the rotation.

I’d really appreciate any advice or suggestions you might have on these issues.

Thanks,
Reza.

[srtee]

There are many problems with your script:

1. Are you sure your material actually assumes triclinic geometry in real life? Pressure fluctuations are much smaller for nanoblocks of the system than for the overall macroscopic system. You should always use an isotropic barostat unless proven otherwise (for example in lipid bilayers where the compressibility parallel to and transverse to the bilayer is clearly different).
2. You’re applying fix npt to a subset of your system (already tricky) and then you don’t use dilate to exclude your rigid molecules. As the documentation says:

Regardless of what atoms are in the fix group (the only atoms which are time integrated), a global pressure or stress tensor is computed for all atoms. Similarly, when the size of the simulation box is changed, all atoms are re-scaled to new positions, unless the keyword dilate is specified with a dilate-group-ID for a group that represents a subset of the atoms.

It’s not hard to imagine how, together with a fully anisotropic barostat, this results in rotation of your methylammonium molecules.

3. You have coupled your methylammonium molecules to a temp/rotate thermometer. This thermometer excludes the centre of mass translations and rotations of the molecules as a single group, rather than per molecule. Please be sure that’s what you actually want.
4. You have coupled the rest of your system to a temp/com thermometer, which is fine, but then you’re also zeroing its momentum every 1000 steps. Firstly this is overkill (if you’re already zeroing the momentum, why not just use the regular compute temp?) More importantly, by resetting the momentum of part of the system, you are introducing spurious and unphysical impulse forces. (More precisely, the total linear momentum of your system should be conserved, so your thermalized group and methylammonium group should develop some random mutual nonzero momentum – whenever you zero the momentum of the thermalized group, you are leaving the methylammonium group momentum untouched and thus “kicking” the methylammonium molecules relative to the thermalized group.)
5. Setting the velocity of all methylammonium particles to zero makes that group start from 0K temperature (and makes your next velocity command redundant). You may get slightly better results if you start the methylammonium molecules from a properly thermalized initial velocity (using the scale keyword if needed).

[RezaNamakian]

Hi Shern,

Thanks for getting back to me!

1. I’m confident about the material existing in that low-symmetry phase based on experimental observations.

2. I tried following the barostating recommendations from the documentation:

“If you wish to perform NPT or NPH dynamics (barostatting), you cannot use both fix npt and the NPT or NPH rigid styles. This is because there can only be one fix which monitors the global pressure and changes the simulation box dimensions…Use fix npt for the group of non-rigid particles. Use the dilate all option so that it will dilate the center-of-mass positions of the rigid bodies as well. Use one of the 4 NVE or 2 NVT rigid styles for the rigid bodies.”

Just to let you know, I need to stick with the triclinic structure for the entire material system. There’s an anticipated phase transition from triclinic to tetragonal in a specific temperature range, which lines up with what has seen in experiments. My goal is to pinpoint the effect of MA rotation during this transition.

3. Thanks for flagging that.
4. You’re on point – I do tend to build new scripts on top of old ones without cleaning up properly.
5. Your explanation totally make sense.

Even after making these changes, I’m still struggling with the same issue. I have my doubts about fix npt, like you mentioned, but I’m still scratching my head trying to figure out another way around this issue.

Reza.

[srtee]

Again, as the documentation states, fix npt scales all atoms outside the group specified by dilate. Rescaling all atom positions in the molecule is not compatible with purely translational motion of your rigid molecules.

Any recommendations in the LAMMPS manual are very general and may not apply to very specialised or unusual simulations.

But why not just use fix nvt and measure the pressure? That way you get rid of a big source of conceptual uncertainty. You can still manipulate the pressure by running simulations at different box sizes and comparing them.

[RezaNamakian]

Just a quick update: I ran the simulation again at 125 K without changing the temperature. This time, I switched from using NPT to NVT and used a really small timestep of 0.01 fs. Despite these adjustments, I’m still noticing the MA molecules starting to rotate after a few timesteps. It seems like there might be some issues with how the “fix rigid” is being integrated.

I’m thinking about different ways to deal with the rotational motion without relying on “fix rigid.” If you’ve got any suggestions or ideas on how to tackle this, I’d really appreciate hearing them.

[RezaNamakian]

I managed to get it working using a mix of “fix rigid/nve” and Langevin thermostating. Surprisingly, the simulation seems to be unaffected by “fix npt,” with a timestep of 1.0 and even “fix momentum” in place. Just one quick question: would it be a good idea to use the Nose-Hoover thermostat and barostat for the non-rigid part, while sticking to nve+ Langevin thermostating for the rigid part?

[srtee]

Good job identifying fix rigid/nve as a solution. I wouldn’t bother with a separate Langevin thermostat – if you are thermostatting the rest of your system at the desired temperature then your rigidified molecules will thermalize to the correct temperature at equilibrium.