Orca

Orca is a free software for all types of QM calculations.
Documentation and tutorials are here: https://www.orcasoftware.de/tutorials_orca

The current version is 5.0.3.
It is installed on Aurora:
module add GCC/11.2.0  OpenMPI/4.1.1 ORCA/5.0.3-static
Then, it is run by
orca input_file >output

Avogadro is a good option to illustrate Orca output


This is a simple sample input file:

! RHF def2-SVP          # Level of theory.
 ! SP                             # Single-point energy calculation; it is default or
! OPT                          # Geometry optimisation
* xyz 0 1                     # Charge, and multiplicity followed by coordinates in Å (can be changed to a.u. by %coords Units bohrs)
 h           1.57623125918282     -1.40098816005358      1.92944524289998
c           3.11804788185000     -1.26044723466300      3.33347664459600
h           4.63927730299500     -2.45475405971100      2.62293967273200
h           3.85126156558200      0.66329382213900      3.40150678020000
h           2.59270405690800     -1.91429242685700      5.21564372964000
*

Alternatively, internal coordinates can be used:
* int 0 2
 O 0 0 0 0.0       0.0 0.0
 H 1 0 0 0.9903 0.0 0.0
 *

Alternatively, an xyz file can be used:
* xyzfile 0 2 hydroxide.xyz

File hydroxide.xyz:
2 # number of atoms; next line blank or title
O 0 0 0
 H 0 0 0.9903

Alternatively, an pdb file can be used:
* pdbfile 0 2 hydroxide.pdb
It seems that the pdb file must contain at least one row with occupancy and B-factors; otherwise, you get strange error messages:
HETATM    1  O     2     2     -20.282   7.868  -0.800  1.00  0.00

Input file is not case sensitive.

Constraints:
%geom
   Constraints           
   {B 0 1 1.25 C}    # Bond length constraint
    end
  end

Cartesian coordinate constraints:
%geom
   Constraints 
    {C 3 C }        # Atom 4 is fixed         
    {C 10:17 C}   # A list of atoms is also allowed
     end
   end

Note that Orca number atoms starting with 0, which may be very confusing in the beginning.

Options

Use different basis sets for different atoms

! B3LYP def2-SVP
%basis
  NewGTO C "def2-TZVP" end # Applies TZVP to all C atoms
  NewGTO 6 "def2-TZVP" end # Applies TZVP to all C atoms
 end 

Unfortunately, such simple input cannot be used for individual atoms (a 
number is interpreted as an element number, not an atom number.

For individual atoms, it needs to be given in the coordinate input, e.g.:

! TPSS D4 DEF2-SV(P)



* xyz -1 1

H       2.384   3.207   1.496

S      -0.935   0.285   2.326 newgto "ma-def2-SV(P)" end

H      -2.837  -2.703   1.492

*



Strangely, it seems that the coordinates then must be in the main input 
file, not in a file reference. You can have it in the pdb file, but it 
is ignored.



However, you may use fragments to change basis sets of certain atoms 
(note that the fragment definition changed between Orca 6.0 and 6.1):

! TPSS D4 DEF2-SV(P)

%PAL NPROCS 24 END

%basis

FragBasis 1 "def2-SV(P)"

FragBasis 2 "ma-def2-SV(P)"

end



* pdbfile -1 1 coord.pdb

%frag

 Definition

  1 {0:25 27:107 113:114 116:146 150:164} end

  2 {26 108:112 115 147:149} end 

 end

end





Use RI and RIJCOSX (exchange)

  !B3LYP def2-SVP def2/J RIJDX RIJCOSX


In fact, this is default for DFT; to turn it off, use !NORI and !NOCOSX

  

Dispersion correction

!D3 or !D4


  

Increase amount printed out (e.g. to print orbitals)

!HF def2-SVP Largeprint



Frequency calculation

!FREQ


!TightOpt or !VeryTightOpt may remove imaginary freqs



Implicit solvent (COSMO is no longer any option)

!CPCM(solvent)

or 

%CPCM


 EPSILON 80


 REFRAC 1.33


END

or

%CPCM SMD TRUE


 SMDSolven "solvent"


END



Relativistic corrections

!B3LYP ZORA ZORA-def2-TVP RIJDX SARC/J



DLPNO CC

!DLPNO-CCSD(T) def2-TZVPP def2-TZVPP/C



Change maximum number of geometry steps

%geom



maxiter 1               



end



Cartesian optimisation

! COpt



SCF options

%scf


MaxIter 300


CNVDIIS 1


CNVSOSCF 1


end





%scf    


        guess moread


        moinp "TS2-opt.gbw"


        maxiter 999


end 





Level shift

%scf


 CNVShift true # default = true


 LShift 0.1 # default = 0.25


 ShiftErr 0.1 #default = 0.0


# more convenient:


 Shift Shift 0.1 ErrOff 0.1 end


end



Output options

%output


print[p_mos] true


print[p_basis] 5


end



Paralell calculations

%pal


nprocs 
20              
 #Maximum number of cores of your machine.


end



or simply pal16, e.g.

!QMMM PBE0 6-31G* D3ZERO Opt pal16



Point charges with file reference

%pointcharges "pointcharges.pc"



File pointcharges.pc

418 # number of point charges


 0.294300            
 6.494988     -4.096926     
-5.179739  # charge and coordinates


 0.164200            
 7.320798     -3.807798     
-6.874823 


 0.164200            
 6.417509     -5.943188     
-4.694079 

...

No special ending

Atomic units



Geometry options

%geom


 MaxIter 200


 TS_search EF


 # ! OptTs 


 TolE 5e-6


 TolRMSG 1e-4


 TolMaxG 3e-4


 TolRMSD 2e-3


 TolMaxD 4e-3


 inhess GFNFF


end



TS optimisation

!QMMM PBE0 6-31G* D3ZERO OptTS SlowConv TightSCF pal16



Scan a reaction coordinate

%geom scan


 B 11 2275 [2.0 1.9 1.8] end 


end  



Dihedral:


%geom


  scan


    D N1 N2 N3 N4 = initial-angle, final-angle, Nsteps


  end


end




Steps of 10–15º are appropriate, not bigger





Read a wavefunction from an old calculation

! MOREAD


%moinp "myorbitals.gbw"

For single-point calculations, if a *.gbw file of the same name is 
found, this file is renamed to *.ges and the orbitals are used as start 
orbitals. If this is not descired, set AUTOSTART false


  

Example for a system that is hard to converge the SCF (nitrogenase)

! TPSS D4 DEF2-SV(P) UKS

! ENGRAD

! NOTRAH


! SlowConv

! TightSCF

%scf

  MaxIter 1000

  DIISMaxEq 40

  DirectResetFreq 1

end

%PAL NPROCS 48 END

*pdbfile -2 1 qm_1_h.pdb




NOTRAH - turns off a newer (failsafe) scf converger

SlowConv - Selects appropriate SCF converger criteria for difficult cases. Most transition metal complexes fall into this category

VerySlowConv - Try this for very difficult cases

%MaxCore 2000



Flip spin

First converge HS state and then flip to BS stat

!UKS DEFGRID3 R2SCAN D4 ZORA

!TightSCF TRAH



%scf

SCFMode Direct

Thresh 1e-15

TCut 1e-15

Flipspin 0 # Here atom number zero is flipped

FinalMs 0.0 # Here the final broken-symmetry state to converge to is selected (here, MS=0)

end



*xyzfile -1 11 fe2h_tzvp.xyz



TD-DFT calculations

! CAM-B3LYP DEF2-SV(P) def2/j  RIJCOSX UKS

! Opt

! CPCM

%cpcm

  epsilon 4.0

end

%tddft

  nroots 100

  iroot 1

  maxdim 7

  maxiter 5000

end

%PAL NPROCS 32 END

*pdbfile 0 1 qm_1_h.pdb



Results are found after ABSORPTION SPECTRUM




Output

File *.in contains the input

*.out the output (very long)

*.engrad energy, gradient and coordinates.


*** FINAL ENERGY EVALUATION AT THE STATIONARY POINT ***


Dispersion correction           -0.077894975


FINAL SINGLE POINT ENERGY     -2474.907119623929


FINAL SINGLE POINT ENERGY (QM/MM)    -2507.427748516620





Scans:


*** OPTIMIZATION RUN DONE ***


RELAXED SURFACE SCAN STEP   6





RELAXED SURFACE SCAN RESULTS


----------------------------





Column   1: NONAME





The Calculated Surface using the 'Actual Energy'


   2.00000000 -2507.39966580


   1.90000000 -2507.42805085


   1.80000000 -2507.42978390





The Calculated Surface using the SCF energy


   2.00000000 -2474.80098184


   1.90000000 -2474.82795119


   1.80000000 -2474.82969327







QM/MM with Orca

Described in the Orca manual under 
		
	
	
		

8.13 ORCA Multiscale Implementation
  
				
			
		


	


You need a prmtop file and a pdb file (and preferably a syst1 file)



1. Make an Orca input file:

!QMMM PBE0 6-31G* D3ZERO





%qmmm


    ORCAFFFilename   "prmtop3.ORCAFF.prms"


    QMAtoms          { 0:33  


                      














 642:653


                      














 2238:2239


                      














 2248:2275


                      














 4716:4718 } end


    ActiveAtoms      { 0:33  


                      














 642:653


                      














 2238:2239


                      














 2248:2275


                      














 4716:4718 } end


end





*pdbfile 0 1 pdb3.pdb 



Made from the syst1 file, but renumbering to start on 0, changing "-" to ":" and removing link atoms.

The link atoms should be removed from QMAtoms list, but not from ActiveAtoms list.



2. Convert prmtop file to Orca format

orca_mm -convff -AMBER prmtop

Gives you the file prmtop3.ORCAFF.prms

It seems to have problems with metals and other two-character elements.



3. Finally start Orca with

orca orca.in >orca.out



This does additive QM/MM and electrostatic embedding.

You may change the charges around the junction atoms (p. 375). Orca 
implements the Charge Shift (CS; default), Redistributed charge and 
dipole (RCD), as well as the Z0, Z1, Z2 and Z3 methods.

%qmmm


   chargeAlteration RCD


end



4. Output (in *.out): 

FINAL SINGLE POINT ENERGY (MM)      -32.522715848742


FINAL SINGLE POINT ENERGY     -2474.898029955347


FINAL SINGLE POINT ENERGY (QM/MM)    -2507.420745804089


*** OPTIMIZATION RUN DONE ***



Protein free

!QMMM PBE0 6-31G* D3ZERO OPT


%PAL NPROCS 24 END


%qmmm


    ORCAFFFilename   "prmtop3.ORCAFF.prms"


    QMAtoms          { 0:33  


                      










 642:653


                      










 2238:2239


                      










 2248:2275


                      










 4716:4718 } end


    ActiveAtoms      { 0:33


                      










 642:653


                      










 2238:2239


                      










 2248:2275


                      










 4716:4718 } end


      ExtendActiveRegion        distance 


      Dist_AtomsAroundOpt       6.


      PrintOptRegionExt         true 


end


*pdbfile 0 1 pdb3.pdb 



If you start from a ComQum calculation, the ActiveAtom list can be 
taken from the belly in the sander.in3 file, but it must be extended 
also with all QMAtoms and also all junction atoms. 





NEB calculations

!QMMM PBE0 6-31G* D3ZERO NEB-TS SlowConv TightSCF pal16


%neb


   Product_PDBFile "orca-prod.pdb"


   TS_PDBFile "TSGuess.pdb" 


   NImages 8


   NEB_TS false 


end


%qmmm


    ORCAFFFilename   "prmtop3.ORCAFF.prms"


    QMAtoms          { 0:33  


                      










 642:653


                      










 2238:2239


                      










 2248:2275


                      










 4716:4718 } end


    ActiveAtoms      { 0:33  


                      










 642:653


                      










 2238:2239


                      










 2248:2275


                      










 4716:4718 } end


end





*pdbfile 0 1 orca-reac.pdb 



Output in orca.out

                   
 ***        THE NEB OPTIMIZATION HAS 
CONVERGED     ***


                        










 PATH 
SUMMARY             










 


  5     4.727   
-2507.39972     
20.33       0.02832   0.00417 
<= CI


             
 INFORMATION ABOUT SADDLE POINT     


Energy                                   










 ....  -2507.39972016 Eh



It seems that Orca requires coordinates of the junction atoms to be the same in all three pdb files (in variance to ComQum).








Useful greps




grep 'FINAL SINGLE' *.out



# Greps for all occurrences of convergence parameters in ORCA as they are displayed after each geom opt step


grep -A8 "Geometry convergence" 







grep "   MAX gradient" 

  

grep "ENERGY (QM/MM)" $1 | tail -n 1 ;

  





neb_energies() {


 
       # Gets the energies of the last step 
of an ORCA .NEB.log file and creates a temporary path vs energy 
file.


        # $1: file to grep for QM/MM energies


 
       tail -n 30 $1 | grep energy | tail -n 1 | sed 
's/energy          : //g'| tr -s ' ' '\n' > 
energies ;


 
       tail -n 30 $1 | grep distance $1 | tail -n 1 | sed
 "s/distance        :   //g" | tr -s ' ' 
'\n' > abscissa ;


        paste abscissa energies > tmp_path_energies ;


        rm -fr abscissa energies ;


}