Memory management issue in the diskless API

[Honza-R]

Describe the bug
I ran into a problem with MOZYME calculations ran through the API via the mopactools python interface. When The calculations of larger molecules fail with a core dump (but water molecule is OK). The same calculation ran normally with MOPAC works, so the issue seems to be specific to the API. I tried to debug it, and found that the calculation actually finishes (the energy is calculated), but then it crashes when de-allocating the arrays.

To Reproduce
Here is a minimal example where normal SCF works and MOZYME fails:

#!/usr/bin/env python3

from mopactools import api
import numpy as np

natom = 76
elements = ['C', 'C', 'C', 'C', 'C', 'C', 'C', 'O', 'C', 'C', 'C', 'O', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H']
coord = np.array([0.96, 10.6998, 7.3661, 2.4796, 10.6822, 7.3561, 3.018, 10.5288, 5.9367, 4.5462, 10.5084, 5.9162, 5.0556, 10.3551, 4.4859, 6.5771, 10.3398, 4.4199, 7.0439, 10.1943, 2.9783, 8.464, 10.201, 2.9569, 8.9714, 10.0543, 1.632, 10.4866, 9.9753, 1.6918, 10.9319, 8.7911, 2.5384, 12.3528, 8.7074, 2.4917, 12.8316, 7.6504, 3.3203, 14.3513, 7.679, 3.3192, 14.8765, 8.9869, 3.9066, 16.3992, 9.0057, 3.9096, 16.9132, 10.3169, 4.4896, 18.4349, 10.3382, 4.4917, 18.9493, 11.6468, 5.0758, 20.4714, 11.6657, 5.076, 20.9903, 12.9734, 5.6584, 22.513, 12.9889, 5.6596, 23.0386, 14.2958, 6.2427, 24.5643, 14.3068, 6.2437, 25.1009, 15.61, 6.8268, 26.62, 15.6215, 6.8273, 0.5884, 10.8107, 8.3895, 0.5551, 9.7694, 6.9552, 0.5744, 11.5339, 6.771, 2.8364, 9.8557, 7.9814, 2.8552, 11.6118, 7.7988, 2.6545, 11.3559, 5.3142, 2.6332, 9.6012, 5.495, 4.9152, 9.6798, 6.5325, 4.9354, 11.4367, 6.3514, 4.6742, 11.1802, 3.8719, 4.6631, 9.4256, 4.0553, 6.9814, 9.5157, 5.0231, 6.9927, 11.2594, 4.8517, 6.6668, 11.0254, 2.3705, 6.6774, 9.2497, 2.56, 8.6643, 10.9137, 1.026, 8.5675, 9.1394, 1.1822, 10.8967, 10.8947, 2.128, 10.9029, 9.8956, 0.6811, 10.5058, 7.8604, 2.1458, 10.6098, 8.9205, 3.5777, 12.4787, 6.6921, 2.9248, 12.4553, 7.7737, 4.343, 14.7176, 7.5764, 2.29, 14.7316, 6.8295, 3.8974, 14.5019, 9.1113, 4.9303, 14.5072, 9.8391, 3.3225, 16.7746, 8.8817, 2.886, 16.7797, 8.1635, 4.4998, 16.5387, 10.4426, 5.5129, 16.5319, 11.1597, 3.8999, 18.8089, 10.2166, 3.4677, 18.817, 9.4942, 5.0788, 18.5758, 11.768, 6.1001, 18.5667, 12.4911, 4.4893, 20.8432, 11.5421, 4.0515, 20.8519, 10.8208, 5.6629, 20.6159, 13.0969, 6.6818, 20.6097, 13.8172, 5.0703, 22.8837, 12.8636, 4.635, 22.8896, 12.142, 6.2458, 22.6646, 14.4189, 7.2663, 22.6604, 15.1402, 5.6549, 24.9337, 14.1779, 5.2188, 24.9363, 13.4578, 6.8306, 24.7359, 15.7376, 7.8526, 24.7307, 16.4596, 6.2416, 27.0154, 14.7979, 7.43, 26.9912, 16.5614, 7.2456, 27.0096, 15.5207, 5.8094])

print("MOPAC without MOZYME - works")
system = api.MopacSystem()
system.natom = natom
system.atom = elements
system.coord = coord
state = api.MopacState()
properties = api.from_data(system, state)
print("heat of formation = ", properties.heat)
print()


print("MOPAC + MOZYME - fails")
system = api.MopacSystem()
system.natom = natom
system.atom = elements
system.coord = coord
mozyme = api.MozymeState()
properties = api.from_data(system, mozyme)
print("heat of formation = ", properties.heat)
print()

Operating system
Linux. I get the same issue (sometimes with different memory-management related error messages) when using libmopac.so from the conda package as well as one I built myself.

[godotalgorithm]

Fixing API issues is a high priority for me, but I'll be busy for the next 2 weeks with moving and switching jobs, so I can't deal with this immediately. MOPAC development will no longer be part of my new job and career path (I couldn't find a job that would allow me to continue developing for MOPAC professionally), but I'll still be maintaining MOPAC for the foreseeable future as an unpaid hobby.

I just looked over the MOZYME initialization in the API, and it looks like I'm not initializing all of the sizing information for the MOZYME API data structure when it is created from scratch in Python. I initialize the primary size variable, numat, to zero and it is safe to infer that the other array size variables are zero when that is zero, but I think this needs to be done explicitly for all of the logic on the Fortran side to behave correctly. This problem isn't necessarily the cause of the bug you are reporting, but it is a reasonably likely candidate.

[Honza-R]

I looked at it a bit more and I tracked the issue down to parsing the bond orders in mopac_api_finalize.F90.

The code was writing beyond the end of the bond_order array in the mozyme part of the mopac_record routine. This is caused by a mismatch between the 1st pass and 2nd pass loops. I suppose the 1st pass is correct, while the 2nd pass was writing the (i|i) term multiple times as it was in the loop over j.

I reworked the 2nd pass loop to match pass 1, and it solved the crashes. I also verified that the kk now index stays in the allocated bounds. I haven't looked at the bond order matrix it produces, but I think the logic of the algorithm is now correct - please check it too. I created a pull request here

Note that there may be other bugs unrelated to this issue - the example I posted above still crashes with an error elsewhere in the code, but this fix helped with multiple other calculations that previously failed.

[Honza-R]

I also identified another issue that caused the example to crash for some molecules: I had not set the system.natom_move variable. Initialising it to the number of atoms in the system reliably resolved the issue.

I did not expect this for a single-point calculation, which should not depend on setting the number of moving atoms. The segfault error message was not very helpful in identifying the issue. Ideally, this could be solved by either fixing the MOPAC code that computes the derivatives to accept the default value of 0, or by ensuring that the API does not use 0 as a default (it would make more sense to make all atoms moving as a default used silently), and providing a descriptive error message when it is deliberately set to 0.

[godotalgorithm]

natom_move is relevant to single-point calculations because you can choose to evaluate forces on a limited set of atoms (move is referring to atoms treated by MOPAC as a dynamical degree of freedom for all purposes, including gradient evaluation). MOPAC is very flexible in turning dynamic variables on or off (each Cartesian coordinate of each atom can be flagged in a MOPAC input file), but I wanted to streamline it in the API because that is too much flexibility for the vast majority of use. Specifically, the current design is mainly catering to molecular docking, which I see as the primary use case of freezing some atoms. If there is enough interest, future major-version updates to the API could change how dynamic variables are chosen (e.g. an atom-based freeze/move flag might be a better compromise).

I'll consider a better choice of defaults as I debug and adjust the Python wrapper to the API. Python might accommodate this with None initializations that are converted to sensible defaults before making the API calls.