Talk:Mod:Hunt Research Group/QMMM MD Aqeuous Cu(II)
This input runs a short 0.1 ps (timestep = 0.001 ps) QM/MM molecular dynamics simulation of an aqueous copper system.
#
# Create the dynamics object
#
dynamics dyn1 coords=cu_18water_qmmm_opt.pun \
theory= hybrid : { coupling=shift qm_region = {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22} \
qm_theory=gaussian : { nproc=8 basis=6-311g g98_mem=80000000 charge=2 mult=2 hamiltonian=b3lyp } \
mm_theory=dl_poly : mm_defs=ff.dat } timestep = 0.001 ensemble = NVE
# initialise random velocities
dyn1 initvel
dyn1 output
#
set count 0
while {$count < 100} {
# Update MM pairlist every 10 steps
if { [ expr $count % 10 ] == 0 } { dyn1 update }
dyn1 force
# output the coordinates to the trajectory every 50 steps
if { [ expr $count % 10 ] == 0 } { dyn1 trajectory }
dyn1 printe
dyn1 step
incr count
}
dyn1 destroy
times
The initial configuration is taken from the QM/MM optimisation. To run this, login to cx1 and submit this simulation into the PBS queue.
NB: When I ran this simulation, ChemShell was throwing weird errors. It only worked when I put the code defining the dynamics object dyn1 into a single line.
The simulation will take a while (~ 1 hour). Once finished, you should see md.out and dynamics.trj in your working directory.
If there were no errors, the end of md.out should look something like this:
Energy and gradient evaluations:
QM energies 100
QM gradients 100
MM energies 100
MM gradients 100
Timing Analysis by Executable
g09 [ 100] 3659.293 99.9%
ChemShell 2.041 0.1%
Total 3661.334
Timing Analysis by Module module # calls time time (exc) percent
dynamics [ 1] 3661.144 0.057 0.0%
hybrid.kill [ 1] 0.001 0.001 0.0%
dl_poly.kill [ 1] 0.000 0.000 0.0%
hybrid.calc [ 100] 3661.024 0.038 0.0%
dl_poly.eandg [ 100] 0.237 0.237 0.0%
gaussian [ 100] 3660.750 3660.342 100.0%
set_matrix_element [ 16600] 0.147 0.147 0.0%
set_matrix_size [ 200] 0.005 0.005 0.0%
get_bq_entry [ 6600] 0.084 0.084 0.0%
get_ecp [ 100] 0.116 0.074 0.0%
get_element_data [ 2200] 0.014 0.014 0.0%
hybrid.update [ 10] 0.015 0.003 0.0%
dl_poly.update [ 10] 0.011 0.011 0.0%
hybrid.init [ 1] 0.047 0.031 0.0%
dl_poly.init [ 1] 0.016 0.016 0.0%
get_cell [ 1] 0.000 0.000 0.0%
connect [ 1] 0.000 0.000 0.0%
get_number_of_bqs [ 101] 0.002 0.002 0.0%
get_number_of_shells [ 3] 0.000 0.000 0.0%
get_connected_atoms [ 22] 0.000 0.000 0.0%
get_atom_entry [ 6622] 0.080 0.080 0.0%
copy_object [ 1] 0.000 0.000 0.0%
check_molecule_connectivity [ 2] 0.000 0.000 0.0%
get_number_of_atoms [ 303] 0.003 0.003 0.0%
Other 0.190 0.0%
Total 3661.334
ChemShell exiting code 0
dynamics.traj contains the trajectory - the coordinates of each atom at each step. A set of coordinate at each time step should have this form:
block=update_coordinates records=55 cu -3.852296 6.526986 4.214061 o -0.108291 0.662057 0.875603 h 1.092434 1.823054 -0.012381 h 0.627664 0.030774 2.426616 o -4.631849 8.838777 1.580035 h -3.431835 10.206983 1.094289 h -5.499046 8.200385 0.033426 o -4.920261 8.782003 6.914051 h -4.111675 10.210162 7.668491 h -6.819169 8.713904 6.867511 o -9.448665 7.685461 6.340187 h -10.770518 8.515366 5.275003 h -10.228182 6.378659 7.420263 o 0.085853 8.206543 4.709663 h 1.766433 8.007057 3.966225 h 0.046867 9.572561 5.956115 o -4.029740 3.586616 1.719834 h -2.593939 2.413401 1.398341 h -5.662705 3.028327 1.047417 o -3.424548 3.720787 6.727423 h -3.250861 4.121120 8.551552 h -3.779189 1.917627 6.378730 o 4.995600 7.367619 3.005797 h 6.008234 6.186735 4.094361 h 6.206435 8.666970 2.315186 o -4.363873 -2.437806 0.960948 h -4.448308 -4.171186 0.167384 h -2.786802 -1.608173 0.373933 o -11.423871 3.700929 8.978197 h -12.947247 3.918002 10.049199 h -10.751475 1.952781 9.341881 o -0.514082 12.112154 8.287439 h 0.199786 11.966829 10.068623 h -0.441351 13.994396 8.003580 o -1.970013 12.709458 0.040809 h -2.568429 14.494590 0.244146 h -0.555362 12.687158 -1.259090 o -3.869360 -1.065867 5.736818 h -4.166092 -1.866178 4.012471 h -4.067672 -2.475658 6.994878 o -13.061802 9.637356 3.877019 h -14.925099 9.634282 4.282584 h -12.774033 10.948868 2.555701 o 3.063029 3.913205 -0.861599 h 4.049648 4.013092 -2.499460 h 3.933026 4.893652 0.441136 o -7.172084 7.174326 -2.295341 h -8.128873 5.649121 -1.793073 h -7.870503 7.719024 -3.988433 o -8.804313 2.661263 0.146495 h -9.216160 1.127254 -0.908414 h -10.320564 3.069537 1.204719 o -3.268579 4.627503 11.687971 h -4.225436 5.978500 12.669773 h -1.970046 3.930527 12.896909
NB: The trajectory cannot be directly visualised in VMD. I am currently writing a script to convert a general ChemShell trajectory file into a VMD-readable one.