[lammps-users] Bead-spring DNA Model

LAMMPS LAMMPS Mailing List Mirror
1

Hi,
I am trying to simulate and find the persistence length of a bead spring model, which is not going so well. I am wondering if lamms is doing strange things, my energy is a bit off. Energy is nicely conserved, except it sometimes massively spikes. I am also getting fluctuation in my e_pair output. I am using a harmonic potential to connect beads, and an angle style harmonic to make the chain stiff. E_pair should always be 0 correct? For example, from log.lammps. I am not sure if this is a problem.:

33500000 300 0 84.061409 217.30362 -1.068686e-05
33600000 300 0 96.196347 229.43856 1.3172717e-05
33700000 300 0 85.88458 219.12679 -2.9232279e-05
33800000 300 -1.7134004 89.953673 221.48248 -6.9664614e-05
33900000 300 0 83.799358 217.04157 -8.1265138e-05
34000000 300 0 75.640836 208.88304 -1.9692338e-06
34100000 300 0 72.048203 205.29041 1.9413903e-05

I get a spike of -1.7 for some reason. Then it will be 0 for a while. If anyone has calculated the persistence length of a simple bead spring model, I would love to talk. I need to get this nailed down soon so I can continue working on the fun parts of my research. My lammps file is here:

2

You have a pair-style defined, so if you get a non-zero pair
energy, some pair(s) of atoms must be interacting. Have you
tried visualizing the dynamics? Maybe the chain(s) are doing
something you don't expect. You can also use commands
like compute pair/local and dump local to dump out all pairwise
forces, so it shouldn't be hard to find what is contributing to
the -1.7.

Steve

3

The nergy is in evdwl, which is VanderWaals. I am not quite sure what this does. I tried turning it off with pair_modify tail no, but it still happens. Some illumination on the topic? Off topic, lammps quits at about 1.4 bil timesteps as if nothing went wrong. Is there a reason?

Thanks,
Josh

4

Again, you have defined a pair_style lj/cut which will compute
Vdwl energy. If you get a non-zero value printed, that's where
it is coming from. There is no particular reason for the code
to stop at 1.4B or any other number of timesteps.

Steve