Dear openbabel support team,

during the attempt to search for the minimum energy rotational conformer of 
some nitrogenated and oxygenated species, like they are important for detailed 
modeling of combustion situations, I'm occuring problems to find correct 
results with the current set of functions I'm aware of within the python 
frontend of openbabel.


The tool I'm working on shall help to identify rotational conformers as a 
preceding step towards a higher level quantum mechanical calculation. The focus 
is set on intermediate to longer chain hydrocarbons with attached functional 
groups. Those might be nitrogenated species with integrated or added nitrogen, 
like nitrites, nitrates or hydrazines, and also oxygenated with peroxy or 
hydroperoxy-groups attached. When it comes to e.g. n-alkylnitrites it starts to 
get difficult to identify the minimum rotational conformer for further 
subsequent calculations because of the large degree of rotational freedom. Such 
a tool would be very helpful for those situations.


Within the development process of the tool I initially followed the suggestions 
of the openbabel manual on single conformer searching like it is described 
here: 
https://open-babel.readthedocs.io/en/latest/3DStructureGen/SingleConformer.html


Unfortunately, on some test cases, e.g. the ROOHs of acetaldehyde (O=CCOO, 
O=C(OO)C) the results of two different force fields (gaff, mmff94) were 
inconsistent and different in the regard of internal hydrogen bonds between the 
carbonyl group and the hydrogen atom of the hydroperoxy group. Similarily, the 
two different force fields lead to different results for n-pentylnitrite 
(Within the gen3d approach). Also for hydrazine (NN) the method is not capable 
to identify both and more importantly the minimum energy conformer at the 
internal bond rotation. Just the anti conformer is identified and not the 
gauche conformer, where the hydrogen bond is formed between the lone-pairs of 
the nitrogen atoms with one of the opposing hydrogens.


Generally, all the available force fields in openbabel and in particular gaff 
and mmff94 are capable to describe the hydrogen bonded structures satisfyingly 
of species containing H, C, N and O. Therefore, the reason for the inconsistent 
results seems to be more statistically within the nature of the weighted rotor 
search method, and actually also the other ones which I tried as well. The 
common denominator of all these methods on first glance seems to be the 
torsional database (torlib.txt), where seemingly no definitions for ...N-N..., 
...O-N..., ...O-O... bonds are available. So I tried to alter it (without any 
knowledge) by adding the following entries:

# hydrazines
[C,H]N([C,H])~N([C,H])[C,H] 1 2 4 5 180.0 60.0 -60.0 120.0 -120.0 0.0 180.0 
30.0 -30.0

# Peroxy, hydroperoxy, nitrites
*O[O,N]* 1 2 3 4 180.0 60.0 -60.0 120.0 -120.0 0.0 180.0 30.0 -30.0

But this did not change anything on the results. Maybe you could give me some 
more tipps or a better description of the smarts and groups in that regard.


Because on my lacking knowledge on all of the details of openbabel, I thought 
it is a good idea to start with a new approach, a brute force analysis of all 
rotational conformers of the internal bonds, loosely following that example:

https://gist.github.com/ghutchis/0d89de9f81157cf1aa9726018d41e8e2


And then I recognized, that not all internal rotational bonds are identified, 
pretty similar like it was previously also stated from Scott McKechnie

https://sourceforge.net/p/openbabel/mailman/message/32125004/


In my case, with the following minimal code snippet:

import openbabel
#smi = 'O=CCOO'
smi = 'NN'

# . openbabel objects
mol = openbabel.OBMol()
conv = openbabel.OBConversion()
build = openbabel.OBBuilder()
ff = openbabel.OBForceField.FindForceField('gaff')

# . initialize molecule
conv.SetInAndOutFormats('smi', 'mdl')
conv.ReadString(mol, smi)
mol.AddHydrogens()
build.Build(mol)

ff.Setup(mol)
ff.SteepestDescent(250)
ff.GetCoordinates(mol)

# set up the pool of rotatable bonds
rl = openbabel.OBRotorList()
rl.Setup(mol)
print("Number of rotatable bonds: {}".format(rl.Size()))

the result for NN gives 0 out of 1 internal rotors and for O=CCOO 2 out of 3 
internal rotors. Seemingly, openbabel does not recognize the ...O-O... and the 
...N-N... coordinate as an internal bond. Additionally I'm presuming, that the 
same is happening with ...N-O... bonds, like in the ...ONO or ...ONO2 groups of 
nitrites and nitrates.


For my research it is quite important to include those kind of rotations in the 
conformational analysis. My suspicion would be, that the inclusion of such 
rotors also solves the problems on the different results with different 
forcefields for these cases. Also for rotational scans with openbabel, 
identifying single rotational bonds in small molecules, like NN, or CC could be 
helpful. Could you please give me suggestions, hints and ideas to persuade 
openbabel to consider this kind of structures/rotations?


Best regards,

Heiko Minwegen
_______________________________________________
OpenBabel-discuss mailing list
OpenBabel-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/openbabel-discuss

Reply via email to