Starting configuration for a molecular simulation

Question

When doing a computer simulation, what’s the best way to prepare a starting configuration to avoid biasing your results?

The book Computer Simulation of Liquids suggests putting the molecules on an fcc lattice, but is that optimal?

I know there are programs like PACKMOL that try to pack the molecules in a more or less random starting configuration.

Or, as I've been told, should you set up the initial configuration using one of the previously mentioned methods (which one?) and run it at high temperature and then cool it down (if so why, how high, and how exactly should you cool it down)?

Finally, how does the answer vary based on the type of system you wish to simulate, e.g., molecular liquid, polymer, protein, ...?

Answer 1

All depends on what are you working on.

When doing a computer simulation, what’s the best way to prepare a starting configuration to avoid biasing your results?

If working with molecules in gas phase (isolated systems):

• you can draw them in your favorite program, save them as 3D model, and then do a simple Molecular Mechanics geometry optimization.
• you can download the file with the structure from an online databases. Some of them already have the structure in 3D, but some has in 2D (smile format, for example). In the later case, you have to convert to 3D and do a geometry optimization.

If you are interested in working with molecules in solid phase (periodic systems):

• search for crystallographic information in online databases.
• obtain the 3D structure (see above), and use specialized software that generate random crystals that uses as inputs the molecules and the crystal information.

The book Computer Simulation of Liquids suggests putting the molecules on an fcc lattice, but is that optimal?

If you are working with isolated system, the answer is no. If you are working with periodic system, does your system has a FCC symmetry? If yes, ok. If not, then, no, it is not optimal.

I know there are programs like PACKMOL that try to pack the molecules in a more or less random starting configuration.

The idea behind PACKMOL is not pack molecules. The idea is to help setting an initial structure file for running molecular dynamics simulations for specific systems.

Or, as I've been told, should you set up the initial configuration using one of the previously mentioned methods (which one?) and run it at high temperature and then cool it down (if so why, how high, and how exactly should you cool it down)?

Again, all depends on what are you working on. For example, if you need the electronic structure of a molecule (molecular orbitals, electron density, charges, electrostatic potential), you don't need to run any simulation like that. If you are studying how a group of molecules crystalize or how a system of magnetic particles (spins) behave, then, yes, it is useful.

Finally, how does the answer vary based on the type of system you wish to simulate, e.g., molecular liquid, polymer, protein, ...?

As wrote above, the correct answer strongly depends on what you want.

Answer 2

In order to perform a molecular dynamics simulation, you need to equilibrate the system before you can get good statistics.

In principle, you can use any starting point, since the simulation should be ergodic; however, this will mean that the equilibration time will be very long.

Placing the molecules on a lattice is a simple way to form the starting point, but it may also be very far from equilibrium: if you want to simulate a mixed liquid, but your molecules don't have uniform shape (e.g. size in one dimension is much bigger in the two other ones), your initial lattice will be quite far from the one you're targeting at equilibrium, which will have an uniform shape.

The idea in PACKMOL is to obtain a randomized starting structure. This will allow for a much faster equilibration, since you can e.g. form the initial structure in a simulation cell that already has a density that is close to the equilibrated one, at the wanted temperature and pressure.

Now, let's assume you have a starting configuration. The next thing to do is to equilibrate it. Like you said, there are many ways to accomplish this. Generally, one runs the initial simulation at high temperature, since this will allow the system to escape any local minima. If the density is already close to correct, the initial equilibration can be run at constant volume. Otherwise, as far as I understand, it is common practice to run the initial equilibration at a large pressure. The idea is again to allow the system to "mix" faster, so that the equilibration takes less CPU time.

What is the right temperature (and pressure)? This probably depends a lot on what you are trying to do. Obviously, there are a few practical issues to keep in mind. You don't want to set the temperature or pressure so high that you need to decrease your time step, since this would beat the purpose of speeding up the equilibration. You also don't want your system to dissociate; however, since molecular MD potentials often don't allow dissociation due to harmonic potentials, this may not be a big issue. The specifics will depend on your system, and what you are trying to accomplish. For instance, if you are simulating a protein in a given a given folding pattern, increasing the temperature and pressure for the protein may be a bad idea. The pressure, in turn, should be high enough so that the system does not evaporate, i.e. that the density does not decrease a lot.

Once you have run the initial equilibration of the system at high temperature and pressure, you can start the production run. However, since the temperature and pressure/volume will change, you should discard the initial part of the simulation from the analysis. Key quantities to look at are the different parts of the energy, density, etc; you should only look at the part of the simulation where these are oscillating around the equilibrated value.