Combine two ReaxFF or Machine-learned FF

[Alireza]

Dear Experts, My system under study consists of a twisted bilayer MoS2 (2700 atoms) which is sandwiched by two hBN (each 722 atoms). I would like to do a optimization in this system within DFTB level of theory. Since there are no suitable parameters for Optimization, I decided to run a reactive force field calculation. I specified first single hBN layer as tags='layer1', the twisted bilayer MoS2 as tags='layer2', and the last hBN layer as tags='layer3'. The calculator and optimization in the script file is as follows:

Code

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------

sw_layer1 = ReaxFF_CHOSMoNiLiBFPN_2021(tags='layer1')
sw_layer2 = ReaxFF_HSMo_2017(tags='layer2')
sw_layer3 = ReaxFF_CHOSMoNiLiBFPN_2021(tags='layer3')

# Combine all 3 potential sets in a single calculator.
calculator = TremoloXCalculator(parameters=[sw_layer1, sw_layer2, sw_layer3])

bulk_configuration.setCalculator(calculator)
bulk_configuration.update()
nlsave('hBN-MoS2-hBN.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
bulk_configuration = OptimizeGeometry(
    bulk_configuration,
    max_forces=0.01*eV/Ang,
    max_stress=0.1*GPa,
    max_steps=400,
    max_step_length=0.2*Ang,
    trajectory_filename='hBN-MoS2-hBN_trajectory.hdf5',
    trajectory_interval=1.0*Minute,
    restart_strategy=RestartFromTrajectory(),
    optimizer_method=LBFGS(),
    enable_optimization_stop_file=True,
)
nlsave('hBN-MoS2-hBN.hdf5', bulk_configuration)
nlprint(bulk_configuration)

But the job ran into this error:

Code

Traceback (most recent call last):
  File "hBN-MoS2-hBN.tags.py", line 3495, in <module>
    sw_layer1 = ReaxFF_CHOSMoNiLiBFPN_2021(tags='layer1')
  File "build/lib/python3.8/site-packages/tremolox/TremoloXReaxFF.py", line 2740, in __init__
  File "build/lib/python3.8/site-packages/tremolox/TremoloXPotentialSet.py", line 24513, in setTags
  File "build/lib/python3.8/site-packages/tremolox/TremoloXPotentialSet.py", line 24430, in actOnlyOnTaggedRegion
  File "build/lib/python3.8/site-packages/tremolox/TremoloXPotentialSet.py", line 9759, in _limitToOneTag
RuntimeError: The ReaxFF potential does not support the usage of tags!
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 14
Traceback (most recent call last):
  File "hBN-MoS2-hBN.tags.py", line 3495, in <module>
    sw_layer1 = ReaxFF_CHOSMoNiLiBFPN_2021(tags='layer1')
  File "build/lib/python3.8/site-packages/tremolox/TremoloXReaxFF.py", line 2740, in __init__
  File "build/lib/python3.8/site-packages/tremolox/TremoloXPotentialSet.py", line 24513, in setTags
  File "build/lib/python3.8/site-packages/tremolox/TremoloXPotentialSet.py", line 24430, in actOnlyOnTaggedRegion
  File "build/lib/python3.8/site-packages/tremolox/TremoloXPotentialSet.py", line 9759, in _limitToOneTag
RuntimeError: The ReaxFF potential does not support the usage of tags!
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 15
slurmstepd: error: task_p_post_term: rmdir(/dev/cpuset/slurm21924163/slurm21924163.4294967294_0) failed Device or resource busy

where the most important part tells

Code

The ReaxFF potential does not support the usage of tags!

I also have plan to use the recent and beautiful feature in qatk, Machine-learned FF. But as the reference I need dftb parameter for Mo and S. For the electronic calculation I am using the parameters by DOI: 10.1021/ct4004959, but these parameters do not contain the repulsion potential. So, optimization is not possible. I need to find a way whether with combining two ReaxFF or using machine-learned FF. Any suggestions are appreciated, Cheers, A

[Julian Schneider]

Hi,

that's an interesting problem.
Unfortunately it is not possible in QATK to combine 2 different ReaxFF potential using tags, and even if that was possible, there would still be the questions how the layers should interact with each other.

As you mention, machine learned FFs, in particular the Moment Tensor Potential (MTP) in QATK, could be a good alternative for such a system. However, there is some effort associated with training a ML FF. In your case, since you there are no repulsive parameters in your DFT model, you'd have to use DFT (probably with dispersion correction) to generate the training data for the ML model.

In your case the training data would probably consist of the following 3 parts:
Training data for MoS2
Training data for BN
Training data for mixed combined BN / MoS2 system.
All of these can proabbly generated using random displacements around the optimized crystal structure, using the functionality described here:https://docs.quantumatk.com/manual/Types/crystalTrainingRandomDisplacements/crystalTrainingRandomDisplacements.html.
In the end all 3 training sets should be combined in the training of a single MomentTensorPotential using this class:
https://docs.quantumatk.com/manual/Types/MomentTensorPotentialTraining/MomentTensorPotentialTraining.html#momenttensorpotentialtraining-c


The cheaper alternative might be to combine 2 Stillinger-Weber potentials, which should work well to descibe for the interactions in each layer, and add Lennard-Jones potentials for the interactions between the layers. Here, the parameters can be taken e.g. from the universal force field (UFF). I have attached a script with an example for such a potential.
Note, there is one subtlety with such a potential: By default the Lennard-Jones potential would act also between the atoms within a layer, and would be added to the Stillinger-Weber interactions, which should already be sufficient alone. This can be avoided to a large extent by setting static bonds (using e.g. the findBonds() method on the configuration), and setting the bonded_mode parameter in the Lennard-Jones potentials to 'mode_14'. In this case the Lennard-Jones potential would only act between atoms that are separated by at least 3 bonds (see the attached example script for details). In this case the effect of the LJ potential within the layer should be minor, but it still makes sense to check e.g. the optimized lattice constants of each material.

Hope this helps.