Moment tensor analysis

The following example takes the name of an XYZ file and tries to obtain the masses for all atoms, and sets all masses to 1 if it fails (printing also a warning). It uses the coordinates and the masses to calculate the moment of inertia tensor. The program performs then principal component analysis on the tensor, and prints the widest eigenvector.

package main

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

func main() {
	mol, err := chem.XyzRead(os.Args[1])
	if err != nil {
		panic(err.Error())
	}
	mass, err := mol.MassCol()
	if err != nil {
		fmt.Println(err, "Will assign mass 1 to all atoms")
		mass = chem.OnesMass(mol.Len())
	}
	coords := mol.Coords[0]
	moment, err := chem.MomentTensor(coords, mass)
	if err != nil {
		panic(err.Error())
	}
	eigvectors, eigvalues, err := chem.Eigenwrap(moment)
	if err != nil {
		panic(err.Error())
	}
	main := eigvectors.GetRowVector(2) //The last is the widest.
	fmt.Printf("Widest eigenvector: %v  Corresponding eigenvalue: %4.1f\n", main, eigvalues[2])
}