10
$\begingroup$

I have been using RDKit to generate 3D coordinates, cleaning-up the structures using a general-purpose force field as follows:

 (1) Adding hydrogens: AllChem.AddHs
 (2) Compute 3D coordinates: AllChem.EmbedMolecule
 (3) Clean-up with force field: AllChem.MMFFOptimizeMolecule

  • But I would like to take into account different solvents so that the 3D structure more closely resembles in vitro activity. For example, an in vitro screen where small molecules are diluted in (1) DMSO or (2) ethanol or (3) water.

  • I know that Psi4 can estimate properties in different solvents. But I have been unable to find a python library that generates an actual 3D structure that's optimized for solvents (e.g. a solvent-specific force field that would be applied after the 3D coordinates have been computed as in the RDKit example above).

  • Any help/suggestions on how to do this, exporting a .mol/.sdf file (preferably using python) would be much appreciated.

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ Good first question, I fear you might not find anything satisfactory beyond running a real DFT calculation with solvent. An expert in forcefields can hopefully prove this incorrect though. $\endgroup$
    – Tristan Maxson
    Mar 20, 2021 at 20:19
  • 1
    $\begingroup$ If the solvent is considered implicit, it is a matter to define it's dielectric permisivity. $\endgroup$
    – Camps ♦
    Mar 21, 2021 at 2:59
  • $\begingroup$ NAMD does have an implicit (GB) solvent option, so optimizing the structure with that might work. However, I am not sure if any other solvent except water are implemented, and even if they are, I am not sure the results would be accurate enough for your purpose. If your molecule is small, you are better off using DFT with PCM or SVPE type implicit solvation (more accurate), as Tristan already mentioned. $\endgroup$
    – S R Maiti
    Mar 21, 2021 at 15:58
  • $\begingroup$ These are all good suggestions. But my concern is that they may not scale, if, for example, it was a high throughput screen or the downstream goal was predictive modeling (QSAR). Here, each molecule that's predicted by the trained model would then have to be passed through this pipeline, which could be millions of molecules. So, accuracy may need to be relaxed. It's more a matter of if I can somehow coarsely account for solvent in computing 3D structure. $\endgroup$
    – jjkow
    Mar 21, 2021 at 16:38
  • 1
    $\begingroup$ Depending on how you want to use this data, you may want to assess if you are able to see any difference, though. The error coming from forcefield parametrization, finite temperature effects, equilibrium between conformers etc can all be much bigger than the effect of solvent on geometry. $\endgroup$
    – Greg
    Mar 24, 2021 at 0:23

1 Answer 1

Reset to default
7
$\begingroup$

From reading your comments, and the question, you want to generate optimized structures of small bioactive molecules in solvents, but it has to be fast. I would definitely recommend using semi-empirical methods, because that would give you the best of both worlds. The calculation will be faster than DFT, and you can use those methods with relatively accurate implicit solvation models.

Here I have covered three softwares, with several semi-empirical methods, and one force-field based method:

MOPAC

MOPAC has loads of semi-empirical methods, prominent among them being PM3 (old), PM6 (newer), PM6-DH+ (empirical dispersion and hydrogen bonding) etc. MOPAC uses the COSMO model for implicit solvation. In the input file, you have to put the keyword EPS=78.6 to specify the dielectric constant of the solvent (default 78.6 for water). The MOPAC manual recommends that you use EF algorithm for optimizing with solvent, instead of the default L-BFGS.

As a test, I optimized an esomeprazole molecule (sdf file obtained from PubChem) on my laptop (Intel core i3, 2.2 GHz 4 core) with PM6/water. It took 8.53 seconds for the optimization to finish (multi-threaded run).

GAMESS

GAMESS can use PM3, AM1 and RM1 semi-empirical methods officially. However, PM6 and PM6-D3H+ (D3 dispersion and hydrogen bond corrected) methods are also implemented for elements until fluorine (unofficially). All the semi-empirical methods are interfaced with SMD solvation model (which is more accurate than COSMO as it also includes non-electrostatic contributions to solvation). See this paper for details of the implementation. For example, with PM3/SMD the RMS error for solvation energy is 2.8 kcal/mol (neutral solute), so the method is quite accurate. To use this you have to specify $BASIS GBASIS=PM3 $END and $PCM SOLVNT=WATER SMD=.T. $END. For the other solvents, use SOLVNT=DMSO or SOLVNT=ETHANOL instead.

I redid the same optimization of esomeprazole on GAMESS and it took 3.8 minutes on 4 cores. This is probably because the implementation of semi-empirical methods in GAMESS is slightly slower itself, also because GAMESS uses different criteria for convergence of optimization and SMD is slightly more demanding.

XTB

XTB is uses a semi-empirical tight-binding approach. There are three parameterizations—GFN0-xTB, GFN1-xTB, GFN2-xTB, each newer than the last. XTB is open-source, so you will have to compile it. It has analytical linearized Poisson-Boltzmann(alpb) implicit solvation. According to the official manual, DMSO and water are available, ethanol is not mentioned. However, from my test, the latest version seems to be able to run with ethanol solvation as well.

XTB is run from command line:

xtb.exe --alpb water --opt -P 4 --gfn 2 esomeprazole.sdf > out.log

This uses GFN2 parameters with water solvation and 4 cores. The calculation took 7.346 second on my laptop.

For the other solvents, use --alpb dmso and --alpb ethanol.

The XTB software also contains a polarizable force field GFN-FF, which can also be run with implicit solvation. I am not sure how accurate it is, but the optimization of esomeprazole takes 0.445 seconds, a significant gain in speed. To use the force field, you have to use --gfnff instead of --gfn 2 in the command line.

Running with Python

I have limited knowledge of Python, so I can't help much here. For MOPAC and GAMESS, you can use the python package ASE. It is supposed to be able to interface to the softwares directly. ASE can also write input files, so use that if the interface does not work. As XTB is used from the command-line and can read sdf files, you can directly use os.subprocess.run(). XTB automatically outputs the optimized geometry in the same format as input (.sdf here).

Not all semi-empirical methods are well-parallelized. MOPAC and XTB uses thread-based parallelization, whereas GAMESS uses an MPI-based one. I assume you are using a cluster for your calculations. It would probably be the most efficient if you ran each calculation on one core, and parallelized the handling of different molecules. You could use something like mpi4py for this. Also check if the software you are using allows multiple instances to be launched at the same time. On my laptop, XTB can do this, but I have not tested GAMESS or MOPAC.

Improve this answer
$\endgroup$
1

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .