ORCA & Chemshell

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)

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

}

}