Polarizable Continuum Model (PCM) in G09 (Part II)

One of the most successful posts this blog has ever published was on certain nuances of the solvation calculations on PCM in G03. However there are some differences in the SCRF modules between G09 and G03 and I here present some of them as well as some tips to work with the new version.

The SCFVAC keyword used to calculate the Gibbs Solvation Energy change is no longer available. It is now replaced by DoVacuum which should be included in the SCRF options as SCRF=(DoVacuum,etc.). However, the absolute solvation energy now requires a gas-phase optimization along with a frequency calculation followed by the same calculations with the SCRF=SMD option in the desired solvent and with the appropriate variables.

Gaussian 03 used to calculate and report non-electrostatic contributions to the solvation energy, however they were not included in the total energy nor during optimization procedures. These non-electrosatic interactions are no longer calculated in the default. In order to include these terms during the SCF procedure, and to have them reported separately, the SCRF=SMD option should be used.

My previous post on PCM mentioned the usage of the options OFac=0.8 & RMin=0.5 as part of the additional input. These ‘magic numbers’ (I hate the term) were used to modify the way by which the overlapping spheres were treated in order to create the surface which in turn defined the cavity. G09 uses a new algorithm to make the overlaps generate a smoother surface. I recommend to use the default values before including ‘magic parameters’.

All the default values which G03 used can be retrieved with the G03Defaults keyword, but it is strongly suggested to use it only for comparison with calculations previously done with the older version.

As with some other so-called ‘white papers’ this post will be further improved as more information arises during my own calculations. Thanks for reading! Please comment/like/share this post, as well as others in the blog, if you found useful the information within.


About joaquinbarroso

Computational and theoretical chemist in his early forties, in love with life, science, baseball, and literature. Science literacy makes us responsible citizens, it is therefore a scientific duty to talk, write, and engage with the general public; as Feynman said, if you find science boring, your learning from the wrong teacher. "Make like a molecule and react!"

Posted on October 24, 2011, in Computational Chemistry, Gaussian, Models, Software, Theoretical Chemistry, White papers and tagged , , , , , , . Bookmark the permalink. 25 Comments.

  1. Dear Sir, may you share the detailed procedure on calculating Gibbs Solvation Energy change ? I want to know how to calculate it? Any reference ?

  2. Dear Prof Barroso, let me say that your blog is rich in info.
    I need an help concerning the solvation energy with G09 in a ONIOM (QM/MM) system:
    how is written in the guide I used the precise key word ( in may case SCRF=(ONIOMPCM=B,DoVacuum) because the system is quite large and the B option is not quite expensive in term of memory).
    The obtained output is like: ONIOM: extrapolated energy = ……

    There is no trace of the Solvatation value.
    On the guide thare are no informations concerning any key to use to obtain the Solvation Energy.

    Thanks a lot for the Help

    Best Regards

  3. Dear Prof Barroso,
    I would like to know how to calculate the delta G value of (water) PCM at WB97XD level using g09.
    I tried with SCRF=(PCM,Read,DoVacuum,SDM) keyword and lastly i used radii=uahf. Job is finished but i could not get direct value of delta G from it.
    This calculation i used for transition metal complex.(Ru metal) This is i require for the solvent corrected free energy profile for studying the mechanism.
    Please do the needful,

    Thanking you,
    Sandhya K S

  4. Dear Prof Barroso,
    As many of them, I’m trying to calculate deta G value of solvation, in water or acetonitrile.
    I tried with scrf=(pcm,solvent=water,dovacuum) and with smd as well, but I don’t understand how do you calculte the delta G value, I just don’t understand how the log file works.
    I need to calculate this value in order to evaluate oxidation potential of organic molecules for my master thesis.
    I would be very thankful if you could help me, it’s driving me crazy for the moment.

    Best Regards,
    Silviu Preda

  5. %chk=ethanol.chk
    #p b3lyp/6-31G(d) test nosymm scrf(pcm,self,skipvac,smd) scf(tight)

    solvation exp -5.01kacal/mol

    0 1
    C -2.99174682 2.84044012 0.00000000
    H -2.63507398 3.34483831 0.87365150
    H -2.63507398 3.34483831 -0.87365150
    H -4.06174682 2.84045331 0.00000000
    C -2.47843110 1.38850797 0.00000000
    H -2.83510335 0.88410999 0.87365187
    H -2.83510454 0.88410957 -0.87365114
    O -1.04843110 1.38849035 -0.00000097
    H -0.72799618 0.48354757 0.00017528

  6. %chk=ethanol.chk
    #p b3lyp/6-31G(d) test nosymm scrf(pcm,self,dovac,smd) scf(tight)

    solvation exp -5.01kacal/mol

    0 1
    C -2.99174682 2.84044012 0.00000000

  7. Luis F. de Figueiredo


    I tried to use the G03Defaults keyword in order to reproduce previous results but I often get an error during a tighter optimization step:
    # B3LYP/6-311++G** Opt=tight int=ultrafinegrid Test geom=allcheck guess=read SCRF=(PCM,solvent=Water,G03D

    and when I get it working, I then try to calculate the vibrational frequencies with
    # B3LYP/6-311++G** Freq Test Geom=allcheck Guess=Tcheck SCRF=(PCM,solvent=Water,G03Defaults)

    For this step I always get the following error message:
    These settings can be overridden by following keywords.
    Analytic 2nd derivatives not available using g03defaults.
    Error termination via Lnk1e in /path_to_gaussian09/g09/l301.exe

    Have you experienced similar issues using the G03Defaults keyword?


  8. Hi, Dr. Barroso

    On G03 print this:

    Variational PCM results
    (a.u.) = -574.656387
    (a.u.) = -574.663975
    Total free energy in solution:
    with all non electrostatic terms (a.u.) = -574.650706
    (Polarized solute)-Solvent (kcal/mol) = -4.76
    Cavitation energy (kcal/mol) = 22.64
    Dispersion energy (kcal/mol) = -15.08
    Repulsion energy (kcal/mol) = 0.77
    Total non electrostatic (kcal/mol) = 8.33

    But in G09 not are print neither used IOPs(3/70).
    Do you know other keywork to g09 displey Variational PCM results?, so g09 skip PCM results


    P.D. puedes contestarme en español

  9. Dear Sir,
    Greetings.I am working with Gaussian 09.I have tried to optimize the rosebengal dye in the solvent acetonitrile with scipcm model.when i working with it , i often get the error result. i have enclosed the input and the output result. please do the needful.
    Thanking you,

    #p opt freq=(raman) ub3lyp scrf=(scipcm,solvent=acetonitrile) pop=nbo polar density=current


    Total “Solvent Accessible Surface Area” of solute cavity = 1.848246E+03 au
    Volume of solute cavity = 4.137515E+03 au
    Total number of cavity surface points = 5881
    Percentage of total area from each origin:
    Org: 1: 1.6% Org: 2: 0.9% Org: 3: 0.5% Org: 4: 1.0%
    Org: 5: 1.5% Org: 6: 1.2% Org: 7: 0.0% Org: 8: 1.1%
    Org: 9: 0.5% Org: 10: 0.9% Org: 11: 2.4% Org: 12: 1.6%
    Org: 13: 1.2% Org: 14: 1.6% Org: 15: 2.8% Org: 16: 7.6%
    Org: 17: 7.0% Org: 18: 5.9% Org: 19: 7.5% Org: 20: 4.1%
    Org: 21: 2.6% Org: 22: 0.1% Org: 23: 1.0% Org: 24: 0.7%
    Org: 25: 1.4% Org: 26: 1.4% Org: 27: 1.4% Org: 28: 5.2%
    Org: 29: 5.9% Org: 30: 6.0% Org: 31: 4.7% Org: 32: 0.9%
    Org: 33: 3.0% Org: 34: 2.1% Org: 35: 4.4% Org: 36: 3.3%
    Org: 37: 4.8%

    WARNING! Serious error in surface integrals.
    Nuclear flux = 455.73 Qnuc = 464.00 Error in int = 0.00
    It is probable that some of the solute is outside the cavity and/or
    parts of the cavity surface cannot be reached from the origin.
    Try more integration points or a different set of integration origins.

    Surface Problems in SciFoc
    Error termination via Lnk1e in C:\G09W\l502.exe at Wed Apr 30 14:10:04 2014.
    Job cpu time: 0 days 0 hours 1 minutes 33.0 seconds.

    • Hi Punitha,

      Try using scrf=(scipcm,solvent=acetonitrile,read) and use the OFac variable described in this very same post and the other one on PCM that I have written before.

      Have a nice day!

      • Hi, Prof. joaquinbarroso,

        Experiencing the same issue.

        Tried using scrf=(scipcm,solvent=acetonitrile,read), modifying OFac and Rmin values, does not help. Receiving the same error.

        Is there anything else that could aid?

        Thank you beforehand,
        Kind regards,
        Dmitry Novikov

      • Dear Dmitry, try using a larger definition of radii, like radii=UAKS or UA0. This post is getting a little outdated so it will depend if you are using G03, G09 or G16 because the SCRF calculation details have changed from one version to the next. I need to make a revision, I guess.
        I hope this helps!

  10. Dear Dr. Barroso,
    What role, if any, does the RSOLV parameter play in cavity definition in G09? I am computing a relaxed potential energy scan which fragments a molecule, and I am following it to large distances, and at a particular distance (3.4 Å in acetonitrile), the cavities become disjoint. Is this distance determined or even affected by the RSOLV parameter? If not, how is it determined? (I know it is the Separa routine in the PCM link, but I don’t know how it works.) I am including nonelectrostatic effects, but I don’t think their inclusion affects how Gaussian defines the cavity…

    Any assistance or comment you could provide would be helpful.

    Thank you very much,
    Ted Lorance

    • Dear Ted,

      The RSOLV parameter, as I’m sure you know, modifies the solvent radius from that which is stored by default for each solvent. The only way a change in RSOLV could affect the cavity surface is if you also use the keyword SURFACE=SAS, which stands for Solvent Accessible Area. As with other models, a sphere is placed on each atom (or group of atoms if using a united topological group of atoms, UAKS or UAK0) overlapping along your solute to define the cavity; and then, with SAS, a solvent molecule sphere is “rolled” over the previous surface; the surface defined by the “rolling” of the solvent sphere is used as the cavity. A bigger solvent sphere will yield a larger cavity albeit with less resemblance to the original solute. If you use a large solv.sphere this might prevent the disjointing of your two resulting cavities (until the limit where the separation between them is slightly smaller than 2xRSOLV*), but also this could lead to unrealistic results by defining an extremely large radius for a solvent molecule.
      If this fragments are getting separated during your scan, it’s very likely that the molecular dipole changes a lot and thus I’d suggest (as maybe a first approximation) to perform a series of Onsager (keyword DIPOLE) calculations (the cavity is the best sphere possible that includes the solute), and I mean a series of isolated calculations because I believe optimizations aren’t possible with this model.
      Inclusion of non-electrostatic terms doesn’t affect the definition of the cavity in PCM models.

      Thank you for a very thought provoking question!

      I hope this helps

      *accounting for overlap of a sphere that could fit between both fragments

      • Dear Dr. Barroso,
        Your previous comments were very helpful. However, the transition state in the disjoint cavity region has turned out to be rate-limiting in our system, so we are now faced with the necessity of making some comment on the accuracy of the our application of the solvation model.
        My concern is that I don’t know for sure exactly what Gaussian09 does when it is computing the solvation properties of a system treated by a single quantum calculation with two cavities in the polarizable continuum model. I strongly suspect that the charges induced by the molecule on the surface of one cavity do not interact with the other cavity or its contents, which seems physically inaccurate to me. Are escaped-charge and similar errors minimized in each cavity or only overall? If there is significant electron density outside of either cavity, what effect does that cause? Most importantly, is there a reason to be more suspicious of disjoint-cavity results in PCM than other PCM results? I would hate to find that a significant portion of the apparent free energy of the TS is due to solvation error caused by the disjoint cavities (one might hope or pretend that OTHER solvation errors would cancel since the energies we are working with are relative).
        I cannot find anything in the literature (or even the Gaussian manual) that refers to disjoint cavities; this only increases my unease and suspicion.

        Thank you very much,
        Ted Lorance

  11. Xiankai Jiang

    Dear Dr. Barroso,
    I have studied your post. However, I meet with difficulties. Could you help me?
    The problem is :

    INPUT(gaussion 09)
    #p sp b3lyp/6-31+g(d,p) geom=connectivity iop(2/16=3,5/13=1)
    scf=(maxcyc=150,xqc) scfcon=5 symm=loose scrf=(Externaliteration,PCM,Read,DoVacuum,SMD)

    Title Card Required

    1 1
    C -1 -2.96822500 -1.86080200 0.06078100
    O -1 -3.97831800 -1.35035900 -0.43498300
    N -1 -1.78967800 -1.21031200 0.23173400

    OUTPUT(something wrong)
    (Enter /home/jxk/g09/l124.exe)
    SC-PCM: Using the SCF density.
    SC-PCM: The equilibrium contribution to the free energy is -0.138704313 a.u.
    SC-PCM: The total energy including the PCM contribution is -729.367155228 a.u.
    SC-PCM: Erms= 0.61885D-06 Emax= 0.30686D-05 Qconv= 0.10000D-06
    SC-PCM: Iterating to achieve self-consistent PCM reaction field, IExtIt= 50.
    SC-PCM: Convergence not achieved within 50 iterations … aborting job.
    Error termination via Lnk1e in /home/jxk/g09/l124.exe at Sun Apr 5 14:48:13 2015.
    Job cpu time: 0 days 3 hours 28 minutes 41.4 seconds.
    File lengths (MBytes): RWF= 97 Int= 0 D2E= 0 Chk= 7 Scr= 1

    what should I do?
    Dr. Barroso, Could you help me?
    Thank you very much,
    Xiankai Jiang

  12. Dear Dr. Barroso,

    I have tried calculate an energy in g09. It calculation finish, but shows this message “Error on total polarization charges = 0.01606”. What is wrong and what can I do to fix it?

    The features that I used:

    #p M062X/QZVP density=current NoSymmetry scrf=(cpcm,solvent=water) scf=(qc,maxcycles=1024).

    Another question. How can I use a dielectric constant=40? Use a solvent with a value closer is the only way?

    Thanks for your time and help.

  13. Dear Doctor Barroso! Perhaps you don’t remember me, I asked you a question regarding the study of transition states and you helped me quite a lot! I thank you again for your help. My tutor has offered me a Magister tesis and I’m quite looking forward it.

    I have a question regarding SMD and DoVacuum. I’ve run the calculation using scrf=(solvent=benzene,smd,dovacuum) but I’m not sure how I can find the value of DeltaG in the output file. Must I put something else in the command?

    Thank you so much for your help!

    Saludos desde este rinconcito de Caracas,

    Maria Fernanda.

  1. Pingback: Delta G of solvation in Gaussian09 « Dr. Joaquin Barroso's Blog

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: