func TestWater(Te *testing.T) {
	//	runtime.GOMAXPROCS(2) ///////////////////////////
	mol, err := XYZFileRead("test/sample.xyz")
	if err != nil {
		Te.Error(err)
	}
	for i := 0; i < 6; i++ {
		s := new(Atom)
		if i == 0 || i == 3 {
			s.Symbol = "O"
		} else {
			s.Symbol = "H"
		}
		mol.AppendAtom(s)
	}
	mol.SetCharge(1)
	mol.SetMulti(1)
	c2 := v3.Zeros(mol.Len())
	v := v3.Zeros(6)
	l, _ := mol.Coords[0].Dims()
	fmt.Println(l, mol.Len())
	c2.Stack(mol.Coords[0], v)
	mol.Coords[0] = c2
	c := mol.Coords[0].VecView(43)
	h1 := mol.Coords[0].VecView(42)
	coords := v3.Zeros(mol.Len())
	coords.Copy(mol.Coords[0])
	w1 := MakeWater(c, h1, 2, Deg2Rad*30, true)
	w2 := MakeWater(c, h1, 2, Deg2Rad*-30, false)
	tmp := v3.Zeros(6)
	tmp.Stack(w1, w2)
	fmt.Println("tmp water", w1, w2, tmp, c, h1)
	coords.SetMatrix(mol.Len()-6, 0, tmp)
	XYZFileWrite("test/WithWater.xyz", coords, mol)
}

//EulerRotateAbout uses Euler angles to rotate the coordinates in coordsorig around by angle
//radians around the axis given by the vector axis. It returns the rotated coordsorig,
//since the original is not affected. It seems more clunky than the RotateAbout, which uses Clifford algebra.
//I leave it for benchmark, mostly, and might remove it later. There is no test for this function!
func EulerRotateAbout(coordsorig, ax1, ax2 *v3.Matrix, angle float64) (*v3.Matrix, error) {
	r, _ := coordsorig.Dims()
	coords := v3.Zeros(r)
	translation := v3.Zeros(ax1.NVecs())
	translation.Copy(ax1)
	axis := v3.Zeros(ax2.NVecs())
	axis.Sub(ax2, ax1) //now it became the rotation axis
	f := func() { coords.SubVec(coordsorig, translation) }
	if err := gnMaybe(gnPanicker(f)); err != nil {
		return nil, CError{err.Error(), []string{"v3.Matrix.Subvec", "EulerRotateAbout"}}

	}
	Zswitch := RotatorToNewZ(axis)
	coords.Mul(coords, Zswitch) //rotated
	Zrot, err := RotatorAroundZ(angle)
	if err != nil {
		return nil, errDecorate(err, "EulerRotateAbout")
	}
	//	Zsr, _ := Zswitch.Dims()
	//	RevZ := v3.Zeros(Zsr)
	RevZ, err := gnInverse(Zswitch)
	if err != nil {
		return nil, errDecorate(err, "EulerRotateAbout")
	}
	coords.Mul(coords, Zrot) //rotated
	coords.Mul(coords, RevZ)
	coords.AddVec(coords, translation)
	return coords, nil
}

//Super determines the best rotation and translations to superimpose the coords in test
//listed in testlst on te atoms of molecule templa, frame frametempla, listed in templalst.
//It applies those rotation and translations to the whole frame frametest of molecule test, in palce.
//testlst and templalst must have the same number of elements. If any of the two slices, or both, are
//nil or have a zero lenght, they will be replaced by a list containing the number of all atoms in the
//corresponding molecule.
func Super(test, templa *v3.Matrix, testlst, templalst []int) (*v3.Matrix, error) {
	//_, testcols := test.Dims()
	//_, templacols := templa.Dims()
	structs := []*v3.Matrix{test, templa}
	lists := [][]int{testlst, templalst}
	//In case one or both lists are nil or have lenght zero.
	for k, v := range lists {
		if v == nil || len(v) == 0 {
			lists[k] = make([]int, structs[k].NVecs(), structs[k].NVecs())
			for k2, _ := range lists[k] {
				lists[k][k2] = k2
			}
		}
	}
	//fmt.Println(lists[0])
	if len(lists[0]) != len(lists[1]) {
		return nil, CError{fmt.Sprintf("Mismatched template and test atom numbers: %d, %d", len(lists[1]), len(lists[0])), []string{"Super"}}
	}
	ctest := v3.Zeros(len(lists[0]))
	ctest.SomeVecs(test, lists[0])
	ctempla := v3.Zeros(len(lists[1]))
	ctempla.SomeVecs(templa, lists[1])
	_, rotation, trans1, trans2, err1 := RotatorTranslatorToSuper(ctest, ctempla)
	if err1 != nil {
		return nil, errDecorate(err1, "Super")
	}
	test.AddVec(test, trans1)
	//	fmt.Println("test1",test, rotation) /////////////77
	test.Mul(test, rotation)
	//	fmt.Println("test2",test) ///////////
	test.AddVec(test, trans2)
	//	fmt.Println("test3",test) ///////
	return test, nil
}

//MassCenter centers in in the center of mass of oref. Mass must be
//A column vector. Returns the centered matrix and the displacement matrix.
func MassCenter(in, oref *v3.Matrix, mass *mat64.Dense) (*v3.Matrix, *v3.Matrix, error) {
	or, _ := oref.Dims()
	ir, _ := in.Dims()
	if mass == nil { //just obtain the geometric center
		tmp := ones(or)
		mass = mat64.NewDense(or, 1, tmp) //gnOnes(or, 1)
	}
	ref := v3.Zeros(or)
	ref.Copy(oref)
	gnOnesvector := gnOnes(1, or)
	f := func() { ref.ScaleByCol(ref, mass) }
	if err := gnMaybe(gnPanicker(f)); err != nil {
		return nil, nil, CError{err.Error(), []string{"v3.Matrix.ScaleByCol", "MassCenter"}}
	}
	ref2 := v3.Zeros(1)
	g := func() { ref2.Mul(gnOnesvector, ref) }
	if err := gnMaybe(gnPanicker(g)); err != nil {
		return nil, nil, CError{err.Error(), []string{"v3.gOnesVector", "MassCenter"}}
	}
	ref2.Scale(1.0/mass.Sum(), ref2)
	returned := v3.Zeros(ir)
	returned.Copy(in)
	returned.SubVec(returned, ref2)
	/*	for i := 0; i < ir; i++ {
			if err := returned.GetRowVector(i).Subtract(ref2); err != nil {
				return nil, nil, err
			}
		}
	*/
	return returned, ref2, nil
}

//Dihedral calculates the dihedral between the points a, b, c, d, where the first plane
//is defined by abc and the second by bcd.
func Dihedral(a, b, c, d *v3.Matrix) float64 {
	all := []*v3.Matrix{a, b, c, d}
	for number, point := range all {
		pr, pc := point.Dims()
		if point == nil {
			panic(PanicMsg(fmt.Sprintf("goChem-Dihedral: Vector %d is nil", number)))
		}
		if pr != 1 || pc != 3 {
			panic(PanicMsg(fmt.Sprintf("goChem-Dihedral: Vector %d has invalid shape", number)))
		}
	}
	//bma=b minus a
	bma := v3.Zeros(1)
	cmb := v3.Zeros(1)
	dmc := v3.Zeros(1)
	bmascaled := v3.Zeros(1)
	bma.Sub(b, a)
	cmb.Sub(c, b)
	dmc.Sub(d, c)
	bmascaled.Scale(cmb.Norm(0), bma)
	first := bmascaled.Dot(cross(cmb, dmc))
	v1 := cross(bma, cmb)
	v2 := cross(cmb, dmc)
	second := v1.Dot(v2)
	dihedral := math.Atan2(first, second)
	return dihedral
}

//This program will align the best plane passing through a set of atoms in a molecule with the XY-plane.
//Usage:
func main() {
	if len(os.Args) < 2 {
		fmt.Printf("Usage:\n%s file.xyz [indexes.dat]\nindexes.dat is a file containing one single line, with all the atoms defining the plane separated by spaces. If it is not given, all the atoms of the molecule will be taken to define the plane.\n", os.Args[0])
		os.Exit(1)
	}
	mol, err := chem.XYZFileRead(os.Args[1])
	if err != nil {
		panic(err.Error())
	}
	var indexes []int
	//if no file with indexes given, will just use all the atoms.
	if len(os.Args) < 3 {
		indexes = make([]int, mol.Len())
		for k, _ := range indexes {
			indexes[k] = k
		}
	} else {
		indexes, err = scu.IndexFileParse(os.Args[2])
		if err != nil {
			panic(err.Error())
		}
	}
	some := v3.Zeros(len(indexes)) //will contain the atoms selected to define the plane.
	some.SomeVecs(mol.Coords[0], indexes)
	//for most rotation things it is good to have the molecule centered on its mean.
	mol.Coords[0], _, _ = chem.MassCenter(mol.Coords[0], some, nil)
	//As we changed the atomic positions, must extract the plane-defining atoms again.
	some.SomeVecs(mol.Coords[0], indexes)
	//The strategy is: Take the normal to the plane of the molecule (os molecular subset), and rotate it until it matches the Z-axis
	//This will mean that the plane of the molecule will now match the XY-plane.
	best, err := chem.BestPlane(some, nil)
	if err != nil {
		panic(err.Error())
	}
	z, _ := v3.NewMatrix([]float64{0, 0, 1})
	zero, _ := v3.NewMatrix([]float64{0, 0, 0})
	fmt.Fprintln(os.Stderr, "Best  Plane", best, z, indexes)
	axis := v3.Zeros(1)
	axis.Cross(best, z)
	fmt.Fprintln(os.Stderr, "axis", axis)
	//The main part of the program, where the rotation actually happens. Note that we rotate the whole
	//molecule, not just the planar subset, this is only used to calculate the rotation angle.
	mol.Coords[0], err = chem.RotateAbout(mol.Coords[0], zero, axis, chem.Angle(best, z))
	if err != nil {
		panic(err.Error())
	}
	//Now we write the rotated result.
	final, err := chem.XYZStringWrite(mol.Coords[0], mol)
	fmt.Print(final)
	fmt.Fprintln(os.Stderr, err)
}

//ScaleBond takes a C-H bond and moves the H (in place) so the distance between them is the one given (bond).
//CAUTION: I have only tested it for the case where the original distance>bond, although I expect it to also work in the other case.
func ScaleBond(C, H *v3.Matrix, bond float64) {
	Odist := v3.Zeros(1)
	Odist.Sub(H, C)
	distance := Odist.Norm(0)
	Odist.Scale((distance-bond)/distance, Odist)
	H.Sub(H, Odist)
}

func TestFrameDCDConc(Te *testing.T) {
	traj, err := New("../test/test.dcd")
	if err != nil {
		Te.Error(err)
	}
	frames := make([]*v3.Matrix, 3, 3)
	for i, _ := range frames {
		frames[i] = v3.Zeros(traj.Len())
	}
	results := make([][]chan *v3.Matrix, 0, 0)
	for i := 0; ; i++ {
		results = append(results, make([]chan *v3.Matrix, 0, len(frames)))
		coordchans, err := traj.NextConc(frames)
		if err != nil {
			if _, ok := err.(chem.LastFrameError); ok && coordchans == nil {
				break
			}
			Te.Error(err)
			break
		}
		for key, channel := range coordchans {
			results[len(results)-1] = append(results[len(results)-1], make(chan *v3.Matrix))
			go SecondRow(channel, results[len(results)-1][key], len(results)-1, key)
		}
		res := len(results) - 1
		for frame, k := range results[res] {
			if k == nil {
				fmt.Println(frame)
				continue
			}
			fmt.Println(res, frame, <-k)
		}
	}
}

func BackBone(stdin *bufio.Reader, options *chemjson.Options, i int) (coords, optcoords *v3.Matrix, optatoms *chem.Topology, list, frozen []int) {
	mol, coordarray, err := chemjson.DecodeMolecule(stdin, options.AtomsPerSel[i], 1)
	if err != nil {
		fmt.Fprint(os.Stderr, err.Marshal())
		log.Fatal(err)
	}
	coords = coordarray[0]

	//chem.PDBWrite("OPTpp.pdb", mol,coords,nil) /////////////////////////////////////
	resid, chain := GetResidueIds(mol)
	fmt.Fprintln(os.Stderr, "resid, chains, atomspersel, i", resid, chain, options.AtomsPerSel[i], i)
	var err2 error
	list, err2 = chem.CutBackRef(mol, []string{chain[0]}, [][]int{resid[1 : len(resid)-1]}) //in each backbone selection the chain should be the same for all residues
	if err != nil {
		panic(err2.Error()) //at least for now
	}
	optcoords = v3.Zeros(len(list))
	optcoords.SomeVecs(coords, list)
	optatoms = chem.NewTopology(nil, 0, 0) //the last 2 options are charge and multiplicity
	optatoms.SomeAtoms(mol, list)
	chem.ScaleBonds(optcoords, optatoms, "NTZ", "HNZ", chem.CHDist)
	chem.ScaleBonds(optcoords, optatoms, "CTZ", "HCZ", chem.CHDist)
	frozen = make([]int, 0, 2*len(list))
	for i := 0; i < optatoms.Len(); i++ {
		curr := optatoms.Atom(i)
		//In the future there could be an option to see whether C and N are fixed
		if curr.Name == "NTZ" || curr.Name == "CTZ" || curr.Name == "C" || curr.Name == "N" {
			frozen = append(frozen, i)
		}
	}
	return coords, optcoords, optatoms, list, frozen
}

func SideChains(stdin *bufio.Reader, options *chemjson.Options) (coords, optcoords *v3.Matrix, optatoms *chem.Topology, list, frozen []int) {
	mol, coordarray, err := chemjson.DecodeMolecule(stdin, options.AtomsPerSel[0], 1)
	if err != nil {
		fmt.Fprint(os.Stderr, err.Marshal())
		log.Fatal(err)
	}
	coords = coordarray[0]
	resid, chains := GetResidueIds(mol)
	//	fmt.Fprintln(os.Stderr,"SIDE! resid, chains", resid, chains)
	toscale := []string{"CA", "HA2", "HA3"}
	if options.BoolOptions[0][1] {
		list = chem.CutAlphaRef(mol, chains, resid)
	} else {
		list = chem.CutBetaRef(mol, chains, resid)
		toscale = []string{"CB", "HB4", "HB4"} //Yes, I am doing this twice for no reason other to have 3 elements in this slice.
	}
	optcoords = v3.Zeros(len(list))
	optcoords.SomeVecs(coords, list)
	optatoms = chem.NewTopology(nil, 0, 0) //the last 2 options are charge and multiplicity
	optatoms.SomeAtoms(mol, list)
	chem.ScaleBonds(optcoords, optatoms, toscale[0], toscale[1], chem.CHDist)
	chem.ScaleBonds(optcoords, optatoms, toscale[0], toscale[2], chem.CHDist)
	frozen = make([]int, 0, 2*len(list))
	for i := 0; i < optatoms.Len(); i++ {
		curr := optatoms.Atom(i)
		if curr.Name == "HA" || curr.Name == "CA" || curr.Name == "CB" {
			frozen = append(frozen, i)
		}
	}
	return coords, optcoords, optatoms, list, frozen
}

//Projection returns the projection of test in ref.
func Projection(test, ref *v3.Matrix) *v3.Matrix {
	rr, _ := ref.Dims()
	Uref := v3.Zeros(rr)
	Uref.Unit(ref)
	scalar := test.Dot(Uref) //math.Abs(la)*math.Cos(angle)
	Uref.Scale(scalar, Uref)
	return Uref
}

//AntiProjection returns a vector in the direction of ref with the magnitude of
//a vector A would have if |test| was the magnitude of its projection
//in the direction of test.
func AntiProjection(test, ref *v3.Matrix) *v3.Matrix {
	rr, _ := ref.Dims()
	testnorm := test.Norm(0)
	Uref := v3.Zeros(rr)
	Uref.Unit(ref)
	scalar := test.Dot(Uref)
	scalar = (testnorm * testnorm) / scalar
	Uref.Scale(scalar, Uref)
	return Uref
}

//CenterOfMass returns the center of mass the atoms represented by the coordinates in geometry
//and the masses in mass, and an error. If mass is nil, it calculates the geometric center
func CenterOfMass(geometry *v3.Matrix, mass *mat64.Dense) (*v3.Matrix, error) {
	if geometry == nil {
		return nil, CError{"goChem: nil matrix to get the center of mass", []string{"CenterOfMass"}}
	}
	gr, _ := geometry.Dims()
	if mass == nil { //just obtain the geometric center
		tmp := ones(gr)
		mass = mat64.NewDense(gr, 1, tmp) //gnOnes(gr, 1)
	}
	tmp2 := ones(gr)
	gnOnesvector := mat64.NewDense(1, gr, tmp2) //gnOnes(1, gr)

	ref := v3.Zeros(gr)
	ref.ScaleByCol(geometry, mass)
	ref2 := v3.Zeros(1)
	ref2.Mul(gnOnesvector, ref)
	ref2.Scale(1.0/mass.Sum(), ref2)
	return ref2, nil
}

func TestProjectionAndAntiProjection(Te *testing.T) {
	A := v3.Zeros(1)
	A.Set(0, 0, 2.0)
	B, _ := v3.NewMatrix([]float64{1, 1, 0})
	C := AntiProjection(A, B)
	D := Projection(B, A)
	fmt.Println("Projection of B on A (D)", D)
	fmt.Println("Anti-projection of A on B (C):", C)
	fmt.Println("Norm of C: ", C.Norm(0), " Norm of A,B: ", A.Norm(0), B.Norm(0), "Norm of D:", D.Norm(0))
}

//TestChangeAxis reads the PDB 2c9v.pdb from the test directory, collects
//The CA and CB of residue D124 of the chain A, and uses Clifford algebra to rotate the
//whole molecule such as the vector defined by these 2 atoms is
//aligned with the Z axis. The new molecule is written
//as 2c9v_aligned.pdb to the test folder.
func TestChangeAxis(Te *testing.T) {
	//runtime.GOMAXPROCS(2) ///////////////////////////
	mol, err := PDBFileRead("test/2c9v.pdb", true)
	if err != nil {
		Te.Error(err)
	}
	PDBFileWrite("test/2c9v-Readtest.pdb", mol.Coords[0], mol, nil)
	//The selection thing
	orient_atoms := [2]int{0, 0}
	for index := 0; index < mol.Len(); index++ {
		atom := mol.Atom(index)
		if atom.Chain == "A" && atom.MolID == 124 {
			if atom.Name == "CA" {
				orient_atoms[0] = index
			} else if atom.Name == "CB" {
				orient_atoms[1] = index
			}
		}
	}
	//Get the axis of rotation
	//ov1:=mol.Coord(orient_atoms[0], 0)
	ov2 := mol.Coord(orient_atoms[1], 0)
	//now we center the thing in the beta carbon of D124
	mol.Coords[0].SubVec(mol.Coords[0], ov2)
	PDBFileWrite("test/2c9v-translated.pdb", mol.Coords[0], mol, nil)
	//Now the rotation
	ov1 := mol.Coord(orient_atoms[0], 0) //make sure we have the correct versions
	ov2 = mol.Coord(orient_atoms[1], 0)  //same
	orient := v3.Zeros(ov2.NVecs())
	orient.Sub(ov2, ov1)
	//	PDBWrite(mol,"test/2c9v-124centered.pdb")
	Z, _ := v3.NewMatrix([]float64{0, 0, 1})
	axis := cross(orient, Z)
	angle := Angle(orient, Z)
	oldcoords := v3.Zeros(mol.Coords[0].NVecs())
	oldcoords.Copy(mol.Coords[0])
	mol.Coords[0] = Rotate(oldcoords, mol.Coords[0], axis, angle)
	if err != nil {
		Te.Error(err)
	}
	PDBFileWrite("test/2c9v-aligned.pdb", mol.Coords[0], mol, nil)
	fmt.Println("bench1")
}

//BestPlaneP takes sorted evecs, according to the eval,s and returns a row vector that is normal to the
//Plane that best contains the molecule. Notice that the function can't possibly check
//that the vectors are sorted. The P at the end of the name is for Performance. If
//That is not an issue it is safer to use the BestPlane function that wraps this one.
func BestPlaneP(evecs *v3.Matrix) (*v3.Matrix, error) {
	evr, evc := evecs.Dims()
	if evr != 3 || evc != 3 {
		return evecs, CError{"goChem: Eigenvectors matrix must be 3x3", []string{"BestPlaneP"}} //maybe this should be a panic
	}
	v1 := evecs.VecView(2)
	v2 := evecs.VecView(1)
	normal := v3.Zeros(1)
	normal.Cross(v1, v2)
	return normal, nil
}

func TestNWChem(Te *testing.T) {
	mol, err := chem.XYZFileRead("../test/ethanol.xyz")
	fmt.Println(mol.Coords[0], len(mol.Coords), "Quiere quedar leyenda, compadre?", err)
	if err != nil {
		Te.Error(err)
	}
	if err := mol.Corrupted(); err != nil {
		Te.Error(err)
	}
	mol.SetCharge(0)
	mol.SetMulti(1)
	calc := new(Calc)
	calc.SCFTightness = 1 //quite tight
	calc.SCFConvHelp = 1
	calc.Job = Job{Opti: true}
	calc.Method = "TPSS"
	calc.Dielectric = 4
	calc.Basis = "def2-SVP"
	calc.HighBasis = "def2-TZVP"
	calc.Grid = 4
	calc.Memory = 1000
	calc.HBAtoms = []int{2}
	calc.HBElements = []string{"O"}
	calc.CConstraints = []int{1}
	calc.SetDefaults()
	nw := NewNWChemHandle()
	orca := NewOrcaHandle()
	nw.SetName("gochem")
	orca.SetName("gochemII")
	atoms := v3.Zeros(mol.Len())
	mol.Next(atoms)
	if err = os.Chdir("../test"); err != nil {
		Te.Error(err)
	}
	err = nw.BuildInput(atoms, mol, calc)
	if err != nil {
		Te.Error(err)
	}
	_ = orca.BuildInput(atoms, mol, calc)
	//The files are already in ./test.
	os.Chdir("../test")
	defer os.Chdir("../qm")
	energy, err := nw.Energy()
	if err != nil {
		Te.Error(err)
	}
	fmt.Println("NWChem Energy: ", energy)
	newg, err := nw.OptimizedGeometry(mol)
	if err != nil {
		Te.Error(err)
	}
	chem.XYZFileWrite("optiNW.xyz", newg, mol)

}

//RotateAbout about rotates the coordinates in coordsorig around by angle radians around the axis
//given by the vector axis. It returns the rotated coordsorig, since the original is not affected.
//Uses Clifford algebra.
func RotateAbout(coordsorig, ax1, ax2 *v3.Matrix, angle float64) (*v3.Matrix, error) {
	coordsLen := coordsorig.NVecs()
	coords := v3.Zeros(coordsLen)
	translation := v3.Zeros(ax1.NVecs())
	translation.Copy(ax1)
	axis := v3.Zeros(ax2.NVecs())
	axis.Sub(ax2, ax1) // the rotation axis
	f := func() { coords.SubVec(coordsorig, translation) }
	if err := gnMaybe(gnPanicker(f)); err != nil {
		return nil, CError{err.Error(), []string{"v3.Matrix.SubVec", "RotateAbout"}}
	}
	Rot := v3.Zeros(coordsLen)
	Rot = Rotate(coords, Rot, axis, angle)
	g := func() { Rot.AddVec(Rot, translation) }
	if err := gnMaybe(gnPanicker(g)); err != nil {
		return nil, CError{err.Error(), []string{"v3.Matrix.AddVec", "RotateAbout"}}

	}
	return Rot, nil
}

//Aligns the main plane of a molecule with the XY-plane.
//Here XYZRead and XYZWrite are tested
func TestPutInXYPlane(Te *testing.T) {
	myxyz, _ := os.Open("test/sample_plane.xyz")
	mol, err := XYZRead(myxyz)
	if err != nil {
		Te.Error(err)
	}
	indexes := []int{0, 1, 2, 3, 23, 22, 21, 20, 25, 44, 39, 40, 41, 42, 61, 60, 59, 58, 63, 5}
	some := v3.Zeros(len(indexes))
	some.SomeVecs(mol.Coords[0], indexes)
	//for most rotation things it is good to have the molecule centered on its mean.
	mol.Coords[0], _, _ = MassCenter(mol.Coords[0], some, nil)
	//The test molecule is not completely planar so we use a subset of atoms that are contained in a plane
	//These are the atoms given in the indexes slice.
	some.SomeVecs(mol.Coords[0], indexes)
	//The strategy is: Take the normal to the plane of the molecule (os molecular subset), and rotate it until it matches the Z-axis
	//This will mean that the plane of the molecule will now match the XY-plane.
	best, err := BestPlane(some, nil)
	if err != nil {
		err2 := err.(Error)
		fmt.Println(err2.Decorate(""))
		Te.Fatal(err)
		//		panic(err.Error())
	}
	z, _ := v3.NewMatrix([]float64{0, 0, 1})
	zero, _ := v3.NewMatrix([]float64{0, 0, 0})
	fmt.Println("Best  Plane", best, z)
	axis := v3.Zeros(1)
	axis.Cross(best, z)
	fmt.Println("axis", axis)
	//The main part of the program, where the rotation actually happens. Note that we rotate the whole
	//molecule, not just the planar subset, this is only used to calculate the rotation angle.
	//	fmt.Println("DATA", mol.Coords[0], zero, axis, Angle(best, z))
	mol.Coords[0], err = RotateAbout(mol.Coords[0], zero, axis, Angle(best, z))
	if err != nil {
		Te.Error(err)
	}
	//	fmt.Println("after!", mol.Coords[0], err)
	//Now we write the rotated result.
	outxyz, _ := os.Create("test/Rotated.xyz") //This is the XYZWrite written file
	XYZWrite(outxyz, mol.Coords[0], mol)
}

/*	for framebunch, j := range results {
		if j == nil {
			break
		}
		for frame, k := range j {
			if k == nil {
				fmt.Println(framebunch, frame)
				continue
			}
			fmt.Println(framebunch, frame, <-k)
		}
	}
}
*/
func SecondRow(channelin, channelout chan *v3.Matrix, current, other int) {
	if channelin != nil {
		temp := <-channelin
		viej := v3.Zeros(1)
		vector := temp.VecView(2)
		viej.Copy(vector)
		fmt.Println("sending througt", channelin, channelout, viej, current, other)
		channelout <- vector
	} else {
		channelout <- nil
	}
	return
}