[lammps-users] Self-Interacting Chain

LAMMPS LAMMPS Mailing List Mirror
1

I’m trying to model a strand of DNA as a chain of pseudoatoms which interact with each other by harmonic bonds and angles and by pairwise volume exclusion and Coulomb potentials. Specifically, each of N pseudoatoms in the chain is bonded to its 1-2 neighbors (N-1 total bonds) and “angled” with its 1-2-3 neighbors (N-2 total angles). Also–and this is the part where I’m having trouble–each pseudoatom interacts with every >1-4 neighbor by means of a long neigh_modify exclude list and a custom “ss/debye/mod” potential which combines soft-sphere repulsion and a screened Coulomb force.

I can’t get the chain to self-interact via bonds and pairwise potentials at the same time; when I turn the bonds on, my Coulomb and VDW energies are always 0; only when I turn the bonds off can I get nonzero pairwise energies.

I created the ss/debye/mod potential by modifying the lj/cut/coul/debye potential. It’s sort of sloppy right now–the soft-sphere energies are labeled VDW energies in the console output, for example–but it does work with non-bonded pairs.

I’ve attached a simple data file and runfile for inspection (I can include the .cpp and .h files for the custom potential if necessary). Is there something simple I’m missing?

Thanks,
Scott Douglas
Georgia Institute of Technology
School of Physics

chainTestData (794 Bytes)

chainTestRun (2.29 KB)

2

Since you're adding your own custom potential, I have no idea.
I would debug this by turning on one option at a time (bonds, angles,
pairwise, exclusoins, etc) and verify that you are getting
the answer and non-zero terms you expect.

Steve

3

Thanks for your quick reply! Unfortunately, I still run into the same problem with the default coul/cut pairstyle. I’ve included a very simple data and runfile; two charged atoms a short distance away from each other. When they are not bonded together, the coulomb energy is nonzero and the bond energy is (naturally) zero. When they are bonded, the coulomb energy is zero and the bond energy is nonzero. I want LAMMPS to calculate both the bond energy and the coulomb energy when the atoms are bonded. How can I go about doing this?

Thanks,
Scott

chainTestData (223 Bytes)

chainTestRun (680 Bytes)

4

Hi Scott,

I think you need to use the "special_bonds" command
http://lammps.sandia.gov/doc/special_bonds.html

Jan-Michael

5

Thank you Jan-Michael,

It turns out that the behavior I wanted was already built into LAMMPS; by default, all pairwise interactions are turned on except for <1-5 neighbors. I’d been confining my DNA molecule to a spherical volume by binding every pseudoatom to a central anchor atom with a custom bondstyle, which meant that all the atoms were at most 1-3 neighbors. Now I’m confining them via a pairwise potential, and everything’s working fine.