Journal of Chemical Theory and Computation Article Figure 4. Convergence of free energy profiles, positions of energy minima and energy barriers. The simulation windows within 1.0−2.5 nm have been simulated for 50 and 100 ns in case of UP and UF simulations of CP and CF simulations, respectively; the rest of each profile is calculated from 30 ns of simulation. Free energy profiles calculated from short simulation times (<5 ns) are biased by high error, because of small data set and nonequilibrium starting structures. Free energy profiles obtained by UP and CP simulations (left) show slow elimination of an artificial minimum (∼2.0 nm) and deepening of the global minimum (∼1.5 nm). Free energy profiles obtained from UF and CF simulation are consistent in coumarin positioning and have deeper energy minima than those obtained from UP and CP simulations. The global minimum energy is considered as reference. UP and UF refer to umbrella simulations with initial structures obtained by pulling and free unbiased simulation, respectively; CP and CF refer to constraint simulation with pulling and free initial structures, respectively. gradually narrowed from 1.38 nm after 5 ns of simulation to 0.51 nm after 40 ns and thereafter remained constant. ΔGwat gradually rose throughout the simulation, to 5.9 ± 0.2 kcal/mol, while ΔGpen grew during the first 11 ns of simulation, until 20 ns of simulation it gradually declined to 2.8 ± 0.1 kcal/mol and then fluctuated around this value. In the last 50 ns of the simulation prolonged to 100 ns the position of the energy minimum remained constant; ΔGpen fluctuated around 2.8 ± 0.1 kcal/mol and ΔGwat continued to rise, to 6.2 ± 0.2 kcal/ mol. The position of the minimum in the UF free energy profile was almost constant (within 1.29−1.35 nm) during the whole simulation time, with the area accessible by thermal motion slowly widening from 0.44 to 0.66 nm. ΔGwat slowly decreased during the first 19 ns of simulation, then very slowly increased, and after 30 ns ΔGwat converged to 6.7 ± 0.1 kcal/mol, while ΔGpen became convergent after 20 ns, fluctuating within 2.6 ± 0.1 kcal/mol. The CF free energy profile showed a minimum position within 1.29−1.49 nm, and the area of thermal motion slowly widened from 0.41 to 0.70 nm. ΔGwat fluctuated around 6.4− 7.0 kcal/mol during the whole simulation time, while ΔGpen decreased during the first 10 ns of simulation and then fluctuated around 2.9−3.3 kcal/mol. The prolonged simulation to 100 ns showed similar trends − a free energy minimum at 1.29 nm, thermal motion within 0.6 nm, a constant ΔGwat value of 7.0 kcal/mol after 80 ns, and ΔGpen already convergent with a final value of 3.1 ± 0.1 kcal/mol. Effect of Coumarin Partial Charges. As assignment of partial charges might introduce another systematic sampling error into free energy calculations, we carried out 10 ns long CF simulations with PRODRG and Mulliken charges (assigned partial charges are listed in Supporting Information Table S1), to assess the extent to which the partial charges affected the free energy profiles. Coumarin with partial charges assigned by PRODRG bore a dipole moment of 9.5 D, assignment of Mulliken partial charges led to 6.0 D, and RESP partial charges resulted in a dipole moment of 4.9 D (Figure 5). The dipole moment based on RESP charges was close to that of coumarin in the gas phase calculated by the hybrid DFT method (B3LYP/cc-pvDZ) of 4.6 D, Mulliken partial charges represent a compromise between the dipole moment in water (represented by continuum dielectrics with εr = 78.39) and heptane (εr = 1.92), which we calculated by the CPCM/ B3LYP/cc-pVDZ method and that resulted in 6.7 and 5.4 D, respectively. Considering these values, the dipole moment stemming from PRODRG charges seemed to be unreliably overestimated, which could systematically bias free energy profiles based on PRODRG charges. The global minimum of the ΔG profile of coumarin bearing RESP partial charges was located at 1.29 nm (CF with a 10 ns sampling window), energy minimum of coumarin bearing Mulliken partial charges was localized at 1.20 nm, and the minimum for PRODRG-charged coumarin was shifted toward the bilayer/water interface, at 1.62 nm. The global free energy minimum for RESP-charged coumarin was also considerably deeper than for Mulliken- charged or PRODRG-charged coumarin (ΔGwat: 7.5, 5.6, and 3.3 kcal/mol, respectively), and the bilayer center penetration barriers of the systems also differed (ΔGpen: 3.1, 4.6, and 10.1 kcal/mol, respectively) (Figure 5). As expected, the energy cost of bilayer center penetration grows with the increasing 1206 dx.doi.org/10.1021/ct2009208 | J. Chem. Theory Comput. 2012, 8, 1200−1211