Local frame missing during potential fitting

[CanisW]

Hi,

I encountered a few molecules where some atoms' multipoles become all zeros after potential fitting and lost their local frame reference atoms. Here're the input/output files for one of them.

Tinker version: 8.10.3

The command: potential 6 combined.xyz -k gdb11_s07-3656.key_2 combined.pot N 0.1

The problematic multipoles: The input gdb11_s07-3656.key_2:

multipole    404    403    405              0.05046
                                      0.36015    0.00000    0.54555
                                      0.15547
                                      0.00000   -0.03102
                                     -0.37181    0.00000   -0.12445
multipole    405    404    406             -0.05961
                                      0.35875    0.00000    0.15007
                                      0.13518
                                      0.00000   -0.93998
                                     -0.49566    0.00000    0.80480

The output combined.key:

multipole   405                         0.00000
                                        0.00000    0.00000    0.00000
                                        0.00000
                                        0.00000    0.00000
                                        0.00000    0.00000    0.00000
multipole   404                         0.00000
                                        0.00000    0.00000    0.00000
                                        0.00000
                                        0.00000    0.00000
                                        0.00000    0.00000    0.00000

All the input / output files: example_mol.zip

Any advice / hints on that? Thanks!

[jayponder]

There is a problem with your XYZ file and/or your parameters. For example, in the XYZ file you attached, atom 5 is an oxygen with atom type 405, and it is bonded only to atom 4 which is a nitrogen with atom type 406. But in your parameters the only multipole parameter you have for atom type 405 is:

multipole 405 404 406 -0.05961 0.35875 0.00000 0.15007 0.13518 0.00000 -0.93998 -0.49566 0.00000 0.80480

This parameter can only be used for an atom type 405 that is bonded to an atom type 404. So this parameter cannot be used for atom 5 in your XYZ file. Thus atom 5 in your XYZ file has no parameters, and all of its multipoles are set zero since no match was found in the parameters you supplied. Thus the result of the POTENTIAL program also has all zeros for that atom. A similar thing occurs for atom 6 in your XYZ file. In fact any multipole values that start at exactly zero remain at zero during POTENTIAL fitting.

You can see that these atoms have no assigned multipole parameters. If you run the Tinker ANALYZE program before trying to run POTENTIAL, you will get a warning that your atoms 5 and 6 have no multipole parameters.

[leucinw]

Hi @jayponder I understand the issue is due to the local frame definition. However, there exists another problem: this tinker XYZ file from poledit is not correct. I believe this is again because of the radii definition in Tinker.

FYI, I tried another approach using babel to generate the MM2 style format, and the connections look fine. See below:

9 MM2 parameters
1  N      2.205271   -0.638209    0.113429     9     2     8
2  C      1.020063   -0.174257   -0.040700     2     7     1     3
3  C      0.476578    1.145136    0.159952     2     2     9     4
4  N     -0.800627    1.387306    0.331632     9     3     5
5  O     -1.649083    0.339350    0.298171     7     6     4
6  N     -1.144953   -0.993835   -0.395125     8     7     5
7  N      0.058678   -1.128189   -0.538897     9     6     2
8  H      2.817003    0.103261    0.477225     5     1
9  H      1.111340    2.026844    0.164467     5     3

[jayponder]

Hi Chengwen, You certainly may be right. But in order to check that, I would need to have the DMA output from either Gaussian or Psi4 to run through poledit. Is the original poster from your lab? And if so can you send the DMA output for me to test.

I know the atomic radii and distance-based bond length method used in Tinker to find the connectivity from only input XYZ coordinates is not perfect.. I'll see if I can fix this case, and maybe make the connectivity assignment more robust in general.

[leucinw]

Yes it's Yanxing's @CanisW poltype job :)

All the files are here: https://biomol.bme.utexas.edu/~liuchw/gdb11_s07-3656.tar

[jayponder]

Sorry for the huge delay. I kind of forgot this, and just took a look again.

This is a really weird molecule! In fact I can't even find it an any databases I've checked. Are you sure it really exists? Regardless, the problem is that POLEDIT is not finding the N-O single bond in the ring. In the structure you got from QM, the bond length is 1.585 Ang, which is incredibly long! The "typical" N-O single bond length is 1.36 Ang. This is why POLEDIT is not finding the bond. The difference of more than 0.2 Ang from "typical" is huge, and it falls outside POLEDIT's built-in tolerances.

I'm tempted to say that even if this structure is "real", that it is a very special case. I'm not sure I want to change POLEDIT just to handle this single case. If you have other molecules where bonds are not assigned correctly let me know.

[CanisW]

Thank you very much for the explanation! Yes, it's a generated molecule.