Abrupt fluctuations in temperature and high positive potential energy

LAMMPS
1

Dear Users,

I am simulating interfacial behaviour of liquids (fluids) on a solid surface. Initially, I ran my simulations keeping the surface rigid using setforce 0 0 0 for the surface atoms(groups), and initializing velocities only to the fluid group. During the simulation, I only considered the fluid -solid cross interaction due to vdw, as well as electrostatic interactions. I kept the interactions within the surface atoms themselves off using neigh_modify exclude molecule/intra. Simulations ran fine, with negative potential energy in NVT, in timestep of 2fs.

As soon as I change my surface to a flexible one from a rigid one using fix spring/self with a high spring constant, along with initializing velocities in all the group of atoms (surface+fluid), I get bond atoms missing error, as well as a very high potential energy. I also tried with timestep 1fs with and without keeping neigh_modify excluded for the surface, and for both the cases, the temperature of the whole system is fluctuating abruptly (for example, if I run simulation at 300K, temperature is touching 400K at some time steps), which distorts the actual dynamics of the system.

#Commands used in the first case#

neighbour 2.0
neigh_modify delay 0 every 1 check yes
neigh_modify exclude molecule/intra surface

timestep 2.0

velocity fluid create 12345 300

fix bal all balance 10000 1.05 shift z 100 1.05
fix sf surface setforce 0 0 0
fix 1 fluid nvt 300 300 100

#Commands used in second case#

neighbour 2.0
neigh_modify delay 0 every 1 check yes
neigh_modify exclude molecule/intra surface (off/on)

timestep 2.0 (1.0)

velocity all create 12345 300

fix bal all balance 10000 1.05 shift z 100 1.05
fix sf surface spring/self 1000
fix 1 all nvt 300 300 100

Is there any remedy to tackle that? Thanks in advance…

With Kind Regards,
Prateek.

2

With a spring constant of that magnitude your wall particles’ harmonic fundamental period will be at most a few femtoseconds (if not substantially less). This means a timestep of 1fs will be far too long for integrating the wall particles’ motion, leading to numerical instability and then to the symptoms you are describing.

Please note that even though this is almost certainly the immediate problem you are facing, there is no guarantee that once you fix this your simulation will be valid. Please make sure that you can either defend your simulation methodology to a reviewer or you are following closely some previously published methodology.

3

Thank you for your feedback on that!! But, what I can do now to fix the problem? Should I decrease the spring constant? Or, anything else?

In literature (a PRL article), the spring constant taken to be ~650-700 kcal/mol.A…I have taken 700 too…this showed the same errors…

With Kind Regards,
Prateek.

4

In cases like this, there rarely is a simple change of a command that can “fix” the problem.
The most common case is that the geometry of the surface is unstable and unrelaxed.

According to your description, you have excluded the force computations inside this part of your system and kept them immobile by setting all forces and velocities to zero. So if your geometry is not consistent with your force field for that solid, it will show the high-energy behavior that you describe. This is something that must be remedied, e.g. by correcting your initial geometry or using a per-equilibrated geometry or changing your simulation protocol where you alternate between having the liquid and the solid immobile (for minimization you only need to use fix setforce 0 0 0 on either group, for MD you may apply time integration only to that part) until either part is properly equilibrated for the desired final conditions.

5

Thank you so much Prof. Axel for your valuable feedback…I have tried with the pre-equilibrated system of liquid on a surface which was kept rigid using setforce 0 0 0 a priori. To observe whether the behaviour of the liquid at the solid-liquid interface changes in changing the surface from a rigid to a flexible one, I introduced fix spring/self on the surface, along with velocities on both solid and liquid. If I don’t use neigh-modify exclude command for the solid, then immediately I get bond-atom missing error with very high potential energy. If I consider neigh-modify exclude for solids, it runs fine at timestep 1fs, with negative potential energy, but the temperature fluctuates abruptly.

I haven’t tried with minimising the surface, but isn’t the structure of the surface prone to get distorted upon minimising it?

With Kind Regards,
Prateek.

6

Well, yes, a thermalized “surface” will in general show some disorder. Beyond that, it is impossible to know what’s “correct” for your system without knowing what your system is. At 300K, a wall of argon should gasify and a wall of gold should stay solid.

Having said that, a spring constant of 650 kcal/mol/square-angstrom is almost certainly too large for a meaningful simulation. An isolated gold atom in such a harmonic trap will vibrate with a fundamental period of about 27 femtoseconds. That is just about accessible with a timestep of 2 fs – but I am guessing you are using much lighter wall atoms whose oscillations will therefore be faster.

7

It will do so in any case, liquid or not and differently in either case, too. But if you don’t make any test whether the geometry/force field is adequate then how can you know it should run correctly without equilibration??

At any rate, this is no longer a discussion about LAMMPS but about how to properly design simulation workflows and how to setup and validate geometries. This is beyond the scope of this forum and something for a discussion with an adviser or tutor from which you are learning the craft of MD simulations. None of us here has a desire to take over that job and teaching these kind of skills over forum messages is extremely inefficient, specifically with somebody that doesn’t respond well to being given outlines and instead needs detailed instructions. In addition, as @srtee noted, we don’t know what your research is and what your specific system is, so how can anybody know what would have to be corrected.