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

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?

• 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. – Camps May 12 '20 at 16:46
• @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? – Kavya Mrudula May 13 '20 at 9:15
• 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? – Camps May 13 '20 at 10:23
• @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. – Kavya Mrudula May 14 '20 at 20:22

## 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.

The 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 vaccum 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 must be 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.

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.

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   !