Thermostat the reactive system composed of atoms and molecules

Hello Everyone!

I want to apply a thermostat to my reactive system composed of Nitrogen molecules and Nitrogen atoms. I want the thermostat to operate only on the translational temperature (center of mass motion) and don’t want the thermostat to ‘mess’ with the internal energies of molecules which I think is happening currently in my simulation using ReaxFF (no matter what thermostat I select),.
What is the correct approach to achieve this?

Thanks a lot,
Gunjan

This is not as simple as you might think as it requires multiple steps to work in unison:

  • you need to have access to the ReaxFF bond order information to determine which atoms have bonds to what other atoms
  • you need to use that information to “cluster” the bonded atom information into “molecules”
  • for each molecule you need to separate the internal degrees of freedom from the center of mass motion. BTW: you didn’t mention rotation. technically, that would be a third category and different from the internal vibrations.
  • with that information you can compute the “internal” and the “translational” temperature (and rotational)
  • after being able to compute it, you may apply a thermostat to those, which would, of course likely be needed to programmed special unless you can make use of the “thermostat bias” functionality in LAMMPS

Most certainly, this will require a significant programming effort.

Thanks, Axel for your reply. Yes, you got that right, I was referring to both rotational and vibrational mode.
I understand the steps. Do you have any suggestions on how best to proceed with this? Are you suggesting interfacing say python with LAMMPS to program this thermostat manually or do you have any other suggestions?
Also, is it possible to tweak the “Tdamp” parameter in say “csvr thermostat” to apply very weak coupling to the heat bath?
“fix ID group-ID temp/csld Tstart Tstop Tdamp seed”
Will this work almost NVE-like for internal energies and still provide some temperature control for translational degrees?

No.

LAMMPS is written in C++. New features need to be implemented in C++, too. You need access to internal data of the ReaxFF pair style. That is not possible from Python currently unless you add an interface. But the computation you need to do is significant, so using Python would be to slow.

You are free to choose a suitable value. Obviously, the larger the time constant the weaker the coupling.

The temp/csld or temp/csvr or any other of those thermostats will work on individual atoms and monitor the temperature of the entire system. So “no”.