Golang Matrix.At方法代码示例

//RotatorToNewY takes a set of coordinates (mol) and a vector (y). It returns
//a rotation matrix that, when applied to mol, will rotate it around the Z axis
//in such a way that the projection of newy in the XY plane will be aligned with
//the Y axis.
func RotatorAroundZToNewY(newy *v3.Matrix) (*v3.Matrix, error) {
	nr, nc := newy.Dims()
	if nc != 3 || nr != 1 {
		return nil, CError{"Wrong newy vector", []string{"RotatorAroundZtoNewY"}}
	}
	if nc != 3 {
		return nil, CError{"Wrong mol vector", []string{"RotatorAroundZtoNewY"}} //this one doesn't seem reachable

	}
	gamma := math.Atan2(newy.At(0, 0), newy.At(0, 1))
	singamma := math.Sin(gamma)
	cosgamma := math.Cos(gamma)
	operator := []float64{cosgamma, singamma, 0,
		-singamma, cosgamma, 0,
		0, 0, 1}
	return v3.NewMatrix(operator)

}

//writePDBLine writes a line in PDB format from the data passed as a parameters. It takes the chain of the previous atom
//and returns the written line, the chain of the just-written atom, and error or nil.
func writePDBLine(atom *Atom, coord *v3.Matrix, bfact float64, chainprev string) (string, string, error) {
	var ter string
	var out string
	if atom.Chain != chainprev {
		ter = fmt.Sprint(out, "TER\n")
		chainprev = atom.Chain
	}
	first := "ATOM"
	if atom.Het {
		first = "HETATM"
	}
	formatstring := "%-6s%5d  %-3s %-4s%1s%4d    %8.3f%8.3f%8.3f%6.2f%6.2f          %2s  \n"
	//4 chars for the atom name are used when hydrogens are included.
	//This has not been tested
	if len(atom.Name) == 4 {
		formatstring = strings.Replace(formatstring, "%-3s ", "%-4s", 1)
	} else if len(atom.Name) > 4 {
		return "", chainprev, CError{"Cant print PDB line", []string{"writePDBLine"}}
	}
	//"%-6s%5d  %-3s %3s %1c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f          %2s  \n"
	out = fmt.Sprintf(formatstring, first, atom.ID, atom.Name, atom.Molname, atom.Chain,
		atom.MolID, coord.At(0, 0), coord.At(0, 1), coord.At(0, 2), atom.Occupancy, bfact, atom.Symbol)
	out = strings.Join([]string{ter, out}, "")
	return out, chainprev, nil
}

//RotatorToNewZ takes a matrix a row vector (newz).
//It returns a linear operator such that, when applied to a matrix mol ( with the operator on the right side)
//it will rotate mol such that the z axis is aligned with newz.
func RotatorToNewZ(newz *v3.Matrix) *v3.Matrix {
	r, c := newz.Dims()
	if c != 3 || r != 1 {
		panic("Wrong newz vector")
	}
	normxy := math.Sqrt(math.Pow(newz.At(0, 0), 2) + math.Pow(newz.At(0, 1), 2))
	theta := math.Atan2(normxy, newz.At(0, 2))      //Around the new y
	phi := math.Atan2(newz.At(0, 1), newz.At(0, 0)) //First around z
	psi := 0.000000000000                           // second around z
	sinphi := math.Sin(phi)
	cosphi := math.Cos(phi)
	sintheta := math.Sin(theta)
	costheta := math.Cos(theta)
	sinpsi := math.Sin(psi)
	cospsi := math.Cos(psi)
	operator := []float64{cosphi*costheta*cospsi - sinphi*sinpsi, -sinphi*cospsi - cosphi*costheta*sinpsi, cosphi * sintheta,
		sinphi*costheta*cospsi + cosphi*sinpsi, -sinphi*costheta*sinpsi + cosphi*cospsi, sintheta * sinphi,
		-sintheta * cospsi, sintheta * sinpsi, costheta}
	finalop, _ := v3.NewMatrix(operator) //we are hardcoding opperator so it must have the right dimensions.
	return finalop

}

//MakeWater Creates a water molecule at distance Angstroms from a2, in a direction that is angle radians from the axis defined by a1 and a2.
//Notice that the exact position of the water is not well defined when angle is not zero. One can always use the RotateAbout
//function to move the molecule to the desired location. If oxygen is true, the oxygen will be pointing to a2. Otherwise,
//one of the hydrogens will.
func MakeWater(a1, a2 *v3.Matrix, distance, angle float64, oxygen bool) *v3.Matrix {
	water := v3.Zeros(3)
	const WaterOHDist = 0.96
	const WaterAngle = 52.25
	const deg2rad = 0.0174533
	w := water.VecView(0) //we first set the O coordinates
	w.Copy(a2)
	w.Sub(w, a1)
	w.Unit(w)
	dist := v3.Zeros(1)
	dist.Sub(a1, a2)
	a1a2dist := dist.Norm(0)
	fmt.Println("ala2dist", a1a2dist, distance) ////////////////7777
	w.Scale(distance+a1a2dist, w)
	w.Add(w, a1)
	for i := 0; i <= 1; i++ {
		o := water.VecView(0)
		w = water.VecView(i + 1)
		w.Copy(o)
		fmt.Println("w1", w) ////////
		w.Sub(w, a2)
		fmt.Println("w12", w) ///////////////
		w.Unit(w)
		fmt.Println("w4", w)
		w.Scale(WaterOHDist+distance, w)
		fmt.Println("w3", w, WaterOHDist, distance)
		o.Sub(o, a2)
		t, _ := v3.NewMatrix([]float64{0, 0, 1})
		upp := v3.Zeros(1)
		upp.Cross(w, t)
		fmt.Println("upp", upp, w, t)
		upp.Add(upp, o)
		upp.Add(upp, a2)
		//water.SetMatrix(3,0,upp)
		w.Add(w, a2)
		o.Add(o, a2)
		sign := 1.0
		if i == 1 {
			sign = -1.0
		}
		temp, _ := RotateAbout(w, o, upp, deg2rad*WaterAngle*sign)
		w.SetMatrix(0, 0, temp)
	}
	var v1, v2 *v3.Matrix
	if angle != 0 {
		v1 = v3.Zeros(1)
		v2 = v3.Zeros(1)
		v1.Sub(a2, a1)
		v2.Copy(v1)
		v2.Set(0, 2, v2.At(0, 2)+1) //a "random" modification. The idea is that its not colinear with v1
		v3 := cross(v1, v2)
		v3.Add(v3, a2)
		water, _ = RotateAbout(water, a2, v3, angle)
	}
	if oxygen {
		return water
	}
	//we move things so an hydrogen points to a2 and modify the distance acordingly.
	e1 := water.VecView(0)
	e2 := water.VecView(1)
	e3 := water.VecView(2)
	if v1 == nil {
		v1 = v3.Zeros(1)
	}
	if v2 == nil {
		v2 = v3.Zeros(1)
	}
	v1.Sub(e2, e1)
	v2.Sub(e3, e1)
	axis := cross(v1, v2)
	axis.Add(axis, e1)
	water, _ = RotateAbout(water, e1, axis, deg2rad*(180-WaterAngle))
	v1.Sub(e1, a2)
	v1.Unit(v1)
	v1.Scale(WaterOHDist, v1)
	water.AddVec(water, v1)
	return water
}

//BuildInput builds an input for NWChem based int the data in atoms, coords and C.
//returns only error.
func (O *NWChemHandle) BuildInput(coords *v3.Matrix, atoms chem.AtomMultiCharger, Q *Calc) error {
	//Only error so far
	if atoms == nil || coords == nil {
		return Error{ErrMissingCharges, NWChem, O.inputname, "", []string{"BuildInput"}, true}
	}
	if Q.Basis == "" {
		log.Printf("no basis set assigned for NWChem calculation, will use the default %s, \n", O.defbasis)
		Q.Basis = O.defbasis
	}
	if Q.Method == "" {
		log.Printf("no method assigned for NWChem calculation, will use the default %s, \n", O.defmethod)
		Q.Method = O.defmethod
		Q.RI = true
	}
	if O.inputname == "" {
		O.inputname = "gochem"
	}
	//The initial guess
	vectors := fmt.Sprintf("output  %s.movecs", O.inputname) //The initial guess
	switch Q.Guess {
	case "":
	case "hcore": //This option is not a great idea, apparently.
		vectors = fmt.Sprintf("input hcore %s", vectors)
	default:
		if !Q.OldMO {
			//If the user gives something in Q.Guess but DOES NOT want an old MO to be used, I assume he/she wants to put whatever
			//is in Q.Guess directly  in the vector keyword. If you want the default put an empty string in Q.Guess.
			vectors = fmt.Sprintf("%s %s", Q.Guess, vectors)
			break
		}
		//I assume the user gave a basis set name in Q.Guess which I can use to project vectors from a previous run.
		moname := getOldMO(O.previousMO)
		if moname == "" {
			break
		}
		if strings.ToLower(Q.Guess) == strings.ToLower(Q.Basis) {
			//Useful if you only change functionals.
			vectors = fmt.Sprintf("input %s %s", moname, vectors)
		} else {
			//This will NOT work if one assigns different basis sets to different atoms.
			vectors = fmt.Sprintf("input project %s %s %s", strings.ToLower(Q.Guess), moname, vectors)
		}
	}
	vectors = "vectors " + vectors

	disp, ok := nwchemDisp[Q.Dispersion]
	if !ok {
		disp = "vdw 3"
	}
	tightness := ""
	switch Q.SCFTightness {
	case 1:
		tightness = "convergence energy 5.000000E-08\n convergence density 5.000000E-09\n convergence gradient 1E-05"
	case 2:
		//NO idea if this will work, or the criteria will be stronger than the criteria for the intergral evaluation
		//and thus the SCF will never converge. Delete when properly tested.
		tightness = "convergence energy 1.000000E-10\n convergence density 5.000000E-11\n convergence gradient 1E-07"
	}

	//For  now, the only convergence help I trust is to run a little HF calculation before and use the orbitals as a guess.
	//It works quite nicely. When the NWChem people get their shit together and fix the bugs with cgmin and RI and cgmin and
	//COSMO, cgmin will be a great option also.
	scfiters := "iterations 60"
	prevscf := ""
	if Q.SCFConvHelp > 0 {
		scfiters = "iterations 200"
		if Q.Guess == "" {
			prevscf = fmt.Sprintf("\nbasis \"3-21g\"\n * library 3-21g\nend\nset \"ao basis\" 3-21g\nscf\n maxiter 200\n vectors output hf.movecs\n %s\nend\ntask scf energy\n\n", strings.ToLower(mopacMultiplicity[atoms.Multi()]))
			vectors = fmt.Sprintf("vectors input project \"3-21g\" hf.movecs output %s.movecs", O.inputname)
		}
	}
	grid, ok := nwchemGrid[Q.Grid]
	if !ok {
		grid = "medium"
	}
	if Q.SCFTightness > 0 { //We need this if we want the SCF to converge.
		grid = "xfine"
	}
	grid = fmt.Sprintf("grid %s", grid)
	var err error

	//Only cartesian constraints supported by now.
	constraints := ""
	if len(Q.CConstraints) > 0 {
		constraints = "constraints\n fix atom"
		for _, v := range Q.CConstraints {
			constraints = fmt.Sprintf("%s %d", constraints, v+1) //goChem numbering starts from 0, apparently NWChem starts from 1, hence the v+1
		}
		constraints = constraints + "\nend"
	}

	cosmo := ""
	if Q.Dielectric > 0 {
		//SmartCosmo in a single-point means that do_gasphase False is used, nothing fancy.
		if Q.Job.Opti || O.smartCosmo {
			cosmo = fmt.Sprintf("cosmo\n dielec %4.1f\n do_gasphase False\nend", Q.Dielectric)
		} else {
			cosmo = fmt.Sprintf("cosmo\n dielec %4.1f\n do_gasphase True\nend", Q.Dielectric)
//.........这里部分代码省略.........

//BuildInput builds an input for Fermions++ based int the data in atoms, coords and C.
//returns only error.
func (O *FermionsHandle) BuildInput(coords *v3.Matrix, atoms chem.AtomMultiCharger, Q *Calc) error {
	//Only error so far

	if atoms == nil || coords == nil {
		return Error{ErrMissingCharges, Fermions, O.inputname, "", []string{"BuildInput"}, true}
	}
	if Q.Basis == "" {
		log.Printf("no basis set assigned for Fermions++ calculation, will used the default %s, \n", O.defbasis)
		Q.Basis = O.defbasis
	}
	if Q.Method == "" {
		log.Printf("no method assigned for Fermions++ calculation, will used the default %s, \n", O.defmethod)
		Q.Method = O.defmethod
	}

	disp, ok := fermionsDisp[strings.ToLower(Q.Dispersion)]
	if !ok {
		disp = "disp_corr  D3"
	}

	grid, ok := fermionsGrid[Q.Grid]
	if !ok {
		grid = "M3"
	}
	grid = fmt.Sprintf("GRID_RAD_TYPE %s", grid)
	var err error

	m := strings.ToLower(Q.Method)
	method, ok := fermionsMethods[m]
	if !ok {
		method = "EXC XC_GGA_X_PBE_R\n ECORR XC_GGA_C_PBE"
	}
	task := "SinglePoint"
	dloptions := ""
	jc := jobChoose{}
	jc.opti = func() {
		task = "DLF_OPTIMIZE"
		dloptions = fmt.Sprintf("*start::dlfind\n JOB std\n method l-bfgs\n trust_radius energy\n dcd %s.dcd\n maxcycle 300\n maxene 200\n coord_type cartesian\n*end\n", O.inputname)
		//Only cartesian constraints supported by now.
		if len(Q.CConstraints) > 0 {
			dloptions = fmt.Sprintf("%s\n*start::dlf_constraints\n", dloptions)
			for _, v := range Q.CConstraints {
				dloptions = fmt.Sprintf("%s cart %d\n", dloptions, v+1) //fortran numbering, starts with 1.
			}
			dloptions = fmt.Sprintf("%s*end\n", dloptions)
		}
	}
	Q.Job.Do(jc)
	cosmo := ""
	if Q.Dielectric > 0 {
		cosmo = fmt.Sprintf("*start::solvate\n pcm_model cpcm\n epsilon %f\n cavity_model bondi\n*end\n", Q.Dielectric)
	}

	//////////////////////////////////////////////////////////////
	//Now lets write the thing.
	//////////////////////////////////////////////////////////////
	file, err := os.Create(fmt.Sprintf("%s.in", O.inputname))
	if err != nil {
		return Error{ErrCantInput, Fermions, O.inputname, err.Error(), []string{"os.Create", "BuildInput"}, true}
	}
	defer file.Close()
	//Start with the geometry part (coords, charge and multiplicity)
	fmt.Fprintf(file, "*start::geo\n")
	fmt.Fprintf(file, "%d %d\n", atoms.Charge(), atoms.Multi())
	for i := 0; i < atoms.Len(); i++ {
		fmt.Fprintf(file, "%-2s  %8.3f%8.3f%8.3f\n", atoms.Atom(i).Symbol, coords.At(i, 0), coords.At(i, 1), coords.At(i, 2))
	}
	fmt.Fprintf(file, "*end\n\n")
	fmt.Fprintf(file, "*start::sys\n")
	fmt.Fprintf(file, " TODO %s\n", task)
	fmt.Fprintf(file, " BASIS %s\n PC pure\n", strings.ToLower(Q.Basis))
	fmt.Fprintf(file, " %s\n", method)
	fmt.Fprintf(file, " %s\n", grid)
	fmt.Fprintf(file, " %s\n", disp)
	fmt.Fprintf(file, " %s\n", O.gpu)
	if !Q.Job.Opti {
		fmt.Fprintf(file, " INFO 2\n")
	}
	fmt.Fprintf(file, "*end\n\n")
	fmt.Fprintf(file, "%s\n", cosmo)
	fmt.Fprintf(file, "%s\n", dloptions)
	return nil
}

//BuildInput builds an input for ORCA based int the data in atoms, coords and C.
//returns only error.
func (O *MopacHandle) BuildInput(coords *v3.Matrix, atoms chem.AtomMultiCharger, Q *Calc) error {
	if strings.Contains(Q.Others, "RI") {
		Q.Others = ""
	}
	//Only error so far
	if atoms == nil || coords == nil {
		return Error{ErrMissingCharges, Mopac, O.inputname, "", []string{"BuildInput"}, true}
	}
	ValidMethods := []string{"PM3", "PM6", "PM7", "AM1"}
	if Q.Method == "" || !isInString(ValidMethods, Q.Method[0:3]) { //not found
		log.Printf("no method assigned for MOPAC calculation, will used the default %s, \n", O.defmethod)
		Q.Method = O.defmethod
	}
	opt := "" //Empty string means optimize
	jc := jobChoose{}
	jc.sp = func() {
		opt = "1SCF"
	}
	jc.opti = func() {}
	Q.Job.Do(jc)
	//If this flag is set we'll look for a suitable MO file.
	//If not found, we'll just use the default ORCA guess
	hfuhf := "RHF"
	if atoms.Multi() != 1 {
		hfuhf = "UHF"
	}
	cosmo := ""
	if Q.Dielectric > 0 {
		cosmo = fmt.Sprintf("EPS=%2.1f RSOLV=1.3 LET DDMIN=0.0", Q.Dielectric) //The DDMIN ensures that the optimization continues when cosmo is used. From the manual I understand that it is OK
	}
	strict := ""
	if Q.SCFTightness > 0 {
		strict = "GNORM=0.01 RELSCF=0.00001"
	}
	multi := mopacMultiplicity[atoms.Multi()]
	charge := fmt.Sprintf("CHARGE=%d", atoms.Charge())
	MainOptions := []string{hfuhf, Q.Method, strict, opt, cosmo, charge, multi, Q.Others, "BONDS AUX\n"}
	mainline := strings.Join(MainOptions, " ")
	//Now lets write the thing
	if O.inputname == "" {
		O.inputname = "input"
	}
	file, err := os.Create(fmt.Sprintf("%s.mop", O.inputname))
	if err != nil {
		return Error{ErrCantInput, Mopac, O.inputname, "", []string{"os.Create", "BuildInput"}, true}

	}
	defer file.Close()
	if _, err = fmt.Fprint(file, "* ===============================\n* Input file for Mopac\n* ===============================\n"); err != nil {
		return Error{ErrCantInput, Mopac, O.inputname, "", []string{"BuildInput"}, true}
		//After this check I just assume the file is ok and dont check again.
	}
	fmt.Fprint(file, mainline)
	fmt.Fprint(file, "\n")
	fmt.Fprint(file, "Mopac file generated by gochem :-)\n")
	//now the coordinates
	for i := 0; i < atoms.Len(); i++ {
		tag := 1
		if isInInt(Q.CConstraints, i) {
			tag = 0
		}
		//	fmt.Println(atoms.Atom(i).Symbol)
		fmt.Fprintf(file, "%-2s  %8.5f %d %8.5f %d %8.5f %d\n", atoms.Atom(i).Symbol, coords.At(i, 0), tag, coords.At(i, 1), tag, coords.At(i, 2), tag)
	}
	fmt.Fprintf(file, "\n")
	return nil
}