Talk:Mod:Hunt Research Group/Starting MD
create the starting simulation box for your MD package
- for gromacs (-g) or lammps (-l) or dlpoly (-d)
fftool/fftool 256 c4c1im.zmat 256 ntf2.zmat -r 3.4 -g
- for gromacs this provides the files:
- field.top (contains the force field)
- config.pdb (contains the atom coordinates)
- run.mdp (default control file for md run)
GROMACS now prepare file to run (energy relax) your simulation
- now prepare your control file we will call this en_min.mdp
- don't use the default from fftool run.mdp this is a md file, not an energy minimization
- create an energy minimization en_min.mdp file en_min.mdp_file
- this uses steepest descent
- we need to generate the file that will run the simulation
- the *.tpr file actually runs the simulation
- this is a binary file and so not readable
gmx grompp -f en_min.mdp -c config.pdb -p field.top
- produces the files
- topol.tpr (binary)
- mdout.pdp (record and explanation of the options)
- run the relaxation
- the -v flag causes the potential energy and max force to be printed each step
gmx mdrun -v
- produces many files!
- confout.gro (structure gromacs format)
- ener.edr (binary energies,temp,pressure, box size, density etc)
- mdout.mdp (run parameters)
- md.log (log file)
- traj.trrr (binary x,v and f)
- use vmd to look at your final structure
- cp confout.gro conf.gro
- vmd: source load.tcl load
- gmx editconf -f confout.gro -o confout.pdb
- head confout.pdb for the cell dimensions
- view cell boundaries "pbc set { 63.001 63.001 63.001}"
- "pbc box_draw"
DLPOLY now prepare file to run (energy relax) your simulation
- for gromacs (-g) or lammps (-l) or dlpoly (-d)
fftool/fftool 256 c4c1im.zmat 256 ntf2.zmat -r 3.4 -d
- now carry out an NVE simulation (this is essential to get a good starting point)
- we have a fixed volume, and we want to relax the individual molecules in other MD programs you would do a minimisation, however we will give it a little bit of energy to allow things to move around a bit. This step is more important for ILs which are less fluid, or for larger solvent molecules where moving a molecule costs a bit more energy. For something like water can go straight to room temperature, but we will go through the motions so that you have experience.
- process
- the minimisation options in DLPOLY don't work very well for ILs
- do T = 5K for 50 steps
- use smaller timestep (e.g. 0.2fs) to make pressure tensor ~E+00
- if use large timestep e.g. 2fs, pressure tensor ~E+06
- do T=400K for 1000 steps
- start from the above relaxed structure by copying the REVCON file into a new directory and calling it CONFIG
- check that the pressure tensor is small ~E-01/-02.
- you need a CONTROL file
- here is our example
- for the shake algorithm use 1*10-6 if you have pdb or inaccurate geometry, if you have xyz 1*10-8 is better. The shake constrains all the H-bonds to this bond length.
dlpoly run temperature 5.15 #pressure 0.001 ensemble nve #ensemble npt hoover 0.2 0.2 integrator leapfrog shake 1.0E-05 quaternion 1.0E-08 steps 5000 equilibration 500 scale 100 print 100 stack 100 stats 100 timestep 0.001 cutoff 16.000 delr 0.100 rvdw 15.000 ewald precision 1.0E-6 traj 1 100 0 job time 1000000.0 close time 1000.0 finish
- start from a grid or crystal structure and melt the solid
- T=1000K for 20ps (20,000 steps)
- use capping, default=1000 but use 2000 kbT per Å, capping constrains the forces on individual atoms to be smaller or equal to a particular value and resets large forces to capping the capping value
- use scaling, every 5 steps, scaling scales all atomic velocities by the same amount to reproduce desired temperature (??)
- test with 100 steps, no scaling or capping of forces, and step ≈1fs, check that temp is not exploding
- if there is a problem, we want to nudge system back to a good position, problem may well be two atoms coming too close together and repelling each other strongly
- shorten the time step 0.1fs
- add capping, default of 1000
Common problems
e.g. box may explode check the potential for the correct units are two cations next to each other - start from crystal structure scale forces using the 'cap (forces)' keyword sets the max forces to what is in the brackets, careful this can be dangerous, run a second equilibration without the cap keyword.
old for dlpoly
- old instructions to convert the packmol to a dlpoly readable file
- use this script to convert the output of packmol into a dlpoly readable file conv_xyz2conf.pl
- dlpoly is highly format driven so check the dlpoly manual for the exact details
- essentially each atom is now numbered with the xyz coordinates on the following line
- notice that in the script specific atomic labels are searched and printed so if you have a different molecule you will need to edit the script to put in the new atom types you need.
- input is the filename you wish to convert excluding the .xyz extension
- for example conv_xyz2conf nacl_1330-spce
- this file needs to be copied to the CONFIG file and appropriately edited
- so nacl_1330_spce.xyz looks like this
3992 Built with Packmol NA 0.001779 9.402974 9.416280 CL -9.999671 -0.580455 -0.000345 OH2 -1.549154 16.628984 6.058143 HO -1.166292 16.884755 6.945835 HO -0.975087 16.999934 5.328183 OH2 -3.418933 0.347746 -2.288811 HO -3.805476 1.262661 -2.172555 HO -3.388480 0.119970 -3.262048 OH2 -15.941200 6.736377 16.160857
- this will produce a file with the same name but extension .cfg
- and you will get something that looks like this
Converted from Packmol XYZ
0 0 0 0.000000
0.000000000000 0.000000000000 0.000000000000 #Note: You are required to fill in the correct numbers here!
0.000000000000 0.000000000000 0.000000000000 #Note: Do not change the spacing!!!!
0.000000000000 0.000000000000 0.000000000000
NA 1
0.001779 9.402974 9.416280
CL 2
-9.999671 -0.580455 -0.000345
OH2 3
-1.549154 16.628984 6.058143
HO 4
-1.166292 16.884755 6.945835
HO 5
-0.975087 16.999934 5.328183
OH2 6
-3.418933 0.347746 -2.288811
HO 7
-3.805476 1.262661 -2.172555
HO 8
-3.388480 0.119970 -3.262048
OH2 9
-15.941200 6.736377 16.160857