Performing QM/MM calculations for organic reactions by ORCA

For QM/MM calculations using ORCA, we require an equilibrated PDB file, a prmtop (topology) file, and an ORCA input file describing the QM/MM calculation. Here, we prepare the PDB and topology files using AMBER. We describe the procedure for the S_N2 reaction between aminoxide (AMO) and methyl iodide (MEI). This reaction is illustrated in the following line.

(O−NH₂)⁻ + CH₃I → [O−H₂N…….CH₃………I]^#− → O−H₂NCH₃ + I⁻

1. Calculate force field parameters of the reactants (AMO and MEI in this example)

We describe the reactants with the general AMBER force field (GAFF) [1]. To this, first, optimize the geometrical structures of the reactants at the B3LYP/6-31G(d) [2–6] level of theory and then calculate electrostatic potentials at the Hartree−Fock/6-31G(d) level of theory with points sampled according to the Merz−Kollman scheme [7]. You can use the Gaussian 16 program package for these [8]. Atomic charges are then fitted to these potentials using the RESP procedure [9], as implemented in the antechamber module of the AMBER program [9], which also assigns GAFF atom types to the molecules.

1.1 Optimize the structure of each reactant using the Gaussian program (you can skip this step if you do not have access to the Gaussian program; see section ‎1.2). The following are the Gaussian input files of AMO and MEI.

Gaussian input file for AMO (AMO.gjf)

%chk=AMO.chk

%nprocshared=28

%mem=10GB

#P B3LYP/6-31g* Opt

AMO

-1 1

N 4.00844372120790 -0.40353050497740 1.20130235025608

H 3.83540033261703 -1.40686035366881 0.97191652808456

H 4.08033668933599 0.06872086683233 0.27345411952844

O 5.11522037614787 -0.23619501622205 1.94219401877628

--Link1--

%chk=AMO.chk

%nprocshared=28

%mem=10GB

#P HF/6-31G* SCF=Tight Geom=AllCheck Guess=Read

Pop=MK IOp(6/33=2,6/41=10,6/42=17)

Gaussian input file for MEI (MEI.gjf)

%chk=MEI.chk

%nprocshared=28

%mem=10GB

#P b3lyp/gen pseudo=read Opt

MEI

0 1

C 0.07702343763513 -2.92823719021080 -5.29252485295213

H -0.61791204772384 -3.09775193244021 -4.46350759513109

H 0.97554093106452 -3.55261519252021 -5.25036273289306

H 0.28327955513878 -1.87045309595898 -5.48646799734132

I -1.12533299542341 -3.70107758083385 -7.27800383832773

C H 0

6-31g*

****

I 0

lanl2dz

****

I 0

lanl2dz

--Link1--

%chk=MEI.chk

%nprocshared=28

%mem=10GB

#P HF/gen pseudo=read SCF=Tight Geom=AllCheck Guess=Read pop=ReadRadii

Pop=MK IOp(6/33=2,6/41=10,6/42=17)

C H 0

6-31g*

****

I 0

lanl2dz

****

I 0

lanl2dz

I 1.98

We use the LANL2DZ [10] basis set for iodine to optimize the MEI structure. The last line in the input file of MEI is the radius of iodine which is used by the program for fitting the charges.

Execute the following commands after preparing the Gaussian input files.

g16 < AMO.gjf > AMO.out

g16 < MEI.gjf > MEI.out

1.2 Extract the reactants' GAFF parameters using the antechamber module of AMBER.

antechamber -fi gout -fo prepi -c resp -i AMO.out -o amo.in -rn AMO -at gaff -pf y

antechamber -fi gout -fo prepi -c resp -i MEI.out -o mei.in -rn MEI -at gaff -pf y

If you do not have access to Gaussian, optimize the structure of the reactants using an MM method in Avogadro [11] (Ctrl+Alt+O) and save the structures as PDB files (AMO.pdb and MEI.pdb). Next, extract the GAFF parameters for the reactants using the antechamber module of AMBER with BCC charges.

antechamber -fi pdb -i AMO.pdb -o amo.in -fo prepi -nc -1 -c bcc -rn AMO -at gaff

antechamber -fi pdb -i MEI.pdb -o mei.in -fo prepi -nc 0 -c bcc -rn MEI -at gaff

The ‘-1’ and ‘0’ in the above commands represent the AMO and MEI charges, respectively.

You may need to prepare some missing parameters using AMBER's parmchk2 module. The following commands should be executed if this is the case.

parmchk2 -i amo.in -f prepi -o amo.dat

parmchk2 -i mei.in -f prepi -o mei.dat

We need the amo.in, amo.dat, mei.in, and mei.dat files from the above commands.

2. Prepare a PDB file containing the reactant state of the reaction

To prepare a PDB file for the reactant state, you need to draw structures of the reactants in one single file. GaussViwe or any other graphical interface such as Avogaro [11] can be used for this. The following is the input PDB file for the reactants (reactants.pdb). It is recommended that hydrogens be removed from the reactants.pdb file and that they are added by the tleap module of the AMBER program (see the next step). It is critical to insert TER between the reactants' coordinates. The order and names of atoms in each reactant must match the parameter files (amo.in or mei.in).

The reactants.pdb file

REMARK input PDB file of the reactant state

ATOM 1 N1 AMO 1 4.008 -0.404 1.201 0.159200

ATOM 2 O1 AMO 1 5.115 -0.236 1.942 0.198400

TER

ATOM 3 C1 MEI 2 0.077 -2.928 -5.293 -0.347900

ATOM 4 I1 MEI 2 -1.125 -3.701 -7.278 0.274700

3. Protonation and solvation of the reactants

Run tleap as:

tleap -s -f tleap.in

where the tleap.in file is as follows

source leaprc.gaff

source leaprc.water.tip3p

loadAmberPrep amo.in

loadAmberPrep mei.in

loadAmberParams amo.dat

loadAmberParams mei.dat

x=loadpdb reactants.pdb

solvateOct x TIP3PBOX 20

saveamberparm x prmtop prmcrd

savepdb x pdbout

quit

The program solvates the solute in an octahedral system with TIP3P [12] water molecules extending the solvent molecules at least 20 angstroms from the reactants.

If you wish to use a solvent other than water, please refer to the Packmol page. It is also necessary to equilibrate the solvated structure in a nonperiodic environment (see below). If you prepared your system using Packmol, do not forget to insert TER between the coordinates of the solvent molecules. The TER can be inserted using pdbanalyser in the following manner, where r1-Methanol.pdb is the name of the PDB file generated by Packmol.

pdbanalyser <<EOF

r1-Methanol.pdb

ter

cubic-box.pdb

EOF

A sample tleap.in for nonstandard solvents is as follows.

source leaprc.gaff

source leaprc.water.tip3p

addPath /home/zafari/Leap

loadAmberPrep eth.in

loadAmberPrep aen.in

#loadAmberParams met.dat

loadAmberParams aen.dat

x=loadpdb cubic-box.pdb

saveamberparm x prmtop prmcrd

savepdb x pdbout

quit

4. Equilibrate the solvated structure with the sander or pmemd module of AMBER as follows

4.1 Equilibrate a periodic system, normally an octahedral box of water molecules

Run MD simulations as follows.

mpirun -bind-to core -np 28 sander.MPI -O -i sander.in00H -o sander.out00H -p prmtop -c prmcrd -r mdrest00H -ref prmcrd

mpirun -bind-to core -np 28 sander.MPI -O -i sander.in00 -o sander.out00 -p prmtop -c mdrest00H -r mdrest00 -ref prmcrd

mpirun -bind-to core -np 28 sander.MPI -O -i sander.in0 -o sander.out0 -p prmtop -c mdrest00 -r mdrest0 -ref prmcrd

mpirun -bind-to core -np 28 sander.MPI -O -i sander.in1 -o sander.out1 -p prmtop -c mdrest0 -r mdrest1 -ref prmcrd

mpirun -bind-to core -np 28 sander.MPI -O -i sander.in2 -o sander.out2 -p prmtop -c mdrest1 -r mdrest2 -ref prmcrd

In the above lines, replace sander with pmemd if you have installed this module. Also, replace '28' with the maximum number of CPU cores of your machine. The followings are the sander.in files.

The sander.in00H file.

Minimization without shake, All hydrogen atoms

&cntrl

irest=0,ntx=1,

nstlim=0,dt=0.0005,

imin=1,maxcyc=1000,drms=0.001,

temp0=300.0,ntt=3,gamma_ln=2.0,

ntc=1,ntf=1,

nsnb=40,cut=8.0,dielc=1.0,

ntpr=100,ntwx=0,ntwv=0,ntwe=0,

ntb=1,ntp=0,taup=0.2,

ipol=0,igb=0,

ncyc=10,ntmin=1,dx0=0.01,

ntr=1,restraint_wt=1000,

restraintmask=" !@H="

&end

The sander.in00 file

Minimization without shake.

&cntrl

irest=0,ntx=1,

nstlim=0,dt=0.0005,

imin=1,maxcyc=1000,drms=0.001,

temp0=300.0,ntt=3,gamma_ln=2.0,

ntc=1,ntf=1,

nsnb=40,cut=8.0,dielc=1.0,

ntpr=100,ntwx=0,ntwv=0,ntwe=0,

ntb=1,ntp=0,taup=0.2,

ipol=0,igb=0,

ncyc=10,ntmin=1,dx0=0.01,

&end

The sander.in0 file

10 ps MD simulation without shake

&cntrl

irest=0,ntx=1,

nstlim=20000,dt=0.0005,

imin=0,maxcyc=100,drms=0.001,

temp0=300.0,ntt=3,gamma_ln=2.0,

ntc=1,ntf=1,

nsnb=40,cut=8.0,dielc=1.0,

ntpr=1000,ntwx=0,ntwv=0,ntwe=0,

ntb=1,ntp=0,taup=0.2,

ipol=0,igb=0,

ncyc=10,ntmin=1,dx0=0.01,

&end

The sander.in1 file

1 ns equilibration with shake and constant volume

&cntrl

irest=1,ntx=5,

nstlim=500000,dt=0.002,

temp0=300.0,ntt=3,gamma_ln=2.0,

ntc=2,ntf=2,

nsnb=15,cut=8.0,dielc=1.0,

ntpr=1000,ntwx=0,ntwv=0,ntwe=0,

ntb=1,ntp=0,taup=0.2,

ipol=0,igb=0,

&end

The sander.in2 file

1 ns simulated annealing at constant pressure

&cntrl

irest=1,ntx=5,

nstlim=500000,dt=0.002,

temp0=300.0,ntt=3,gamma_ln=2.0,

ntc=2,ntf=2,

nsnb=15,cut=8.0,dielc=1.0,

ntpr=1000,ntwx=0,ntwv=0,ntwe=0,

ntb=2,ntp=1,pres0=1.0,taup=1.0,

ipol=0,igb=0,

nmropt=1,

&end

&wt type='TEMP0',istep1= 0,istep2=160000,value1=370.00,value2=370.00 &end

&wt type='TEMP0',istep1=160001,istep2=500000,value1=370.00,value2= 0.00 &end

&wt type='TAUTP',istep1= 0,istep2=200000,value1= 0.20,value2= 0.20 &end

&wt type='TAUTP',istep1=200001,istep2=320000,value1= 1.00,value2= 1.00 &end

&wt type='TAUTP',istep1=320001,istep2=400000,value1= 0.50,value2= 0.50 &end

&wt type='TAUTP',istep1=400001,istep2=500000,value1= 0.05,value2= 0.05 &end

&wt type='END' &end

4.2 Equilibrate a non-periodic system, normally a PDB structure from Packmol

Run MD simulations as follows.

mpirun -bind-to core -np 4 sander.MPI -O -i minH.in -o minH.out -p prmtop -c prmcrd -r minH.mdrest -ref prmcrd

ambpdb -p prmtop -c minH > em.pdb

mpirun -bind-to core -np 4 sander.MPI -O -i heat1.in -o heat1.out -p prmtop -c minH.mdrest -r heat1.mdrest

ambpdb -p prmtop -c heat1.mdrest > heat1.pdb

mpirun -bind-to core -np 4 sander.MPI -O -i heat2.in -o heat2.out -p prmtop -c heat1.mdrest -r heat2.mdrest

ambpdb -p prmtop -c heat2.mdrest > heat2.pdb

mpirun -bind-to core -np 4 sander.MPI -O -i heat3.in -o heat3.out -p prmtop -c heat2.mdrest -r heat3.mdrest

ambpdb -p prmtop -c heat3.mdrest > heat3.pdb

mpirun -bind-to core -np 4 sander.MPI -O -i heat4.in -o heat4.out -p prmtop -c heat3.mdrest -r heat4.mdrest

ambpdb -p prmtop -c heat4.mdrest > heat4.pdb

mpirun -bind-to core -np 4 sander.MPI -O -i heat5.in -o heat5.out -p prmtop -c heat4.mdrest -r heat5.mdrest

ambpdb -p prmtop -c heat5.mdrest > heat5.pdb

mpirun -bind-to core -np 4 sander.MPI -O -i heat6.in -o heat6.out -p prmtop -c heat5.mdrest -r heat6.mdrest

ambpdb -p prmtop -c heat6.mdrest > heat6.pdb

mpirun -bind-to core -np 4 sander.MPI -O -i equil.in -o equil1.out-p prmtop -c heat6.mdrest -r equil.mdrest -ref prmcrd -x equil1.nc

ambpdb -p prmtop -c equil.mdrest > equil.pdb

For a non-periodic system, you can only use the sander module. Also, replace '4' with the maximum number of CPU cores of your machine. The followings are the sander.in files.

The minH.in file.

Optimizing_H-atoms

&cntrl

irest=0,ntx=1,

nstlim=0,dt=0.0001,

imin=1,maxcyc=10000,drms=0.001,

temp0=300.0,ntt=1,tautp=0.2,

ntc=1,ntf=1,

nsnb=40,cut=999.0,dielc=1.0,

ntpr=100,ntwx=0,ntwv=0,ntwe=0,ntwr=10,

ntb=0,ntp=0,taup=0.2,

ipol=0,igb=0,

ncyc=10,ntmin=1,dx0=0.01,

ibelly=1,bellymask="@H= "

&end

The heat1.in file

Stage 1 heating of TC5b 0 to 50K

&cntrl

imin=0, irest=0, ntx=1,

nstlim=10000, dt=0.0005,

ntc=2, ntf=2,

ntt=1, tautp=0.2,

tempi=0.0, temp0=50.0,

ntpr=50, ntwx=50,

ntb=0, igb=0,

cut=999,rgbmax=999.

The heat2.in file

Stage 1 heating of TC5b 50K to 100K

&cntrl

imin=0, irest=1, ntx=5,

nstlim=10000, dt=0.0005,

ntc=2, ntf=2,

ntt=1, tautp=0.2,

tempi=50.0, temp0=100.0,

ntpr=50, ntwx=50,

ntb=0, igb=0,

cut=999,rgbmax=999.

The heat3.in file

Stage 1 heating of TC5b 100 to 150K

&cntrl

imin=0, irest=1, ntx=5,

nstlim=10000, dt=0.0005,

ntc=2, ntf=2,

ntt=1, tautp=0.2,

tempi=100.0, temp0=150.0,

ntpr=50, ntwx=50,

ntb=0, igb=0,

cut=999,rgbmax=999.

The heat4.in file

Stage 1 heating of TC5b 150 to 200K

&cntrl

imin=0, irest=1, ntx=5,

nstlim=10000, dt=0.0005,

ntc=2, ntf=2,

ntt=1, tautp=0.2,

tempi=150.0, temp0=200.0,

ntpr=50, ntwx=50,

ntb=0, igb=0,

cut=999,rgbmax=999.

The heat5.in file

Stage 1 heating of TC5b 200 to 250K

&cntrl

imin=0, irest=1, ntx=5,

nstlim=10000, dt=0.0005,

ntc=2, ntf=2,

ntt=1, tautp=0.2,

tempi=200.0, temp0=250.0,

ntpr=50, ntwx=50,

ntb=0, igb=0,

cut=999,rgbmax=999.

The heat6.in file

Stage 1 heating of TC5b 250 to 300K

&cntrl

imin=0, irest=1, ntx=5,

nstlim=10000, dt=0.0005,

ntc=2, ntf=2,

ntt=1, tautp=0.2,

tempi=250.0, temp0=300.0,

ntpr=50, ntwx=50,

ntb=0, igb=0,

cut=999,rgbmax=999.

The equil.in file

Stage 2 equilibration 1 0-2ns

&cntrl

irest=1, ntx=5,

nstlim=1000000, dt=0.002,

ntc=2, ntf=2,

ntt=1, tautp=0.2,

tempi=325.0, temp0=325.0,

ntpr=500, ntwx=500,

ntb=0, igb=0,

cut=999,rgbmax=999,

ntr=1,restraint_wt=0.1,

restraintmask='!@H=',

&end

If the above procedure does not work, you can download the following sander.in and bash files and equilibrate your structures using them.

https://drive.google.com/drive/folders/1DnOsXWLk1yx-iLsqylo7jhhEGwa3GVOO?usp=sharing

5. convert the equilibrated system to a PDB file

The final coordinate of the MD simulations of the octahedral box of water molecules is stored in the mdrest2 file, and we have to convert it to a PDB file. However, the final structure of the simulations of the system from Packmol is stored in the equil.pdb file and does not need the following steps.

5.1 Execute the following command. It wraps the water molecules in the mdrest2 file into another file called mdrest2-wrap.

cpptraj prmtop ptraj.in

The ptraj.in file is as follows, in which "1" is the residue number of the central residue in the pdbout file.

The ptraj.in file

trajin mdrest2

trajout mdrest2-wrap restart

image origin center familiar com :1

You can obtain the residue number of the central residue by running the command "cap center" in changepdb. A graphical program such as RasMol can also be used to find the central residue.

5.2 Convert the mdrest2-wrap file to PDB format by ambpdb or changeparm

execute

ambpdb -p prmtop -c mdrest2-wrap > equilibrated_system.pdb

changeparm <<EOF

prmtop

equilibrated_system.pdb

mdrest2-wrap

EOF

6. Convert the topology file to the ORCA format

orca_mm -convff -AMBER prmtop

7. Set up QM/MM calculations

To perform QM/MM calculations by ORCA we need to assign the QM system, the active region and the extension shell. The QM system is treated by a DFT method, the active region includes the QM system and a sphere of atoms around the QM system which is optimized by MM methods. Finally, the extension shell is those atoms around the active region whose positions are fixed during optimization however their point charges are included in the QM/MM calculations. The rest of the system is ignored in the calculations.

In the QM system, there are a limited number of atoms that can be assigned by the user. Nevertheless, the active region can contain hundreds of atoms, and it is tedious to assign such a system by hand. To accomplish this, we developed pdbanalyser (a homemade program) to set up QM/MM calculations through ORCA. You need only a text file containing the atom numbers of the atoms you intend to include in the QM system (s1 by default) and a PDB file containing the coordinates of the system (equilibrated_system.pdb). The atom numbers in the s1 file correspond to those in the PDB file. This program converts the numbers into ORCA format (subtracting each number by one since atom numbering in ORCA begins at 0 instead of 1). The following is the s1 file for our example. In this instance, it is simply indicating that the QM system includes atoms 1-9. It is important to note that ORCA automatically identifies junction atoms and the user does not need to specify which atoms are junctions.

The s1 file.

1-9

Run the pdbanalyser program and enter the name of your PDB file to set up the QM/MM calculations. The program reads the PDB file and asks for a command. Execute the qmmm command. Next, the program asks for a file containing atom numbers from the PDB file that you wish to include in the QM system (the s1 file). Afterward, it asks for the radius of the active region surrounding the QM system (the default value is 6 angstroms) and includes all atoms or groups within the desired distance in the region. Further, it asks for the radius of the extension shell around the active region, the level of theory you wish to use in the QM part of the QM/MM calculation, the number of CPU cores, charge, and multiplicity of the QM region. In the end, a file (qmmm.inp by default) is generated, which can be used to perform QM/MM calculations.

Finally, run the following command to perform the QM/MM calculations.

orca qmmm.inp > qmmm.out

A sample qmmm.inp file is shown in the following.

# ORCA's input file for performing QM/MM calculations.

# This file is generated by a home-made Python program (pdbanalyser) by Mehdi and Maryam.

# The program is for free. However, cite it as follows, if you use it in your project.

# DOI: ................

# The working directory is: D:\Dropbox\python\uni\z-PDBanalyser

# The setup time and date are 18:24:26 and 2023-04-26, respectively.

! RIJCOSX TPSSh def2-SVP def2/J D3BJ TIGHTSCF Opt

!QMMM

%qmmm

ORCAFFFilename "prmtop.ORCAFF.prms"

# Excluding junctions, 12 atoms are included in the QM system, as follows.

QMatoms {386 388:398 } end

Use_QM_InfoFromPDB false

# The thickness of the active region is 6.0 angstroms (from the QM region). A total of 242 atoms are included in the active region, as follows.

ActiveAtoms {234:235 238 329:331 335 337:341 344:346 348:362 364:365 367 375:402 407 420:422 424:433 435:436 439:440 442:447 2072 2075:2080 2082:2100 2102:2105 2109 2122 2355:2357 2364:2366 2466:2468 2487:2489 2502:2504 2517:2519 2532:2534 2579 2598:2600 2631:2633 2667 2682:2684 2703:2705 2814:2816 2871:2873 3222:3224 3687:3689 4132 4140:4142 4365:4367 4605:4607 4693 5186 5210 5817:5819 5861 6138 6140 6255:6257 6507:6509 6706:6707 6900:6902 7236:7238 7692 7694 7818:7820 7866:7868 8902 9054:9056 9561:9563 9933:9935 10377:10379 10404:10406 11094:11096 11442:11444 11784:11786 11970:11972 12726:12728 13326 13597 13761:13762} end

Use_Active_InfoFromPDB false

# Non-active atoms that have a distance of less than the sum of their VDW radii plus Dist AtomsAroundOpt will be included in the extension shell. This is the default option of ORCA.

end

% pal

nprocs

end

*pdbfile 0 1 out_pdb_ORCA_QMMM.pdb

8. References

[1] J. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, Development and testing of a general Amber force field, J. Comput. Chem. 25 (2004) 1157–1174. https://doi.org/10.1002/jcc.20035.

[2] C. Lee, W. Yang, R. Parr, G, Development of the Colic-Salvetti correlation-energy into a functional of the electron density, Am. Phys. Soc. 37 (1988) 785–789.

[3] A.D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A. 38 (1988) 3098–3100. https://doi.org/10.1103/PhysRevA.38.3098.

[4] A.D. Becke, Density-functional thermochemistry. III. The role of exact exchange, J. Chem. Phys. 98 (1993) 5648–5652. https://doi.org/10.1063/1.464913.

[5] G.A. Petersson, M.A. Al‐Laham, A complete basis set model chemistry. II. Open‐shell systems and the total energies of the first‐row atoms, J. Chem. Phys. 94 (1991) 6081–6090. https://doi.org/10.1063/1.460447.

[6] G.A. Petersson, A. Bennett, T.G. Tensfeldt, M.A. Al-Laham, W.A. Shirley, J. Mantzaris, A complete basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first-row elements, J. Chem. Phys. 89 (1988) 2193–2218. https://doi.org/10.1063/1.455064.

[7] B.H. Besler, K.M. Merz, P.A. Kollman, Atomic charges derived from semiempirical methods, J. Comput. Chem. 11 (1990) 431–439. https://doi.org/10.1002/jcc.540110404.

[8] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, G.A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V Marenich, J. Bloino, B.G. Janesko, R. Gomperts, B. Mennucci, H.P. Hratchian, J. V Ortiz, A.F. Izmaylov, J.L. Sonnenberg, Williams, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V.G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M.J. Bearpark, J.J. Heyd, E.N. Brothers, K.N. Kudin, V.N. Staroverov, T.A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A.P. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, J.M. Millam, M. Klene, C. Adamo, R. Cammi, J.W. Ochterski, R.L. Martin, K. Morokuma, O. Farkas, J.B. Foresman, D.J. Fox, Gaussian 16, (2016).

[9] C.I. Bayly, P. Cieplak, W. Cornell, P.A. Kollman, A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model, J. Phys. Chem. 97 (1993) 10269–10280. https://doi.org/10.1021/j100142a004.

[10] P.J. Hay, W.R. Wadt, Ab initio effective core potentials for molecular calculations. Potentials for the transition metal atoms Sc to Hg, J. Chem. Phys. 82 (1985) 270. https://doi.org/10.1063/1.448799.

[11] M.D. Hanwell, D.E. Curtis, D.C. Lonie, T. Vandermeersch, E. Zurek, G.R. Hutchison, Avogadro: an advanced semantic chemical editor, visualization, and analysis platform, J. Cheminformatics 2012 41. 4 (2012) 1–17. https://doi.org/10.1186/1758-2946-4-17.

[12] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys. 79 (1983) 926–935. https://doi.org/10.1063/1.445869.