GFN-FF energy output showing nonphysical trend

[hillaper]

Describe the bug In the output of GFN-FF molecular dynamics for a sphere of 500 water molecules (and other systems) the total energy and the potential energy in xtb.trj file show a nonphysical trend. Water is confined using an external logfermi potential. The same trend appears in both NVE and NVT ensembles. I attached a plot of potential energy vs time of both ensembles for clarification.

To Reproduce Input files for different ensembles differ only by the obvious "nvt" parameter in xcontrol file.

1. happens with input: xcontrol.txt, water500.xyz

2. start xtb with: xtb --gfnff -I xcontrol.txt water.xyz --md (or with --omd)

3. run xtb with your options and the --verbose flag

4. output showing the error out file and xtb.trj file

Input and output files are attached. The coordinate file is in .txt format as .xyz could not be uploaded.

Expected behaviour Potential energy is expected to eventually settle down and fluctuate around some (local or global) minima after a decent equilibration time without any type of trend.

Additional context Seems like the trend is almost periodic, but with closer inspection it appears not. I have not seen the same trend in GFN0 or GFN2 simulations, but those haven't advanced quite long enough for a good comparison.

xcontrol.txt water500.txt

[awvwgk]

@sespic could you have a look at this?

[sespic]

The input structure contains single oxygen and hydrogen atoms, that are not part of a water molecule. GFN-FF is a dissociative FF, but it is not able to form bonds. Hence the initial set up of the topology does not work and a warning is printed: warning: no bond partners for atom xx. In general we always recommend to start MD simulations from optimized structures (--omd). Here in this case, a GFN2-xTB geometry optimization would be necessary before the GFN-FF MD simulation.

[hillaper]

I've done the geometry optimization before running GFN-FF MD for some other systems (including water), but NOT with GFN2. They showed the same behaviour.

I shall do the optimization with GFN2 now before moving on to GFN-FF MD, and see what happens.

Thank you for a very quick and clear response!

[sespic]

Did you perform the geometry optimization with GFN-FF? For your input structure water500.txt it can't work, because the initial topology setup is wrong. Try to optimize the input structure with GFN2-xTB. Then start the MD simulation with GFN-FF.

[hillaper]

Okay, that definitely sounds reasonable. And yes, the geometry optimization was performed with GFN-FF using the flags --gfnff and --omd. I will perform the optimization with GFN2 now, and and then move on to GFN-FF MD.

I'll report back in a couple days.

[hillaper]

Hello again!

So I did the GFN2-xTB geometry optimization on my original input structure and then ran GFN-FF MD starting from the optimized structure (xtbopt_water.txt). There were no warnings concerning bond partners in the output file and the following lines seemed to be okay as well:

      #atoms :   1500
      #bonds :   1000
      #angl  :   500
      #tors  :   0
      #nmol  :   500
      #optfrag :   500

After plotting the potential energy, it seems like the same type of behaviour still occurs (just not as often as earlier). I decided to check if this had something to do with the logfermi confinement potential, so I repeated the simulation but this time without confining my system. Without confinement the potential energy seems to be acting as expected.

Is it possible that this has something to do with the potential energy printout in xtb.trj? I tried to look for clues in the source code, but couldn't nail anything down. I'll attach my inputs and a picture for clarification again (xcontrol_conf.txt for confined system and xcontrol_noconf.txt for not confined system).

xcontrol_conf.txt xcontrol_noconf.txt xtbopt_water.txt

[hillaper]

Hello,

When you have the time, any type of further assistance on this or guidance towards the files that need to be changed in the source code would be greatly appreciated. I still haven't been able to produce energies without this trend. @sespic @awvwgk

Kind Regards, hillaper

[sespic]

In GFN-FF the H-bond list is updated if a certain RMSD threshold between the initial and current structure is exceeded. Could you please check if the jumps appear, when the H-bond list is updated? Besides that, if it works without the constraining potential, why do you want to apply it then?

[hillaper]

Thank you for the respond! I haven't yet had the opportunity to check if the H-bond list is updated when the jumps appear (not a trivial task for me as I'm not that experienced in Fortran). Decreasing the amount of atoms constrained by the confinement potential, however, decreases the amplitude of these jumps, so that would imply that the jumps have something to do with the constraint.

Ultimately the problem is about knowing what are the terms included in the energy output of xtb.trj and the general output file. I would be particularly interested in the unbiased energy of the system, which does not seem to be the default output (unless I'm mistaken). Same interest goes for some other cases in which I'm using metadynamics, as it could provide a negative image of the underlying PES with the applied Gaussians. When you have the time, can you point out which files are involved and which to edit for a more detailed potential energy output. Thanks!

The constraining potential is necessary for some parts of other systems I'm studying, which experience the same behaviour, so it is necessary to have unfortunately.