Cgdna structure wholly unstable during equilibration

LAMMPS LAMMPS General Discussion
1

Structural Instability in DNA Origami Switch Simulations Using LAMMPS

I’m encountering persistent structural instability in my LAMMPS simulations of a DNA origami switch. Despite numerous modifications to both my input script (using oxDNA2) and initial topology in oxView, the structure remains unstable. Here’s a summary of my troubleshooting efforts:

oxView Modifications:

  1. Shape: Alternated between hexagonal and square (square caused segmentation faults; hexagonal advised).
  2. Scaffold sequence: Tested ‘RANDOM’ (usual choice) and M13mp18; instability persisted.
  3. Default base: Changing from T to N (A/C/G/T) resulted in segmentation faults.

For each structure, I consistently:

  • Rotate the top beam 90 degrees
  • Translate it downward to form a cross-beam
  • Output top and config files
  • Convert ‘oxDNA → LAMMPS’
  • Use the dat.lammps file as input for LAMMPS read_data

LAMMPS Input Script Modifications:

  1. Fix commands: Changed from nve to nve/asphere; added angmom 3.33 parameter.
  2. Environmental parameters: Adjusted rhos, NaCl concentration, and temperature.
  • Higher rhos produced larger, spikier structures.
  • Note: oxDNA team advised rhos computations are unreliable below ~0.1.
  1. Current attempt: Replacing ‘seqdep’ with ‘seqav’.
  2. Experimented with altering timestep in combination with other modifications.

I suspect my core issue lies in failing to properly manage energy minimization during the equilibration run of this coarse-grained structure. However, omitting energy minimization altogether would be methodologically unsound, also going against the wishes of my advisor (even if doing so is typical of relaxing CG-DNA structs).

Query: How can I address this structural instability while maintaining proper energy minimization? Are there additional parameters or approaches I should consider?

Any insights or suggestions from the community would be greatly appreciated. I’m happy to provide additional details if needed.

Script:

log examplelog.lammps

variable T      equal 0.1	#0.1
variable rhos   equal 0.2 	#0.175

units lj
dimension 3
newton on
boundary p p p 

atom_style hybrid bond ellipsoid oxdna
atom_modify sort 0 1.0

neighbor 1.4 bin
neigh_modify every 1 delay 0 check yes

# Read the data file
read_data output10.dat.lammps.txt 

# Set masses (already set in data file)
#set atom * mass 3.1575

# Define nucleotide groups
group all type 1 4

#oxDNA bond style
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
special_bonds lj 0 1 1

# oxDNA2 pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk     seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond   seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna2/hbond   seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxdna2/hbond   seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff * * oxdna2/xstk    47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh      ${T} ${rhos} 0.815


#Energy Minimization settings
min_style cg
minimize 1.0e-6 1.0e-8 1000 100000


#fix cmds (time integration)
velocity create 0.2
fix 1 all nve/asphere
fix 2 all langevin ${T} ${T} 0.01 12345 angmom 3.333

# Run settings
timestep 0.002
thermo 100
#dump out all xyz 10000 oxdna2o10.lammpstrj 
dump out all xtc 10000 oxdna2o10.xtc
dump_modify out unwrap yes
run 1000000 #1m

Last thread has initial topology attached, which I should add has initial velocities set at zero via oxView: How to anneal a cg-dna for equilibration run

2

(1) What is “rho”?

(2) Have you tried simulating a short duplex (10-15 base pairs) to start with? This way, you can run some very fast simulations to determine if your settings can adequately reproduce oxDNA behaviour for very simple systems to start with, or if there is some fundamental difficulty.

3

‘rho’ is how the oxdna2 authors denote salt ion concentration. Sorry, I should have used parenthesis in my explanation, not commas. I’ll give this a shot but most of their example if not all do not utilize traditional conjugate gradient descent. Will report back later.