API Features

[gbarnes128]

Background Several years ago Jimmy shared mopac2012 object files along with a few select .F90 source files with Bill Hase and me. I used these files to create a tight integration with VENUS as well as my own in-house direct dynamics code.

As has been discussed in this project, creating an API for mopac is challenging due to the inherent file-based input. The technique I implemented creates a dummy mopac input file ('MOPAC input') for initialization purposes only, turns off what output can be easily turned off, and directs the rest to /dev/null. Calculation results can then be obtained by using mopac subroutines and copying data between user and mopac variables. While this ad hoc approach is not truly an API, it may be of use to the community.

PR Features This pull request implements the same approach as that used for mopac2012, but it is more extensive as most of the source was previously not available to me. (BTW, thank you to all of the contributors to this project - open source mopac is an exciting development!)

Several logical flags have been added to allow for the same code to run calculations on systems of different sizes and/or charge states.

It includes interface routines that initialize mopac for a given system, return the SCF energy, return the SCF energy+gradient, return partial charges, and return bond orders. There is also a program included that illustrates how to use the new features.

Status

• [x] Ready for discussion
• [ ] Ready for merge

[godotalgorithm]

I welcome contributions like this, and these features might be of immediate use in implementing a first version of the Python interface to MOPAC, as I will be looking to suppress superfluous outputs.

Unfortunately, I'm about to leave on vacation for 3 weeks, so I can't review this PR right now and no one else is reviewing PRs for MOPAC at this time. This will have a high priority when I get back.

[gbarnes128]

Thanks for the heads up and enjoy your vacation! I start my new position at Illinois State University on Aug 1, so I may be a bit slow to respond around that time.

I had the same thoughts regarding a possible python interface as well.

[thegodone]

any update ?

[gbarnes128]

Please let me know if I can assist.

[godotalgorithm]

I wasn't clear on who was taking initiative on the discussion implied by the "ready for discussion" radio button in the PR and/or how long you were going to be distracted by starting a new faculty position. I was also distracted by an invited J. Chem. Phys. paper that I've been frantically trying to finish for the last several months.

It's probably a good time to resolve this now. What do we need to discuss, and what actions remain before this PR is ready for merging? Is this ready for a full code review, or are there specific portions of the PR that I should focus on and modify first?

From a brief examination, it seems like the example program is serving in place of any documentation of this interface. Is the example illustrating all of the working functionality, or should it be expanded in some way? I see that you are using targeted reset flags to perform the same reset functionality that is usually performed by globally incrementing the numcal variable - have you been thorough about this, or should I take another pass at looking for all reset points in the code and adding similar reset flags?

[gbarnes128]

Sorry about the confusion on who would take the lead - this is my first PR and contribution to an open-source project. We both have had a lot going on!

I would be happy to write up more formal documentation if it would be useful - is there a format/style you would like to follow or an example? The example program does illustrate all of the features that are currently present.

Because of the history of this approach - namely, only having access to a handful of source code files - using the targeted reset flags made sense to 1) avoid having unintended consequences and 2) touch as few mopac variables as possible.

I would have to test if incrementing the numcal variable would work similarly. One advantage of the targeted reset flags is that specific aspects of the calculation can be reset without reinitializing the entire calculation.

I believe the reset flags have been implemented globally, but I have not tested all calculation types. Mozyme, in particular, has not been tested. Having a second set of eyes could be useful.

If the targeted reset flags are acceptable, this should not impact how Mopac functions outside of the "API" and has no current merge conflicts.

[godotalgorithm]

While I will eventually request that documentation be provided with the submission of new features in MOPAC, I don't yet have an appropriate place to put such documentation. Ideally, it will go here when the new documentation website is ready for deployment, but it is still a work in progress and not really ready for contributions.

I need to examine what you are doing more closely during the code review, but I'm happy to accept your approach of using reset flags. While it is making a somewhat complicated system even more complicated, it still seems like an improvement in the absence of refactoring the control flow logic in MOPAC, which is unlikely to ever happen. There is another open interfacing request [#52], and I might be able to use your interface to fulfill that request as well.

From this brief discussion, it seems like this PR is ready for a code review, and I will proceed with that soon. I can probably make any pre-merge adjustments myself, but there is a small chance that I will request modifications from you before merging.

[godotalgorithm]

Can you explain the API recommendations regarding GEO-OK and NOSYM? I assume that certain observed behaviors led to these recommendations, and I would like to understand that better, especially for future documentation purposes.

[gbarnes128]

If GEO-OK is not included, calculations are halted if two atoms are closer than 0.1 Ang. For high-energy collision systems, this can occur and would disrupt a direct dynamics simulation.

NOSYM is included to guard against two subsequent calls for a mopac energy/gradient having a different point group for the same molecule due to a change in geometry between calls.

[godotalgorithm]

I'm sorry about the very extended delay, but I've finally finished a first draft of a diskless/stateless MOPAC API that was inspired by this PR. Once testing is complete, I expect that the new API will supersede this API, and it should be able to do everything that this API can do and more. I'm interested in feedback, and I'm happy to include features that I might have overlooked.

I used the same strategy as this API to eliminate disk-based output by piping the output to a dummy location (e.g. /dev/null). I have also eliminated all disk-based input by creating a new, minimal initialization process for the set of features exposed by the API. To avoid confusion, the set of options is very limited and controlled by the input data type, while the raw keyword line is hidden by the API. The goal is to provide simple access to MOPAC's core features to a larger audience that is otherwise unfamiliar with MOPAC and its idiosyncrasies.

[thegodone]

I see that you merge the branch any documentation on how to use the API ?

[godotalgorithm]

At the moment, the best way to understand the new API would be to look at src/interface/mopac_api.F90 for a summary of the data structures and functions, and tests/mopac_api_test.F90 for an example of a program that calls the API. To access and use the mopac_api module, a Fortran program needs to link with the libmopac shared library and a BLAS/LAPACK library.

Basically, you put all input data into the mopac_system derived type, and all of the accessible property outputs from a calculation are returned in the mopac_properties derived type. The mopac_state/mozyme_state derived types save the electronic ground state at the end of the calculation to initialize the next calculation if you are running a sequence of 1SCF calculations to drive molecular dynamics. Uninitialized mopac_state/mozyme_state correspond to starting a calculation from the default atomic or Lewis-structure initial guesses for the electronic state. The operations supported right now are 1SCF, geometry relaxation, and vibrational calculations, using either MOPAC or MOZYME. I'm satisfied with this API design, but I'm still willing to consider changes upon reasonable request. I'm trying to keep the API as simple as I can.

I still need to test robustness of the API and add error handling (integer error codes with a list of meanings). For example, I'm not presently using GEO-OK on the internal keyword line of the API, and I need to test close approaches of atoms to see if I need to add it or find a different work-around.

Before the release, I will write a thin wrapper for this API in Python, where the different derived types are mapped to Python classes and the Fortran API calls are mapped to Python functions. The plan is to have MOPAC usable in Python through this API after installation through Conda Forge or PyPI, but that will be post-release.

I'll put some temporary documentation for the API somewhere before release (maybe the GitHub Wiki for MOPAC, which is empty right now). The next development activity after the API release will be a new MOPAC website, which is where the API documentation will eventually go. I'm probably going to scrap the Sphinx-based test website and use Hugo.

[thegodone]

this is excellent, let me know when I can test the python layer of integration. Really great improvement are coming!