Nucleic acids and MD

[gnarw]

I'm trying to run molecular dynamics on a small double stranded DNA (pdb ID: 3ot0) using xtb, ver. 6.5.1
I'm using a potential wall to keep the DNA strands from dissociating. I used the ASE tool to shift the coordinates of the DNA so it's centered relative to the centre of mass.
The DNA has two modified bases that each have an unpaired electron, hence I'm using --uhf 2.
I use the following:
xtb DNA.xyz --gfnff --input DNA_md.inp --gbsa water --chrg -18 --uhf 2 --omd

The DNA_md.inp has the following lines:
$md
temp=295 # in K
time=50 # in ps
dump=50 # in fs
step=2 # in fs
hmass=4
shake=2
$end

$wall
potential=logfermi
ellipsoid: 120.47630861971221, 105.46915010347814, 101.93127904476376, all
$end

I have also run this same DNA in crest with the -nci option.
I am using the ellipsoid parameters from the crest output file.

The MD run stops early with the message:
MD is unstable, emergency exit

Do you have any suggestions how to successfully run MD on DNA using xtb ?

[cplett]

Hi,
One possibility might be that the wall potential is too small. You can check this by looking at the potential energy of the wall potential written as add. restraining in the xtb. output If this value is large, the wall potential is so small that the potential energy of the system heavily increases and the MD becomes unstable. Then you should try to increase the size of the wall potential manually until your MD becomes stable.
Another possibility is that your input geometry has some unreasonable defects/structures, e.g., missing hydrogen atoms, atoms that are in too close proximity, or wrong charges.
If both are not the case, you might try to change the settings. For example, you could try to reduce the step length even further. Or, if the MD performs a few steps until it crashes, you could check, if something happens that leads to the emergency exit. This could be prevented by additional geometric constraints.
You can also provide the xtb output and your input geometries so that we can help you specifically.

[gnarw]

I checked the potential energy of the wall potential, the add.restraining.
That value is 0.00000 Eh.
I'm not sure what "large" value for the potential wall energy is.
Here are more detailed information on the input files and xtb results for the experiments I have done.

1. I removed all solvent atoms and added missing hydrogens to the initial DNA X-ray structure (PDB ID 3ot0).

2. Used xtb and gfnff to obtain a geometry optimised structure of the DNA. Saved as 3ot0edit_gfnff_opt.xyz
   xtb 3ot0edit.xyz --gfnff --gbsa water --chrg -18 --uhf 2 --opt
   This completes without problems.

3. Used the ASE tool to shift the coordinates of the DNA structure
   Output stored as 3ot0edit_gfnff_opt_shifted.xyz

4. Ran MD with no potential wall
   xtb 3ot0edit_gfnff_opt_shifted.xyz --gfnff --input 3ot0edit_mdNoWall.inp --gbsa water --chrg -18 --uhf 2 --md
   The output "3ot0edit_mdNoWall.out" says MD is unstable
   add. restraining is 0 Eh

5. Ran MD with logfermi potential and ellipsoid parameters from CREST
   xtb 3ot0edit_gfnff_opt_shifted.xyz --gfnff --input 3ot0edit_mdWall.inp --gbsa water --chrg -18 --uhf 2 --md
   The output "3ot0edit_mdWall.out" says MD is unstable
   add.restraining is 0 Eh

6. From the "3ot0edit_mdWall.out" I can see that the ellipsoid radii are fairly large compared to the DNA
   ellipsoidal wallpotenial with radii 63.7533229 Å 55.8118759 Å 53.9397150 Å
   I measured the approximate length and width of the DNA structure in pymol and set the radii according to those measurements.
   New radii are:
   ellipsoidal wallpotenial with radii 17.9920268 Å 11.1127225 Å 11.1127225 Å
   Running:
   xtb 3ot0edit_gfnff_opt_shifted.xyz --gfnff --input 3ot0edit_mdWall2.inp --gbsa water --chrg -18 --uhf 2 --md
   The output "3ot0edit_mdWall2.out says MD is unstable
   add.restraining is now 4.34 Eh

All MD input and output files are attached.
Hope this information helps figuring out what I could be doing wrong.

3ot0edit_mdNoWall.out.txt
3ot0edit_mdWall.out.txt
3ot0edit_gfnff_opt_shifted.xyz.txt
3ot0edit_mdNoWall.inp.txt
3ot0edit_mdWall.inp.txt
3ot0edit_mdWall2.inp.txt
3ot0edit_mdWall2.out.txt

[cplett]

As the system is highly charged and rather large, it is quite dynamic. Thus, a time step of 2 fs seems to be too large. For me, the MD worked by setting step=0.5 in the additional input file. The wall potential should not be a problem as it seems to be large enough.

[gnarw]

This is still not working for me.
I changed the time step to 0.5 (step=0.5) in the input file. I ran the MD using both the large and smaller ellipsoid wall potential, as previously mentioned.
In both cases the MD run exits with "MD is unstable, emergent exit"
The MD outputs are attached.

3ot0edit_md.out.txt
3ot0edit_md2.out.txt

[cplett]

It seems that the time step is still too large for a longer MD.
I tried step=0.1 and the simulation seems to be stable now.

[gnarw]

I tried setting the time step to 0.1 fs and the MD length was 50 ps.
The run still exits. So I tried a step size of 0.01 fs. Now the MD run for
much longer but still exits before finished.

3ot0edit_md.out.txt

[cplett]

Another possibility to increase the stability of the MD is to increase the hmass. For example, setting hmass=5 already makes the MD much more stable. Maybe this in combination with the small time step will make the MD stable.

[gnarw]

I set the hmass to 5. Still the MD run exits because it's unstable. I even tried to reduce the time step to 0.01 fs but the MD is still unstable.
So I tried to remove the negative charge of the DNA by adding hydrogens to the negatively charged oxygen atoms on the phosphate backbone.
Now the MD run completes without error. Even with a step size of 0.5 ps.
I also tried running MD on the original structure (i.e. not adding hydrogens to the negatively charged oxygen atoms) but not setting the --chrg flag to -18.
Here the MD also finishes without warning.
The question is. can I use the structure as is and simply treat the system as uncharged by setting --chrg 0 ?

[cplett]

This seems to be some kind of problem that GFN-FF has with the high charge. I would not suggest setting the charge of a charged structure to 0 as this might produce some unphysical behavior during the MD. The best way to solve this issue would be to add some counter ions (like sodium cations) to your negatively charged moieties so that your total charge will be naturally 0. You may also try to set shake=0, as this can also increase the stability.

[github-actions]

This issue is set to auto-close in 7 days if no further updates would be detected.