Simple Water Calculation

Geometry Optimization

Lets first make a new directory called water_opt using mkdir water_opt. Go into the directory with cd water_opt. This is the directory that we will setup and run our geometry optimization calculation for the simple water molecule.

Setting Up the Calculation

Geometry

We need a coordinate file that tells TURBOMOLE what atoms are included in our calculation and some starting positions. For water, we have 1 oxygen and 2 hydrogens. We will create a simple .xyz file that holds our coordinates. We will guess the coordinates very simply (though a better guess will result in a faster optimization). Further, programs such as molden, chemcraft, avogadro, chimera, etc, are capable to generating these .xyz files, hence it is not likely that you will be editing these files by hand for large systems. Regardless, save the following as coord.xyz

3
Initial Water Geometry
O 0 0 0 
H 1 0 0
H 0 1 0

The first line indicates the number of atoms in the system (3 in our case). The second line is a comment line. You can place anything that you want here, but there needs to be a line here! Then our coordinates are in the form element x y z where the x,y,z are in angstroms (note: this differs from TURBOMOLE which uses bohrs (atomic units) by default). If you view this file in your favorite chemical model viewer (molden for example) you should a 90° angle between the 2 hydrogens, whereas we would expect the angle to be about 104.5° (based on the oxygen being tetrahedral with 2 lone pairs). Hence, we should wish the optimizer should result in a angle of about 104.5° in the converged geometry.

We can now convert the coord.xyz file to a coord file (which is the native version TURBOMOLE uses) using TURBOMOLE's x2t command (xyz 2 turbomole). If you type x2t coord.xyz you should see

$coord
    0.00000000000000      0.00000000000000      0.00000000000000       o
    1.88972613288564      0.00000000000000      0.00000000000000       h
    0.00000000000000      1.88972613288564      0.00000000000000       h
$end

it is similar to the xyz file, except the elements are in the last column, and we are now in Bohr's rather than angstroms. Also note, that the output was to stdout. To capture the output into the coord file, we use the redirection symbol > on the x2t command to push the content into the file: x2t coord.xyz > coord. You should see a new file called coord with the above content. At this point, all we need to do is setup the parameters for our calculation.

definput.json

Create a new file called definput.json with the following contents

{
	"geometry" : {
		"cartesians" : true,
		"idef" : {
			"idef_on" : false,
			"freeze_stretch" : [
				"4,54",
				"3,24"
			],
            "freeze_dihedral": [
                "1,2,3,4"
            ]
		},
		"ired" : false,
		"iaut" : {
			"iaut_on" : false,
			"bonds" : [
				"1,5",
				"10,8"
			]
		}
	},

	"basis" : {
		"all" : "def2-SVP"
	},

	"charge" : 0,

	"open_shell" : {
		"open_shell_on" : false,
		"unpaired" : 1
	},

	"dft" : {
		"dft_on" : true,
		"func" : "b3-lyp",
		"grid" : "m4"
	},

	"rij" : true,
	"marij" : true,
	"dsp" : true,
	"stp" : {
		"itvc" : 0,
		"trad": 0.1
	},

	"scf" : {
		"iter" : 300,
		"conv" : 6
	},

	"cosmo" : null,

	"freeze_atoms" : [],
	"calculation" : "geo",
	"geo_iterations" : 200,
	"weight" : true,
	"gcart" : null,
	"denconv" : null
}

Note that charge is set to 0 (since water is neutral), our functional is b3-lyp, basis set is def2-SVP, and the calculation is set to "geo". open_shell_on is false since this is a closed shell system (no unpaired electrons) and cosmo is null since we are performing the calculation in the gas phase. Calling setupturbomole.py now will result in the following being printed out

	[Define]          ==>> Ended Normally
	[setup turbomole] ==>> SUCCESS

At this point, a bunch of new files should be present (basis, auxbasis, control, etc). The control file is the actual TURBOMOLE file that contains all of the settings for your calculations and direct manipulation of this file will result in changes to your settings. Certain properties, such as excited state calculations, or writing out data is done by directly editing the file. Though, for our purposes it should not be necessary. At this point we are readying to run the calculation!

Running the Calculation

We can easily run the calculation with runturbomole.py -n 1. This will run the calculation with 1 processor in the current directory. The job should finish fairly quickly. If you inspect the directory, you should see

GEO_OPT_CONVERGED  basis      coord       definput.json  hessapprox  jobex.out
MFILE              control    coord.xyz   energy         job.last    mos
auxbasis           converged  define.out  gradient       job.start   statistics

The GEO_OPT_CONVERGED file lets us know the geometry is converged. The energy file contains the energy at each step of the geometry optimization, and the gradient contains the coordinates and forces on each atom at each step of the optimization. If we run t2x > traj.xyz will can convert this data back into a .xyz file which can then be opened in molden. The coord file also contains the optimized geometry:

$coord
    0.09420457462053      0.09420457403034     -0.00000000000000  o
    1.90942595850592     -0.11390440008893      0.00000000000000  h
   -0.11390440024085      1.90942595894425      0.00000000000000  h
$user-defined bonds
$end