How to anneal a cg-dna for equilibration run

LAMMPS LAMMPS General Discussion
1

Hello,

I tried to solve my problem ( Part of structure I've simulated has poor fidelity ) of an unstable, ‘wonky’ structure by adding fix nve/asphere as well as angular momentum value to my fix langevin parameter following the fix nve/asphere. Later, I realized that these fix parameters are not invoked during energy minimization, which I require for this run.

From what I understand, I must perform simulated annealing within the input script such that I can minimize and relax the structure.

I also understand that there is only one way to do this minimization + equilibration: by adding ‘velocity create’ with a relatively high temperature value, followed by a third fix known as ‘fix viscous’.

If there is any immediate misunderstanding in the reasoning for the annealing procedure, please let me know. Below is my current input script followed by my intended modifications. Also is attached the data file of my structure of interest.


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)
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

After the minimization command and before/during the fix commands, I was thinking of adding my ‘annealing’ procedure here:

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

Other concerns:

  • How high should I set the initial ‘velocity create’ temperature, and does this mean my fix langevin must change?
  • Is ‘fix viscous’ necessary if I already have fix langevin?**
  • Is annealing the right approach to circumvent the ‘fix’ commands not invoked during EM warning?

**[from the ‘fix viscous’ docs]: Brownian dynamics models also typically include a randomized force term to thermostat the system at a chosen temperature. The fix langevin command does this. It has the same viscous damping term as fix viscous and adds a random force to each atom. The random force term is proportional to the square root of the chosen thermostatting temperature. Thus if you use fix langevin with a target T=0, its random force term is zero, and you are essentially performing the same operation as fix viscous. Also note that the gamma of fix viscous is related to the damping parameter of fix langevin, however the former is specified in units of force/velocity and the latter in units of time, so that it can more easily be used as a thermostat.

output10.dat.lammps.txt (6.9 MB)

2

Follow the instructions at:

https://lorenzo-rovigatti.github.io/oxDNA/relaxation.html

First Monte Carlo then MD with oxDNA, followed by tacoxdna to convert to Lammps input script

That’s what I did based on the advice of Oliver Henrich, the author of CG-DNA package.

1 Like
3

I’ve just read a thoughtful post on edtech, by Dan Meyer, about one of OpenAI’s co-founders starting an edtech company. There’s a quote from that article about compilers that applies just as well to molecular dynamics, and LAMMPS in particular, and so I will modify it and quote it here:

A satisfying aspect of learning [molecular dynamics] is that [LAMMPS] never responds to your [input script] with just a :white_check_mark: or an :x:. If you have an error, you receive information about which line of code had what kind of error. If you don’t have an error, you’ll see something happen. It may not be what you wanted to see happen, though, which is generally a very interesting and productive moment of learning.

Every single question you have asked in the first post is a question you can answer yourself by experimentation, assuming you have two things. The first is enough time, but that should not be an issue, especially with oxDNA being an incredibly efficient coarse-grained model.

The second is either incredibly simple but tedious or incredibly complicated but fast: you need to understand what “should” happen based on the commands you input. (If you do not know what to expect when you input one command instead of another, then how could you know what commands are “right” and what are not?)

The incredibly simple but tedious way is, as renowned physicist Shakira said, Try Everything. (And keep careful records.) If you try every combination of input commands and observe how the output changes, you will understand how each command works.

The incredibly complicated but fast way is to write down a lot of mathematics. LAMMPS “just” solves partial differential equations that physically model how particles move over time. You can write down those equations and analytically work out how certain system quantities will evolve as a result. You can write down the equations of motion for

  • particles with fix langevin
  • particles with fix brownian
  • particles with both

and work out how the total average energy changes over time.

Of course, most of us do not work at Google or Microsoft, so we do not have enough compute to Try Everything, and we did not receive PhDs in mathematical physics in Russia, so we do not have enough brainpower to Derive Everything. So we combine both approaches: do just enough math to form an initial hypothesis, then do just enough simulations to test that hypothesis, and repeat until we have somewhat publishable results.

Now this is a very long-winded post where I could’ve simply answered your first questions with “… I don’t know, try?” But I’m writing all this to emphasise that you must achieve your own personal understanding of what exactly you have simulated, and how, and why.

After all, it’s not obvious to me why your DNA initially misfolded. I’m not sure it’s obvious to you either. But it’s even less obvious to me that simulated annealing (which usually means you heat a system and then cool it down slowly) is what you want, given that 95% of your initial structure looks okay. Why not just whack the misbehaving 5% into the correct initial structure by directly editing the data file? If you’ve already done that, and the oxDNA potential is spontaneously undoing your hard work, then maybe that unfolding is physically feasible? How would you know? How would you rule it out?

Those are the questions you really need to answer. Which fix to choose is much more straightforward once you understand what you are actually hoping to do to begin with.

2 Likes
4

Thanks. And no, I’m not sure myself nor am I certain the problem with the structure is standalone to a singular beam anymore-- the problem is the full structure shortly during equilibrating. I will write a new forum post that considers yours’ and community’s reasoning guidelines better. Sorry!