Golang Matrix.NVecs方法代码示例

//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
}

func EncodeCoords(coords *v3.Matrix, enc *json.Encoder) *Error {
	c := new(Coords)
	t := make([]float64, 3, 3)
	for i := 0; i < coords.NVecs(); i++ {
		c.Coords = coords.Row(t, i)
		if err := enc.Encode(c); err != nil {
			return NewError("postprocess", "chemjson.EncodeCoords", err)
		}
	}
	return nil
}

//SelCone, Given a set of cartesian points in sellist, obtains a vector "plane" normal to the best plane passing through the points.
//It selects atoms from the set A that are inside a cone in the direction of "plane" that starts from the geometric center of the cartesian points,
//and has an angle of angle (radians), up to a distance distance. The cone is approximated by a set of radius-increasing cilinders with height thickness.
//If one starts from one given point, 2 cgnOnes, one in each direction, are possible. If whatcone is 0, both cgnOnes are considered.
//if whatcone<0, only the cone opposite to the plane vector direction. If whatcone>0, only the cone in the plane vector direction.
//the 'initial' argument  allows the construction of a truncate cone with a radius of initial.
func SelCone(B, selection *v3.Matrix, angle, distance, thickness, initial float64, whatcone int) []int {
	A := v3.Zeros(B.NVecs())
	A.Copy(B) //We will be altering the input so its better to work with a copy.
	ar, _ := A.Dims()
	selected := make([]int, 0, 3)
	neverselected := make([]int, 0, 30000)     //waters that are too far to ever be selected
	nevercutoff := distance / math.Cos(angle)  //cutoff to be added to neverselected
	A, _, err := MassCenter(A, selection, nil) //Centrate A in the geometric center of the selection, Its easier for the following calculations
	if err != nil {
		panic(PanicMsg(err.Error()))
	}
	selection, _, _ = MassCenter(selection, selection, nil) //Centrate the selection as well
	plane, err := BestPlane(selection, nil)                 //I have NO idea which direction will this vector point. We might need its negative.
	if err != nil {
		panic(PanicMsg(err.Error()))
	}
	for i := thickness / 2; i <= distance; i += thickness {
		maxdist := math.Tan(angle)*i + initial //this should give me the radius of the cone at this point
		for j := 0; j < ar; j++ {
			if isInInt(selected, j) || isInInt(neverselected, j) { //we dont scan things that we have already selected, or are too far
				continue
			}
			atom := A.VecView(j)
			proj := Projection(atom, plane)
			norm := proj.Norm(0)
			//Now at what side of the plane is the atom?
			angle := Angle(atom, plane)
			if whatcone > 0 {
				if angle > math.Pi/2 {
					continue
				}
			} else if whatcone < 0 {
				if angle < math.Pi/2 {
					continue
				}
			}
			if norm > i+(thickness/2.0) || norm < (i-thickness/2.0) {
				continue
			}
			proj.Sub(proj, atom)
			projnorm := proj.Norm(0)
			if projnorm <= maxdist {
				selected = append(selected, j)
			}
			if projnorm >= nevercutoff {
				neverselected = append(neverselected, j)
			}
		}
	}
	return selected
}

//XYZStringWrite writes the mol Ref and the Coord coordinates in an XYZ-formatted string.
func XYZStringWrite(Coords *v3.Matrix, mol Atomer) (string, error) {
	var out string
	if mol.Len() != Coords.NVecs() {
		return "", CError{"Ref and Coords dont have the same number of atoms", []string{"XYZStringWrite"}}
	}
	c := make([]float64, 3, 3)
	out = fmt.Sprintf("%-4d\n\n", mol.Len())
	//towrite := Coords.Arrays() //An array of array with the data in the matrix
	for i := 0; i < mol.Len(); i++ {
		//c := towrite[i] //coordinates for the corresponding atoms
		c = Coords.Row(c, i)
		temp := fmt.Sprintf("%-2s  %12.6f%12.6f%12.6f \n", mol.Atom(i).Symbol, c[0], c[1], c[2])
		out = strings.Join([]string{out, temp}, "")
	}
	return out, nil
}

//RMSD returns the RSMD (root of the mean square deviation) for the sets of cartesian
//coordinates in test and template.
func RMSD(test, template *v3.Matrix) (float64, error) {
	//This is a VERY naive implementation.
	tmr, tmc := template.Dims()
	tsr, tsc := test.Dims()
	if tmr != tsr || tmc != 3 || tsc != 3 {
		return 0, fmt.Errorf("Ill formed matrices for RMSD calculation")
	}
	tr := tmr
	ctempla := v3.Zeros(template.NVecs())
	ctempla.Copy(template)
	//the maybe thing might not be needed since we check the dimensions before.
	f := func() { ctempla.Sub(ctempla, test) }
	if err := gnMaybe(gnPanicker(f)); err != nil {
		return 0, CError{err.Error(), []string{"v3.Matrix.Sub", "RMSD"}}
	}
	var RMSD float64
	for i := 0; i < template.NVecs(); i++ {
		temp := ctempla.VecView(i)
		RMSD += math.Pow(temp.Norm(0), 2)
	}
	RMSD = RMSD / float64(tr)
	RMSD = math.Sqrt(RMSD)
	return RMSD, nil
}

//XYZStringWrite writes the mol Ref and the Coord coordinates in an XYZ-formatted string.
func XYZWrite(out io.Writer, Coords *v3.Matrix, mol Atomer) error {
	iowriterError := func(err error) error {
		return CError{"Failed to write in io.Writer" + err.Error(), []string{"io.Writer.Write", "XYZWrite"}}
	}
	if mol.Len() != Coords.NVecs() {
		return CError{"Ref and Coords dont have the same number of atoms", []string{"XYZWrite"}}
	}
	c := make([]float64, 3, 3)
	_, err := out.Write([]byte(fmt.Sprintf("%-4d\n\n", mol.Len())))
	if err != nil {
		return iowriterError(err)
	}
	//towrite := Coords.Arrays() //An array of array with the data in the matrix
	for i := 0; i < mol.Len(); i++ {
		//c := towrite[i] //coordinates for the corresponding atoms
		c = Coords.Row(c, i)
		temp := fmt.Sprintf("%-2s  %12.6f%12.6f%12.6f \n", mol.Atom(i).Symbol, c[0], c[1], c[2])
		_, err := out.Write([]byte(temp))
		if err != nil {
			return iowriterError(err)
		}
	}
	return nil
}

//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
}

func main() {
	//This is the part that collects all the data from PyMOL, with all  the proper error checking.
	stdin := bufio.NewReader(os.Stdin)
	options, err := chemjson.DecodeOptions(stdin)
	if err != nil {
		fmt.Fprint(os.Stderr, err.Marshal())
		log.Fatal(err)
	}
	mainName := options.SelNames[0]
	if len(options.AtomsPerSel) > 1 {
		for _, v := range options.SelNames[1:] {
			mainName = mainName + "__" + v //inefficient but there should never be THAT many selections.
		}
	}
	dielectric := options.FloatOptions[0][0]
	charge := options.IntOptions[0][0]
	multi := options.IntOptions[0][1]
	qmprogram := options.StringOptions[0][0]
	method := options.StringOptions[0][1]
	calctype := options.StringOptions[0][2]
	var osidemol *chem.Topology
	var osidecoords, sidecoords *v3.Matrix
	var sidelist, sidefrozen []int
	selindex := 0
	total := 0
	selections := len(options.AtomsPerSel)
	if options.BoolOptions[0][0] { //sidechain selections exist
		sidecoords, osidecoords, osidemol, sidelist, sidefrozen = SideChains(stdin, options)
		selections--
		total += osidemol.Len()
		selindex++
	}
	fmt.Fprint(os.Stderr, selections)
	obbmol := make([]*chem.Topology, selections, selections)
	obbcoords := make([]*v3.Matrix, selections, selections)
	bbcoords := make([]*v3.Matrix, selections, selections)
	bblist := make([][]int, selections, selections)
	bbfrozen := make([][]int, selections, selections)
	for i := 0; i < selections; i++ {
		bbcoords[i], obbcoords[i], obbmol[i], bblist[i], bbfrozen[i] = BackBone(stdin, options, selindex)
		total += obbmol[i].Len()
		selindex++
		fmt.Fprint(os.Stderr, "chetumanga")
	}
	//Now we put the juit together
	bigC := v3.Zeros(total)
	bigA := chem.NewTopology([]*chem.Atom{}, 0, 0)
	bigFroz := make([]int, 0, total)
	setoffset := 0
	if options.BoolOptions[0][0] {
		bigC.SetMatrix(0, 0, osidecoords)
		setoffset += osidecoords.NVecs()
		bigA = chem.MergeAtomers(bigA, osidemol)
		//	bigA = osidemol
		bigFroz = append(bigFroz, sidefrozen...)
	}
	for k, v := range obbcoords {
		bigC.SetMatrix(setoffset, 0, v)
		bigA = chem.MergeAtomers(bigA, obbmol[k])
		tmpfroz := SliceOffset(bbfrozen[k], setoffset)
		bigFroz = append(bigFroz, tmpfroz...)
		setoffset += v.NVecs()

	}
	bigA.SetCharge(charge)
	bigA.SetMulti(multi)
	chem.PDBFileWrite(mainName+"toOPT.pdb", bigC, bigA, nil) /////////////////////////////////////
	chem.XYZFileWrite(mainName+"toOPT.xyz", bigC, bigA)      /////////////////////////////////////
	//Ok, we have now one big matrix and one big atom set, now the optimization

	calc := new(qm.Calc)
	if calctype == "Optimization" {
		calc.Optimize = true
	}
	calc.RI = true //some options, including this one, are meaningless for MOPAC
	calc.CConstraints = bigFroz
	calc.Dielectric = dielectric
	calc.SCFTightness = 1
	calc.Dispersion = "D3"
	calc.Method = "TPSS"
	if method == "Cheap" {
		calc.BSSE = "gcp"
		if qmprogram == "ORCA" {
			calc.Method = "HF-3c"
			calc.RI = false
		} else if qmprogram == "MOPAC2012" {
			calc.Method = "PM6-D3H4 NOMM MOZYME"
		} else {
			calc.Basis = "def2-SVP"
		}
	} else {
		calc.Basis = "def2-TZVP"
	}
	//We will use the default methods and basis sets of each program. In the case of MOPAC, that is currently PM6-D3H4.
	var QM qm.Handle
	switch qmprogram {
	case "ORCA":
		orca := qm.NewOrcaHandle()
		orca.SetnCPU(runtime.NumCPU())
		QM = qm.Handle(orca)
//.........这里部分代码省略.........