I am running a B3LYP-D3/6-31+G geometry optimization of a protein-ligand interaction. Well, not quite, I have simplified the protein to four amino acids surrounding the active site, but I keep the entire ligand (170 atoms total).

I am wondering if I should optmize the geometry in solution, i.e. will including solvation effects make my results (energies) be more realistic? My mind is split.

On the one hand, I am modelling a biological system which usually takes place in aqueous solution. On the other hand, my system is so simplified, that it no longer has much root in reality. By "placing my molecules in water", the water molecules will interact with parts of my system that otherwise would not interact as strongly (there would be much less water inside the protein than outside it, I mean). So, I feel there is no guarantee that my results will improve, and I am not able to find any data indicating one or the other. If I included the entire protein, then the case would be different, perhaps, since at least then that protein would be "solvated correctly".

So, given the simplified system, will solvation improve the resulting geometry, and thus the interaction energy between the protein residues and the ligand, or not?

  • 2
    $\begingroup$ Solvation corrections at this level of theory will not do you much good. Unfortunately B3LYP/6-31+G is probably far too simple for that system in any case. Calibration is necessary. I am obviously not familiar with your project, so any advice is flawed to begin with. Have a look at this question and you might see how difficult it can be to even selct the appropriate functional, not to speak of basis set or any other shortcomings. $\endgroup$ Dec 16, 2015 at 15:20
  • $\begingroup$ I should add that I use the B3LYP-D3 method. Still, you think solvation will not benefit much? Your answer in the linked question was insightful, thank you for the link. $\endgroup$
    – Yoda
    Dec 16, 2015 at 15:37
  • $\begingroup$ In case you are interested in binding to protein, the method doesn't mean much, a fair amount of research is done using molecular mechanics of all things. What is more important is to approximate solvation using actual solvent molecules and some stochaistic or systematic method to find the global minimum. $\endgroup$
    – permeakra
    Dec 16, 2015 at 15:54
  • $\begingroup$ I don't think there is a definitive answer to this. The best thing is to do it both ways to see if it makes a difference. See for example this paper dx.doi.org/10.1002/jcc.21458 $\endgroup$
    – Jan Jensen
    Dec 17, 2015 at 10:23

1 Answer 1


This is difficult to answer without full knowledge of the system, but some tips:

1) The overall energy of the system is arguably irrelevant. Rather, the relative difference between different conformations, and perhaps, different ligands is of utility.
2) The level of theory and basis set your using is quite low. If you have the time and the resources, scale it up. For biomolecules I've found the minnesota functionals to be most accurate (if insisting on DFT) but I still believe MP2 is king.
3) Referring back to #1, if you solvate, then your comparison should too be solvated. Perhaps your approach would be to run a solvated and unsolvated model for both your ligand, and perhaps, another known ligand. This could too help convince you of the utility of the solvation model.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.