View source for Mod:Hunt Research Group/Aug09QtoRuth

Summary of discussions with Ruth about classical MD simulation of ILs (12-Aug-09, Cambridge) (Notes by Ling)

1. Temp in NVE

NVE when we set the temp, this sets the initial temp (kinetic energy), but this doesn't mean the simulation will end at the set temp. Total energy = kinetic energy + potential energy (also called 'configurational energy'). Normally we don't have the potential energy corresponds to what would have at the set temp, so the simulation would often finish higher or lower than the set temp. 

2. RESTART

Q: is the jump in temp from "equilibration" to "thermostat" an issue?

A: No.

Note: If 'RESTART' is used, the simulation would restart from the velocity of the last run. If 'RESTART' is not used, it will restart from a random velocity just keeping the kinetic energy corresponds to the set temp. Don't use 'RESTART' except in special circumstances. Look at output for first step of thermostat. 

From the manual 2.20 P86 about how to restart:

"You can also extend a simulation beyond its initial allocation of timesteps, provided you still have the REVCON and REVIVE files. These should be copied to the CONFIG and REVOLD files respectively and the directive timesteps adjusted in the CONTROL file to the new total number of steps required for the simulation. For example if you wish to extend a 10000 step simulation by a further 5000 steps use the directive timesteps 15000 in the CONTROL file and include the *restart* directive. You can use the *restart scale* directive if you want to reset the temperature at the restart, but note that this also resets all internal accumulators (timestep included) to zero."

3. Thermostat:

Nose-Hoover: small oscillations in temperature. To do: check the time step in sample CONTROL files. Try larger timestep.

Berendsen: no oscillations in temperature. 

The longer the timestep parameter in the thermostat, the gentler the thermostat.  

4. How to compute MSD of whole molecule?

Ref: Ruth's immsd.f

Note: when |x1-x0|>box/2, delta x = x1-x0, apply pbc to delta x -> delta x new, and x1 new = x0 + delta x new

To do: hack DL-POLY and write my own MSD code based on HISTORY file (fortran core + C wrapper). 

Structure of code: 

5. How to compute diffusion coefficients 

Note: have to run for a long long time to get the correct slope (three stages). The diffusion coefficients in the current literature might be getting the slope of the second stage.

Ref: Liu06. The self-diffusion coefficient of component i can be calculated by the Einstein relation (Ref: Mor02) D_i = 1/6lim_{t→infinity}d/dt<[r_i(t) - r_i(0)]²>

where the quantity in <> is the ensemble-averaged mean square displacement (MSD) of the species i, and r_i is the vector coordinate of the centre of mass. The self-diffusion coefficients were obtained from the slope from 20-100ps and applying the above equation.

6. how do we get the kinetic and potential energy out of DLPOLY files?

Kinetic energy = temperature

Potential energy = configuration energy

7. Forcefield

Set the ring and attached atoms 'rigid' (e.g. the numbering is 1-10), then I need to delete from 'constraints', 'bonds', 'angles' and 'dihedrals' all the lines that only involve numbers 1-10, but leave the rest. I have already constrained C-H distances. (Note: don't forget to delete the three constraints of C-H bonds attached to the ring or it would result in the error 'SHAKE algorithm fail to converge'!) 

8. Viscosity

to calculate viscosity independently from diffusion, we can print pressure tensor to the HISTORY file (note just choose certain atoms or the file size would be huge). 

'shear viscosity': from off-diagonal elements of pressure tensor. When people talk about 'viscosity' in the literature, they normally mean 'shear viscosity' (imagine sliding layers). 

'bulk viscosity': from diagonal elements of pressure tensor. (Imagine constant volume, elongate x, shorten y, z). 

viscosity and diffusion should be related by the Nernst-Einstein equation, but for IL relation may not hold.

People in the literature normally use Stoke-Einstein relations to get viscosity from diffusion rather than independently. Computing viscosity independently can be way-off, e.g. Ruth's former student got negative viscosity.

9. 3D probability density map (probability distribution of anions around bmim, I can later develop this program to alkyl chains rotating around the ring). 

a. use im3ddir.f to get .pdens file. This program prepares 3d histograms from HISTORY file. Format of .pdens file (1st line: no. of grid points (e.g. 41), 2nd line: size of box (in A), 3rd line: coordinates of bottom left point of the box)

b. use ATEN

- load cation coordinates, e.g. dmim.xyz
- File -> open grids
- grid window: change colour or change cut-off (e.g. to 12A, like isosurfaces)

I have download this software and read a document about it (DOI 10.1002/jcc.21359).

Note: Several common functional forms for intramolecular and intermolecular energies/forces are implemented, suitable for the redefinition of many general forcefields within Aten, e.g., OPLS-AA
or AMBER, or specific forcefields such as Canongia Lopes ionic liquids parameters [Lop04]

FYI: can download it from 
http://www.projectaten.net/tiki-list_file_gallery.php?galleryId=3

Choose:
aten-latest-mac.dmg	Latest Macintosh installer	Source revision 987	13.42 MB