Geometry manipulations and QM input building

The Following program reads an XYZ file, which name has been passed as an argument. It rotates a bunch of atoms about an axis defined by two more atoms. For simplicity, the rotated atoms and the ones defining the axis have been hardcoded in the program, so they have to be modified in the code for the system of interest.

The program rotates the atoms by several angles (also hardcoded for simplicity), writing a new XYZ and preparing an ORCA single-point energy calculation for each angle of rotation. Note that in order to run the energy calculations you need the ORCA software version 2.9, which you can obtain from their developers here.

Note that the calculation settings are defined once and used as many times as needed. In setting the inputs, it is possible to go with the defaults to a rather large degree. In this case we have made some changes, and made the density functional explicit (the default is a revPBE/def2-SVP single-point calculation using the RI approximation and Grimme et al. D3 dispersion correction).

The reference for the methods and the basis sets used in this example can be found in the ORCA manual.

Note: Gochem can now also interact with the ab initio program Turbomole and the semiempirical program Mopac.

package main

import (
	"fmt"
	"github.com/rmera/gochem"
	"math"
	"os"
)

func main() {
	mol, err := chem.XyzRead(os.Args[1])
	if err != nil {
		panic(err.Error())
	}
	mol.SetCharge(1)
	mol.SetUnpaired(1) //I just hardcoded these.
	//Setting up the calculations
	//We will mosstly go with the defaults.
	calc := new(chem.QMCalc)
	calc.SCFTightness = 2 //rather demanding.
	calc.Method = "revPBE"
	calc.Dielectric = 4
	calc.Basis = "def2-TZVPP"
	calc.RI = true //RI approximation.
	calc.Disperssion = "D3" // D3 dispersion correction (link in the explanation above).
	orca := chem.MakeOrcaRunner()
	//Now we play with a bond and make orca inputs
	// to calculate a SP energy for each geometry.
	//*******************************************
	//I just hard-coded these ones, they make sense
	//for my test system but you will have to change them for yours.
	axis1 := mol.Coords[0].GetRowVector(32) //the 2 atoms defining the rotation axis
	axis2 := mol.Coords[0].GetRowVector(36)
	torotate_indexes := []int{31, 33, 34, 35, 63, 64, 65}
	torotate, _ := mol.SomeCoords(torotate_indexes)                     //The atoms that will rotate
	angles := []float64{0, 0.2 * math.Pi, 0.4 * math.Pi, 0.5 * math.Pi} //The rotation angles in radians.
	for _, angle := range angles {
		rotated, err := chem.RotateAbout(torotate, axis1, axis2, angle)
		if err != nil {
			panic(err.Error())
		}
		//now we put the rotated coords in the molecule
		mol.SetCoords(rotated, torotate_indexes, 0)
		orca.SetName(fmt.Sprintf("angle%1.1f", angle))
		if err := orca.BuildInput(mol, mol.Coords[0], calc); err != nil {
			panic(err.Error())
		}
		chem.XyzWrite(mol, mol.Coords[0], fmt.Sprintf("angle%1.1f.xyz", angle))
	}
}

##Disclaimer Gochem only provides a convenient way to produce input for the ORCA program. Obtaining the software, complying with its license terms (as with any other third party software's license terms), and citing the corresponding publications for the methods used, are exclusively your (the final user's) responsibilities.