Some questions during using Kinbot

[wangdu817]

Hi, There are some questions during my Kinbot using

1. During reaction search, some fake TS possessing small negative frequency are always found, I think these cases can be filtered by limit the negative frequency higher than -100 or some value to avoid invalid IRC and conformer calculations, hence, I modified the reaction_generator.py as:

if status == 0 and freq[0] < -100. and freq[1] > 0.: self.species.reac_ts_done[index] = 1 elif status == 0 and 0>freq[0] >= -100.: logger.info('\tRxn search failed for {}, imaginary freq is lower than -100.' .format(obj.instance_name)) self.species.reac_ts_done[index] = -999

and it works well for my system, but I am not sure whether other place also should be modified accordingly?

2. During Kinbot running, a error was occurred, is it caused by my modification?

30-Mar-23 10:40:35-ERROR: Traceback (most recent call last): File "/share/home/wd/anaconda3/envs/kinbot-new/bin/kinbot", line 8, in sys.exit(main()) File "/share/home/wd/Kinbot_309/kinbot/kb.py", line 168, in main rg.generate() File "/share/home/wd/Kinbot_309/kinbot/reaction_generator.py", line 493, in generate pr_opt.do_optimization() File "/share/home/wd/Kinbot_309/kinbot/optimize.py", line 319, in do_optimization self.species.kinbot_freqs, self.species.reduced_freqs = frequencies.get_frequencies(self.species, File "/share/home/wd/Kinbot_309/kinbot/frequencies.py", line 35, in get_frequencies hess /= np.sqrt(np.outer(masses, masses)) ValueError: operands could not be broadcast together with shapes (0,) (39,39)

3, Recently I try to use compound method in Guass, such as G4(SP) to calculate the L3 energy, I found codes in pes.py for Gauss SP energy reading, but it seems kinbot does not contain the code to generate input and submit files for Gauss L3 calculation. Besides, the code seems can not read the energy correctly from the log file. (for example, the single point keyword is "G4 Energy")

[juditzador]

Hi @wangdu817, we are already allowing negative frequencies to be accepted. What we do, as described in our recent paper, https://pubs.acs.org/doi/full/10.1021/acs.jpca.2c06558, is: "Some of the stable species in the C5H5 network exhibit a single small imaginary frequency. Such a behavior is hard to avoid, and, while large imaginary frequencies mean that the minimization was unsuccessful, small values simply reflect the flatness of the potential in the direction of a given normal mode that is hard to capture accurately numerically. The problem becomes more prevalent for larger, floppier molecules or when using looser convergence criteria. In KinBot, we accept such outcomes if there is just one imaginary frequency and it is smaller in magnitude than 50 cm–1. For the purposes of the master equation calculation we simply take the absolute value of the normal-mode frequency and provide a warning." As you can see, our threshold, however, is 50 cm-1. The error you are experiencing is most likely because you didn't modify this value, and I think your < sign should be the opposite. The simplest approach is that we will put in a keyword, which defaults to 50, but will allow the user to set a different threshold. This should not take us too long, we could get this done even today after a bit of discussion. At the same time, you should consider running your calculations with tighter convergence criteria, this might be useful, 100 cm-1 sounds like a large value to me - but I cannot know how big or floppy your structures are.

[juditzador]

We will investigate the L3 calculations with composite methods and will get back to you soon.

[juditzador]

@wangdu817 , I added a new keyword called imagfreq_threshold. It's a positive number, and in your case you want to set it to 100. The default is 50.

[wangdu817]

Thanks, professor @juditzador What I meant to say in Q1 is that it is possible to consider some structures with a small imaginary frequency(such as lower than -100) in L1 transition state search are failed and avoid the following IRC, conformer search and L2 opt, because these structures are apparently not correct and the very small negative frequency can be a basis for estimation. For example, in my system, a TS structure is obtained from the "1042482923172701110102_r12_insertion_R_2_1_10" search, Gauss was normally terminated and the structure has a imaginary frequency, but it is obviously not a true transition state.

I have check many similar cases with small imaginary frequency , all of them are not ture TS.

[juditzador]

Hi @wangdu817 , I see, I misunderstood your problem, nevermind. However, I am not sure about some details here.

1. I do not think that we want to exclude a saddle just because the imaginary frequency is smaller (in absolute value) than 100i cm-1. There are reactions where this happens and are valid. The example you show of course has a much smaller imaginary frequency, and it might be possible to implement such a cutoff, but having it in a general way might still be difficult. IRCs are not that expensive in the big picture of things.
2. I know that some of the saddles are simply conformational saddles or fragments drifting apart, like your example. But I do not understand how they then get to conformational search and L2. If the saddle is not reactive (do not connect the original reactant to a new species), they will be filtered out in the IRC step. No conformational search or L2 is done for these species. Are you sure the conformation search and L2 are done for your system??

[juditzador]

May I ask what your threshold energy is in your calculations? Here, you are essentially making :CH2 and and an N-centered radical, this should be really high in energy, and should be excluded already on energetic grounds, which in turn even block the IRC. Of course it depends on your application, but glancing at the filename, you seem to be working on some peroxy chemistry.

[wangdu817]

Hi @juditzador, you are right, most non-reactive saddles are filtered out by IRC, only minor cases, conformation search and L2 were done, such as shown below, The threshold energy is 100 kcal, almost 70 kcal higher than the bimolecular reactants, because it was for combustion application instead of autooxidation.

[juditzador]

Hi @wangdu817 , if you set the threshold to 100 kcal/mol, then you will have to deal with the plethora of pathways. If this reaction shown above was picked up by the IRC, it means that it did produce some sort of a reaction, I assume that this is a bond breaking process. Then it all comes down to the level of theory you are using, if it shows a barrier, then KinBot cannot second guess it. Moreover, the small imaginary frequency signals a very tiny barrier, which means that in your energy range you could expect this reaction to be likely barrierless. I am not sure what you mean by "it was for combustion application instead of autooxidation". Autoxidation is part of combustion. These peroxy radicals are typically ~35 kcal/mol deep (shallower if begin with resonantly stabilized species, and deeper if the ROO itself becomes resonantly stabilized, see https://www.science.org/doi/pdf/10.1126/science.aaa1495 and the picture below). Studying these high energy channels don't make sense for most systems, because both for entropic and energetic reasons the lower energy channels will dominate, including the one leading to the reactants. 100 kcal/mol can break many bonds, while peroxy radicals and their isomers are full of weak bonds, which should dominate reactivity - even at high T.

[wangdu817]

Thank you, professor @juditzador, I lean a lot from your patient replies! I chose a high threshold above the ROO because in a similar reaction system, the pathway of QOOH isomerization is nearly 50 kal/mol and other pathways higher than 60 kcal/mol were also reported. After your guidance, I think may be a lower barrier of 40 or 50 kcal/mol may be more appropriate? Which already includes all potential QOOH intermediates. https://www.sciencedirect.com/science/article/abs/pii/S1540748920303783?via%3Dihub

[juditzador]

The reference you provide is quite extensive, I think these are the gray pathways that are going to of no relevance. It is sometimes useful to show the pathways that are not relevant though to tell the audience that you looked... Nevertheless, as you point out, even those pathways are within at most 70 kcal/mol of the ROO radical. I would say that doing 50 or maybe 60 kcal/mol is more than enough to ensure that you found everything that matters. Of course, you need to use your judgement and know what questions you are trying to answer. The nice thing about KinBot is that you can do a search with a lower cutoff, say 50, and then still see in the kinbot.log files which pathways were abandoned based on high energies. And any time you can restart the calculation with a raised or decreased threshold. It makes more sense to gradually increase the threshold, rather than go with a high one right at the beginning if you are exploring new systems. Otherwise the computer is going to spend a lot of time on high-energy, less important channels, and it'll take longer for the human to understand what the most important pathways are.