Dear users,
I am trying to simulate the heat transport of graphene-like materials. Before simulation, both ends of the material are fixed using the fix setforce command. The system is first balanced in NVT and then switched to NVE. However, the total energy of the system increased linearly (in NVE). Simultaneously, the temperature is going up linearly. Although the rate of rise is relatively slow, the change is obvious after a long time of accumulation （such as 100ps). And with the simulation continues, the energy will keep going up.
I’ve done a lot of tests (including changing timestep and removing the fix setforce, as well as checking the potential parameter and so on), and I’m pretty sure this nonconservation is caused by the fix setforce command. Is there a reasonable way to solve this problem.
Thank you.
Your description is too vague and thus nobody can verify your claims and observations or give meaningful advice.
Hi void
Although you were asking a question about energy conservation, you also mentioned that you were using “fix setforce” to immobilize a group of atoms. Although I am not certain why your energy is increasing (see below), I do have a comment regarding how to immobilize atoms. Generally speaking, “fix setforce” is not a good way to immobilize atoms. Several alternate methods to freeze atoms were discussed in this thread:
https://sourceforge.net/p/lammps/mailman/message/36801562/
cheers
andrew
P.S. Incidentally, are you sure you zeroed out the forces? If you use fix setforce with non-zero values, then this would cause acceleration (causing the kinetic energy to increase, although not linearly). Also, you should always visualize your system to make sure it is doing what you think it is doing. (For example, during the simulation do the two ends that you immobilized remain fixed?)
Thank you for your reply. I feel I need to describe my problem in more detail. I mistakenly thought that this was a common problem, so I didn’t present my script and the result information before.
In fact, the non-conservation never occurs when only AIREBO potential is used. In this case, the interlayer interaction is described by 12-6 LJ. Here, the interlayer potential is replaced by the latest version of ILP because it describes the interlayer interaction more accurately. However, once ILP is used, nonconservation occurs. If the time required by the simulation is very short （such as 2ps）, this non-conservation is likely to be ignored. However, in my simulation of heat transport, the time required is up to several nanoseconds. The non-conservation cannot be ignored, which can also be seen from the linear increase of temperature. On the other hand, when I remove the fix setforce command, the nonconservation disappears. I’m confused what’s the root cause of the problem. (Lammp is the latest version).
#my script
#The box size is large enough and the model is located in the center of the box
#box=[-100 100; -100 100; -50 50]
#model=[-50 50; -50 50; -3.4 0]
units metal
boundary f f f
atom_style bond
neighbor 3.0 bin
read_data Bilayer_Graphene100x100A.data
pair_style hybrid/overlay rebo ilp/graphene/hbn 16.0 1
pair_coeff * * rebo CH.rebo C C
pair_coeff 1 2 none
pair_coeff * * ilp/graphene/hbn BNCH.ILP C C
mass 1 12.01
mass 2 12.01
timestep 0.0002
neigh_modify every 1 delay 0 check yes
#define boundary: left and right
region left block -50 -45 INF INF INF INF
group left region left
region right block 45 50 INF INF INF INF
group right region right
group boundary union left right
group mobile subtract all boundary
velocity boundary set 0.0 0.0 0.0
fix 1 boundary setforce 0.0 0.0 0.0
velocity mobile create 300 123456
fix NVT mobile nvt temp 300 300 0.02
thermo_style custom step temp press etotal
thermo 100
run 30000
unfix NVT
fix NVE mobile nve
run 100000
same as the input you posted before, this does not follow the force field settings from the example.
also, does your data file follow the molecule ID requirement of the inter-layer potential?
axel.
Thank you for your reply again. This problem still arises when I use the same force field parameter as in the examples. I changed these parameters because I wanted to shorten the test time, after all, the calculation speed of long range interaction is too slow. I found the reason after repeated tests. It’s not the force field settings but my script lacks the command “processors * * 1”. Although I don’t quite understand this command, and I rarely use it before. But it does ensure the conservation of energy.
best regards.