Constant Potential Energy

Hello,

I am trying to run a simulation that couples two small chains together via a spring. I need to output the total potential energy of the simulation (E_spring + E_system), but when I do, it seems to be nearly constant throughout the simulation regardless of the length of the spring. I have included my input script here:

############################################
############# General Settings #############
############################################

Set units to Lennard-Jones units

units lj

Specify atom type

atom_style hybrid charge bond

Create box

boundary p p p

############################################
############ Set Atom Parameters ###########
############################################

Read data file

read_data readfile_A_2_15-mer.in

Set parameters for FENE style bond

bond_style fene
bond_coeff * 30.0 2.0 1.0 1.0
special_bonds fene

Set parameters for electrostatic interactions

pair_style lj/cut/coul/debye 0.0002 2.5 7500.0
pair_coeff * * 0.256 1.0

Set parameters for angular interaction / bond stiffness

dielectric 1

Set parameters for near neighbor list construction

neigh_modify every 1 delay 0 check no

Set groups

group chain1 molecule 1
group chain2 molecule 2

############################################

Compute COMs and Set Output Variable

############################################

Compute COMs

compute com1 chain1 com
compute com2 chain2 com

Set variable for distance between COMs

variable R_coms equal “sqrt(((c_com2[1]-c_com1[1])^2)+((c_com2[2]-c_com1[2])^2)+((c_com2[3]-c_com1[3])^2))”

Set spring

fix coms_spring chain1 spring couple chain2 3 2.89 2.89 2.89 0
fix_modify coms_spring energy yes

############################################
######### General Output Settings ##########
############################################

Output thermodynamic properties

thermo_style custom step temp emol pe etotal press vol pxx pyy pzz pxy pxz pyz v_R_coms
thermo 1000

Output particle trajectories for analysis

dump dump_vals all custom 20000 dump_dyn.lammpstrj id type x y z ix iy iz q mol
dump_modify dump_vals sort id

Output restarts for archive

restart 500000 archive.*.restart
restart 10000 15-mer_d5.restart 15-mer_d5.restart

############################################
############## Set Ensemble ###############
############################################

Set fixes which act on all particles, set the ensemble

fix 1 all nve
fix 2 all langevin 2 2 0.1 699483

############################################
############# Timer Settings ##############
############################################

Set timer

Time format is HH:MM:SS

timer timeout 44:30:00 every 500

############################################
########### Energy Minimization ###########
############################################

Initial energy minimization

min_style quickmin
minimize 0.0001 0.0001 100000 10000000

############################################
####### Equilibration and Production #######
############################################

fix load_bal all balance 100 1.1 shift xyz 10 1.01

timestep 0.01
run 500000 upto

write_data A_2_15-mer_run.data
write_restart A_2_15-mer_run.restart

Any help figuring out what went wrong would be greatly appreciated! Thanks in advance!

Please note that your quoted input file is very hard to read since you didn’t quote it correctly and thus the typesetting engine of the forum interpreted comment and other special characters as formatting instructions. Please see the forum guidelines about some suggestions.

The first thing you should do, is to output the energy of the fix spring instance directly. It is computed as a global scalar, so you can just add f_coms_spring to the custom thermo style. Then you can see how large the energy contributed by the fix is compared to your total energy.