Embedded Files
ORCA can be used as a QM code in either QM or QM/MM geometry optimizations or molecular dynamics simulations in Chemshell. The current ORCA interface in Chemshell is documented on: http://www.cse.scitech.ac.uk/ccg/software/chemshell/manual/orca.html
When using ORCA with Chemshell, then Chemshell handles the coordinates, geometry optimization, molecular dynamics, QM-MM interactions etc. and only asks ORCA to calculate energies and gradients for a particular set of coordinates.
As the current ORCA-Chemshell interface requires pre-programmed Chemshell keywords that are used to create ORCA inputfiles, it can not take advantages of all the newest ORCA features. A rewritten ORCA interface is available here: orca-chemsh.tcl
Downloading this file and sourcing it in a Chemshell inputfile directly updates the ORCA interface in Chemshell. This new interface allows you to specify exactly in the Chemshell inputfile what SimpleInput and Blockinput keywords you want to be in the automatically generated ORCA inputfile.
QM geometry optimization with Chemshell and ORCA (new interface)
QM geometry optimization with Chemshell and ORCA (new interface)
An example inputfile using the newer interface:
test.chm:
c_create coords=test.c {
coordinates angstrom
Li 0 0 0
Li 0 0 4
}
# Sourcing the new ORCA interface
source orca-chemsh.tcl
# ORCA settings for simple input line
set orcasimpleinput " ! B3LYP RIJCOSX def2-SVP def2/J Grid4 FinalGrid5 D3BJ"
# ORCA block settings are specified here.
set orcablocks "
%scf
MaxIter 500
end
%pal
nprocs 6
end
"
#######
dl-find coords=test.c \
result=result.c coordinates=cartesian \
maxcycle=700 maxene=3000 tolerance=0.00045 \
theory= orca : [ list \
executable=/full/path/to/orca \
orcasimpleinput= $orcasimpleinput \
orcablocks= $orcablocks \
charge=0 \
mult=1 ] \
The Chemshell input above will create an ORCA inputfile that looks like this:
! B3LYP RIJCOSX def2-SVP def2/J Grid4 FinalGrid5 D3BJ Engrad
%scf
MaxIter 500
end
%pal
nprocs 6
end
%scf
AutoStart false #This part will change in the next optimization step.
end
%coords
CTyp xyz
Charge 0
Mult 1
Units bohrs
coords
Li 0.0000000000 0.0000000000 0.0000000000
Li 0.0000000000 0.0000000000 7.5589066540
end
end
Nudged elastic band (NEB) reaction path optimization with Chemshell and ORCA
If the intent is only to find the transition state then the NEB method is rather slow. Yet, since the input is only the reactant and product it can be convenient, especially for complex transition states that are hard to find by simple surface scans. Since the NEB method uses only gradients, its strength is in finding reaction paths and transition states for large molecular systems where calculation of the Hessian becomes too expensive.
c_create coords=react.c {
coordinates angstrom
C -1.74467 -0.436835 -0.00184383
Cl 0.0874915 -0.438907 0.00877167
H -2.09273 -1.22942 0.662529
H -2.07667 -0.617326 -1.02476
H -2.08697 0.537709 0.349213
F -4.40505 -0.450349 0.0900013
}
c_create coords=prod.c {
coordinates angstrom
C -2.50231 -0.439974 0.0202477
Cl 0.457516 -0.434744 -0.0210778
H -2.12542 -1.22786 0.677941
H -2.15121 -0.619108 -0.999564
H -2.13381 0.529453 0.366348
F -3.86338 -0.442901 0.040015
}
# Sourcing the new ORCA interface
source orca-chemsh.tcl
# ORCA settings for simple input line
set orcasimpleinput " ! PM3 TIGHTSCF "
# ORCA block settings
set orcablocks "
%scf MaxIter 500
end
%pal
nprocs 1
end
"
#######
dl-find neb=free coords= react.c coords2=prod.c \
result= result.c \
maxcycle=700 \
theory=orca : [ list \
executable=/full/path/to/orca \
charge=-1 \
mult=1 \
orcasimpleinput= $orcasimpleinput \
orcablocks= $orcablocks ] \
Dimer method TS optimization with Chemshell and ORCA
The dimer method searches for the transition state when the input geometries are given, ideally both close to the real TS. Since it only uses gradients it may be convenient for systems where calculation of the Hessian is too expensive.
# 007.xyz
c_create coords=dimer1.c {
coordinates angstrom
C -1.89114901629189 -0.43869113724225 0.01103837233128
Cl 0.06001816965482 -0.43550634752739 -0.01479978615730
H -2.11271779443058 -1.25066004555130 0.69629294199360
H -2.13576409736978 -0.62416113261704 -1.03000199025385
H -2.11954392095269 0.55776737346318 0.37558184980154
F -4.11944334060987 -0.44387871052520 0.04579861228474
}
# 008.xyz
c_create coords=dimer2.c {
coordinates angstrom
C -1.99283303076949 -0.43906204115263 0.01309547131865
Cl 0.08654162182067 -0.43514846620583 -0.01623959206210
H -2.09531444966205 -1.26265127699215 0.70686814906647
H -2.12108831451492 -0.62708960496606 -1.04443065188875
H -2.10329090994441 0.57197871671963 0.38153156676382
F -4.09261491692978 -0.44315732740297 0.04308505680191
}
# Sourcing the new ORCA interface
source orca-chemsh.tcl
# ORCA settings for simple input line
set orcasimpleinput " ! PM3 TIGHTSCF "
# ORCA block settings
set orcablocks "
%scf MaxIter 500
end
%pal
nprocs 1
end
"
#######
dl-find dimer=true coords= dimer1.c coords2=dimer2.c \
result= dimerfromscan-result.c \
maxcycle=700 \
theory=orca : [ list \
executable=/full/path/to/orca \
charge=-1 \
mult=1 \
orcasimpleinput= $orcasimpleinput \
orcablocks= $orcablocks ] \
QM/MM geometry optimization with Chemshell and ORCA
QM/MM geometry optimizations in Chemshell are documented in the Chemshell manual. Using ORCA as the QM code in QM/MM calculations just means that you have to modify the qm_theory part. As always in QM/MM, setting up the system and the MM region is the complicated part.
Below is an example input file (see Chemshell manual for details regarding the interface and MM settings):
qm-mm-test.chm:
source orca-chemsh.tcl
####################
# ORCA QM REGION SETTINGS
####################
# ORCA settings for simple input line
set orcasimpleinput " ! BP86 def2-SVP def2/J TIGHTSCF "
# ORCA block settings
set orcablocks "
%scf MaxIter 500
end
%pal
nprocs 6
end
"
#######
dl-find coords= ${sys_name_id}.c \
result= ${dir}/result.c \
residues= $residues \
active_atoms= $act \
maxcycle=700 \
theory=hybrid : [ list \
groups = $groups \
coupling=shift \
cutoff=1000.0 \
atom_charges= $charges \
qm_region= $qmatoms \
debug=no \
qm_theory= orca : [ list \
executable=/full/path/to/orca \
charge=0 \
mult=1 \
orcasimpleinput= $orcasimpleinput \
orcablocks= $orcablocks ] \
mm_theory=dl_poly : [ list \
list_option=none \
debug=no \
exact_srf=yes \
use_pairlist=no \
mxlist=16000 \
cutoff=1000 \
mxexcl=2000 \
debug_memory=no \
scale14 = { 1.0 1.0 } \
atom_types= $types \
use_charmm_psf=yes \
charmm_psf_file= ${sys_name_id}.psf \
charmm_parameter_file= $prm \
charmm_mass_file= $top ] ]
DFT Molecular Dynamics Simulation in Chemshell using ORCA as QM code
Chemshell can perform molecular dynamics simulations at the MM, QM or QM/MM level. To perform QM MD simulations using ORCA as QM code one only needs to specify ORCA as QM code, ORCA will calculate the atomic forces which the Chemshell MD module will use to calculate trajectories. Chemshell MD scripts are a little complicated but allows a lot of flexibility in changing the way the calculation is carried out.
md-test.chm:
c_create coords=initial.c {
coordinates angstrom
C -1.318564168 1.167866561 -11.695235956
...
...
}
source orca-chemsh.tcl
######
# ORCA settings for simple input line
# Here using the semiempirical PM3 method
set orcasimpleinput " ! PM3 TIGHTSCF"
# ORCA block settings
set orcablocks "
%scf MaxIter 500
end
"
#######
# Optional read-in of masses. Useful to read-in deuterium masses for hydrogen.
# Remember to add masses keyword to the dynamics module if you use this.
#source masses
########################################
# SIMULATION PARAMS
########################################
set TITLE simulation
set NSTEP 500000
# F=start new
set RESTART F
set RANDOMIZE T
set WRITETRA 1
set WRITERSTRT 1
set UPDATE 25
set TODO "MD"
#######################################
# Heating up steps. Increasing temperature at selected steps.
set tempa 10
set stepb 10
set tempb 200
#
set stepc 50
set tempc 250
#
set stepd 100
set tempd 400
#
set stepe 200
set tempe 600
#
#set finalstep 4000
#set finaltemp 300
#set finaltaut 0.02
dynamics dyn1 coords=initial.c \
timestep = 0.001 \
temperature = 10 \
ensemble = NVT \
nosehoover=4 \
taut=0.02 \
theory= orca : [ list \
executable=/full/path/to/orca \
orcasimpleinput= $orcasimpleinput \
orcablocks= $orcablocks \
charge=0 \
mult=1 ]
#
# MD Loop
#
# Organize restart, protocol, etc.
if { $RESTART == "T" } {
dyn1 load
# Get the last step no. from protocol (if it exists)
if { [ file exists $TITLE.prot ] && [ file writable $TITLE.prot ] } {
set fprot [ open "|tail -1 $TITLE.prot" r ]
gets $fprot line
close $fprot
scan $line %d last
set first [ expr $last + 1 ]
set fmode a
} else {
set first 1
set fmode w
}
} else {
# No restart
set first 1
set fmode w
}
exec rm -f dynamics_tra.xyz
set fprot [ open $TITLE.prot $fmode ]
puts $fprot \
{ #Step time[ps] E_pot[au] E_kin[au] E_tot[au] T[K] F_c[au] Fric }
# Remove existing EXIT file
if [ file exists EXIT ] { exec rm EXIT }
# Initial velocity distribution
if { $RANDOMIZE == "T" } { dyn1 initvel }
########################################
# DYNAMICS PROCEDURES
########################################
# Take a step
proc mdstep { } {
global istep fprot WRITETRA WRITERSTRT UPDATE stepb stepc stepd stepe tempb tempc tempd tempe finalstep finaltemp finaltaut
if { ! ($istep % $UPDATE) } { dyn1 update }
dyn1 force
dyn1 step
dyn1 printe
if { ! ($istep % $WRITERSTRT) } { dyn1 dump }
if { ! ($istep % $WRITETRA) } { writetra full }
# if { ! ($istep % 10) } { dyn1 output }
if { ($istep == $stepb) } { dyn1 configure temperature = $tempb ensemble=NVT }
if { ($istep == $stepc) } { dyn1 configure temperature = $tempc ensemble=NVT }
if { ($istep == $stepd) } { dyn1 configure temperature = $tempd ensemble=NVT }
if { ($istep == $stepe) } { dyn1 configure temperature = $tempe ensemble=NVT }
# if { ($istep == $finalstep) } { dyn1 configure temperature = $finaltemp ensemble=NVT taut=$finaltaut }
puts $fprot [ format "%8d %14.4f %12.6f %12.6f %12.6f %14.2f %13.5e %12.3e " \
$istep \
[ dyn1 get time ] \
[ dyn1 get potential_energy ] \
[ dyn1 get kinetic_energy ] \
[ dyn1 get total_energy ] \
[ dyn1 get temperature ] \
[ dyn1 get constraint_force ] \
[ dyn1 get friction ] ]
flush $fprot
return 0
}
# Write trajectory
proc writetra { {type full} } {
switch $type {
full {
# The standard ChemShell trajectory
dyn1 trajectory
puts "--- ChemShell trajectory written."
#write_xyz coords=hybrid.mndo.coords file=temp1.xyz
#exec cat temp1.xyz >> dynamics_tra_mndo.xyz
puts "--- MM xyz trajectory written."
}
none { }
default {
puts "*** Invalid trajectory type: $type."
puts "*** No trajectory written."
return 1
}
}
return 0
}
# Closing down
proc mdstop { } {
global fprot
dyn1 dump
close $fprot
dyn1 destroy
times
return OK
}
########################################
# SIMULATION LOOP
########################################
for { set istep $first } { $istep < $first + $NSTEP } { incr istep } {
# Try to catch errors and abort gracefully (1 if error occurred)
#set errorflag [ catch {
mdstep
#} ]
#if { $errorflag } {
# puts "*** Error during MD step. Aborting..."
# mdstop
# }
# Allow graceful termination via EXIT file (1 if file exists)
set exitflag [ file exists EXIT ]
if { $exitflag } {
puts "*** EXIT file detected. Terminating..."
mdstop
}
}
Page updated
Google Sites
Report abuse