Matter Modeling Stack Exchange is a question and answer site for materials modelers and data scientists. It only takes a minute to sign up.
Sign up to join this community
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?
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.
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".
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.
To directly answer your question now, Force-fields take a lot of work. If you are making one from scratch you need to
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.
-- 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.
