Data & Codes used in: Quantum Mechanical Derived (VdW-DFT) Transferable Lennard-Jones and Morse Potentials to Model Cysteine and Alkanethiol Adsorption on Au(111)

Here we provide the data, codes, and outline the procedure to reproduce the results presented in the paper "Quantum Mechanical Derived (VdW-DFT) Transferable Lennard-Jones and Morse Potentials to Model Cysteine and Alkanethiol Adsorption on Au(111)" by E. Ventura-Macias, P. M. Martinez, Rubén Pérez, and J. G. Vilhena.

The following sections contain a detailed description of the data and codes. At the end of this README, you will find instructions on how to generate the classical force-field parameters (Morse and Lennard-Jones) from the potential energy surfaces (PES) computed at the DFT level.

The procedure is general and applies to any given pair of molecule and surface. The provided codes will allow you to swiftly generate the PES at the DFT level, fit the Lennard-Jones and Morse potentials, and test them in a LAMMPS MD simulation.

For a thorough explanation of the procedure and the relevance of these results, please refer to the original publication (ARTICLE_DOI).

If you find this helpful, please consider citing the article (ARTICLE_DOI).

Data structure

The data is organized in the following way:

DFT

The equilibrium adsorption geometry of methanethiol (MTH), propanethiol (PTH), and cysteine (CYS) for the Au-mol configuration with PBE+DFT-D3.

Potential Energy Surface (PES) computed at the DFT level (Figure 3 of the main manuscript).

The scripts used to generate the PES.

MD

Fitting code and general instructions on how to use it.

General input scripts used to generate MD data within LAMMPS.

DFT Data

The DFT data in DFT.zip is organized in the following way:

DFT_PES/

This folder contains the DFT potential energy surfaces (PES) for the interaction of the sulfur atom of methanethiolate (mth), propanethiolate (pth), and cysteine (cys) with the Au(111) surface and the necessary scripts to reproduce it.

results/

The PES are given in one csv file per molecule named as mol + _PES_S-Au111.csv, where mol is the molecule name. Columns are as follows:

Label
Definition

i
calculation number

site
adsorption site

z
distance of the S atom to the surface

deltaz
distance difference from the minimum energy position of the S atom

pbed3
PBE+D3 binding energy

pbe
PBE component of the binding energy

d3
DFT-D3 component of the binding energy

mth/ | pth/ | cys/

Each folder contains the CONTCAR (VASP) file for the optimized geometry of the molecule adsorbed on the Au(111) surface with PBE+D3.

setup_grid.py

Python script to set up the POSCAR files for the PES calculations.

read_results.py

Python script to read the results of the PES calculations.

sub_array.sh

Bash script to submit the PES calculations to an SLURM-based cluster.

INCAR | KPOINTS | surf.CONTCAR

VASP input files for the PES calculations.

MD fitting

The MD data in MD.zip is are organized as:

Fitting/

optimize_Morse.py

Contains the Python script used for the fitting procedure of Morse potential.

optimize_LJ.py

Contains the Python script used for the fitting procedure of LJ potential.

Histogram/

in.test

LAMMPS input script to extract an XY file for the position of the S atom in an NVT simulation.

coord.data

Au surface for LAMMPS simulations.

SCH3.data

SCH3 molecule for LAMMPS simulations.

sheng.eam

EAM potential file in case Au dynamics are wished to be included.

Single_point_scan/

in.scan

LAMMPS input script to perform single-point energy scan of a given molecule.

coord.data

Au surface for LAMMPS simulations.

SCH3.data

SCH3 molecule for LAMMPS simulations.

launch.sh

Launches the single point scan.

plot_scan.py

Molecules

LAMMPS sample geometries for the 3 molecules.

Steps to reproduce the results

The PES calculations were performed using the VASP code and two Python scripts for setup and results parsing.

Requirements

VASP (tested with version 5.4.4)

PAW_PBE pseudopotentials set version 5.4

Python3 with packages:

ASE (Atomic Simulation Environment)

Numpy

Pandas

matplotlib (optional)

Steps

Each molecule has its own directory with the necessary files to reproduce the results. The following steps are for the methanethiolate molecule (mth).

Set up the grid of points for the PES calculations by running the setup_grid.py script. It will create a subfolder run/ inside the molecule's directory with the POSCAR files for each point in the grid.

python setup_grid.py mth

Create the corresponding POTCAR file and place it in the molecule's directory.

Change the sub_array.sh script to match the number of calculations in the array numbers and the MOL variable.

#SBATCH --array=1-number of calculations

MOL=mth

Submit the calculations to a SLURM-based cluster by running the sub_array.sh script from the DFT_PES directory.

sbatch sub_array.sh

After the calculations are finished, run the read_results.py script to parse the results and generate the PES csv files. It takes the arguments --surf and --mol to specify the surface and molecule PBE+D3 and D3 total energies.

python read_results.py mth --surf -128.5228 -17.7147 --mol -22.6026 -0.0083

Reference values for the PBE+D3 and D3 surface and molecule energies

Surface

Molecule

Molecule
PBE+D3
D3
PBE+D3
D3

mth
-128.5228
-17.7147
-22.6026
-0.0083

pth
-128.5228
-17.7147
-55.8615
-0.0086

cys
-128.5228
-17.7147
-74.0338
-0.1611

Fitting procedure

Both the .xyz and .csv files should be located at the same folder as the script. Then, simply run the code (mth is used as an example):

python optimize_Morse.py mth

or

python optimize_LJ.py mth

The script will print the optimized parameters: [De, re, α] or [ϵ, σ] for Morse or LJ respectively. It will also plot a fitting plot and a birdview of the resulting PES.

Using the potential. Histogram example

The code will run for the optimized mth Morse parameters and extract a occ.lammpstrj containing the (x,y) positions of the S atom throughout the NVT simulation. Note that only Au-S interaction is included. Start as:

lmp -in in.test

This is easily adaptable to other routines or molecules and is thought to be a generic LAMMPS starting input.

Single-point energy scan

Go to the 2C) folder and launch the scan with:

./launch.sh

This will create a folder named fine_scan containing 128 folders. Each folder is assigned to an (x,y) position. Then, inside each folder, a single-point energy evaluation is performed at various Z heights around the absolute minima.

The in.scan file should be modified accordingly with the appropriate potentials. It is set to perform the scan with the optimized Morse potential by default.

The output is gathered in the E_readout.dat folder with the following structure (all energies in eV):

Total Energy
Intramolecular energy
Au-mol vdW interaction energy
Au-S interaction energy

-133.62
0.123112
-0.2403
-1.31765

-133.721
0.123112
-0.28-401
-1.37853

Therefore, the total adsorption energy will be the sum of the last two columns.

Energies are ordered in increasing Z for the same (x,y) point. That is, the first 12 lines correspond to 12 heights of the starting (x,y) coordinate, the next 12 lines to heights at the second (x,y) configuration and so on.

The python script plot_scan.py may be used to plot the results. The E_readout.dat file and .csv must be in the same folder.

python plot_scan.py mth