PyConSolv: A Python Package for Conformer Generation of (Metal-Containing) Systems in Explicit Solvent

We introduce PyConSolv, a freely available Python package that automates the generation of conformers of metal- and nonmetal-containing complexes in explicit solvent, through classical molecular dynamics simulations. Using a streamlined workflow and interfacing with widely used computational chemistry software, PyConSolv is an all-in-one tool for the generation of conformers in any solvent. Input requirements are minimal; only the geometry of the structure and the desired solvent in xyz (XMOL) format are needed. The package can also account for charged systems, by including arbitrary counterions in the simulation. A bonded model parametrization is performed automatically, utilizing AmberTools, ORCA, and Multiwfn software packages. PyConSolv provides a selection of preparametrized solvents and counterions for use in classical molecular dynamics simulations. We show the applicability of our package on a number of (transition-metal-containing) systems. The software is provided open source and free of charge.


Multiwfn charge settings
Multiwfn 3.8 1 is used to compute the charges for the parametrizations, as follows.
All atoms are used for fitting 4 layers, with the scaling factors: 1.4,1.6,1.8,2.0.Automatic radii are utilized for fitting, missing radii are taken from UFF and scaled by 1/1.2.For fitting, the tightness parameter is 0.1, the restraint strengths for the two stages are set to 0.0005 and 0.0010, respectively.
The maximum number of iterations is 50 and the convergence threshold is set to 10 -6 .Charge equivalence constraints are enabled for CH 2 and CH 3 groups.No charge constraint settings are used and the connectivity is determined automatically.

CREST settings
CREST 2 was used for generating conformations as comparisons for the chosen test cases.For this, the GFN2-xTB 3 method was chosen, using ALPB 4 as the solvent model, with an appropriate solvent for each system.The energy threshold, RMSD and conformer pair energy differences were set to 6 kcal/mol, 0.125Å and 0.05 kcal/mol, as detailed in the CREST documentation.For the case of the Mo-catalyst, the energy threshold was raised to 50 kcal/mol, because a low number of conformers was obtained.

Solvent Parametrization
The solvents were parametrized using RESP charges and antechamber, after which they were equilibrated and a short, 10ns, classical molecular dynamics (cMD) run was performed to examine the density.

Timings
The timing for the complete run of PyConSolv for each the testcase system is detailed in Table S2.The calculations were performed on a Ryzen R9 5900x CPU, using 12 threads and a total system memory of 64 GB.The simulations were performed using the GPU accelerated implementation of pmemd 6 on an RTX 3060 GPU.

Amber MD settings
To perform the classical molecular dynamics simulations, the Amber20 6 package was used, in combination with AmberTools23 6 .For the system equilibration, a multi-step approach with repeated heating and cooling was chosen, as follows: 1.The entire simulation box is restrained using cartesian restraints (restraint weight = 1000 kcal mol -1 Å -2 ), with the exception of the hydrogen atoms, which undergo an energy minimization for 500 cycles.
3. The system is equilibrated under constant volume to achieve a temperature of 300K, using Langevin dynamics 7 with the collision frequency of 2.0 ps -1 , a timestep of 0.001 ps and the SHAKE 8 algorithm for hydrogen bonds.
4. The solute heavy atoms remain restrained (restraint weight = 1000 kcal mol -1 Å -2 ) and the whole system is equilibrated at constant pressure (1 atm), to achieve a reasonable system density and eliminate the voids created by the automated solvent placement in tleap.This is achieved using Langevin dynamics with the collision frequency of 2.0 ps -1 , a timestep of 0.001 ps and the SHAKE algorithm for hydrogen bonds.

5.
The system is subsequently cooled under constant volume to achieve a temperature of 100K, using Langevin dynamics with the collision frequency of 2.0 ps -1 , a timestep of 0.001 ps and the SHAKE algorithm for hydrogen bonds.
8. A final constant pressure (1 atm) equilibration takes place, using Langevin dynamics with the collision frequency of 2.0 ps -1 , a timestep of 0.001 ps and the SHAKE algorithm for hydrogen bonds.
The simulations then take place under constant pressure (1 atm, Berendsen barostat 9 ) with isotropic position scaling, a pressure relaxation time of 2.0 ps with Langevin dynamics, SHAKE for the hydrogen atoms and a timestep of 0.001 ps.
The size of the simulation box is automatically detected based on the size of the system and a cubic shape is chosen in all instances.These settings resulted in the following number of particles for the simulations: Cu(I)-Calix [8]arene in chloroform: 966; molybdenum N-heterocyclic carbene (NHC) catalyst in dichloromethane: 844; hydrogenobyric acid in water: 5529; methylcobalamine in water: 6766.In the case of the Molybdenum catalyst, it is apparent that CREST 2 achieves a greater sampling for the conformational space, thanks to its metadynamics based approach.PyConSolv can achieve further sampling through prolonging the classical molecular dynamics simulations or employing enhanced sampling techniques.Because of to the use of explicit solvation, unphysical intramolecular hydrogen bonding is minimized, when conformers are generated using PyConSolv, as shown in Figure S3.

Package Usage
The package was installed following the instructions on the github page(https://github.com/PodewitzLab/PyConSolv).
A new folder for the chosen system was created and the structure was saved as an XMOL file (input.xyz).This serves as the input for the pyconsolv command.
The parametrization was started using the following command: "pyconsolv input.xyz-m MMM -b BBB -s SSS -d DDD", where MMM represents the DFT functional, BBB represents the basis set, SSS represents the solvent and DDD represents the dispersion correction model.
After the geometry optimization and frequency calculation are complete, a pop-up window appears, where the charges for each of the detected fragments are assigned (see Figure S4).
Following the charge assignment, the parameters are created, the system solvated and the equilibration is automatically started.
The simulation run is started in the newly created simulation folder, using the "run_simulation.sh"script, generated by PyConSolv.
Following the 100ns cMD production run, the analysis was performed using "pyconsolv -a sim-01", choosing kmeans as the clustering method and default values (10 mininum points per cluster, 10 clusters total).
For a more detailed guide, please consult the user manual provided with PyConSolv, which can be found on the github repository.

Figure S1
Figure S1 Structures of parametrized counterions; color coding: C atoms in grey, B atoms in cyan,

Case 2 :
Figure S2 Comparison of conformers generated with CREST and PyConSolv.

Case 4 :Figure S3 .
Figure S3.Conformers generated with PyConSolv, compared to a conformer with significant

Figure S4
Figure S4Fragment charge assignment dialogue.Here the user must provide the charges for each

Table S1 .
Density of boxes with parametrized solvent, after 10ns cMD.Densities taken from CRC Handbook of Chemistry and Physics, 95 th Edition. 5 * denotes gas phase, while the calculated density is in the condensed phase.

Table S2 .
Timings for electronic structure calculations and simulation steps for PyConSolv and CREST.