Geometry manipulations and QM input building

The Following program reads an XYZ file (hardcoded). 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 5.X, 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 programs Turbomole and NWChem, as with the semiempirical programs xtb and Mopac.

package main

import (
	"fmt"
	"math"
	"os"

	chem "github.com/rmera/gochem"
	"github.com/rmera/gochem/qm"
	v3 "github.com/rmera/gochem/v3"
)

// This is a throw-away mini program, so the data is hardcoded.
func main() {
	mol, err := chem.XYZFileRead("../FILES/sample.xyz")
	if err != nil {
		panic(err.Error())
	}
	mol.SetCharge(1)
	mol.SetMulti(1)
	//Setting up the calculations
	//We will mosstly go with the defaults.
	calc := new(qm.Calc)
	calc.SCFTightness = 1 //rather demanding.
	calc.Method = "TPSS"
	calc.Dielectric = 80 //COSMO with epsilon=80. Just delete this line to avoid COSMO usage.
	calc.Basis = "def2-TZVP"
	calc.RI = true //RI approximation (also called charge-density fitting, in some programs).
	calc.Dispersion = "D3"
	orca := qm.NewOrcaHandle() //Change this line to use MOPAC2012 or Turbomole
	//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].VecView(6) //the 2 atoms defining the rotation axis
	axis2 := mol.Coords[0].VecView(7)
	torotate_indexes := []int{8, 9, 10, 11, 70, 83, 69}
	torotate := v3.Zeros(len(torotate_indexes))
	torotate.SomeVecs(mol.Coords[0], 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.Coords[0].SetVecs(rotated, torotate_indexes)
		orca.SetName(fmt.Sprintf("angle%1.1f", angle))
		//We first write the QM input and then an XYZ witht he non optimized geometry.
		oridir, err := os.Getwd()
		if err != nil {
			panic(err.Error())
		}
		os.Chdir("../FILES/results")
		if err := orca.BuildInput(mol.Coords[0], mol, calc); err != nil {
			panic(err.Error())
		}
		chem.XYZFileWrite(fmt.Sprintf("angle%1.1f.xyz", angle), mol.Coords[0], mol)
		os.Chdir(oridir)
	}
}

Disclaimer

goChem only provides a convenient way to produce inputs for the ORCA program. As with any other third-party program, obtaining the software, complying with its license terms, and citing the corresponding publications for the methods used, are exclusively your (the final user's) responsibilities.