GFN-FF MD run stability problems

[TadKalQM]

Dear developers of xTB, I am asking for a help with the unstable MD run. I am trying to simulate the peptide composed of some non canonical amino acids, 25 residues long (close to 500 atoms), for 10 ns. There are no metals in the structure, only one selenium atoms and few phosphates. The overall charge of the peptide is -10.

the xtb was invoked like this:

export OMP_STACKSIZE=16G
export OMP_NUM_THREADS=8

xtb structure.xyz --omd --gfnff --input md_prod.inp --chrg -10 --alpb water

where "structure.xyz" was already optimized on the gfn2 level in advance, and it is further optimized on GFN-FF level using --omd flag. the input file md_prod.inp contain the default input with recommended settings for GFN-FF:

$md
   temp=300 # in K
   time= 10000.0  # in ps
   dump= 2000.0  # in fs
   step=  2.0  # in fs
   velo=false
   nvt =true
   hmass=4
   shake=0
   sccacc=2.0
$end

So everything should be fine in my eyes. However, the simulation runs for few hundreds of ps at best, and then suddenly crash like this:

     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
  -87.8815208562450        12779.9691723493
 MD is unstable, emergency exit
 but still taking it as converged!
 average properties
 Epot               :  -88.7944534681385
 Epot (accurate SCC):  -88.7947429662098
 Ekin               :  0.650226318120320
 Etot               :  -88.1442271500182
 T                  :   300.183370269401
 normal exit of md()
########################################################################
[WARNING] Runtime exception occurred
-1- gfnff_setup: Could not read topology file.
########################################################################

------------------------------------------------------------------------
 * finished run on 2024/08/21 at 22:03:11.857
------------------------------------------------------------------------

In the last geometry, one methyl group just exploded - one hydrogen is flying 7 A form the rest of the peptide, and the remaining two and a carbon are unphysically close to the surrouding amino acids.

What am I doing wrong? I tried to change the time step to 1 fs, but it crashed anyway, should I further reduce it? It keeps happening, and it happens also for other similar peptides, so I assume that its not a structure related problem. Should I try to turn on shake to 1 or 2 (instead of 0 recommended for GFN-FF) ? Should I perform some equilibration before the run, with for example nvt=false? The documentation do not mention any so I assumed that the preoptimization using --omd should be enough. Should I somehow change the thermostat? I run xtb 6.7.1, I tried this also on xtb 6.6.1, and it crashed again. I assume I am doing something fundamentally wrong, but couldn't find what it is.

In some attempts (turned on/off the parellelization, different time step) I also got the "thermostating problem" message (here if this helps):

     iter. time                 ...        0 min,  0.036 sec
     iter. time                 ...        0 min,  0.036 sec
  -84.1104101908045        1281871.65493222     
 MD is unstable, emergency exit 
 but still taking it as converged!
 average properties 
 Epot               :  -88.7192583876595     
 Epot (accurate SCC):  -88.7201250053026     
 Ekin               :  0.687133589469596     
 Etot               :  -88.0321247981899     
 T                  :   317.221974817274     
 **thermostating problem**
 normal exit of md()
########################################################################
[WARNING] Runtime exception occurred
-1- gfnff_setup: Could not read topology file.
########################################################################

Thank you very much!

Sincerely, Tadeas

P.S. I am a little bit confused by the output, mostly consisting of almost empty lines with three dots like this:

    iter. time                 ...        0 min,  0.038 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
 586400 1172.80    -88.79424   0.6217  300.  287.   -88.16830
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.041 sec
     iter. time                 ...        0 min,  0.041 sec
     iter. time                 ...        0 min,  0.041 sec
     iter. time                 ...        0 min,  0.041 sec
     iter. time                 ...        0 min,  0.041 sec
     iter. time                 ...        0 min,  0.040 sec

or

   iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
block <Epot> / <T> :     -88.81599  300.     drift: -0.37D-02   Tbath : 300.
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec
     iter. time                 ...        0 min,  0.039 sec

Is this normal?

[cplett]

Hi, Thank you very much for the detailed report. Reducing the step to 1.0 fs is already a good idea to make it more stable. You could additionally try to increase the hmass to 5 or 6 and test different shake modes. However, the systems might simply be unstable with a GFN-FF MD. Another option is to add some cations to reduce the negative net charge. If these suggestions don’t resolve the issue, please feel free to share your input structure so that we can test the system on our end.