# Question about fix langevin: the energy cumulative energy of "sink" and "source" is not equal

Dear all
I’m trying to calculate the interfacial thermal conductivity between LiF and polymer(PEO) by NEMD method. The simulation model is showed as below. It contains group source, gourp sink, the main sample, two fixed layers, two vacuum layers at both ends. Source gourp and sink group have the same number of atoms which is 768.

I controlled the heat source temperature to 355K and the heat sink to 255K by using fix langevin.The specific settings about this part is illustrated below. And before that, I have used the npt ensemble to stabilize the system at 305K. .

But when I tried to calculate the heat flux using the data producted from “fix heat” and “fix cold”, I found they differ from each other a lot. The slope of energy with respect to time is showed below.

In general, their gap should be small. So what’s the reason for that big difference between the energy accumulated in source and sink ?
Any comment on this issue is highly appreciated.

255 K is -20°C. My first guess would be that your cold region is freezing the surrounding PEO and that you are recording the continuous removal of enthalpy of freezing at that interface.

Notice how the heat outflow there begins immediately with a consistent gradient, whereas the heat inflow from hot region appears to only start with some lag. This may be consistent with an initial enthalpic “pulse” from the onset of freezing.

I don’t have any other guesses (and it IS a guess that you will have to assess further by yourself unless you provide further data).

This is likely not something that somebody can just tell from looking at your input (if it was that easy, we would not need to do simulations in the first place), but here are a couple of questions that may help you look for explanations.

Does your simulation conserve energy without the hot/cold region thermostats?

Is there a converged and linear temperature gradient across the system?

Thank you for your reply! I check the temperature of PEO and the lowest point of PEO temperature is about 236K. Is that mean PEO is frozen? And I don’t know much about the phenomenon of “enthalpic pulse” . But I will check it out. Thanks again!

It’s impossible for me to tell from here whether your PEO is frozen, or if that is even the source of the problem. Remember that the real-life freezing point of a substance is not guaranteed to be the freezing point that the force field exhibits in a simulation.

Simply downloading and visualising a trajectory will tell you whether the polymers are locally frozen or not.

“Enthalpic pulse” is a made-up phrase. It essentially means that if you had a polymer equilibrated at a temperature at which it is fluid, and it is suddenly cooled to freezing, the enthalpy of freezing must go somewhere. And if the cold source thermostat is already removing heat at its maximum rate, then the rest of the PEO must heat up while the steady state is being reached… And during that time the hot thermostat will have to add less heat to the system.

I am purely guessing the most interesting possible explanation. Whether it is valid or not, I am not even sure, and you will have to double check your own simulation results.

Dear Mr. Srtee:
I’m sorry for not getting back to you in a timely manner. Your idea is enlightening. I conducted the simulation again. In terms of the possible freezing phenomenon. I checked the trajectory as you said. It seems the polymer is not frozen. The gif below is a screenshot of the fix langevin process from 1.0ns to 2.0ns.

``````reset_timestep  0
timestep        \${TS}
thermo          1000
thermo_style    custom step temp vol press pe ke etotal

fix             e_output all ave/time 100 10 1000 v_energy file e_output.dat
fix             nvt all nvt temp 255 255 \${tdamp}
run             500000
unfix           nvt

write_restart   nvt.restart
write_data      nvt.data

# define regions
variable        z_l equal zlo
variable        z_h equal zhi
variable        z1  equal \${z_l}
variable        z2  equal \${z_h}
variable        c equal \${z1}+8
variable        h equal \${z2}-8
variable        fix_left equal \${z1}+2
variable        fix_right equal \${z2}-2

region          lfixed block INF INF INF INF \${z1} \${fix_left} units box
region          rfixed block INF INF INF INF \${fix_right} \${z2} units box
region          fixed union 2 lfixed rfixed

region          cold   block INF INF INF INF \${fix_left} \${c} units box
region          heat   block INF INF INF INF \${h} \${fix_right} units box
region          main   block INF INF INF INF \${fix_left} \${fix_right} units box

group           fixed  region  fixed
group           heat   region  heat
group           cold   region  cold
group           main   region  main
compute         KE      main ke/atom
variable        KB      equal 8.617333262145e-5*23.0609 # convert to kcal/mol
variable        TEMP    atom c_KE/\${KB}/1.5

compute         T_heat  heat temp/region heat
compute         T_cold  cold temp/region cold
compute         T_main  main temp/region main
compute         T_fixed fixed temp/region fixed

# fix boundary
velocity        fixed set      0.0 0.0 0.0
fix             fixed fixed   nve/noforce

change_box      all z delta -10 10 units box
write_data      vaccum.data

# fix nvt
fix             nvt main nvt temp \${T} \${T} \${tdamp}
run             500000
unfix           nvt

write_restart   nvt_main.restart
write_data      nvt_main.data

# fix langevin to creat temperature gradient
fix             nve     main nve
fix             HEAT    heat langevin \${Thi} \${Thi} \${tdamp} 766144 tally yes
fix             COLD    cold langevin \${Tlo} \${Tlo} \${tdamp} 545812 tally yes
run             2000000 # run 2 ns for steady structure
unfix           e_output
write_restart   heat.restart
write_data      heat.data

# collect data
reset_timestep  0
fix             HEAT    heat langevin \${Thi} \${Thi} \${tdamp} 766144 tally yes
fix             COLD    cold langevin \${Tlo} \${Tlo} \${tdamp} 545812 tally yes

fix             1 all ave/time 10 100 1000 c_T_heat c_T_cold c_T_main c_T_fixed file temp.dat
fix             2 all ave/time 10 100 1000 f_HEAT f_COLD file heat.dat

compute         T1 main chunk/atom bin/1d z lower 0.02 units reduced
fix             T1 main ave/chunk 10 50000 500000 T1 v_TEMP temp file 50chunk_500000.dat
compute         T2 main chunk/atom bin/1d z lower 0.01 units reduced
fix             T2 main ave/chunk 10 50000 500000 T2 v_TEMP temp file 100chunk_500000.dat

dump            heat all custom 50000 heat.dump id mol type x y z
run             1000000

write_restart   product.restart
write_data      product.data
``````

The total energy of system is illustrated below:

During 1.0ns and 3.0ns, I used langevin to first get a stable structure before I collect the energy of sink and source. But I still get the same results as yesterday. Therefore, How do I determine the heat flow through the interface if the heat sink and the heat source energy are still not conserved?

Thanks for providing more information. The script portion you have pasted is still missing important data (such as the values of the timestep and equilibration temperature).

The total energy graph is the most self-explanatory piece of additional information. After the initial equilibration for the first 0.5ns, the total energy of the system jumps to a new average value from 0.5 to 1.0 ns, and then has the same average but larger fluctuations from 1.0 to 1.5ns, after which there is another jump in energy and still larger, slower energy fluctuations.

I am not sure how those correspond to either your changes in simulation parameters, or to some change in the condition of the PEO at particular times or temperatures. Beyond that, I don’t have more specific information or insight. For that, you should refer to the scientific literature for papers about similar studies.

The timestep is set as 1fs, and equilibration temperature is set as 305K. As the total energy graph illustrated, for the first 0.5ns, I relax the “all group” using the nvt ensemble at 305K, and for the second 0.5ns, I do the exactly same for the “main group”, where the main group is the part of the all group after the fixed layer is removed. This main explain the first jump at 0.5ns. And Starting at 1.0ns, I started taking the langevin hot bath as shown in the pasted script. The second jump at 1.0ns maybe is a sign that system is re-equilibrating. The total simulation script is pasted as below:

``````read_restart    equilibrium2.restart

special_bonds lj/coul   0.0 0.0 0.5

kspace_style    pppm 0.0001
# oplsaa coeffs for PEO and LiTFSI provided by moltemplate
pair_style hybrid lj/cut/coul/long 10.0 10.0  buck/coul/long 10.0 lj/cut 10.0

pair_coeff 1 1 lj/cut/coul/long 0.14 2.9        # -O-
pair_coeff 2 2 lj/cut/coul/long 0.066 3.5       # CH3OR C
pair_coeff 3 3 lj/cut/coul/long 0.066 3.5       # -CH2OR O
pair_coeff 4 4 lj/cut/coul/long 0.03 2.5        # H-COR H
#pair_coeff 5 5 0.71 3.05 # F-
pair_coeff 6 6 lj/cut/coul/long 0.0005 2.87     # LiTFSI Li
#pair_coeff 7 7 0.0005 4.07 #Li+
pair_coeff 8 8 lj/cut/coul/long 0.17 3.25       # TFSI N
pair_coeff 9 9 lj/cut/coul/long 0.25 3.55       # TFSI S
pair_coeff 10 10 lj/cut/coul/long 0.21 2.96     # TFSI O
pair_coeff 11 11 lj/cut/coul/long 0.066 3.5     # TFSI C
pair_coeff 12 12 lj/cut/coul/long 0.053 2.95    # TFSI F

pair_modify mix geometric

# buckingham coeff for LiF
pair_coeff 7 7 buck/coul/long 2771.9202  0.28571 3.689744 # li+ -- li+
pair_coeff 5 5 buck/coul/long 11415.146  0.30395 514.2581 # f- --- f-
pair_coeff 7 5 buck/coul/long 7125.8181  0.28409 41.27901 # li+ -- f-

# lj/cut coeff for nonbonded interactions
pair_coeff 7 1 lj/cut 0.038729833 2.609360068   # Li-O
pair_coeff 7 2 lj/cut 0.051234754 2.737075322   # Li-C
pair_coeff 7 3 lj/cut 0.038729833 2.609360068   # Li-O
pair_coeff 7 4 lj/cut 0.033166248 2.369453297   # Li-H
pair_coeff 7 6 lj/cut 0.025       2.183592758   # Li-Li
pair_coeff 7 8 lj/cut 0.041533119 2.668336103   # Li-N
pair_coeff 7 9 lj/cut 0.082764727 2.801700833   # Li-S
pair_coeff 7 10 lj/cut 0.038729833 2.609360068  # Li-O
pair_coeff 7 11 lj/cut 0.051234754 2.737075322  # Li-C
pair_coeff 7 12 lj/cut 0.035355339 2.558161645  # Li-F

pair_coeff 5 1 lj/cut 0.054772256 3.056964179   # F-O
pair_coeff 5 2 lj/cut 0.072456884 3.206587439   # F-C
pair_coeff 5 3 lj/cut 0.054772256 3.056964179   # F-O
pair_coeff 5 4 lj/cut 0.046904158 2.775904309   # F-H
pair_coeff 5 6 lj/cut 0.035355339 2.558161645   # F-Li
pair_coeff 5 8 lj/cut 0.058736701 3.12605684    # F-N
pair_coeff 5 9 lj/cut 0.117046999 3.282298673   # F-S
pair_coeff 5 10 lj/cut 0.054772256 3.056964179  # F-O
pair_coeff 5 11 lj/cut 0.072456884 3.206587439  # F-C
pair_coeff 5 12 lj/cut 0.05        2.996983288  # F-F

bond_coeff 1 441.92 1.323
bond_coeff 2 268.0 1.529
bond_coeff 3 320.0 1.41
bond_coeff 4 340.0 1.09
bond_coeff 5 233.03 1.818
bond_coeff 6 637.07 1.437
bond_coeff 7 374.88 1.570
angle_coeff 1 93.33 107.1
angle_coeff 2 50.0 109.5
angle_coeff 3 33.0 107.8
angle_coeff 4 35.0 109.5
angle_coeff 5 37.5 110.7
angle_coeff 6 82.93 111.7
angle_coeff 7 60.0 109.5
angle_coeff 8 80.19 125.6
angle_coeff 9 103.97 102.6
angle_coeff 10 115.80 118.5
angle_coeff 11 91.30 103.0
angle_coeff 12 94.29 113.6
dihedral_coeff 1 -0.55 0.0 0.0 0.0
dihedral_coeff 2 0.0 0.0 0.468 0.0
dihedral_coeff 3 0.0 0.0 0.3 0.0
dihedral_coeff 4 0.65 -0.25 0.67 0.0
dihedral_coeff 5 0.0 0.0 0.76 0.0
dihedral_coeff 6 0.000 0.000 -0.004 0.0
dihedral_coeff 7 7.833 -2.490 -0.764 0.0
dihedral_coeff 8 0.0 0.0 0.347 0.0
dihedral_coeff 9 0.0 0.0 0.316 0.0
variable        T       equal 255
variable        Thi     equal 305
variable        Tlo     equal 205
variable        TS      equal 1
variable        tdamp   equal 100
variable        pdamp   equal 1000
variable        energy  equal etotal

reset_timestep  0
timestep        \${TS}
thermo          1000
thermo_style    custom step temp vol press pe ke etotal

fix             e_output all ave/time 100 10 1000 v_energy file e_output.dat
fix             nvt all nvt temp 255 255 \${tdamp}
run             500000
unfix           nvt

write_restart   nvt.restart
write_data      nvt.data

# define regions
variable        z_l equal zlo
variable        z_h equal zhi
variable        z1  equal \${z_l}
variable        z2  equal \${z_h}
variable        c equal \${z1}+8
variable        h equal \${z2}-8
variable        fix_left equal \${z1}+2
variable        fix_right equal \${z2}-2

region          lfixed block INF INF INF INF \${z1} \${fix_left} units box
region          rfixed block INF INF INF INF \${fix_right} \${z2} units box
region          fixed union 2 lfixed rfixed

region          cold   block INF INF INF INF \${fix_left} \${c} units box
region          heat   block INF INF INF INF \${h} \${fix_right} units box
region          main   block INF INF INF INF \${fix_left} \${fix_right} units box

group           fixed  region  fixed
group           heat   region  heat
group           cold   region  cold
group           main   region  main

group           LiF type 5 7
group           PEO type 1 2 3 4
group           PEO_LiTFSI type 1 2 3 4 6 8 9 10 11 12
# compute
compute         KE      main ke/atom
variable        KB      equal 8.617333262145e-5*23.0609 # convert to kcal/mol
variable        TEMP    atom c_KE/\${KB}/1.5

compute         T_heat  heat temp/region heat
compute         T_cold  cold temp/region cold
compute         T_main  main temp/region main
compute         T_fixed fixed temp/region fixed

# fix boundary
velocity        fixed set      0.0 0.0 0.0
fix             fixed fixed   nve/noforce

change_box      all z delta -10 10 units box
write_data      vaccum.data

# fix nvt
fix             nvt main nvt temp \${T} \${T} \${tdamp}
run             500000
unfix           nvt
write_restart   nvt_main.restart
write_data      nvt_main.data

# fix langevin to creat temperature gradient
fix             nve     main nve
fix             HEAT    heat langevin \${Thi} \${Thi} \${tdamp} 766144 tally yes
fix             COLD    cold langevin \${Tlo} \${Tlo} \${tdamp} 545812 tally yes
run             2000000 # run 2 ns for steady structure

unfix           e_output
write_restart   heat.restart
write_data      heat.data

# collect data
reset_timestep  0
fix             HEAT    heat langevin \${Thi} \${Thi} \${tdamp} 766144 tally yes
fix             COLD    cold langevin \${Tlo} \${Tlo} \${tdamp} 545812 tally yes

fix             1 all ave/time 10 100 1000 c_T_heat c_T_cold c_T_main c_T_fixed file temp.dat
fix             2 all ave/time 10 100 1000 f_HEAT f_COLD file heat.dat

compute         T1 main chunk/atom bin/1d z lower 0.02 units reduced
fix             T1 main ave/chunk 10 50000 500000 T1 v_TEMP temp file 50chunk_500000.dat
compute         T2 main chunk/atom bin/1d z lower 0.01 units reduced
fix             T2 main ave/chunk 10 50000 500000 T2 v_TEMP temp file 100chunk_500000.dat

dump            heat all custom 50000 heat.dump id mol type x y z
run             1000000

write_restart   product.restart
write_data      product.data

``````

Please feel free if you need more information. Now, I will fix langevin again at a timestep of 0.5fs to check if that’s the key to this issue.