Expanding the range of force fields available for ONIOM calculations

Subscriber access provided by Queen Mary, University of London

Computational Chemistry

Expanding the range of force fields available for ONIOM calculations: the SICTWO interface W. M. Chamil Sameera, and Feliu Maseras J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00332 • Publication Date (Web): 23 Aug 2018 Downloaded from http://pubs.acs.org on August 26, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Expanding the Range of Force Fields Available for ONIOM Calculations: The SICTWO Interface W. M. C. Sameera, * Feliu Maseras * a,b

a

a,c

Institute for Chemical Research of Catalonia (ICIQ), The Barcelona Institute of Science and

Technology, Avgda. Països Catalans 16, 43007, Tarragona, Spain. b

Institute of Low Temperature Science, Hokkaido University, Sapporo, Hokkaido, 060-0819,

Japan. c

Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain

E-mail: [email protected], [email protected]

KEYWORDS QM/MM methods, ONIOM, DFT, molecular mechanics, homogeneous catalysis, enantioselectivity. ABSTRACT The ONIOM scheme is one of the most popular QM/MM approaches, but its extended application has been so far hindered by the limited availability of force fields in most practical implementations. This paper describes a simple software code to overcome this limitation, and its application to three representative chemical problems. The "Shell Interface for Combining Tinker With ONIOM" (SICTWO) program gives access to all force fields available in the Tinker molecular mechanics program from a Gaussian09 or Gaussian16 calculation. The first application presented is the geometry optimization of dithienobicyclo[4.4.1]-undeca-3,8-diene-11-one ethylene glycol ketal. A variety of force fields were tested

ACS Paragon Plus Environment

1

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 26

for the MM part, and molecular structure with ONIOM(B3LYP:OPLS-AA) method shows a good agreement with the experimental results. The second problem is related with enantioselectivity of the vanadium-catalyzed asymmetric oxidation of 1,2-bis(tert-butyl)disulfide. ONIOM(B3LYP:MM3) results, obtained with SICTWO, are consistent with those of a previous IMOMM(B3LYP:MM3) study. In the third application, we have combined AMOEBA09 polarizable force field with ONIOM to study a water molecule binding on ice. Calculated ONIOM(ωB97X-D:AMOEBA09) binding energies are consistent with the ωB97X-D results. These studies show SICTWO as an easy-to-use tool that can further expand the use of the multiscale ONIOM method within computational chemistry and computational biology.


2

Page 3 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION Development of computational methods for the study of complex molecular systems has been a strong interest in chemistry. In 2013, Karplus, Levitt, and Warshel won the Nobel prize in chemistry for the development of computational methods using quantum mechanics and classical mechanics.

1,2,3

Their pioneering work was the basis for the modern hybrid

quantum mechanics/molecular mechanics (QM/MM) methods, and a number of research groups have made significant contributions for QM/MM method developments.

4,5,6,7,8,9,10

QM/MM methods are implemented in commercially available computational chemistry

QM/MM

software, and successfully applied for modeling organic, organometallic, biochemical, and solid-state systems.

9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28

QM&treatment [describe&chemical&processes] –&transition&metal¢re& –&reaction¢re/active&site& MM&treatment [describe&environmental&effects] (a) (b) W.&M.&C.&Sameera,&F.&Maseras,& Figure 1. (a) Partitioning a molecular system into QM–&bulky&groups/ligands and MM regions (the black ring is WIREs&Comput.&Mol.&Sci.&Wiley&VCH,&2012.& the boundary region) (b) QM/MM molecular model of–&enzyme&structure Tp Cu catalyst (‘ball and stick’ W.&M.&C.&Sameera,&F.&Maseras,& –&explicit&solvent&molecules Phys.&Chem.&Chem.&Phys.&2011,13,&10520. model is the QM region and ‘wireframe’ is the MM region). (One&of&PCCP&most&read&article&in&2011) tBu,Me

25

QM/MM methods partition the molecular system into two regions, where the

Tuesday, April 3, 12

electronically important region is described by a QM approach, while the remainder part is described by a MM method (Figure 1). The total QM/MM energy of a given molecular system can be expressed as; 𝐸!"#$% QM: MM = 𝐸!" QM + 𝐸!! MM + 𝐸!"# QM: MM

….. (1)


3


Page 4 of 26

In this formula, labels in subscript describe the computational method, and labels in parenthesis is the partition of a molecular system. The interaction term, 𝐸!"# QM: MM , includes both bonded (stretching, bending, torsional) and non-bonded interactions between the QM and MM regions (E.g. van der Waals, electrostatic). These weak interactions can be described by QM or MM methods. The particular form of this 𝐸!"# QM: MM term is not defined a priori, and different descriptions give rise to different QM/MM methods. In the simplest case scenario, interactions between QM and MM region are calculated by MM, the so-called mechanical embedding (ME). Another approach is to include the partial charges of the MM region in the QM Hamiltonian, and admits charge distribution of the MM region to interact with that of the QM region. This approach is called electronic embedding (EE).

26

More accurate electrostatic interactions between the QM and MM regions may be calculated by using the EE. An interesting concept within QM/MM methods is the subtractive approach,

17,27

where

the energy of the MM region [𝐸!! MM ] and the QM/MM interaction term [𝐸!"# QM: MM ] can be approached by subtracting internal energy of the ‘model system’ and the ‘real system’ [i.e. 𝐸!! QM: MM − 𝐸!! QM ]. Therefore, the subtractive methods do not need an additional coupling Hamiltonian to represent the interactions between QM and MM regions, which makes them easier to implement. The most popular subtractive QM/MM method is the ONIOM (our Own N-layer Integrated molecular Orbital molecular Mechanics),

6,8,11,15,17

and this

is an evolution of the previously proposed IMOMM (Integrated Molecular Orbital Molecular Mechanics) and IMOMO (Integrated Molecular Orbital Molecular Orbital) methods. 6

Apart from the interactions between the QM and MM regions, the second important aspect in QM/MM methods is the boundary region when there are bonds across the QM/MM partition. In the ONIOM approach, link atoms (LA) are used, where the dangling bonds from the QM region are capped with hydrogen atoms. This LA is placed on the line between the link atom connection (LAC) and substituted atom(s) (the link atom host, LAH). Distance


4

Page 5 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


between LAC and LA is determined by scaling the LAC-LAH distance with a constant factor that depends on the atom type. By applying these simple concepts of subtractive approach and link atoms, the total energy (E), gradient (G), and Hessian (H) of a two-layer ONIOM(QM:MM) method can be defined rigorously as; 𝐸!"#!$ QM: MM = 𝐸!" QM + 𝐸!! QM: MM − 𝐸!! QM

…(2)

𝐺!"#!$ QM: MM = 𝐺!" QM × 𝐽 + 𝐺!! QM: MM − 𝐺!! QM × 𝐽

…(3)

𝐻!"#!$ QM: MM = 𝐽! × 𝐻!" QM × 𝐽 + 𝐻!! QM: MM − 𝐽! × 𝐻!! QM × 𝐽 …(4) where 𝐽 =

!(!"#$ !""#$.) !(!"#$% !""#$.)

The Jacobian, J, projects the gradients of the link atoms onto the LAH and LAC coordinates, and is the identity matrix for atoms which are not in the boundary region. The ONIOM formalism is thus quite simple if QM and MM values are available for energy, gradients, and Hessian. ONIOM approach has been successfully applied by many researchers, but further applications of ONIOM have been hindered by the relatively scarce availability of force fields in the common implementations. For example, the ONIOM(QM:MM) implementation in the commonly used Gaussian09 program has only 28

three force field parameter sets, namely UFF, AMBER and Dreiding. Therefore, development of an ONIOM(QM:MM) implementation to support various polarizable and non-polarizable force fields become critical. In this work, we present the Shell Interface for Combining Tinker with ONIOM (SICTWO) program. The ONIOM(QM:MM) implementation in SICTWO uses the Gaussian09 standardized interface to obtain MM energy and derivatives from Tinker MM program. It is a subroutine (sictwo) that called from the Gaussian09 program, and the same 29,30

subroutine calls the Tinker subroutines. Importantly, this interface does not require modifications to either the Gaussian09 or Tinker original codes, thus being portable and adaptable to future program modifications. All the results reported in the manuscript have been indeed obtained with Gaussian09, but we have confirmed that the interface is also operative with Gaussian16. This interface opens widely the range of force fields for ONIOM,


5


Page 6 of 26

as Tinker supports different versions of a variety of popular force fields: MM2, MM3, CHARMM, OPLS, Merck Molecular Force Field (MMFF), Liam Dang's polarizable model, AMOEBA polarizable atomic multipole force field, AMBER, etc. Therefore, SICTWO allows a substantial extension of the applications ONIOM(QM:MM) method. In this paper, we describe algorithmic and software implementation in SICTWO, as well as a selected set of applications to prove its performance. We published an initial application of this interface with the AMOEBA09 force field to the binding properties of small radicals on crystalline water ice.

31

RESULTS AND DISCUSSION Implementation of the SICTWO interface Figure 2 summarizes the operation mode of the SICTWO interface for the connection between the QM program (Gaussian09) and the MM program (Tinker). The use of SICTWO requires access to working versions of the Gaussian09 and Tinker. No modifications are required in either the Gaussian09 or Tinker codes. The current version of the modules of SICTWO is provided in the Supporting Information (Section 1). Further updates will be kept in a website. All modules of SICTWO program are in the “sictwo_files” directory, and this is kept in

the

working

directory

of

user.

SICTWO

is

activated

by

using

the

external="sictwo_files/sictwo" keyword in the Gaussian09 input file. This keyword requests Gaussian09 to run SICTWO interface that facilitates the use of Tinker for low-level calculations in ONIOM. In the presence of the external keyword, Gaussin09 creates a text file with the EIn extension, and this file consists of Cartesian coordinates, atom types, connectivity details of current structure. The modules of SICTWO then: (i) transforms the information in the EIn file into a format suitable for the Tinker program, (ii) converts MM atom types for Tinker, (iii) calls different Tinker modules to perform MM calculation, and (iv) transforms the Tinker output into a format suitable for Gaussian09, a file with an EOu


6

Page 7 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


extension. Then, Gaussian09 proceeds with its normal operations (single point, geometry optimization, frequency calculation, ...), incorporating the MM energy and derivatives computed with the Tinker program.

Figure 2. Flux diagram for the role of the SICTWO interface between the Gaussian09 and Tinker programs. All the information transfer between Gaussian09, Tinker, and the SICTWO interface takes place through text files containing the different pieces of information. This ensures an easy portability of the code, and an easy adaptation to future evolutions of either the QM or the MM codes, as far as no fundamental changes are made to the program structure. It also suggests a rather simple extension to other codes that may be able to carry out ONIOM calculations. Most of the operation of SICTWO takes place as a black box for the user, which does not need to be aware of the temporary files that are deleted at the end of the calculation. Most of the communication between user and program takes place through the Gaussian09 input and output files. The Gaussian09 input file contains the usual information, such as method, type of calculation, Cartesian coordinates, etc. User must pay special attention to three aspects: (i) the keyword opt=nomicro is mandatory, as the present


7


Page 8 of 26

implementation does not support microiterations, (ii) geometry=connectivity is strongly encouraged, as should be the case for most MM calculations where better control of the connectivity being used, (iii) special care must be taken with the definition of the atom types, which is included in the mm.type file of SICTWO. The GaussView interface can be in this sense very useful to provide at least an initial version of the atom types (e.g. UFF) and connectivity, which can be checked by user. Besides the Gaussian09 input and output files, there are two additional sources of communication between user and the program. The first of them is the inp.par file that contains information to be used by the Tinker program, most importantly the location of the file with the force field parameters. Also, user can define missing potential parameters in the inp.par file. The same inp.par file can in principle be used for a whole set of calculations on a given system. The modules that constitute the SICTWO interface are the following (all scripts can be found in the Supporting information, Section 1):

sictwo (or sictwo_g09d01) modules: This is the master module of SICTWO, and is a bash script that receives control from the QM code to perform several tasks; (i) identify the file name of EIn files of Gaussian09, (ii) execute make_input and make_out modules of SICTWO (see below for more details), (iii) execute analyze, testgrad, and testhess modules of Tinker. These Tinker modules provide MM energy and derivatives of both real (large) and model (small) systems for ONIOM calculation. In particular, SICTWO obtains MM energy and MM dipole moment from the analyze module of Tinker, testgrad provides MM gradients, and MM force constants can be obtained by the testhess module.

make_input (or make_input_g09d01 ) module: This is a perl script that refers Gaussian09 EIn files to prepare input files for Tinker program with xyz extension. Tinker input file contains the number of atoms in the model, atom name, atomic X-, Y-, and Zcoordinates, atom types, and the connectivity details. MM atom type conversions between Gaussian09 and Tinker are included in mm.type file of SICTWO.


8

Page 9 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


make_output module: This perl script collects MM energy and derivatives from Tinker MM calculations, and converts them the Gaussian09 standard text file with the EOu extension. This file contains energy and dipole moment, gradients, polarizability, dipole derivatives, and force constants. Some tuning by the user may be required in the mm.type and inp.par files when dealing with unusual atom types or connectivity. This is however limited to the first calculation on a given system, as the same files can be used in all calculations of that system, including full reaction profiles and isomeric modifications.

Applications (a) Optimization of RESVAN with different force fields We have previously reported the performance of QM/MM methods vs DFT methods for the description of systems where dispersion is expected to be important. We have shown that the QM/MM methods are better in this case because they consist of dispersion in the MM part. One of the criteria we used in that study was the prediction of the structure of 32

dithienobicyclo-[4.4.1]-undeca-3,8-diene-11-one ethylene glycol ketal system (CSD entry: RESVAN, Figure 3). The structure of this molecule is very sensitive to intramolecular dispersion correlations of thiophene rings. The reported X-ray structure showed the sulfursulfur distance of 4.29 Å, and estimated elongation of 0.10 Å due to crystal packing. The 33

34

same molecule was used by Hobza and Grimme for a computational study, where the 35

calculated sulfur-sulfur distance fit with the experimental value at the ab-initio levels: CCSD(T)/6-31G* [r(S–S) = 4.26 Å], MP3/CBS(B) [r(S–S) = 4.22 Å], and SCS-MP2/def2TZVP [r(S–S) = 4.12 Å].


9


Page 10 of 26

Figure 3. QM/MM molecular model of RESVAN system (‘ball and stick’ model is the QM region and ‘wireframe’ model is the MM region. In our previous study, we included one of the two thiophene rings in the MM region 32

(see Figure 3), and optimized the structure at the ONIOM(B3LYP:UFF) level. When we wanted to evaluate other force fields, we had to perform single-point calculations on the structures

resulting

from

a

relaxed

potential

energy

surface

scan

with

the

ONIOM(B3LYP:UFF) method (sulfur-sulfur distance was the reaction coordinate). By using the SICTWO interface, full optimization with other force fields is straightforward. We have performed ONIOM(B3LYP:MM) structure optimizations using the MM2, MM3, MMFF, 36

37

38

and OPLS-AA force fields for the MM description. The mm.type and inp.par files used for 39

this example are provided in the Supporting Information (Section 2). The def2-TZVPP basis 40

sets were used for the QM description. Results are summarized in Table 1 and section 2.2 in Supporting Information, which can be viewed as an updated version of Table 1 in our previous publication. There are minor differences in some of the numbers due to changes in 32

the basis set in the QM calculations, and also to implementation differences between MacroModel and Tinker in the MM calculations. Table 1. Experimental and calculated sulfur–sulfur distance of RESVAN molecule.

X-ray Exp. gas phase (estimated) 33

34

r (S-S)/ Å 4.29 4.19±0.05


10

Page 11 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ONIOM (QM :M M ) optimization

structure

SICTWO ONIOM(B3LYP:MM2) ONIOM(B3LYP:MM3) ONIOM(B3LYP:MMFF) ONIOM(B3LYP:OPLS-AA)

3.97 4.00 4.45 4.17

Gaussian09 ONIOM(B3LYP:UFF)

3.99

M M structure optimization

MM2 MM3 MMFF OPLS-AA UFF

3.95 3.92 4.40 4.09 3.80

The conclusions from Table 1 are the same from the previous work. Dispersion is 32

important to obtain a good structural description, and this can be achieved by ONIOM(B3LYP:MM) when the force field provides a good description, as is the case with OPLS-AA [r(S-S) = 4.17 Å]. This is due to the fact that we have incleded one of two thiophene rings in the MM region, and therefore the interactions between two thiophene rings were calculated by MM. This result emphasized that a careful choice of the force field for ONIOM(QM:MM) calculations give accurate molecular structures at a low-computational cost. In our ONIOM(QM:MM) calculations, computational cost of MM calculations is negligible compared to QM calculations. This example enhances the applications of SICTWO, as changing between the different force fields available in Tinker is straightforward. This only requires minor tweaking in the atom type conversions and defining missing potential parameters.

(b) Computing enantioselectivity with ONIOM(B3LYP:MM3) Our group has been studying enantioselective homogenous catalysis for some time with either IMOMM or ONIOM approaches.

41,42,43,44,45

We want to check whether SICTWO is


11


Page 12 of 26

able to reproduce the results we obtained with IMOMM(B3LYP:MM3) on vanadium catalyzed oxidation of sulphides.

41

R3

*

R2 N

L =

OH OH tBu

tBu

R1

S*

S tBuS

tBuS

VO(acac) 2, L, H 2O 2, rt

O

"R"

Figure 4. Oxidation of 1,2-bis(tert-butyl)-disulfide by a vanadium-catalyst. The reaction under study is the vanadium-catalyzed asymmetric oxidation of 1,2bis(tert-butyl)-disulfide by Cogan et al (Figure 4). This approach was originally developed 46

by Bolm et al. In this method, 1:1 mixture of vanadium(IV)-oxo complex and a Schiff base 47

gives V (O)(L)(OOH) active species, where L isthe Schiff base. The vanadium(V) species V

was characterized by V NMR experiments. In this reaction, hydrogen peroxide acts as the 51

oxidant. Starting from the active vanadium(V) species, oxygen atom of the -OOH group transfers to disulfide. Experimental studies were performed with different substituents for R1, R2, and R3 of the Schiff base, and they have a significant effect on the enantioselectivity. Following our previous work, we have used four different ligands that produce different enantiomeric excess. We label these ligands as system A: R1 = R2 = R3 = tBu; system B: R1 = R3 = tBu, R2 = H; system C: R1 = R2 = tBu, R3 = iPr; system D: R1 = R2 = H, R3 = tBu.


12

Page 13 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. (a) Fully optimized structure of the most stable isomer of V (O)(L)(OOH) active V

intermediate (‘ball and stick’ model is the QM region and ‘wireframe’ model is the MM region), (b) the conformations of the transition state used in this study. For this study, we have applied ONIOM(B3LYP:MM3) approach as implemented in SICTWO. The LanL2DZ basis set for vanadium,

48,49,50

the LanL2DZ basis set and a d

polarization function for sulfur, 6-31G(d) basis sets for oxygen and nitrogen, and 6-31G for carbon and hydrogen were used.

51,52,53

In the MM3 calculations, van der Waals parameters for

vanadium have been obtained from the UFF force field. We have introduced hydrogen bond 54

parameters for the hydrogen atoms of methyl groups and the anionic oxygen atom that transferred to sulfide. The QM/MM partition is indicated in Figure 5a. The mm.type and 41

inp.par files for this example can be found in the Supporting Information (Section 3.1). As is often the case for systems involving enantioselectivity, there are important aspects concerning the isomeric and conformational complexity of the molecules. These two issues affect the discussion below. The first of them is the conformational complexity of the transition state that we have highlighted in Figure 5b. Further, there are two freely rotating single bonds (labeled as x and y), and their rotation produces different local minima of close energy, which have to be explicitly considered. The second issue is that there are two energetically similar diastereomeric forms of the catalyst. The catalyst is generated in situ,


13


Page 14 of 26

and the vanadium becomes a stereogenic center that can be in two different configurations. Several possible conformations of the catalyst and transition states is outlined in the Supporting Information (Section 3.2 and 3.3). We focus here only in the key results regarding the enantiomeric excess of the reaction (Table 2). Table 2. Relative energies (ΔΔE, kcal mol ) of 2aa -1

x1y1z1

, 2aa

x1y1z2

, 2bs

x1y1z2

, and 2ba

x1y1z1

transition states, and experimental and calculated enantiomeric excess at 298.14 K. Values in parenthesis represent the calculated ee% from IMOMM(B3LYP:MM3). A

B

C

D

2aa

x1y1z1

0.0

0.0

0.0

0.0

2aa

x1y1z2

2.2

2.1

2.2

0.3

2bs

x1y1z2

1.2

1.4

0.7

0.2

2ba

x1y1z1

3.6

3.7

2.6

0.5

ee % (calculated)

77 (90)

83 (90)

54 (74)

17 (21)

ee % (experimental)

82

83

60

46

The overall agreement between the computed enantiomeric excess with IMOMM(B3LYP:MM3) and ONIOM(B3LYP:MM3) is very good, and the agreement between both of them and the experimental values is also satisfactory. The minor difference between the two computational values can be explained by minor changes in QM/MM partition, MM3 version, and treatment of the link atoms. The agreement in the ee values is particularly encouraging if one takes into account the exponential dependence of the enantiomeric ratio with respect to the energy difference. The experimental trends in enantiomeric excess are reproduced in the ONIOM(B3LYP:MM3) calculations, and the qualitative explanations outlined in the previous publication are confirmed. 2aa

x1y1z1

is always

the lowest energy transition state, leading to the R enantiomer as the major product. The 2bs

x1y1z2

transition state is the most significant contributor to the S enantiomer. The

discrimination between these two transition states can be attributed to the substituent at the R

3

position, which rules the relative stability of the catalyst isomers. Substituent R rules instead 1


14

Page 15 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the discrimination within one given catalyst, and its nature contributes to the smaller enantiomeric excess for ligand D. The SICTWO interface thus allows the use of the MM3 force field in ONIOM to calculate enantioselectivity. The two obvious alternative approaches have serious shortcomings: ONIOM(B3LYP:UFF) produces bad results, and IMOMM(B3LYP:MM3) was associated with a particular version of Gaussian98, and a very elaborate construction of the input file. In contrast, ONIOM(QM:MM3) implementation in SICTWO requires only the appropriate choice of atom types and missing potential parameters.

(c) ONIOM calculations using the AMOEBA09 polarizable force field In QM/MM computations of large molecular systems, there is a growing interest in computational chemistry field to use polarizable force fields in MM calculations. A number of polarizable force fields, AMOEBA09

55-57

for instance, are under development. However,

their applications are restricted to specialized researchers. In this direction, we herein present the ONIOM(QM:AMOEBA09) approach, which is a simple and accurate approach to study interstellar medium molecular species on water ice.


15


Page 16 of 26

Top view

Side view

Model A

Model B

Figure 6. (a) QM/MM ice cluster models. Each model consists of three water layers. We have applied the ONIOM(QM:AMOEBA09) approach to compute binding energies of a H O molecule on water ice. For this purpose, two large ice cluster models, A 2

and B, were used (Figure 6). Model A holds 162 H O molecules (48 molecules in the QM 2

region and 114 molecules in the MM region), while model B consists of 156 H O molecules 2

(44 molecules in the QM region and 112 molecules are in the MM region). Depending on the number of dangling-H (d-H) and dangling-O (d-O) atoms, four binding sites for model A (A1-A4) and four binding sites for model B (B1-B4) can be formed (See section 4.2 in Supporting Information). The binding site on ice is described by ωB97X-D /def2-TZVP 58

method, while AMOEBA09 polarizable force field was used for the remaining part. ONIOM(ωB97X-D/def2-TZVP:AMOEBA09) optimized ice structures with a H O molecule 2

on the binding site and calculated binding energy are shown in Figure 7.


16

Page 17 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.17 2.00

2.03

1.94

1.90

A2 (0.99)

A1 (0.84)

2.08

2.11

1.98 2.03

2.27 1.92

B1 (0.86)

2.19 1.90

2.35

2.66

1.97

2.78

A3 (0.95)

2.11 2.11

A3 (0.56) 2.07

2.10

1.93

2.12

B2 (0.90)

2.70

B3 (0.74)

B4 (0.71)

Figure 7. ONIOM(ωB97X-D/def2-TZVP:AMBER) optimized structures with a H O 2

molecule on the binding site and calculated binding energies (eV). Only the binding pockets are shown for simplicity. Bond distances (Å) are shown in blue color. Calculated binding energies of H O on ice are is in the range of 0.56 – 99 eV, and H O 2

2

binding is relatively stronger at all binding sites, and is qualitatively in agreement with the experimental

results.

59

Calculated

binding

energies

from

ONIOM(ωB97X-D/def2-

TZVP:AMOEBA09) and ωB97X-D/def2-TZVP (Section 4.3, Supporting Information) found good agreement, where the maximum absolute discrepancy between ONIOM(ωB97XD/def2-TZVP:AMOEBA09) and ωB97X-D/def2-TZVP binding energies is only 0.01 eV. Therefore, we argue that size of the QM region in our QM/MM ice model is good enough to compute binding energies accurately. We have approached to the same conclusion from calculated dipole moments (Section 4.4, Supporting Information). CONCLUSIONS SICTWO is a user-friendly interface that works with the commercial versions of Gaussian09/Gaussian16 and Tinker programs. The ONIOM(QM:MM) approach in SICTWO supports for various polarizable and non-polarizable force fields. Performance of the interface was demonstrated in three practical examples. First, we have proved how different


17


Page 18 of 26

force fields can be used as the MM part in the ONIOM(B3LYP:MM) geometry optimization of a molecular system involving two thiophene units. Then, we have shown how ONIOM(B3LYP:MM3)

rationalized

the

enantioselectivity

of

vanadium-catalyzed

asymmetric oxidation of a disulfide. Finally, we have applied ONIOM(QM:AMOEBA09) approach to calculate binding energy of H O on water ice. These examples evidence that 2

SICTWO is a simple and a very useful computational approach for modeling a wide spectrum of chemical and biological processes in the future.

ASSOCIATED CONTENT Supporting Information. Two files. The first one contains: (1) SICTWO modules, (2) RESVAN system (atom type conversions, defined potential parameters, calculated sulfur– sulfur distance), (3) vanadium-catalyzed asymmetric oxidation of 1,2-bis(tert-butyl)-disulfide [atom type conversions, defined potential parameters, isomers of V (O)(L)(OOH) species, V

transition state structures, ONIOM extrapolated energies of the calculated structures, (4) water ice systems (atom type conversions, description about the dangling atoms, comparison of binding energies, comparison of dipole moments, ONIOM extrapolated energies of the computed structures (5) Cartesian coordinates of the computed structures. The second file contains the manual for the use of the interface.

AUTHOR INFORMATION Corresponding Author W. M. C. Sameera Present Address Institute of Low Temperature Science, Hokkaido University, Sapporo, Hokkaido, 060-0819, Japan.


18

Page 19 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Funding Sources FM would like to acknowledge financial support from the CERCA Programme/Generalitat de Catalunya and MINECO (Project CTQ2017-87792-R). WMCS would like to acknowledge the Japan Society for the Promotion of Science grant (JSPS, No: P14334).

ACKNOWLEDGMENT Prof. Keiji Morokuma (Fukui Institute for Fundamental Chemistry) for his inspiration and mentorship. Dr. Marcus Lundberg (Uppsala University), and Dr. Tsutomu Kawatsu (University of Tokyo), Dr. Charles Goehry (ICIQ), and Dr. Ignacio Funes-Ardoiz (ICIQ) for useful discussions.

References (1) Lifson, S.; Warshel, A. Consistent Force Field for Calculations of Conformations, Vibrational spectra, and Enthalpies of Cycloalkane and n-alkane Molecules. J. Chem. Phys. 1968, 49, 5116-5129. (2) Honig, B.; Karplus, M. Implications of Torsional Potential of Retinal Isomers for Visual Excitation. Nature, 1971, 229, 558-560. (3) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions - Dielectric, Electrostatic and Steric Stabilization of Carbonium-ion in Reaction of Lysozyme. J. Mol.

Biol. 1976, 103, 227-249. (4) Singh, U. C.; Kollman, P. A. Combined Ab-initio Quantum Mechanical and Molecular Mechanical Method for Carrying out Simulations on Complex Molecular Systems -


19


Page 20 of 26

Applications to the CH Cl + Cl Exchange Reaction and Gas-phase Protonation of Polyethers. -

3

J. Comput. Chem. 1986, 7, 718-730. (5) Field, M. J.; Bash, P. A.; Karplus, M. A Combined Quantum Mechanical and Molecular Mechanical Potential for Molecular Dynamics Simulations. J. Comput. Chem. 1990, 11, 700-733. (6) Maseras, F.; Morokuma, K. IMOMM - A New Integrated Ab-Initio Plus Molecular Mechanics Geometry Optimization Scheme of Equilibrium Structures and Transition States.

J. Comput. Chem. 1995, 16, 1170-1179. (7) Gao, J. L. Hybrid Quantum and Molecular Mechanical Simulations: An Alternative Avenue to Solvent Effects in Organic Chemistry. Acc. Chem. Res. 1996, 29, 298-305. (8) Svensson, M.; Humbel, S.; Froese, R. D. J.; Matsubara, T.; Sieber, S.; Morokuma, K. ONIOM: A Multilayered Integrated MO+MM Method for Geometry Optimizations and Single Point Energy Predictions. A Test for Diels-Alder Reactions and Pt(P(t-Bu) ) ) + H 3

2

2

Oxidative Addition. J. Phys. Chem. 1996, 100, 19357-19363. (9) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem. Int.

Ed. Engl. 2009, 48, 1198-1229. (10) Lin, H.; Zhang, Y.; Pezeshki, S.; Truhlar, D. G.; QMMM-version 1.3.8; University of Minnesota, Minneapolis, 2009, based on Frisch, M. J. et. al. Gaussian; Gaussian, Inc.; Pittsburgh PA, 2003, based on TINKER 4.2 by Ponder, J. W. TINKER–version 4.2, Washington University, St. Louis, MO, 2004, and based on Rodgers, J. M.; Lynch, B. J.; Fast, P. L.; Chuang, Y.-Y.; Pu, J.; Zhao, Y.; Truhlar, D. G.; MULTILEVEL-version 3.1/G03; University of Minnesota, Minneapolis, 2002. (11) Vreven, T.; Morokuma, K. Hybrid Methods: ONIOM(QM:MM) and QM/MM. Annu.

Rep. Comput. Chem. 2006, 2, 35–51. (12) Maseras, F. The IMOMM Method Opens the Way for the Accurate Calculation of "Real" Transition Metal Complexes. Chem. Commun. 2000, 1821-1827. (13) Ujaque, G.; Maseras, F. Applications of Hybrid DFT/Molecular Mechanics to Homogeneous Catalysis. Struct. Bonding, 2004, 112, 117-149.


20

Page 21 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(14) Bo, C.; Maseras, F. QM/MM Methods in Inorganic Chemistry. Dalton Trans, 2008, 2911-2919. (15) Chung, L. W.; Hirao, H.; Li, X.; Morokuma, K. The ONIOM Method: Its Foundation and Applications to Metalloenzymes and Photobiology. WIREs Comput. Mol. Sci. 2012, 2, 327-350. (16) Sameera, W. M. C.; Maseras, F. Transition Metal Catalysis by Density Functional Theory and Density Functional Theory/Molecular Mechanics. WIREs Comput. Mol. Sci. 2012, 2, 375-385. (17) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z. F.; Liu, F. Y.; Li, H.-B.; Ding, L. N.; Morokuma, K. The ONIOM Method and its Applications. Chem. Rev. 2015, 115, 5678-5796. (18) Jover, J.; Maseras, F. QM/MM Calculations on Selectivity in Homogenenous Catalysis.

Struct. Bond. 2016, 167, 59-79. (19) Wang, B.; Wang, X.; Wang, W.; Liu, F. Exploring the Mechanism of Fluorescence Quenching and Aggregation-Induced Emission of a Phenylethylene Derivative by QM (CASSCF and TDDFT) and ONIOM (QM:MM) Calculations, J. Phys. Chem. C, 2016, 120, 21850–21857. (20) Biancardi, A.; Barnes, J.; Caricato, M. Point charge embedding for ONIOM excited states calculations, J. Chem. Phys. 2016, 145, 224109-224118. (21) Menger, M. F. S. J.; Caprasecca, S.; Mennucci, B. Excited-State Gradients in Polarizable QM/MM Models: An Induced Dipole Formulation, J. Chem. Theory Comput. 2017, 13, 3778–3786. (22) Roßbach, S.; Ochsenfeld, C. J. Influence of Coupling and Embedding Schemes on QM Size Convergence in QM/MM Approaches for the Example of a Proton Transfer in DNA, J.

Chem. Theory Comput. 2017, 13, 1102–1107. (23) Fernandes, H. S.; Ramos, M. J.; Cerqueira, N. M. F. S. A. molUP: A VMD Plugin to Handle QM and ONIOM Calculations Using the Gaussian Software, J. Comput. Chem. 2018, 39, 1344–1353.


21


Page 22 of 26

(24) Ramozzi, R.; Sameera, W. M. C.; Morokuma, K. Predicting Reaction Pathways from

Reactants, Applied Theoretical Chemistry, World scientific publishing, 2018, PP. 321-349. (25) Muñoz-Molina, J. M.; Sameera, W. M. C.; Álvarez, E.; Maseras. F.; Belderrain, T. R.; Pérez, P. J. Mechanistic and Computational studies of the Atom Transfer Radical Addition of CCl to Styrene Catalyzed by Copper Homoscorpionate Complexes, Inorg. Chem. 2011, 50, 4

2458-2467. (26) Vreven, T.; Byun, K. S.; Komaromi, I.; Dapprich, S.; Montgomery, J. A.; Jr, Morokuma. K.; Frisch, M. J. Combining Quantum Mechanics Methods with Molecular Mechanics Methods in ONIOM, J. Chem. Theory. Comput. 2006, 2, 815-826. (27) Cao, L. L.; Ryde, U. On the Difference Between Additive and Subtractive QM/MM Calculations. Front. Chem. 2018, 6, 1-15. (28) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian, Inc., Wallingford CT, 2009. (29) Ponder, J. W.; Richards, F. M.; An Efficient Newton‐like Method for Molecular Mechanics Energy Minimization of Large Molecules. J. Comput. Chem. 1987, 8, 10161024. (30) Kundrot, C. E.; Ponder, J. W.; Richards, F. M. Algorithms for Calculating Excluded Volume and its Derivatives as a Function of Molecular Conformation and Their Use in Energy Minimization. J. Comput. Chem. 1991, 12, 402-409.


22

Page 23 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(31) Sameera, W. M. C.; Senevirathne, B.; Andersson, S.; Maseras, F.; Nyman, G. ONIOM(QM:AMOEBA09) Study on Binding Energies and Binding Preference of OH, HCO and CH Radicals on Hexagonal Water Ice (I ). J. Phys. Chem. C. 2017, 121, 15223-15232. 3

h

(32) Sameera, W. M. C.; Maseras, F. Quantum Mechanics/Molecular Mechanics Methods can be More Accurate than Full Quantum Mechanics in Systems Involving Dispersion Correlations. Phys. Chem. Chem. Phys. 2011, 13, 10520–10526. (33) Knoblock, K. M.; Silvestri. C. J.; Collard, D. M. Stacked Conjugated Oligomers as Molecular Models to Examine Interchain Interactions in Conjugated Materials, J. Am. Chem.

Soc. 2006, 128, 13680-13681. (34) Moellmann, J.; Grimme, S. Importance of London Dispersion Effects for the Packing of Molecular Crystals: A Case Study for Intramolecular Stacking in a bis-thiophene Derivative.

Phys. Chem. Chem. Phys. 2010, 12, 8500-8504. (35) Pluhackova, K.; Grimme, S.; Hobza, P. On the Importance of Electron Correlation Effects for the Intramolecular Stacking Geometry of a bis-thiophene Derivative. J. Phys.

Chem. A, 2008, 112, 12469–12474. (36) Allinger, N. L. Conformational Analysis. 130. MM2 - Hydrocarbon Force-Field Utilizing V1 and V2 Torsional Terms. J. Am. Chem. Soc. 1977, 99, 8127-8134. (37) Allinger, N. L.; Yuh Y. H.; Lii, J.-H. Molecular Techanics - the MM3 Force Field for Hydrocarbons. 1. J. Am. Chem. Soc. 1989, 111, 8551-8566. (38) Halgren, T. A. Merck Molecular Force Field. I. Basis, Form, Scope, Parameterization, and Performance of MMFF94. J. Comput. Chem. 1995, 17,490-519. (39) Jorgensen, W. L.; Maxwell D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 117, 11225-11236. (40) Weigend, F.; Ahlrichs, R. Balanced Basis Bets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys.

Chem. Chem. Phys. 2005, 7, 3297-3305.


23


Page 24 of 26

(41) Balcells, D.; Maseras, F.; Ujaque, G. Computational Rationalization of the Dependence of the Enantioselectivity on the Nature of the Catalyst in the Vanadium-Catalyzed Oxidation of Sulfides by Hydrogen Peroxide, J. Am. Chem. Soc. 2005, 127, 3624-3634. (42) Ujaque, G.; Maseras, F.; Lledós, A. Theoretical Study on the Origin of Enantioselectivity in the bis(dihydroquinidine)-3,6-pyridazine Osmium tetroxide-Catalyzed Dihydroxylation of Styrene. J. Am. Chem. Soc. 1999, 121, 1317-1323 (43) Carbó, J. J.; Maseras, F.; Bo, C.; van Leeuwen, P. W. N. M. Unraveling the Origin of Regioselectivity in Rhodium Diphosphine-Catalyzed Hydroformylation. A DFT QM/MM Study, J. Am. Chem. Soc. 2001, 123, 7630-7637. (44) Fernández-Pérez, H.; Donald, S. M. A.; Munslow, I. J.; Benet-Buchhoz, J.; Maseras, F.; Vidal-Ferran, A. Highly Modular P-OP Ligands for Asymmetric Hydrogenation: Synthesis, Catalytic Activity and Mechanism. Chem. Eur. J. 2010, 16, 6495-6508. (45) Liu, C. H.; Besora, M.; Maseras, F. Computational Characterization of the Origin of Selectivity in Cycloaddition Reactions Catalyzed by Phosphoric Acid Derivatives. Chem.

Asian J. 2016, 11, 411-416. (46) Cogan, A. D.; Liu, G.; Kim, K.; Backes, B. J.; Ellman, J. A. Catalytic Asymmetric Oxidation of tert-butyl Disulfide. Synthesis of tert-butanesulfinamides, tert-butyl Sulfoxides, and tert-butanesulfinimines, J. Am. Chem. Soc. 1998, 120, 8011-8019. (47) Bolm, C.; Bienewald, F. Asymmetric Sulfide Oxidation with Vanadium Catalysts and H O Angew. Chem., Int. Ed. Engl. 1995, 34, 2640-2642. 2

2

(48) Hay, P. J.; Wadt, W. R. Ab-initio Effective Core Potentials for Molecular Calculations. Potentials for the Transition Metal Atoms Sc to Hg, J. Chem. Phys. 1985, 82, 270-283. (49) Wadt, W. R., Hay, P. J. Ab-initio Effective Core Potentials for Molecular Calculations. Potentials for Main Group Elements Na to Bi. J. Chem. Phys. 1985, 82, 284-298. (50) Hay P. J.; Wadt, W. R. Ab-initio Effective Core Potentials for Molecular Calculations. Potentials for K to Au Including the Outermost Core Orbitals, J. Chem. Phys. 1985, 82, 299310.


24

Page 25 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(51) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257-2261. (52) Hariharan P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta, 1973, 28, 213-222. (53) Gordon, M. S. The Isomers of Silacyclopropane. Chem. Phys. 1980, 76, 163-168. (54) Rappé, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A. ; Skiff, W. M. UFF, A full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J.

Am. Chem. Soc. 1992, 114, 10024-10035. (55) Ponder, J. W.; Case, D. A. Force Fields for Protein Simulations Adv. Protein

Chem. 2003, 66, 27-85. (56) Ren, P. Y.; Ponder, J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation, J. Phys. Chem. B 2003, 107, 5933-5947. (57) Ren, P. Y.; Ponder, J. W. Consistent Treatment of Inter- and Intramolecular Polarization in Molecular Mechanics Calculations, J. Comput. Chem. 2002, 23, 1497-1506. (58) Chai, J.-D.; Head-Gordon, M. Long-range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections, Phys. Chem. Chem. Phys. 2008, 10, 66156620. (59) Fraser, H. J.; Collings, M. P.; McCoustra, M. R. S.; Williams, D. A. Thermal Desorption of Water Ice in the Interstellar Medium, MNRS, 2001, 327, 1165-1172


25


Page 26 of 26

For Table of Contents use only:

Expanding the range of force fields available for ONIOM calculations: the SICTWO interface W. M. C. Sameera and Feliu Maseras

26 ACS Paragon Plus Environment

Expanding the range of force fields available for ONIOM calculations

Recommend Documents