Source of "grad seems to be bogus" error?

[jasondboettger]

Hi folks!

Newbie just getting started with QCxMS here. I've tried to set up a test calculation on PFBS (https://pubchem.ncbi.nlm.nih.gov/compound/Perfluorobutanesulfonic-acid) in negative mode by CID. By dint of being a perfluorinated compound with only a single proton, I thought it'd be a quick example where I wouldn't have to calculate multiple deprotomers, and so could just grab the structure from the sdf there and delete a proton. Here's a chunk of QCxMS's output during the setup phase:

Generating uniform velocity distribution, T=     500.0
 Groundstate eTemp set to   298.149993896484     
 
 initializing GFN2-xTB ...
 initialization successful!
 energy =      -23.34560
 charge =             -1
 mult   =              1
 etemp  =          298.1
 M trajactory, equilibration ...       10625

             E N T E R I N G   M D   M O D U L E

 grad seems to be bogus!
 QC calc failed!
 Gradient fails in MD 1

So it looks like QCxMS was able to calculate a single-point energy successfully, implying to me that my molecule's structure is at least reasonable.

But the grad seems to be bogus! error would I assume imply that the MD module doesn't like the starting energy gradient. Digging through iniqm.f90 I find:

     gn=gnorm(nuc,grad)
     if(gn.lt.1.d-8.or.gn.gt.20.0) then
        write(*,*)'grad seems to be bogus!'
        E=0
        return
     endif

So it looks like this error gets thrown if my gradient is too steep (>20) or possibly too flat/negative? (<1e-8) Am I reading that right? And if so, is there a way to fix (optimize) the starting structure? A starting opt within QCxMS would be ideal, but I'm guessing not possible based on the input parameters I see documented. I'm guessing I'll want to start using CREST anyway, but if generating deprotomers with CREST isn't likely to solve this problem, I'd like to know in advance.

Thanks for your help and perspective!

[JayTheDog]

Hey, welcome to the QCxMS community 👋
so, from my first intuition, I would think that the starting structure might be off. Check if you actually use the 3D structure and this one makes sense (bohr/angström, etc). You could optimize it then by using e.g. xtb or some other QC package like ORCA on GFN2-xTB level (DFT doesn't make sense as QCxMS is using GFN-xTB anyways). QCxMS doesn't do it inherently. Using CREST for conformer search could help getting a better starting guess, but I would not think that it is relevant here.
Then I would check the input file. Are you sure that you are using the correct charge? How does it look?
I would not expect QCxMS to have a problem with such a structure in the first place, so there has to be something off somewhere.

[tobigithub]

@jasondboettger can you share the input file and coordinate file that you used. Tnx

[tobigithub]

@jasondboettger This worked for me, observe there is no hydrogen in the molecule anymore, its de-protonated, directly coming from crest, you can use the crest-runall.sh script to create protomers, deprotomers, tautomers, and openbabel to convert PubChem 3D conformers to tmol (turbomole) or XYZ coordinates.

coord file

$coord
    3.35479475697029     -1.55565353042221     -0.63804602316408      s
   -1.66754326488468      3.43448137789686      2.00710849640390      f
    0.38402074815879      4.12758714495269     -1.45058432795988      f
    3.89267733673916      3.08298969627511      1.52711389418737      f
    1.74837834947823      0.37948200603048      3.70437338096344      f
   -4.05743754885826      2.11569911524996     -1.92841616546046      f
   -1.27022177539320     -0.33684939517807     -3.43845856277377      f
   -5.80077549628295     -0.84737346100393      0.89141827326023      f
   -3.67526369391416     -3.73656531922595     -0.91201283511508      f
   -2.38299116772925     -2.04618553048241      2.58408321048583      f
    5.52387203802844     -2.48147933905698      0.70352616288441      o
    3.82940897707305     -0.38695916664747     -3.04625755493702      o
    1.08924173365026     -3.08594421655317     -0.47493708321273      o
   -0.24086861262847      2.27316381780967      0.18351427956088      c
    2.18094144786630      1.18823444314507      1.32030723326833      c
   -2.12382056646193      0.52163600806371     -1.23789821014719      c
   -3.41267798205039     -1.60333577369542      0.35535942074246      c
$end

and then the input can be simplified or changed, but this worked:

# INPUT FILE
xtb2
cid
elab 20
ecom 20
maxcoll 10
noesi
charge -1
ntraj 64

results and output, I did not check the quality of the spectrum, just if it runs OK, which it does

 re-call of important run parameters:

 QC Program            : xTB
 QC Level              : GFN2-xTB

 total traj.   (ntraj) :   64
 time steps    (tstep) :    0.50 fs
 Initial temp. (tinit) :  500.00 K
 M+ Ion charge(charge) :   -1

----------- CID settings -----------
 Collision Gas         : Ar
 E (LAB)               :   20.00 eV


 internal Energy        :   0.93 eV

and trajectory example.

[jasondboettger]

@tobigithub @JayTheDog Thanks so much for taking the time! Through your help I've realized the issue comes from my unfamiliarity with the Turbomole coord format; all of my interatomic distances are ~0.53x yours, of course because I tried to specify my coordinates in Angstroms instead of Bohrs. I'll be able to test this today, but I expect fixing this should do it.

[jasondboettger]

Thanks again! Using the $coord angs keyword at the start of the coordinate block fixed the issue. Spectrum looks fairly reasonable for a PFAS just using default input parameters: