How to get fragment charge or partial charge

[YasinEl]

Hello and thank you for this package!

I was wondering how QCxMS determines the charge of the produced CID fragments (i.e. when a fragmentation of a [M+H]+ occurs in a given trajectory generally only one of the fragments should be charged)? Also is there a way to get partial charges or electron densities for the individual atoms?

Thank you for the help,
Yasin

[JayTheDog]

That's correct, QCxMS computes the ionization potentials (IPs) of both fragments, and assigns the charge to the fragment with the lower IP (as this one is more likely to be charged). If the IPs are close, the charge is distributed between the fragments based on a Boltzmann distribution.
The charge is printed in the output and the .res file for each fragmentation step. However, because we do not assume partial charge in the GFN-xTB calculations, we assume that the fragments will retain an entire charge for the resuming MD calculation.

[YasinEl]

Thank you, this was what I was looking for! But could it be that this approach fails sometimes? For example, I have a case like the following where a hydrogen-fragment is apparently assigned no charge which seems unlikely. Or is there something I am misunderstanding?

 ============================ 
 Starting Collision Dynamics: 
 ============================ 
 
 initializing GFN2-xTB ...
 initialization successful!
 energy =      -44.13676
 charge =              1
 mult   =              1
 etemp  =         5000.0
 initializing S=0 state using HS-UHF.
 
 Est. no. steps (collision):  540    Distance (A):    19.747423
 
   step  time [fs]    Epot       Ekin       Etot      eTemp    Avg. Temp.
    100   6765.    -1200.4249   19.9889   -24.1259     5000.    1918.399
    200   6815.    -1200.7755   20.2213   -23.9064     5000.    1922.391
    300   6865.    -1202.1471   21.5925   -22.5856     5000.    1943.802
    400   6915.    -1199.8078   19.2218   -24.8704     5000.    1966.716
 FRAGMENTATION occured!
STOP      after   1417 steps
    500   6965.    -1200.3693   19.8696   -24.2431     5000.    1996.925
 
COLLISION after    480 steps
STOP      after   1330 steps
 
    600   7015.    -1201.9581   21.2774   -22.8937     5000.    2142.759
    700   7065.    -1199.9757   19.3094   -24.7889     5000.    2151.020
    800   7115.    -1200.7738   20.1704   -23.9572     5000.    2121.133
    900   7165.    -1201.1449   20.6822   -23.4591     5000.    2069.124
   1000   7215.    -1200.5783   20.0175   -24.1030     5000.    2062.520
   1100   7265.    -1200.3383   19.9336   -24.1780     5000.    2035.645
   1200   7315.    -1201.2061   20.6148   -23.5287     5000.    2024.516
   1300   7365.    -1201.3788   20.7282   -23.4217     5000.    2019.171
 The collision MD finished in        1330  steps
##################################################
 Energetics after collision: 
 
 Internal Energy fragment # 1 :       0.558283 eV
 Internal Temp fragment # 1 :    4520.303692 eV
 Internal Energy fragment # 2 :       0.040705 eV
 Internal Temp fragment # 2 :    8569.182493 eV
 Kinetic Energy  (COM)         :       0.607340 eV
##################################################
 
 E X I T - the CID module is finished
 
 fragment assigment list:111111111111111111111111121
 computing average fragment structures ...
  inter fragment distances (Angst.)

           1         2
 
    1   0.00000
    2   2.71418   0.00000

 computing IPs with             GFN2-xTB at (K)   2014.
 fragment  1 E(N)=    -38.4311  E(I)=    -38.0159             IP/EA(eV)=   11.30
 fragment  2 E(N)=     -0.3935  E(I)=      0.2281             IP/EA(eV)=   16.92

 wall time for IP (s)       0.1

 H ATOM
 
 Summed charges per fragment
           1  0.999999999999991     
           2  8.791524012799498E-015
 
 
Charge of ..
- Fragment # 1          1.00000000
- Fragment # 2          0.00000000
 
  mass                formula                      q pop   spin    |q IPB|  diss time (ps)
 M=193.24             H11C14N1    27   4   1   1   0.743   0.000   1.000    0.000 ~
 M=1.01                     H1    27   4   1   2   0.231   0.000   0.000    0.000
 
-- No of overall fragmentations:   1 --

[YasinEl]

might be a bit of a special case though because it seems like the hydrogen reacted with the collision gas

[JayTheDog]

Well, yes, it can fail, especially for the single hydrogen fragmentation (because this is not reasonable). However, in this case you are right, you have a reaction between the dinitrogen and the hydrogenatom, which is a reaction that can occur (topic of issue #36 ).
Nevertheless, IP calculations can be wrong, especially in this case as only the Hydrogen and not N2H is calculated. In other close cases, using a higher level method for IP calculations can help, but the process of charge distribution during fragmentation is not an easy thing to compute...

[YasinEl]

I think I have found one of the cases where the IP is close between fragments, which made me realize that I still don't entirely understand how charges are handled here. We have two fragments where the assigned charges are

Charge of ..
- Fragment # 1          0.36221358
- Fragment # 2          0.63778642

From all the cases I saw the fragment with the higher charge (in this case Fragment # 2) retains the charge and is further simulated in the next steps. However, here the opposite happens and Fragment # 1 gets the charge. Above you said that where the IP is close the charge is assigned via Boltzmann. Could it be that the result of the Boltzmann is not shown in the .out file, or am I missing something? Is the only way for me to find out which fragment gets the charge to map to the statistical charge at the beginning of the next trajectory (which would only work if there is a "next trajectory":

================================================================================
trajectory    85      collision   4 / 6
--------------------------------------------------------------------------------
 
 total charge          :   1
 statistical charge    :  0.36221
================================================================================

Below I have put more of the respective .out file:

================================================================================
trajectory    85      collision   3 / 9
--------------------------------------------------------------------------------
 
 total charge          :   1
 statistical charge    :  1.00000
================================================================================
 ----- Ion properties -------------------------
 Number of atoms              :      13
 Mass collision Ion           :     109.109460 Da
 Mass Collsion Gas            :      14.007000 Da
  
  ----- Collision Energetics  -----------------
 Kinetic Energy (LAB)         :      30.962772 eV
 Kinetic Energy (COM)         :       3.522645 eV
 Velocity now                 :    7400.045195 m/s
  
  ----- Collision calculations ----------------
 Radius Ion                   :       0.45E-09 m
 Radius Gas Atom              :       0.19E-09 m
 Cross-section                :       0.13E-17 m²
                                ---------------
 Mean free path               :       0.023863 m
 Calculated collisions        :       8.381056
                                ===============
 Number of collisions varied  :       9
  
                                ---------------
  
 ============================ 
 Starting Collision Dynamics: 
 ============================ 
 
 initializing GFN2-xTB ...
 initialization successful!
 energy =      -27.92159
 charge =              1
 mult   =              1
 etemp  =         5000.0
 initializing S=0 state using HS-UHF.
 
 Est. no. steps (collision):  260    Distance (A):    19.240119
 
   step  time [fs]    Epot       Ekin       Etot      eTemp    Avg. Temp.
    100   2626.     -759.1606   35.4028     7.5042     5000.    3149.389
    200   2676.     -759.1743   35.3988     7.4997     5000.    3026.507
    300   2726.     -758.4797   34.5207     6.6471     5000.    3257.109
 
COLLISION after    260 steps
STOP      after   1110 steps
 
 FRAGMENTATION occured!
STOP      after    851 steps
    400   2776.     -757.7593   34.1828     6.3357     5000.    4740.319
    500   2826.     -757.9435   34.2962     6.4423     5000.    4710.951
    600   2876.     -758.0488   34.3055     6.4477     5000.    4642.229
    700   2926.     -757.2186   33.2328     5.4055     5000.    4648.975
    800   2976.     -757.2439   33.4997     5.6715     5000.    4633.654
 The collision MD finished in         851  steps
##################################################
 Energetics after collision: 
 
 Internal Energy fragment # 1 :       0.883523 eV
 Internal Temp fragment # 1 :   20666.268298 eV
 Internal Energy fragment # 2 :       0.169404 eV
 Internal Temp fragment # 2 :    8915.621912 eV
 Kinetic Energy  (COM)         :       2.454489 eV
##################################################
 
 E X I T - the CID module is finished
 
 fragment assigment list:1111221111221
 computing average fragment structures ...
  inter fragment distances (Angst.)

           1         2
 
    1   0.00000
    2   4.89149   0.00000

 computing IPs with             GFN2-xTB at (K)   4638.
 fragment  1 E(N)=    -16.5336  E(I)=    -16.0994             IP/EA(eV)=   11.82
 fragment  2 E(N)=     -5.8377  E(I)=     -5.4118             IP/EA(eV)=   11.59

 wall time for IP (s)       0.0

 
 Summed charges per fragment
           1  0.362213583059648     
           2  0.637786416940352     
 
 
Charge of ..
- Fragment # 1          0.36221358
- Fragment # 2          0.63778642
 
  mass                formula                      q pop   spin    |q IPB|  diss time (ps)
 M=81.08                H3C3N3    85   3   1   1   0.653   0.000   0.362    0.000 ~
 M=28.03                H2C1N1    85   3   1   2   0.275   0.000   0.638    0.000
 
-- No of overall fragmentations:   2 --
================================================================================
 - Entering Mean-Free-Path simulation - 

--------------------------------------------------------------------------------
Run #  85      collision   3 / 9 MFP traj.    1
 total charge                 :   1
 statistical charge           :     0.3622
 
================================================================================

 
 initializing GFN2-xTB ...
 initialization successful!
 energy =      -16.09078
 charge =              1
 mult   =              2
 etemp  =         5000.0

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

Eimp (eV) =    0.0     tauIC (fs) =     0.      nstep =    1000

avcycle =   50    more =  250

step   time [fs]    Epot       Ekin       Etot    error  #F   eTemp   frag. T
    100   3050.    -16.10760   0.1715   -15.93613  0.0068 1    5000.    4086.
    200   3100.    -16.11772   0.1869   -15.93086  0.0009 1    5000.    4167.
    300   3150.    -16.10895   0.1808   -15.92811 -0.0020 1    5000.    4175.
    400   3200.    -16.12800   0.1948   -15.93323  0.0030 1    5000.    4185.
    500   3250.    -16.09860   0.1690   -15.92963 -0.0006 1    5000.    4209.
    600   3300.    -16.12681   0.1975   -15.92935 -0.0010 1    5000.    4206.
    700   3350.    -16.10842   0.1814   -15.92701 -0.0033 1    5000.    4229.
    800   3400.    -16.14169   0.2089   -15.93276  0.0024 1    5000.    4229.
    900   3450.    -16.12375   0.1872   -15.93659  0.0062 1    5000.    4242.
   1000   3500.    -16.09486   0.1644   -15.93047  0.0001 1    5000.    4238.
   1000   3501.    -16.09505   0.1644   -15.93065  0.0001 1    5000.    4238.

        E X I T   M D  because of tmax has been reached
 fragment assigment list:111111111
 computing average fragment structures ...
 
 Summed charges per fragment
           1  0.362213583059648     
 
  mass                formula                      q pop   spin    |q IPB|  diss time (ps)
 M=81.08                H3C3N3    85   3   2   1   1.000   0.000   0.362    0.000
 
 
  Cross-section (M+) (m²)      :   1.089873010551957E-018
  Mean free path (m)           :   2.879085605749278E-002
  chamber length (m)           :   0.200000002980232     
 
 --------------------------------------------------
 -- NOTHING HAPPENED - Next collision simulation --
 --------------------------------------------------
  
  
 
================================================================================
trajectory    85      collision   4 / 6
--------------------------------------------------------------------------------
 
 total charge          :   1
 statistical charge    :  0.36221
================================================================================
 ----- Ion properties -------------------------
 Number of atoms              :       9
 Mass collision Ion           :      81.076140 Da
 Mass Collsion Gas            :      14.007000 Da
  
  ----- Collision Energetics  -----------------
 Kinetic Energy (LAB)         :      19.623613 eV
 Kinetic Energy (COM)         :       2.890817 eV
 Velocity now                 :    6834.214570 m/s
  
  ----- Collision calculations ----------------
 Radius Ion                   :       0.40E-09 m
 Radius Gas Atom              :       0.19E-09 m
 Cross-section                :       0.11E-17 m²
                                ---------------
 Mean free path               :       0.028791 m
 Calculated collisions        :       6.946650
                                ===============
 Number of collisions varied  :       6
  
                                ---------------
  
 ============================ 
 Starting Collision Dynamics: 
 ============================ 
 
 initializing GFN2-xTB ...
 initialization successful!
 energy =      -21.87006
 charge =              1
 mult   =              2
 etemp  =         5000.0
 
 Est. no. steps (collision):  180    Distance (A):    17.000000
 
   step  time [fs]    Epot       Ekin       Etot      eTemp    Avg. Temp.
    100   3551.     -596.0738   24.9695     3.0642     5000.    4368.230
    200   3601.     -594.4909   23.5069     1.6597     5000.    4716.274
    300   3651.     -592.8677   21.6516    -0.1358     5000.    4940.859
 
COLLISION after    260 steps
STOP      after   1110 steps
 
 FRAGMENTATION occured!
STOP      after    348 steps
 The collision MD finished in         348  steps
##################################################
 Energetics after collision: 
 
 Internal Energy fragment # 1 :       0.344837 eV
 Internal Temp fragment # 1 :   18148.485481 eV
 Internal Energy fragment # 2 :       0.167635 eV
 Internal Temp fragment # 2 :    7058.001277 eV
 Kinetic Energy  (COM)         :       1.164214 eV
##################################################
 
 E X I T - the CID module is finished
 
 fragment assigment list:112222112
 computing average fragment structures ...
  inter fragment distances (Angst.)

           1         2
 
    1   0.00000
    2   4.03930   0.00000

 computing IPs with             GFN2-xTB at (K)   6423.
 fragment  1 E(N)=     -5.8837  E(I)=     -5.4083             IP/EA(eV)=   12.94
 fragment  2 E(N)=    -10.4434  E(I)=     -9.9767             IP/EA(eV)=   12.70

 wall time for IP (s)       0.0

 
 Summed charges per fragment
           1  0.393909506072703     
           2  0.606090493927297     
 
 
Charge of ..
- Fragment # 1          0.14267937
- Fragment # 2          0.21953421
 
  mass                formula                      q pop   spin    |q IPB|  diss time (ps)
 M=28.03                H2C1N1    85   4   1   1   0.416   0.645   0.143    0.000
 M=53.04                H1C2N2    85   4   1   2   0.419   0.492   0.220    0.000 ~
 
-- No of overall fragmentations:   3 --
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      run not continued because of small fragment(s)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 

 wall time (min)      2.83
Sat Jul 22 15:17:07 PDT 2023
                      --- QCxMS V5.2.1Sep 20 18:00:00 CEST 2022  ---
 normal termination of QCxMS

[JayTheDog]

Hey, yes, this is a good example of very close IPs. The Boltzmann-weigthing is applied on both fragments, so that not only a single one obtains the entire charge, but both get "some" of the charge. This means that you will see both fragments in your spectrum, with a different intensity. The reason behind this is statistics, as both fragments are likely. But the one with the higher statistical charge will occur more frequently, so that you will see this fragment with a higher intensity.
However, in this case, QCxMS doesn't choose the fragment with the larger charge for the continuous simulation, but as the charges are this close, we also consider the size of the molecule. Here, the higher-charged fragment would be underneath the threshold of the molecular size that QCxMS would consider still fragment-worthy (small molecules are very unlikely to fragment a lot) and the program would stop. However, in reality, the larger fragment with still considerable size would fragment consecutively, which is what you can see here. Fragment H3C3N3 still has the potential to undergo a next fragmentation event, which is correctly simulated here.

So this is the reason, the first fragment gets "the charge" (which is only =1 because of the way the QC calculations can be conducted), but more importantly, it will be simulated further.
Does this explain your question?

[YasinEl]

Thank you this makes sense!