15
$\begingroup$

Before performing Molecular Dynamics on a co-polymer that doesn't have a ready-to-use structure anywhere, the monomer fragment was built in Avogadro software with 0 net charge, B3LYP theory and 6-31G(d) basis set with a multiplicity of 1, being the parameters used for an optimization in Gaussian.

What are the parameters that I should look out for, in the Gaussian output, to know that the structure has been optimised for use in MD?

$\endgroup$
4
  • 2
    $\begingroup$ I think that there is no need to pre-optimize the polymer with Gaussian in order to do a Molecular Dynamics. It will be enough that the structure has no structural problems (clashing atoms, bad dihedrals, etc.). Molecular Mechanics (using free software like OpenBabel) or even semi-empirical method will do the task faster. This is because in the Molecular Dynamics, your system will change its conformation. $\endgroup$
    – Camps
    May 12, 2020 at 16:46
  • $\begingroup$ @I.Camps Thank you. Is this applicable for the structure of monomer as well? What about the charges on atoms that need to be optimised? Because I'm assuming that if monomers are optimised using Gaussian, then only structural problems can be checked while building a polymer. If geometric stability is all that is required, why do we have PDB databases with optimised structures? $\endgroup$ May 13, 2020 at 9:15
  • $\begingroup$ Normally, the forcefields will assign their own charges according the type of atoms you have. If you do a previous calculation and, even when you save the optimized structure in a format that include the charges (not every format do that), you have to manually set your MD to use those charges. About PDB databases with optimized structures, could you tell me an example? $\endgroup$
    – Camps
    May 13, 2020 at 10:23
  • $\begingroup$ @I.Camps Thank you for the clarification. I think the problem comes with molecules that do not have a forcefield and need to tune the existing ones for the use. The site called Automated Topology Builder takes a PDB file and builds topology based on optimization at the B3LYP/6-31G* level of theory + Hessian. This is also a repository for such kinds of optimised molecules. $\endgroup$ May 14, 2020 at 20:22

2 Answers 2

17
$\begingroup$

Ah, yes, the fun of force-field building.

For the answer to a simple geometry optimization, see I. Camps response below.

Skip to the end if you want what is a more general answer to building an entire force-field. Read the whole thing if you want some insight into force-fields, particularly partial charges.

First, doing an electronic structure calculation (in your example you suggest using B3LYP/6-31G(d)) unsurprisingly gives you an electron density. Perhaps interestingly, it does not give you partial charges.

Partial charges are not real

There is no way to actually calculate a partial charge since they are not real observables. However, given an electron density, you can apply a protocol, to calculate partial charges that represent the given electron density according to some principle.

For the record B3LYP/6-31G(d) is typically used with Restrained Electrostatic Potential (RESP) partial charges, however, more frequently, HF/6-31G* is used. Further, this is the default method for the older AMBER forcefields, as well as the General AMBER Force-Field, although other Force-fields use it.

It is my opinion that RESP using HF or B3LYP with 6-31G* (additional note, 6-31G* is the same as 6-31G(d)) is a terrible choice.

If you are using RESP, then your electronic structure calculation only calculates the electron density, and another program must calculate the RESP partial charges. Gaussian will output ESP (not restrained) partial charges, that usually, are identical.

One problem with RESP when using B3LYP/6-31G(d), is that it adds polarization effects implicitly. What I mean is that you just did a calculation for the electron density of the molecule in a vacuum, and calculated partial charges from this density. However, this density by no means should be useful if you put your molecule in a liquid. In a liquid, the neighboring molecules will polarize each other, and a molecules electron density will distort from the density in a vacuum. However, B3LYP/6-31G* is so terribly bad its vacuum calculations actually somewhat polarize the molecule. Unfortunately, it does it almost randomly, ranging from underpolarizing, to overpolarizing by 35% Paper to read. RESP using HF or B3LYP with 6-31G(d) is so bad, it is actually by fluke useful, in an average way, but messes up alot.

RESP is from the 1990's. It used to take time to calculate electron densities back then, so they made a parameterized model AM1-BCC which uses a even worse semiempirically derived electron density to calculate Mulliken partial charges, and then adds corrections to them in an attempt to match RESP using HF/6-31G*. This is faster to calculate, but again, not that great, it often fails and really, you can't expect it to work outside of its training set.

Important point here though, RESP partial charges are derived by starting with an initial guess for the partial charges, and then calculating the electrostatic potential of the molecule using the partial charges. You also calculate it from the QM software. You then change the partial charges until the forcefields electrostatic potential is as close as possible to the QM's "correct" electrostatic potential. Other methods exist however, and Mulliken is one, where for brevity, you "partition" the electron density to each atom. This is a much more intuitive method that doesn't not rely on opptimization. However, Mulliken charges are usually only qualitatively useful.

More advanced methods such as the Hirshfeld familes CM5, DDEC6 and MBIS exist and in my opinion are far superior to RESP. I do not say superior to AM1-BCC because AM1-BCC is nothing but parameters, it is not a "method".

Why don't we use these better methods?

Hopefully we will in the future, unfortunately partial charges are only part of a force-field, and if you change the charges, you should really change the other parts too.

This is most clear with Lennard-Jones parameters, which is typically fit to experimental data for the given partial charge method. It is thus, alot of work to refit all of the LJ parameters (and torsions as well) for a new charge method. But, this is slowly happening.

Building a Force-Field in Practice

To directly answer your question now, Force-fields take a lot of work. If you are making one from scratch you need to

  1. Generate partial charges
  2. Calculate bond parameters
  3. calculate angle parameters
  4. Calculate torsion parameters
  5. Calculate (fit to data) Lennard-Joes parameters

This is a ridiculous amount of work. And, Gaussian does not do this at a push of a button. Each item is progressively more work and its own task. If you do this however, I would strongly recommend ditching the B3LYp/6-31G(d) electronic structure calculation and the presumably RESP or AM1-BCC partial charges that follow. Use better methods. i.e., use $\omega$B97X-D/DEF2TZVP for instance for the electronic structure calculation, and then use CM5/DDEC6 or MBIS for the partial charges. I recommend MBIS.

But, what I recommend even more is to not make your own force-field if this is your first time. There are many ready to go toolkits that will make one for you.

For polymers I would recommend either GAFF or OPLS.

Here is a link to GAFF. It will use AM1-BCC be default, but the other parameters are fit to it, so it will be decent for you. Follow their instructions.

Or, OPLS has a easy to use webserver that will whip up a force-field for you for almost any program you want to use found here. They use empirically scaled CM1 partial charges - a faster(presumably), older version of CM5. It has some options, but basically, you just upload a structure in a pdb or mol2 file, and it will give you back an entire force-field for your molecule.

$\endgroup$
5
  • $\begingroup$ Would it be even better to compute RESP (instead of e.g. CM5) charges on top of ωB97X-D/DEF2TZVP wavefunctions? In my understanding, RESP charges are explicitly designed to reproduce the correct electrostatic potential near the molecule by actually fitting against the electrostatic potential, while CM5 charges are not generated by this way, so RESP charges should be more theoretically sound, as long as the underlying wavefunction is accurate enough. $\endgroup$
    – wzkchem5
    Jul 17, 2022 at 11:25
  • $\begingroup$ @wzkchem5 the short coming of RESP isn't so much the quality of the electron density being passed to it, as the RESP methodologies ability to translate that to partial charges. In theory, it should be quite good. In practice, it is a numerical solution with all the problems of a global minima search, potentially with a very flat minima. $\endgroup$
    – B. Kelly
    Jul 18, 2022 at 11:51
  • $\begingroup$ This looks really interesting aip.scitation.org/doi/abs/10.1063/5.0089466 $\endgroup$
    – B. Kelly
    Jul 18, 2022 at 11:56
  • $\begingroup$ Thanks for the comment. Indeed, if there are atoms that are completely buried by other atoms and simultaneously highly charged, RESP does not know if the charge is carried by the buried atoms or the surface atoms, so it tends to assign close-to-zero charges to the buried atoms. $\endgroup$
    – wzkchem5
    Jul 18, 2022 at 19:13
  • 1
    $\begingroup$ @wzkchem5 The R in RESP does try to address the issue of buried atoms. RESP quite often does a good job, but when it doesn't, it really can go sideways. $\endgroup$
    – B. Kelly
    Jul 19, 2022 at 12:36
10
$\begingroup$

All Gaussian jobs when completed without any error, will have, at the end of the output file, a phrase/thought/quote from someone famous. It is kind off an easteregg. My latest Gaussian successful job finished with this quote:

THERE IS NO SCIENCE WITHOUT FANCY, NOR ART WITHOUT FACTS.
     -- VLADIMIR NABAKOV

Anyway, answering your second question:

What are the parameters that I should look out for, in the Gaussian output, to know that the structure has been optimised for use in MD?

If a successfully optimization has being reached, you will have the direct information below in your output file.

Just a remainder: the full optimization will depends on two parameters, the number of steps in the geometry moves and the threshold on the energy for the self consistent field (SCF).

     Item               Value     Threshold  Converged?
     Maximum Force            0.000022     0.000450     YES
     RMS     Force            0.000003     0.000300     YES
     Maximum Displacement     0.001224     0.001800     YES
     RMS     Displacement     0.000305     0.001200     YES
 Predicted change in Energy=-4.458789D-09
 Optimization completed.
    -- Stationary point found.
                           ----------------------------
                           !   Optimized Parameters   !

For your first question:

What are the parameters that need to be calculated to design the structure of a molecule?

The answer from @CharlieCrown is pretty much what you need to do/know.

$\endgroup$
1
  • 1
    $\begingroup$ This is actually the proper answer. When I went to sleep, I realized the OP probably meant a simple geometry optimization, which, to be sure the analytical derivative is zero, requires a frequency calc. I was thrown off by the use of such a small basis set. I do not recommend ever using 6-31G(d) $\endgroup$
    – B. Kelly
    May 30, 2020 at 19:02

You must log in to answer this question.

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