TURBOMOLE

1. TURBOMOLE AdemTekin Yrd. Doc. Dr. Program Package for abinitio Electronic Structure Calculations 15.06.2012 Istanbul

2. Attention!! Input and lsf scripts can be downloaded from: http://www.be.itu.edu.tr/~adem.tekin/turbomole/filelist.txt http://www.be.itu.edu.tr/~adem.tekin/turbomole/water.xyz http://www.be.itu.edu.tr/~adem.tekin/molpro/filelist.txt http://www.be.itu.edu.tr/~adem.tekin/molpro/ex1.com Turbomole 15.06.2012

3. Activate turbomole and molpro in your shell • Go to your home directory • Open your bash shell script: vi .bashrc • Add the following lines for TURBOMOLE: export TURBODIR=/RS/progs/TURBOMOLE PATH=$PATH:$TURBODIR/scripts PATH=$PATH:$TURBODIR/bin/`sysname` export PARA_ARCH=MPI • Add the following lines for MOLPRO: export INTEL_LICENSE_FILE=/RS/progs/intel/licenses/RS/progs/intel/impi/3.1/bin64/mpivars.sh export molprodir=/RS/progs/molpro/mpp/2009/1-20/x86_64 PATH=$PATH:$molprodir/bin export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/RS/progs/molpro/mpp/2009/1-20/x86_64/lib export PARA_ARCH=MPI • Add the following lines for VESTA: export PATH=$PATH:/RS/progs/VESTA-x86_64 • Log out and log in again then check: export $TURBODIR export $molprodir Turbomole 15.06.2012

4. Outline • Introduction • How to create input (define tool) • Single point calculations at the Hartree-Fock, DFT and RIDFT levels • Geometry optimization of water dimer at the RIMP2 level • Geometry optimization of water dimer at the SCS-MP2 level • Geometry optimization of water dimer at the CP-SCS-MP2 level Turbomole 15.06.2012

5. Turbomole • Developed by Reinhart Ahlrichs (University of Karlsruhe) 1989-2007 • Since 2007 TURBOMOLE GmbH founded by Ahlrichs is responsible • for more information please visit to: • http://www.turbomole.com/ • http://www.cosmologic.de • A forum site is available under http://www.turbo-forum.com Turbomole 15.06.2012

6. Turbomoleusage philosophy The usage of TURBOMOLE has been adapted to the way a UNIX user is working: • Command line driven • Many different programs, each one specialized for methods and/or properties • Scripts are used to combine the functionalities and manage workflows • Input can be changed by text editors • Output processed by standard UNIX tools (editors, grep, awk, ...) Turbomole 15.06.2012

7. Turbomole usage philosophy Turbomole 15.06.2012

8. Turbomole modules Turbomole 15.06.2012

9. Turbomole modules DEFINE interactive input generator which creates the input file control. UFF performs a geometry optimization at a force field level. DSCF for semi-direct SCF (self-consistent field) and DFT calculations. GRAD require a successful DSCF run and calculates the gradient of the energy with respect to nuclear coordinates. RIDFT and perform DFT calculations – as DSCF and DGRAD – within the RI-J approximation. RDGRAD MPGRAD requires a well converged SCF run and performs closed-shell RHF or UHF calculations yielding single point MP2 energies and, if desired, the corresponding gradient. RIMP2 calculates MP2 energies and gradients for RHF and UHF wavefunctions significantly more efficient then MPGRAD. RICC2 calculates electronic excitation energies, transition moments and properties of excited states at the CIS, CIS(D), ADC(2) and CC2 level using a closed shell RHF or a UHF SCF reference function. Turbomole 15.06.2012

10. Turbomole modules RELAX requires a gradient run – by GRAD, RDGRAD, RIMP2 or MPGRAD and proposes a new structure based on the gradient and the approximated force constants. STATPT performs a structure optimization using the “Trust Radius Image Minimization” algorithm. It can be used to find minima or transition structure. FROG executes one molecular dynamics step. AOFORCE require a well converged SCF or DFT run and performs an analytic calculation of force constants, vibrational frequencies and IR intensities. ESCF requires a well converged SCF or DFT run and calculates time dependent and dielectric properties. EGRAD computes gradients and first-order properties of excited states. MPSHIFT computes NMR chemical shieldingns for all atoms of the molecules at the SCF, DFT or MP2 level. FREEH calculates thermodynamic functions from molecular data in a control file: an AOFORCE or a NUMFORCE run is a necessary prerequisite. Turbomole 15.06.2012

11. How to create the input • Build definehas some limited features to build a structure from scratch, but that is neither convenient nor intuitive, except for some small cases. → Use an external builder (molden, ECCE, Hyperchem, Insight, Arguslab, MAPS, Accelrys DS Visualizer, ...there are many!) and convert the structure to xyz format. • Convert Use the TURBOMOLE script x2tto convert xyz files to TURBOMOLE coordinates: x2tstruct.xyz>coord Turbomole 15.06.2012

12. Define defineis an interactive input generator which creates the input file control: • Supports most basis sets in use, especially the only fully atom optimized consistent basis sets of SVP, TZV, and QZV quality available for the atoms H–Rn (including lanthanides and actinides up to Lawrencium for SVP and TZV) • Determines the molecular symmetry • Determines the internal coordinates, allowing efficient geometry optimization • Allows to perform a geometry optimization at a force field level to pre-optimize the geometry and to calculate a Cartesian Hessian matrix • Sets the keywords necessary for single point calculations and geometry optimizations for a variety of methods • Manipulate geometries of molecules Turbomole 15.06.2012

13. Define There are four major menus in define: • Molecular Geometry Menu • Atomic Attribute Definition Menu • Occupation number & Molecular Orbital Definition Menu • General Menu Navigation within define: 1. To finish the menu and proceed to the next one, enter: * 2. To go back to the previous menu, enter: & 3. If you are in one of the four main menus, hitting <Enter> will just print out the menu you are in 4. To quit define immediately (and prematurely), enter: qq Turbomole 15.06.2012

14. Define To call define just type define in the shell : $define Turbomole • First you will be asked whether you want to define a new system or not • If not press the ENTER button 15.06.2012

15. Define Then you will be asked for a n input title: Turbomole • either give a name • or press the ENTER button to proceed 15.06.2012

16. Define – molecular geometry menu Turbomole 15.06.2012

17. Define – atomic attribute definition menu Suppose we created the coord file via: x2tstruct.xyz>coord acoordto load the geometry desyto detect the symmetry, if you expect a higher symmetry, repeat with increased tolerance desy 0.1 ired to use internal redundant coordinates Then pressing * will lead us to the Atomic Attribute Definition Menu Select only these for a geometry optimization based on cartesian coordinates Turbomole 15.06.2012

18. Define – atomic attribute definition menu The default basis set is def-SV(P) You can assign the same basis set for all atoms in the system by: ballaug-cc-pVTZ (Dunning’s basis set) If you want to use a different basis set for any atomic species (let’s say for hydrogen): b“h”aug-cc-pVDZ bl use to see the current basis set assignments TURBOMOLE has a new set of basis sets abbreviated by def2: def2-SV(P), def2-SVP, def2-TZVP, def2-TZVPP, def2-QZVP, def2-QZVPP bi information about improved TURBOMOLE basis sets (At least you must use def2-TZVP for a quantative analysis) Turbomole 15.06.2012

19. Define – atomic attribute definition menu SV(P) or def-SV(P) for routine SCF or DFT. Quality is about 6-31G*. TZVP or def-TZVP for accurate SCF or DFT. Quality is slightly better than 6-311G**. TZVPP or def-TZVPP for MP2 or close to basis set limit SCF or DFT. Comparable to 6-311G(2df). QZVP and QZVPP for highly correlated treatments; quadruple zeta + 3d2f1g or 4d2f1g (beyond Ne), 3p2d1f for H. Turbomole 15.06.2012

20. Define – occupation number & orbital definition menu This menu serves for the definition of occupation numbers (and thereby the charge of the molecule) and the start orbitals for the SCF calculations. For well-behaved molecules, The procedure is very easy. Just choose the option: eht to perform an extended Huckel calculation. Turbomole 15.06.2012

21. Define – occupation number & orbital definition menu Turbomole • If you do not want them type n • You can accept the default parameters by typing y, then you will be asked for the charge of the system (default is 0) 15.06.2012

22. Define – occupation number & orbital definition menu Turbomole • to accept the current occupation type y, then you will enter to the GENERAL MENU 15.06.2012

23. Define – general menu Here, you can choose additional parameters for your calculation. Turbomole 15.06.2012

24. Define – general menu - scf Turbomole 15.06.2012

25. Define – general menu – mp2 Turbomole • Do not forget to freeze (frozen core approximation) • Do not forget to add cbas basis: the basis set selected in atomic attribute definition menu • You can change memory (default is 60 MB) • You can change convergence, use 0.1E-07 15.06.2012

26. Define – general menu – cc2 • CC2 is an approximated coupled cluster singles and doubles (CCSD) model • The CC2 total energy is of second-order Moller-Plesset perturbation theory (MP2) quality • The scaling of CC2 is N5 whereas the scaling of CCSD is N6, where N is the number of orbitals Turbomole • Similar to MP2 submenu • Use always ri option • Spin component scaled (SCS) MP2 is hidden in the ricc2 submenu 15.06.2012 For more info about CC2 check Chem. Phys. Lett. 243 409 (1995)

27. Define – general menu – cc2 – ricc2 Turbomole 15.06.2012

28. Define – general menu – dft • switch on DFT • select functional • select gridsize Turbomole Type on to switch on DFT (default functional is b-p) (default gridsize is m3) 15.06.2012

29. Define – general menu – dft Turbomole 15.06.2012

30. S22 test of interaction energies [JCP 132 (2010) 144104] Mean Absolute Deviation Turbomole 15.06.2012

31. Empirical Dispersion Correction for DFT Calculations • DFT functionals fails for intermolecular interaction energy calculations (They do not result any binding between two monomers) • This is due to insufficient treatment of electron correlation (which is so important to describe dispersion forces) • Stefan Grimme suggested dispersion correction to ordinary DFT energies • When using the dispersion correction, the total energy is given by Turbomole • Edisp is an empirical dispersion correction given by Nat is the number of atoms C6ijdenotes the dispersion coefficient for atom pair ij s6 global scaling factor depends on the selected DF Rij is the interatomic distance between i and j 15.06.2012

32. Define – general menu – dft Turbomole Turbomole offers grids 1 (coarse) to 7 (nest), and the multiple grids m3 to m5. The latter employ a coarser grid during SCF iterations, and grid 3 to grid 5 in the final SCF iteration and the gradient evaluation. Default is grid m3, for clusters with more than 50 atoms use m4. 15.06.2012

33. Define – general menu – RI-J (resolution of identity) • In the RI approximation two-electron integrals are approximated by three-center expansions • RI procedure is recommended for non-hybrid functionals (which speeds up calculations by a factor of 10 at least) Turbomole • activate ri with on • assign jbas • then call jobex script with –ri option (jobex -ri) 15.06.2012

34. Define – general menu – MARI-J • RI-J calculations can be done even more effciently with the Multipole Accelerated RI-J (MARI-J ) option, especially for larger molecules. • Just rely on the defaults. Turbomole 15.06.2012

35. Define – general menu – stp – transition state search • Stationary points are places on the potential energy surface (PES) with a zero gradient. • Two types of stationary points are of special importance to chemists. These are minima (reactants, products, intermediates) and first-order saddle points (transition states) • The type of stationary point optimization depends on the value of irtvec. By default irtvec is set to 0, which implies a structure minimization. A value irtvec > 0 implies a transition state optimization using the eigenvalue-following TRIM (Trust Radius Image Minimization) algorithm . • Transition state search is invoked by jobex –trans Turbomole 15.06.2012

36. Vibrational normal modes • After having performed a single point energy calculation or a geometry optimization, just start the program aoforce(analytic force constant calculation): aoforce > force.out • All details will be written to force.out • The keyword $vibrational spectrum in the controlfile contains a list of the modes: • NumForceenables numerical force constant calculations for all levels of theory with a gradient implemented. To run NumForce just issue the command NumForce –ri –level mp2 -central Turbomole 15.06.2012

37. Structure optimizations using the jobex script • The shell script Jobex controls and executes automatic optimizations of molecular geometry parameters. • It will cycle through the direct SCF, gradient and force relaxation programs and stop if either the maximum number of cycles is reached or the convergence criteria are fulfilled. • Jobex has many options: -energy integerconverge total energy up to 10-integerHartree (default:6) -gcartintegerconverge maximum norm of cartesian gradient up to 10-integera.u. (def:3) -c integer perform up to integer cycles (default:20) -dscfbegin with a direct SCF step -grad begin with a gradient step -statptbegin with a force relaxation step -relax use the RELAX program for force relaxation -trans perform transition state search -level leveldefine the optimization level, level=scf, mp2, cc2 or uff -riuse RI modules -ex perform excited state geometry optimization using EGRAD -mda molecular dynmaics run Turbomole 15.06.2012

38. Turbomole output • There will be an output written to file job.start. • The convergence is signalled by the file converged; otherwise, you should find the file not.converged within your working directory. • If Jobex finds a file named stop or STOP in the working directory, Jobex will stop after the present step has terminated. You can create stop by the command touch stop. • The output of the last complete cycle is written to file job.last, while the output of the running cycle is collected within the file job.<cycle>, where <cycle> is the index of the cycle. • Look at the file gradient: it contains all geometries and gradients at each step of optimization: grep cycle gradient • To see the final geometry call t2x –c > final.xyzt2x > final_all.xyz (to see all trajectory) • At the command line, use dist, bend or tors to get bond lengths, angles, dihedral angles Turbomole 15.06.2012

39. Keywords in the control file • The file control is the input file for TURBOMOLE. • Keywords start with a ‘$’, e.g. $title. $title NO2 c2v UKS SVP $operating system unix $symmetry c2v $coord file=coord $intdef file=coord $atoms n 1 \ basis =n def-SVP o 2-3 \ basis =o def-SVP … $energy file=energy $grad file=grad $forceapprox file=force $lock off $dft functional b-p gridsize m3 $last step define $end Turbomole 15.06.2012

40. Submitting Turbomole jobs on UHEM • single job submission #!/bin/bash #BSUB -m karadeniz_temp1 # ege or anadolu #BSUB -a intelmpi # bukismidegistirmeyin !!! #BSUB -J triazine # isinizitanimlayacakbirisim #BSUB -o triazine.out # bukismidegistirmeyin !!! #BSUB -e triazine.err # bukismidegistirmeyin !!! #BSUB -q workshop# kuyrukismi #BSUB -P workshop # projeismi #BSUB -R "span[ptile=1]" # bukismidegistirmeyin !!! #BSUB -n 1 # kullanilacakolanislemcisayisi echo 'Starting time:' date jobex -ri -level cc2 -energy 7 -gcart 4 -c 1000 echo 'Ending time:' date Turbomole 15.06.2012

41. Submitting Turbomole jobs on UHEM • parallel job submission #!/bin/bash #BSUB -m karadeniz_temp1 #BSUB -m "b074 b075" # submit the job to the b074 and b075 nodes #BSUB -a intelmpi # bukismidegistirmeyin !!! #BSUB -J triazine # isinizitanimlayacakbirisim #BSUB -o triazine.out # bukismidegistirmeyin !!! #BSUB -e triazine.err # bukismidegistirmeyin !!! #BSUB -q workshop# kuyrukismi • #BSUB -P workshop # projeismi #BSUB -R "span[ptile=8]" # bukismidegistirmeyin !!! #BSUB -n 16 # kullanilacakolanislemcisayisi exportHOSTS_FILE=/AKDENIZ/HOME005/users/tekin/hostsfile exportPARNODES=16 echo 'Starting time:' date jobex -ri -level cc2 -energy 7 -gcart 4 -c 1000 echo 'Ending time:' date check availability of nodes with bhostscommand b074 b074 b074 b074 b074 b074 b074 b074 b075 b075 b075 b075 b075 b075 b075 b075 Turbomole 15.06.2012

42. Turbomole examples Turbomole • Single point calculations at the Hartree-Fock, DFT and RIDFT levels • Geometry optimization at the RIMP2 level • Geometry optimization at the SCS-MP2 level • Geometry optimization at the CP-SCS-MP2 level 15.06.2012

43. Single point calculations – HartreeFock, DFT and RIDFT energy calculation of benzene • 1. Create an empty directory, e.g., hf_def2_tzvp_benzene • 2. Call define • 3. Hit <enter> • 4. Give a title: benzene, HF calculation 5. Load benzene from the structure library ($TURBODIR/structures) a !benzene 6. define will scan through the structure library and ask you if you want to take the proposed molecule Turbomole 15.06.2012

44. Single point calculations – HartreeFock, DFT and RIDFT energy calculation of benzene 7. Just say y to accept the structure 8. Switch on the symmetry detection be entering: desy 9. Next step is to provide internal coordinates. (That is not needed for single point calculations.) ired 10. Now we are done with this menu and proceed to the next one by entering * 11. This is the atomic attributes menu. Choose your basis set: b all def2-TZVP 12. Now we are done with this menu and proceed to the next one by entering * 13. Determine the starting molecular orbitals with an extended Huckel guess: eht 14. Then you will be asked a lot of questions. Just accept the defaults. 15. Now you are in general menu. 16. The default method is Hartree-Fock (HF) calculation, and therefore we do not have to change anything here. Simply enter: * 17. To run a HF single point energy calculation, call dscf > dscf.out Turbomole 15.06.2012

45. Single point calculations – HartreeFock, DFT and RIDFT energy calculation of benzene 18. After the run, look in the output (dscf.out) and/or the energy file. The energy file should give: $energy SCF SCFKIN SCFPOT 1 -230.7848712224 230.6217054044 -461.4065766268 $end 19. Look in the output (dscf.out) for the dipole and quadrupole moments and whatever you like to know. 20. Get more details about the orbitals, the HOMO-LUMO gap and occupation: eiger 21. Let’s run a DFT calculation just after the HF calculation . 22. Call define and just accept all the defaults, then you will be in general menu. 23. Choose dft 24. Choose on to switch on DFT (Now you have chosen the default functional (b-p) and gridsize m3) 25. Quit from define with * 26. Run dscfagain. 27. Look again to the energy file: $energy SCF SCFKIN SCFPOT 1 -230.7848712224 230.6217054044 -461.4065766268 2 -232.3365711535 231.4048783475 -463.7414495010 $end Turbomole 15.06.2012

46. Single point calculations – HartreeFock, DFT and RIDFT energy calculation of benzene Turbomole 28. Now, let’s switch on RI approximation following previous steps enter to the general menu 29. enter riand then on. 30. Run ridft program: ridft > ridft.out $energy SCF SCFKIN SCFPOT 1 -230.7848712224 230.6217054044 -461.4065766268 2 -232.3365711535 231.4048783475 -463.7414495010 3 -232.3368258829 231.4106397547 -463.7474656376 $end 15.06.2012

47. Do not be confused with auxiliary basis 1. Coulomb-Fitting ($jbas). Basis sets declared in this data-group are needed for RI-DFT calculations, where only the Coulomb part of the Fock operator is approximated. They can be kept rather small. 2. Coulomb- and Exchange-Fitting ($jkbas). If additionally the exchange part of the Fock operator needs to be approximated (RI-SCF calculations), a larger fitting basis is needed. It also can be used for RI-DFT runs (if the accuracy of the $jbas auxiliary basis is not satisfactory for you). 3. Correlation-Fitting ($cbas). These basis sets are needed for RI-MP2 and RI-CC2 calculations. They contain basis functions with high angular momentum. Turbomole 15.06.2012

48. Geometry optimizations – RIMP2, SCS-MP2 and CP-SCS-MP2 calculations of water dimer 1. Open a new directory, e.g., rimp2_avdz_water 2. Create water.xyz file. Here are the coordinates: 6 O 0.0640488861 -1.5542001572 0.0000000000 H -0.0275603110 -0.6109761929 0.0000000000 H -0.8291328873 -1.8994223998 0.0000000000 O -0.0539423286 1.4751039841 0.0000000000 H 0.3402097464 1.8881478835 0.7547208224 H 0.3402097464 1.8881478835 -0.7547208224 3. Convert xyz file to turbomole coord file: x2t water.xyz > coord 4. Call define 5. Hit <Enter> 6. Give a title 7. Read the coordinates: a coord 8. desy 9. ired Turbomole 15.06.2012

49. Geometry optimizations – RIMP2, SCS-MP2 and CP-SCS-MP2 calculations of water dimer 10. Type * to go to the next menu 11. Choose basis set: b all aug-cc-pVDZ 12. Type * to go to the next menu 13. Activate eht and accept the defaults 14. Now you are in the general menu 15. First select scf 16. Change scf option conv and then type 8 and go back to the general menu by pressing to ENTER 17. Now, select mp2 18. Activate freeze and go back to the mp2 submenu by pressing to * 19. Activate cbas: Automatically the auxiliary basis set will be loaded. 20. Press * to go back to the mp2 submenu. 21. Select the denconv and enter 0.1E-07. 22. Press * to go back to the mp2 submenu. 23. Press * to quit define. 24. Run the jobex script with the following parameters: jobex -ri -level mp2 -energy 7 -gcart 4 -c 1000 25. Check the convergence by grep cycle gradient 26. After the convergence convert the turbomole geometry to xyz: t2x –c > rimp2_opt.xyzor t2x > final_all_trajectory.xyz

50. Geometry optimizations – RIMP2, SCS-MP2 and CP-SCS-MP2 calculations of water dimer 1. Open a new folder, e.g., scsmp2_avdz_water 2. Copy the rimp2 xyz output here: cp ../rimp2_avdz_water/rimp2_opt.xyz . 3. Convert the xyz file to turbomole file: x2t rimp2_opt.xyz > coord 4. Call define 5. a coord 6. desy 7. ired 8. * 9. b all aug-cc-pVDZ 10. * 11. Activate eht and accept the defaults 12. Now you are in the general menu 13. First select scf 14. Change scf option conv and then type 8 and go back to the general menu by pressing to ENTER 15. Now, select cc2 16. Activate freeze and go back to the cc2 submenu by pressing to * 17. Activate cbas: Automatically the auxiliary basis set will be loaded. 18. Press * to go back to the cc2 submenu. 19. Select the denconv and enter 0.1E-07.