Error 88 in DL_POLY_4 with Larger QM Region in QM/MM (legend Array Exceeded)

Py-ChemShell
1

Dear DL_POLY Support Team,

I hope you are well.

I am currently using ChemShell (v23.0.3) with DL_POLY_4 (compiled from the dl-poly-5.1.0-pre release) for QM/MM simulations using Orca/dl_poly. I’m encountering the following error in _dl_poly.out when I increase the QM region beyond 42 atoms:

DL_POLY_4 terminated due to error    88  
error - legend array exceeded in build_book_intra

Smaller QM regions work without issue, but increasing the number of QM atoms beyond this threshold leads to termination.

I would be grateful if you could advise whether there are additional parameters or structural changes I need to apply in the code to support a larger QM region.

Thank you very much for your time and support.

Thanks,
Hafiz

2

Dear Hafiz,

Sorry to hear you are running into a problem. You can look up dl_poly error codes in the manual to get more information (Appendix C: DL_POLY_5 Error Messages & User Action — DL_POLY 5.3.0 documentation). For this error you get:

Message 88: error - legend array exceeded in build_book_intra

The second dimension of a legend array has been exceeded.

Action:

If you have an intra-molecular (like) interaction present in abundance in your model that you suspect is driving this out of bound error increase its legend bound value, mxfinteraction, at the end of scan_field, recompile and resubmit. If the error persists contact authors.

I believe that this parameter is increased by ChemShell when dl_poly is compiled by ChemShell, so was dl_poly compiled by ChemShell for your calculation or are you using a standalone dl_poly compilation?

(To compile dl_poly though ChemShell, you should set the option ./setup --dl_poly when compiling ChemShell)

All the best,
Tom

3

Hi Tom,

Thanks for your response.

Yes, I compiled DL_POLY through ChemShell previously using the ./setup --dl_poly option as you suggested. After encountering the “too many atom types in FIELD (scan_field)” error, I attempted to increase the maximum allowed atom types by modifying the ffield.F90 file — changing:

integer, parameter :: mmk = 1000

to

integer, parameter :: mmk = 2000

However, even after making this change, recompiling, and rerunning the calculation, I’m still getting the same error. It seems that the updated parameter is not being picked up in the compiled executable.

Do you know if there are additional steps required to ensure the modified DL_POLY code is actually used in the ChemShell build? Or is there possibly a cached build or an alternative source of the parameter overriding this?

Thanks,
Hafiz

4

Hi Hafiz,

I think in this case, the option --dl_poly will redownload the dl_poly source code and the change you have made locally will not have an effect.

Could you try replacing the file chemsh/interfaces/dl_poly/hack.py.in with the following version to modify the patch ChemShell is applying?

#  Copyright (C) 2022 The authors of Py-ChemShell
#
#  This file is part of Py-ChemShell.
#
#  Py-ChemShell is free software: you can redistribute it and/or modify
#  it under the terms of the GNU Lesser General Public License as
#  published by the Free Software Foundation, either version 3 of the
#  License, or (at your option) any later version.
#
#  Py-ChemShell is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with Py-ChemShell.  If not, see
#  <http://www.gnu.org/licenses/>.

# you.lu@ukri.stfc.org
# Jan 2022


def hack():
    '''Increase the value of mxb in ffield.f90'''

    from os     import path
    from shutil import copy
    
    CHEMSH_ARCH = '@CHEMSH_ARCH@'
    
    CWD = path.abspath(path.join(path.dirname(__file__)))
    DL_POLY_SRC_DIR = path.join(CWD, CHEMSH_ARCH, 'src', 'source')
    FILENAME = path.join(DL_POLY_SRC_DIR, 'ffield.F90')
    
    if path.isfile(FILENAME):
        BAK = FILENAME+'.BAK'
        if not path.isfile(BAK):
            with open(FILENAME, 'r') as fp:
                lines = fp.readlines()
            # backup file
            copy(FILENAME, BAK)
        else:
            with open(BAK, 'r') as fp:
                lines = fp.readlines()
    
        # replace string
        for i, line in enumerate(lines):
            if 'mxb = 6' in line:
                lines[i] = lines[i].replace('mxb = 6', 'mxb = 14')
            if 'mmk = 1000' in line:
                lines[i] = lines[i].replace('mmk = 1000', 'mmk = 2000')
    
        # write out a new file
        with open(FILENAME, 'w') as fp:
            for line in lines:
                fp.write(line)
    
        print("ChemShell has hacked DL_POLY's ffield.F90 by increasing mxb")

    else:

        print("\n\n### not found!\n\n")


if __name__ == '__main__':

    hack()

Then you will need to recompile, making sure to have the options ./setup --dl_poly --recompile dl_poly

Just let me know if it doesn’t work.

All the best,
Tom

5

Hi Tom,

I believe the problem has already been fixed. Earlier, I was only changing the mmk value and not the mxb value. Now that I have updated both parameters, the calculation is running successfully with more than 50 atoms.

Integer, Parameter :: mmk = 2000, mxb = 100

I am currently struggling to search for the transition state (TS) of my system. I want to scan the distance between two atoms (i.e., atom 10 and atom 12) by reducing the distance from 5.6 Å to 1.0 Å. I am using ORCA and DL_POLY. If you have any advice or sample scripts related to TS research, I would really appreciate it.

Thanks,
Hafiz

6

Hi Hafiz,

Glad to hear you have it working now. It might be worth us increasing the value of mmk in ChemShell then, does your simulation run with mxb = 14 and mmk = 2000?

With respect to transition state searches, have you seen the tutorial on this? Transition state optimisation — Py-ChemShell 23.0.3 documentation

The options for transition state are documented in the manual in this section Geometry Optimisation — Py-ChemShell 23.0.3 documentation

There are also some test examples in tests/dl_find that you can look at.

It is always helpful for us to get feedback on our documentation, so please do ask if anything is unclear or if this doesn’t include the information you want.

All the best,
Tom

7

Thanks for the follow-up.

No I have not tried mxb = 14 with mmk = 2000, but I can test that configuration as well to see if it’s sufficient. But after setting mxb = 100 and mmk = 2000, the simulation is running smoothly now.

Regarding transition state searches, I’ve gone through the tutorial and the relevant sections of the manual. The information there is definitely helpful, particularly for smaller systems. However, for my current system, I’m trying to follow a similar approach to what I typically use with ORCA: performing a scan calculation from the reactant structure toward the product, in hopes of locating a good guess for the TS, which I can then refine through a proper TS optimisation.

I’m not entirely sure if this scan-based strategy can be effectively implemented in ChemShell, since it seems more geared toward cases where both reactant and product geometries are known. I’d appreciate any guidance or confirmation on whether this type of scan-and-refine workflow is feasible within ChemShell.

I wanted to do this
The DL-find optimizer for reaction path calculations via relaxed potential energy scans, with a step size of 0.1Å along the reaction coordinate.

Using something like this
opt = Opt(
theory=qmmm,
coordinates=‘dlc’,
constraints=[(‘scan’, ‘distance’, 0, 1, 1.0, 3.0, 21)],
algorithm=‘lbfgs’,
trust_radius=‘const’,
tolerance=0.0003
)
opt.run()

Thanks,
Hafiz

8

Hi Hafiz,

If I understood you correctly, you want to reproduce this functionality of Orca in ChemShell: 4.2. Surface Scans - ORCA 6.1 Manual?

As it happens, this functionality to perform a constrained minimisation in Py-ChemShell has recently been developed and is being tested internally at the moment.

All the best,
Tom

9

Hi Tom,

Yes, that’s exactly what I was referring to. Great to hear that this feature has recently been developed in Py-ChemShell. I’d be happy to follow its progress or help test it if possible.

Thanks,
Hafiz