SciCodePile/SciCode-Domain-Code

1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Computational Biochemistry","rmera/gochem","atomicdata.go",".go","3454","157","/*3 * atomicdata.go, part of gochem.4 *5 *6 * Copyright 2021 Raul Mera <rmera{at}chemDOThelsinkiDOTfi>7 *8 * This program is free software; you can redistribute it and/or modify9 * it under the terms of the GNU Lesser General Public License as10 * published by the Free Software Foundation; either version 2.1 of the11 * License, or (at your option) any later version.12 *13 * This program is distributed in the hope that it will be useful,14 * but WITHOUT ANY WARRANTY; without even the implied warranty of15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the16 * GNU General Public License for more details.17 *18 * You should have received a copy of the GNU Lesser General19 * Public License along with this program.  If not, see20 * <http://www.gnu.org/licenses/>.21 *22 *23 * goChem is currently developed at the Universidad de Santiago de Chile24 * (USACH)25 *26 */27 28package chem29 30// A map for assigning mass to elements.31// Note that just common ""bio-elements"" are present32var symbolMass = map[string]float64{33	""H"":  1.0,34	""C"":  12.01,35	""O"":  16.00,36	""N"":  14.01,37	""P"":  30.97,38	""S"":  32.06,39	""Se"": 78.96,40	""K"":  39.1,41	""Ca"": 40.08,42	""Mg"": 24.30,43	""Cl"": 35.45,44	""Na"": 22.99,45	""Cu"": 63.55,46	""Zn"": 65.38,47	""Co"": 58.93,48	""Fe"": 55.84,49	""Mn"": 54.94,50	""Cr"": 51.996,51	""Si"": 28.08,52	""Be"": 9.012,53	""F"":  18.998,54	""Br"": 79.904,55	""I"":  126.90,56}57 58// A map for assigning covalent radii to elements59// Values from Cordero et al., 2008 (DOI:10.1039/B801115J)60// Note that just common ""bio-elements"" are present61var symbolCovrad = map[string]float64{62	""H"":  0.4,  // 0.31 I altered this one. Since H always has only one bond, it doesn't matter if I set a longer radius, the extra bonds will get eliminated later.63	""C"":  0.76, //the sp3 radius64	""O"":  0.66,65	""N"":  0.71,66	""P"":  1.07,67	""S"":  1.05,68	""Se"": 1.2,69	""K"":  2.03,70	""Ca"": 1.76,71	""Mg"": 1.41,72	""Cl"": 1.02,73	""Na"": 1.66,74	""Cu"": 1.32,75	""Zn"": 1.22,76	""Co"": 1.5,  // hs77	""Fe"": 1.52, //hs78	""Mn"": 1.61, //hs79	""Cr"": 1.39,80	""Si"": 1.11,81	""Be"": 0.96,82	""F"":  0.57,83	""Br"": 1.2,84	""I"":  1.39,85}86 87// A map for assigning van der Waals radii to elements88// Values from 10.1021/j100785a001 and 10.1021/jp811155689// metal radii from 10.1023/A:101162572880390// Note that just common ""bio-elements"" are present91var symbolVdwrad = map[string]float64{92	""H"":  1.10, // 0.31 I altered this one. Since H always has only one bond, it doesn't matter if I set a longer radius, the extra bonds will get eliminated later.93	""C"":  1.70, //the sp3 radius94	""O"":  1.52,95	""N"":  1.55,96	""P"":  1.80,97	""S"":  1.80,98	""Se"": 1.90,99	""K"":  2.75,100	""Ca"": 2.31,101	""Mg"": 1.73,102	""Cl"": 1.75,103	""Na"": 2.27,104	""Cu"": 2.00,105	""Zn"": 2.02,106	""Co"": 1.95,107	""Fe"": 1.96,108	""Mn"": 1.96,109	""Cr"": 1.97,110	""Si"": 2.10,111	""Be"": 1.53,112	""F"":  1.47,113	""Br"": 1.83,114	""I"":  1.98,115}116 117// A map for checking that atoms don't118// have too many bonds. A value of 0 means119// undefined, i.e. that this atom shouldn't120// be checked for max bonds. I decided not to define it121var symbolMaxBonds = map[string]int{122	""H"":  1, //this is the only one truly important.123	""C"":  4,124	""O"":  2,125	""N"":  0, //undefined126	""P"":  0,127	""S"":  0,128	""Se"": 0,129	""Be"": 0,130	""F"":  1,131	""Br"": 1,132	""I"":  1,133}134 135var Three2OneLetter = map[string]string{136	""SER"": ""S"",137	""THR"": ""T"",138	""ASN"": ""N"",139	""GLN"": ""Q"",140	""SEC"": ""U"", //Selenocysteine!141	""CYS"": ""C"",142	""GLY"": ""G"",143	""PRO"": ""P"",144	""ALA"": ""A"",145	""VAL"": ""V"",146	""ILE"": ""I"",147	""LEU"": ""L"",148	""MET"": ""M"",149	""PHE"": ""F"",150	""TYR"": ""Y"",151	""TRP"": ""W"",152	""ARG"": ""R"",153	""HIS"": ""H"",154	""LYS"": ""K"",155	""ASP"": ""D"",156	""GLU"": ""E"",157}158","Go"
159"Computational Biochemistry","rmera/gochem","handy.go",".go","29609","915","/*160 * handy.go, part of gochem.161 *162 *163 * Copyright 2012 Raul Mera <rmera{at}chemDOThelsinkiDOTfi>164 *165 * This program is free software; you can redistribute it and/or modify166 * it under the terms of the GNU Lesser General Public License as167 * published by the Free Software Foundation; either version 2.1 of the168 * License, or (at your option) any later version.169 *170 * This program is distributed in the hope that it will be useful,171 * but WITHOUT ANY WARRANTY; without even the implied warranty of172 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the173 * GNU General Public License for more details.174 *175 * You should have received a copy of the GNU Lesser General176 * Public License along with this program.  If not, see177 * <http://www.gnu.org/licenses/>.178 *179 *180 * Gochem is developed at the laboratory for instruction in Swedish, Department of Chemistry,181 * University of Helsinki, Finland.182 *183 *184 */185 186package chem187 188import (189	""fmt""190	""math""191	""strings""192 193	v3 ""github.com/rmera/gochem/v3""194)195 196// NegateIndexes, given a set of indexes and the length of a molecule, produces197// a set of all the indexes _not_ in the original set.198func NegateIndexes(indexes []int, length int) []int {199	ret := make([]int, 0, length-len(indexes))200	for i := 0; i < length; i++ {201		if !isInInt(indexes, i) {202			ret = append(ret, i)203 204		}205	}206	return ret207}208 209const allchains = ""*ABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789""210 211// FixGromacsPDB fixes the problem that Gromacs PDBs have when there are more than 10000 residues212// Gromacs simply restarts the numbering. Since solvents (where this is likely to happen) don't have213// chain ID in Gromacs, it's impossible to distinguish between the water 1 and the water 10001. FixGromacsPDB214// Adds a chain ID to the newly restrated residue that is the letter/symbol coming after the last seen chain ID215// in the constant allchains defined in this file. The current implementation does nothing if a chain ID is already216// defined, even if it is wrong (if 9999 and the following 0 residue have the same chain).217func FixGromacsPDB(mol Atomer) {218	//	fmt.Println(""FIXING!"")219	previd := 999999220	const pdbmaxresidue = 9999221	lastchain := ""*""222	j := 1223	for i := 0; i < mol.Len(); i++ {224		at := mol.Atom(i)225		if at.MolID > pdbmaxresidue {226			at.MolID = j227			if j == pdbmaxresidue {228				j = 0229			}230			j++231		}232		if at.Chain == "" "" {233			if previd > at.MolID {234				index := strings.Index(allchains, lastchain) + 1235				//	fmt.Println(""new chain index:"", index)236				lastchain = string(allchains[index])237			}238			at.Chain = lastchain239			//		fmt.Println(lastchain) /////////240		} else {241			lastchain = at.Chain242		}243		previd = at.MolID244 245		//a fix for Martini Waters246		if at.MolName == ""WN"" || at.MolName == ""WN "" || at.MolName == "" WN"" {247			at.MolName = ""WNN""248		}249 250	}251}252 253// Molecules2Atoms gets a selection list from a list of residues.254// It select all the atoms that form part of the residues in the list.255// It doesnt return errors. If a residue is out of range, no atom will256// be returned for it. Atoms are also required to be part of one of the chains257// specified in chains, but a nil ""chains"" can be given to select all chains.258func Molecules2Atoms(mol Atomer, residues []int, chains []string) []int {259	atlist := make([]int, 0, len(residues)*3)260	for key := 0; key < mol.Len(); key++ {261		at := mol.Atom(key)262		if isInInt(residues, at.MolID) && (isInString(chains, at.Chain) || len(chains) == 0) {263			atlist = append(atlist, key)264		}265	}266	return atlist267 268}269 270// EasyShape takes a matrix of coordinates, a value for epsilon (a number close to271// zero, the closer, the more272// strict the orthogonality requriements are) and an (optative) masser and returns273// two shape indicators based on the elipsoid of inertia (or it massless equivalent)274// a linear and circular distortion indicators, as percentages, and an error or275// nil (in that order). If you give a negative number as epsilon, the default276// (quite strict) will be used.277func EasyShape(coords *v3.Matrix, epsilon float64, mol ...Masser) (float64, float64, error) {278	var masses []float64279	var err2 error280	var err error281	if len(mol) == 0 {282		masses = nil283	} else {284		masses, err = mol[0].Masses()285		if err != nil {286			masses = nil287			err2 = err288		}289	}290	moment, err := MomentTensor(coords, masses)291	if err != nil {292		return -1, -1, err293	}294	rhos, err := Rhos(moment, epsilon)295	if err != nil {296		return -1, -1, err297	}298	linear, circular, err := RhoShapeIndexes(rhos)299	if err != nil {300		return -1, -1, err301	}302	return linear, circular, err2303}304 305// MolIDNameChain2Index takes a molID (residue number), atom name, chain index and a molecule Atomer.306// it returns the index associated with the atom in question in the Ref. The function returns also an error (if failure of warning)307// or nil (if succses and no warnings). Note that this function is not efficient to call several times to retrieve many atoms.308func MolIDNameChain2Index(mol Atomer, molID int, name, chain string) (int, error) {309	var ret int = -1310	var err error311	if mol == nil {312		return -1, CError{""goChem: Given a nil chem.Atomer"", []string{""MolIDNameChain2Index""}}313	}314	for i := 0; i != mol.Len(); i++ {315		a := mol.Atom(i)316		if a.Name == """" && err == nil {317			err = CError{""Warning: The Atoms does not seem to contain PDB-type information"", []string{""MolIDNameChain2Index""}} //We set this error but will still keep running the function in case the data is present later in the molecule.318		}319		if a.MolID == molID && a.Name == name && a.Chain == chain {320			ret = i321			break322		}323 324	}325	if ret == -1 {326		var p string327		if err != nil {328			p = err.Error()329		}330		err = CError{fmt.Sprintf(""%s, No atomic index found in the Atomer given for the given MolID, atom name and chain. %s %d"", p, chain, molID), []string{""MolIDNameChain2Index""}}331	}332	return ret, err333}334 335// OnesMass returns a column matrix with lenght rosw.336// This matrix can be used as a dummy mass matrix337// for geometric calculations.338func OnesMass(lenght int) *v3.Matrix {339	return v3.Dense2Matrix(gnOnes(lenght, 1))340}341 342// Super determines the best rotation and translations to superimpose the coords in test343// considering only the atoms present in the slices of int slices indexes.344// The first indexes slices will be assumed to contain test indexes and the second, template indexes.345// If you give only one, it will be assumed to correspond to test, if test has more atoms than346// elements on the indexes set, or templa, otherwise. If no indexes are given, all atoms on each system347// will be superimposed. The number of atoms superimposed on both systems must be equal.348// Super modifies the test matrix, but template and indexes are not touched.349func Super(test, templa *v3.Matrix, indexes ...[]int) (*v3.Matrix, error) {350	var ctest *v3.Matrix351	var ctempla *v3.Matrix352	if len(indexes) == 0 || indexes[0] == nil || len(indexes[0]) == 0 { //If you put the date in the SECOND slice, you are just messing with me.353		ctest = test354		ctempla = templa355	} else if len(indexes) == 1 {356		if test.NVecs() > len(indexes[0]) {357			ctest = v3.Zeros(len(indexes[0]))358			ctest.SomeVecs(test, indexes[0])359			ctempla = templa360		} else if templa.NVecs() > len(indexes[0]) {361			ctempla = v3.Zeros(len(indexes[0]))362			ctempla.SomeVecs(templa, indexes[0])363		} else {364			return nil, fmt.Errorf(""chem.Super: Indexes don't match molecules"")365		}366	} else {367		ctest = v3.Zeros(len(indexes[0]))368		ctest.SomeVecs(test, indexes[0])369		ctempla = v3.Zeros(len(indexes[1]))370		ctempla.SomeVecs(templa, indexes[1])371	}372 373	if ctest.NVecs() != ctempla.NVecs() {374		return nil, fmt.Errorf(""chem.Super: Ill formed coordinates for Superposition"")375	}376 377	_, rotation, trans1, trans2, err1 := RotatorTranslatorToSuper(ctest, ctempla)378	if err1 != nil {379		return nil, errDecorate(err1, ""Super"")380	}381	test.AddVec(test, trans1)382	//	fmt.Println(""test1"",test, rotation) /////////////77383	test.Mul(test, rotation)384	//	fmt.Println(""test2"",test) ///////////385	test.AddVec(test, trans2)386	//	fmt.Println(""test3"",test) ///////387	return test, nil388}389 390// RotateAbout about rotates the coordinates in coordsorig around by angle radians around the axis391// given by the vector axis. It returns the rotated coordsorig, since the original is not affected.392// Uses Clifford algebra.393func RotateAbout(coordsorig, ax1, ax2 *v3.Matrix, angle float64) (*v3.Matrix, error) {394	coordsLen := coordsorig.NVecs()395	coords := v3.Zeros(coordsLen)396	translation := v3.Zeros(ax1.NVecs())397	translation.Copy(ax1)398	axis := v3.Zeros(ax2.NVecs())399	axis.Sub(ax2, ax1) // the rotation axis400	f := func() { coords.SubVec(coordsorig, translation) }401	if err := gnMaybe(gnPanicker(f)); err != nil {402		return nil, CError{err.Error(), []string{""v3.Matrix.SubVec"", ""RotateAbout""}}403	}404	Rot := v3.Zeros(coordsLen)405	Rot = Rotate(coords, Rot, axis, angle)406	g := func() { Rot.AddVec(Rot, translation) }407	if err := gnMaybe(gnPanicker(g)); err != nil {408		return nil, CError{err.Error(), []string{""v3.Matrix.AddVec"", ""RotateAbout""}}409 410	}411	return Rot, nil412}413 414// MatchAxes returns a rotated version of mol such that the vectors ax1 and ax2 are superimposed415// mol is, by default, centered on center (which should be the origin of both ax1 and ax2)416// unless recenter is given and true, in which case, it is not modified.417func MatchAxes(mol, ax1, ax2, center *v3.Matrix, recenter ...bool) (*v3.Matrix, error) {418	r2, c2 := ax2.Dims()419	if c2 != 3 || r2 != 1 {420		panic(""Wrong ax2 vector"")421	}422 423	r1, c1 := ax1.Dims()424	if c1 != 3 || r1 != 1 {425		panic(""Wrong ax1 vector"")426	}427	//let's center all in 'center'428	ax1.Sub(ax1, center)429	ax2.Sub(ax2, center)430	mol.Sub(mol, center)431	normal := v3.Zeros(1)432	normal.Cross(ax1, ax2) //this vector should be the axis of rotation.433	angle := Angle(ax1, ax2)434	zero := v3.Zeros(1)435	rot, err := RotateAbout(mol, normal, zero, angle)436	if err != nil {437		return nil, err438	}439	//we now undo the centering for all of our vectors440	ax1.Add(ax1, center)441	ax2.Add(ax2, center)442	if len(recenter) > 0 && recenter[0] {443		mol.Add(mol, center)444	}445	rot.Add(rot, center)446	return rot, nil447}448 449// EulerRotateAbout uses Euler angles to rotate the coordinates in coordsorig around by angle450// radians around the axis given by the vector axis. It returns the rotated coordsorig,451// since the original is not affected. It seems more clunky than the RotateAbout, which uses Clifford algebra.452// I leave it for benchmark, mostly, and might remove it later. There is no test for this function!453func EulerRotateAbout(coordsorig, ax1, ax2 *v3.Matrix, angle float64) (*v3.Matrix, error) {454	r, _ := coordsorig.Dims()455	coords := v3.Zeros(r)456	translation := v3.Zeros(ax1.NVecs())457	translation.Copy(ax1)458	axis := v3.Zeros(ax2.NVecs())459	axis.Sub(ax2, ax1) //now it became the rotation axis460	f := func() { coords.SubVec(coordsorig, translation) }461	if err := gnMaybe(gnPanicker(f)); err != nil {462		return nil, CError{err.Error(), []string{""v3.Matrix.Subvec"", ""EulerRotateAbout""}}463 464	}465	Zswitch := RotatorToNewZ(axis)466	coords.Mul(coords, Zswitch) //rotated467	Zrot, err := RotatorAroundZ(angle)468	if err != nil {469		return nil, errDecorate(err, ""EulerRotateAbout"")470	}471	//	Zsr, _ := Zswitch.Dims()472	//	RevZ := v3.Zeros(Zsr)473	RevZ, err := gnInverse(Zswitch, nil)474	if err != nil {475		return nil, errDecorate(err, ""EulerRotateAbout"")476	}477	coords.Mul(coords, Zrot) //rotated478	coords.Mul(coords, RevZ)479	coords.AddVec(coords, translation)480	return coords, nil481}482 483// Corrupted is a convenience function to check that a reference and a trajectory have the same number of atoms484func Corrupted(X Traj, R Atomer) error {485	if X.Len() != R.Len() {486		return CError{""Mismatched number of atoms/coordinates"", []string{""Corrupted""}}487	}488	return nil489}490 491//Some internal convenience functions.492 493// isInInt is a helper for the RamaList function,494// returns true if test is in container, false otherwise.495func isInInt(container []int, test int) bool {496	if container == nil {497		return false498	}499	for _, i := range container {500		if test == i {501			return true502		}503	}504	return false505}506 507// Same as the previous, but with strings.508func isInString(container []string, test string) bool {509	if container == nil {510		return false511	}512	for _, i := range container {513		if test == i {514			return true515		}516	}517	return false518}519 520// ""caps"" the atom tocap (which must be part of mol) with an H atom.521// The reference atom, if not nil, is employed to define the new bond as opposite to522func CapWithH(mol *Molecule, tocap int, position *v3.Matrix, atomIndex, bondIndex int) *Molecule {523	if atomIndex < 0 {524		for i := 0; i < mol.Len(); i++ {525			in := mol.Atom(i).Index()526			if in > atomIndex {527				atomIndex = in528			}529		}530		atomIndex++ //one more than the last one531	}532	if bondIndex < 0 {533		for _, v := range mol.Bonds {534			in := v.Index535			if in > bondIndex {536				bondIndex = in537			}538		}539		bondIndex++540	}541	for i := 0; i < mol.Len(); i++ {542		if mol.Atom(i).Index() == tocap {543			tocap = i544			break545		}546	}547	cp := &Atom{Symbol: ""H"", index: atomIndex, ID: mol.Len(), Name: ""H""}548	//We'll use a _very_ naive way of determining the position of the new atom, of not given.549	//Hopefully it's good enough that it can be fixed by optimization.550	if position == nil {551		r := mol.Coords[0].VecView(tocap)552		position = v3.Zeros(0)553		position.Set(0, 0, r.At(0, 0))554		position.Set(0, 1, r.At(0, 1))555		position.Set(0, 2, r.At(0, 2)+CHDist)556	}557	bond := &Bond{At1: mol.Atom(tocap), At2: cp, Dist: CHDist, Index: bondIndex}558	cp.Bonds = []*Bond{bond}559	mol.Atom(tocap).Bonds = append(mol.Atom(tocap).Bonds, bond)560	mol.Atoms = append(mol.Atoms, cp)561	mol.Bonds = append(mol.Bonds, bond)562	r := v3.Zeros(1)563	r.Copy(mol.Coords[0].VecView(tocap))564	ScaleBond(r, position, CHDist) //Check that this works.565	ncoords := v3.Zeros(mol.Len())566	ncoords.StackVec(mol.Coords[0], position)567	mol.Coords[0] = ncoords568	return mol569}570 571// ScaleBond  moves the atom at2 (in place) so the distance between it and a1 is the one given (newdist).572// 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.573func ScaleBond(a1, a2 *v3.Matrix, newdist float64) {574	Odist := v3.Zeros(1)575	Odist.Sub(a1, a2)576	distance := Odist.Norm(2)577	//	println(""dists"", distance, newdist) /////////////////////////578	scaling := math.Abs(distance-newdist) / distance579	if distance < newdist {580		scaling = 1 / scaling581	}582	Odist.Scale(scaling, Odist)583	a2b := v3.Zeros(1)584	a2b.Copy(a2)585	a2.Add(a2b, Odist)586	//DEBUG587	//	Odist.Sub(a1, a2)588	//	distance = Odist.Norm(2)589	//	println(""distsfinal"", distance, newdist) /////////////////////////590 591}592 593// 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.594// Notice that the exact position of the water is not well defined when angle is not zero. One can always use the RotateAbout595// function to move the molecule to the desired location. If oxygen is true, the oxygen will be pointing to a2. Otherwise,596// one of the hydrogens will.597func MakeWater(a1, a2 *v3.Matrix, distance, angle float64, oxygen bool) *v3.Matrix {598	water := v3.Zeros(3)599	const WaterOHDist = 0.96600	const WaterAngle = 52.25601	const deg2rad = 0.0174533602	w := water.VecView(0) //we first set the O coordinates603	w.Copy(a2)604	w.Sub(w, a1)605	w.Unit(w)606	dist := v3.Zeros(1)607	dist.Sub(a1, a2)608	a1a2dist := dist.Norm(2)609	//	fmt.Println(""ala2dist"", a1a2dist, distance) ////////////////7777610	w.Scale(distance+a1a2dist, w)611	w.Add(w, a1)612	for i := 0; i <= 1; i++ {613		o := water.VecView(0)614		w = water.VecView(i + 1)615		w.Copy(o)616		//		fmt.Println(""w1"", w) ////////617		w.Sub(w, a2)618		//		fmt.Println(""w12"", w) ///////////////619		w.Unit(w)620		//		fmt.Println(""w4"", w)621		w.Scale(WaterOHDist+distance, w)622		//		fmt.Println(""w3"", w, WaterOHDist, distance)623		o.Sub(o, a2)624		t, _ := v3.NewMatrix([]float64{0, 0, 1})625		upp := v3.Zeros(1)626		upp.Cross(w, t)627		//		fmt.Println(""upp"", upp, w, t)628		upp.Add(upp, o)629		upp.Add(upp, a2)630		//water.SetMatrix(3,0,upp)631		w.Add(w, a2)632		o.Add(o, a2)633		sign := 1.0634		if i == 1 {635			sign = -1.0636		}637		temp, _ := RotateAbout(w, o, upp, deg2rad*WaterAngle*sign)638		w.SetMatrix(0, 0, temp)639	}640	var v1, v2 *v3.Matrix641	if angle != 0 {642		v1 = v3.Zeros(1)643		v2 = v3.Zeros(1)644		v1.Sub(a2, a1)645		v2.Copy(v1)646		v2.Set(0, 2, v2.At(0, 2)+1) //a ""random"" modification. The idea is that its not colinear with v1647		v3 := cross(v1, v2)648		v3.Add(v3, a2)649		water, _ = RotateAbout(water, a2, v3, angle)650	}651	if oxygen {652		return water653	}654	//we move things so an hydrogen points to a2 and modify the distance acordingly.655	e1 := water.VecView(0)656	e2 := water.VecView(1)657	e3 := water.VecView(2)658	if v1 == nil {659		v1 = v3.Zeros(1)660	}661	if v2 == nil {662		v2 = v3.Zeros(1)663	}664	v1.Sub(e2, e1)665	v2.Sub(e3, e1)666	axis := cross(v1, v2)667	axis.Add(axis, e1)668	water, _ = RotateAbout(water, e1, axis, deg2rad*(180-WaterAngle))669	v1.Sub(e1, a2)670	v1.Unit(v1)671	v1.Scale(WaterOHDist, v1)672	water.AddVec(water, v1)673	return water674}675 676// FixNumbering will put the internal numbering+1 in the atoms and residue fields, so they match the current residues/atoms677// in the molecule678func FixNumbering(r Atomer) {679	resid := 0680	prevres := -1681	for i := 0; i < r.Len(); i++ {682		at := r.Atom(i)683		at.ID = i + 1684		if prevres != at.MolID {685			prevres = at.MolID686			resid++687		}688		at.MolID = resid689	}690}691 692// CutBackRef takes a list of lists of residues and selects693// from r all atoms in each the list list[i] and belonging to the chain chain[i].694// It caps the N and C terminal695// of each list with -COH for the N terminal and NH2 for C terminal.696// the residues on each sublist should be contiguous to each other.697// for instance, {6,7,8} is a valid sublist, {6,8,9} is not.698// This is NOT currently checked by the function!. It returns the list of kept atoms699func CutBackRef(r Atomer, chains []string, list [][]int) ([]int, error) {700	//i:=r.Len()701	if len(chains) != len(list) {702		return nil, CError{fmt.Sprintf(""Mismatched chains (%d) and list (%d) slices"", len(chains), len(list)), []string{""CutBackRef""}}703	}704	var ret []int //This will be filled with the atoms that are kept, and will be returned.705	for k, v := range list {706		nter := v[0]707		cter := v[len(v)-1]708		nresname := """"709		cresname := """"710		for j := 0; j < r.Len(); j++ {711			if r.Atom(j).MolID == nter && r.Atom(j).Chain == chains[k] {712				nresname = r.Atom(j).MolName713				break714			}715		}716		if nresname == """" {717			//we will protest if the Nter is not found. If Cter is not found we will just718			//cut at the real Cter719			return nil, CError{fmt.Sprintf(""list %d contains residue numbers out of boundaries"", k), []string{""CutBackRef""}}720 721		}722		for j := 0; j < r.Len(); j++ {723			curr := r.Atom(j)724			if curr.Chain != chains[k] {725				continue726			}727			if curr.MolID == cter {728				cresname = curr.MolName729			}730			if curr.MolID == nter-1 {731				makeNcap(curr, nresname)732			}733			if curr.MolID == cter+1 {734				makeCcap(curr, cresname)735			}736		}737	}738	for _, i := range list {739		t := Molecules2Atoms(r, i, chains)740		//	fmt.Println(""t"", len(t))741		ret = append(ret, t...)742	}743	//	j:=0744	//	for i:=0;;i++{745	//		index:=i-j746	//		if index>=r.Len(){747	//			break748	//		}749	//		if !isInInt(ret, i){750	//			r.DelAtom(index)751	//			j++752	//		}753	//	}754	return ret, nil755}756 757func makeNcap(at *Atom, resname string) {758	if !isInString([]string{""C"", ""O"", ""CA""}, at.Name) {759		return760	}761	at.MolID = at.MolID + 1762	at.MolName = resname763	if at.Name == ""C"" {764		at.Name = ""CTZ""765	}766	if at.Name == ""CA"" {767		at.Name = ""HCZ""768		at.Symbol = ""H""769	}770}771 772func makeCcap(at *Atom, resname string) {773	if !isInString([]string{""N"", ""H"", ""CA""}, at.Name) {774		return775	}776	at.MolID = at.MolID - 1777	at.MolName = resname778	if at.Name == ""N"" {779		at.Name = ""NTZ""780	}781	if at.Name == ""CA"" {782		at.Name = ""HNZ""783		at.Symbol = ""H""784	}785}786 787/*788//Takes a list of lists of residues and produces a new set of coordinates789//whitout any atom not in the lists or not from the chain chain. It caps the N and C terminal790//of each list with -COH for the N terminal and NH2 for C terminal.791//the residues on each sublist should be contiguous to each other.792//for instance, {6,7,8} is a valid sublist, {6,8,9} is not.793//This is NOT currently checked by the function!794//In addition, the Ref provided should have already been processed by795//CutBackRef, which is not checked either.796func CutBackCoords(r Ref, coords *v3.Matrix, chain string, list [][]int) (*v3.Matrix, error) {797	//this is actually a really silly function. So far I dont check for errors, but I keep the return balue798	//In case I do later.799	var biglist []int800	for _, i := range list {801		smallist := Molecules2Atoms(r, i, []string{chain})802		biglist = append(biglist, smallist...)803	}804	NewVecs := v3.Zeros(len(biglist), 3)805	NewVecs.SomeVecs(coords, biglist)806	return NewVecs, nil807 808}809*/810 811// CutLateralRef will return a list with the atom indexes of the lateral chains of the residues in list812// for each of these residues it will change the alpha carbon to oxygen and change the residue number of the rest813// of the backbone to -1.814func CutBetaRef(r Atomer, chain []string, list []int) []int {815	//	pairs := make([][]int,1,10)816	//	pairs[0]=make([]int,0,2)817	for i := 0; i < r.Len(); i++ {818		curr := r.Atom(i)819		if isInInt(list, curr.MolID) && isInString(chain, curr.Chain) {820			if curr.Name == ""CB"" {821				//			pairs[len(pairs)-1][1]=i //I am assuming that CA will show before CB in the PDB, which is rather weak822				//		paairs=append(pairs,make([]int,1,2))823			}824			if curr.Name == ""CA"" {825				curr.Name = ""HB4""826				curr.Symbol = ""H""827				//		pairs[len(pairs)-1]=append(pairs[len(pairs)-1],i)828			} else if isInString([]string{""C"", ""H"", ""HA"", ""O"", ""N""}, curr.Name) { //change the res number of the backbone so it is not considered829				curr.MolID = -1830			}831 832		}833	}834	newlist := Molecules2Atoms(r, list, chain)835	return newlist836}837 838// CutAlphaRef will return a list with the atoms in the residues indicated by list, in the chains given.839// The carbonyl carbon and amide nitrogen for each residue will be transformer into hydrogens. The MolID of the840// other backbone atoms will be set to -1 so they are no longer considered.841func CutAlphaRef(r Atomer, chain []string, list []int) []int {842	for i := 0; i < r.Len(); i++ {843		curr := r.Atom(i)844		if isInInt(list, curr.MolID) && isInString(chain, curr.Chain) {845			if curr.Name == ""C"" {846				curr.Name = ""HA2""847				curr.Symbol = ""H""848			} else if curr.Name == ""N"" {849				curr.Name = ""HA3""850				curr.Symbol = ""H""851			} else if isInString([]string{""H"", ""O""}, curr.Name) { //change the res number of the backbone so it is not considered852				curr.MolID = -1853			}854 855		}856	}857	newlist := Molecules2Atoms(r, list, chain)858	return newlist859}860 861// TagAtomsByName will tag all atoms with a given name in a given list of atoms.862// return the number of tagged atoms863func TagAtomsByName(r Atomer, name string, list []int) int {864	tag := 0865	for i := 0; i < r.Len(); i++ {866		curr := r.Atom(i)867		if isInInt(list, i) && curr.Name == name {868			curr.Tag = 1869			tag++870		}871	}872	return tag873}874 875// ScaleBonds scales all bonds between atoms in the same residue with names n1, n2 to a final lenght finallengt, by moving the atoms n2.876// the o<ration is executed in place.877func ScaleBonds(coords *v3.Matrix, mol Atomer, n1, n2 string, finallenght float64) {878	for i := 0; i < mol.Len(); i++ {879		c1 := mol.Atom(i)880		if c1.Name != n1 {881			continue882		}883		for j := 0; j < mol.Len(); j++ {884			c2 := mol.Atom(j)885			if c1.MolID == c2.MolID && c1.Name == n1 && c2.Name == n2 {886				A := coords.VecView(i)887				B := coords.VecView(j)888				ScaleBond(A, B, finallenght)889			}890		}891	}892}893 894// Merges A and B in a single topology which is returned895func MergeAtomers(A, B Atomer) *Topology {896	al := A.Len()897	l := al + B.Len()898	full := make([]*Atom, l, l)899	for k, _ := range full {900		if k < al {901			full[k] = A.Atom(k)902		} else {903			full[k] = B.Atom(k - al)904		}905	}906	a, aok := A.(AtomMultiCharger)907	b, bok := B.(AtomMultiCharger)908	var charge, multi int909	if aok && bok {910		charge = a.Charge() + b.Charge()911		multi = (a.Multi() - 1 + b.Multi()) //Not TOO sure about this.912	} else {913		multi = 1914	}915	return NewTopology(charge, multi, full)916}917 918// SelCone, Given a set of cartesian points in sellist, obtains a vector ""plane"" normal to the best plane passing through the points.919// 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,920// 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.921// If one starts from one given point, 2 cgnOnes, one in each direction, are possible. If whatcone is 0, both cgnOnes are considered.922// if whatcone<0, only the cone opposite to the plane vector direction. If whatcone>0, only the cone in the plane vector direction.923// the 'initial' argument  allows the construction of a truncate cone with a radius of initial.924func SelCone(B, selection *v3.Matrix, angle, distance, thickness, initial float64, whatcone int) []int {925	A := v3.Zeros(B.NVecs())926	A.Copy(B) //We will be altering the input so its better to work with a copy.927	ar, _ := A.Dims()928	selected := make([]int, 0, 3)929	neverselected := make([]int, 0, 30000)     //waters that are too far to ever be selected930	nevercutoff := distance / math.Cos(angle)  //cutoff to be added to neverselected931	A, _, err := MassCenter(A, selection, nil) //Centrate A in the geometric center of the selection, Its easier for the following calculations932	if err != nil {933		panic(PanicMsg(err.Error()))934	}935	selection, _, _ = MassCenter(selection, selection, nil) //Centrate the selection as well936	plane, err := BestPlane(selection, nil)                 //I have NO idea which direction will this vector point. We might need its negative.937	if err != nil {938		panic(PanicMsg(err.Error()))939	}940	for i := thickness / 2; i <= distance; i += thickness {941		maxdist := math.Tan(angle)*i + initial //this should give me the radius of the cone at this point942		for j := 0; j < ar; j++ {943			if isInInt(selected, j) || isInInt(neverselected, j) { //we dont scan things that we have already selected, or are too far944				continue945			}946			atom := A.VecView(j)947			proj := Projection(atom, plane)948			norm := proj.Norm(2)949			//Now at what side of the plane is the atom?950			angle := Angle(atom, plane)951			if whatcone > 0 {952				if angle > math.Pi/2 {953					continue954				}955			} else if whatcone < 0 {956				if angle < math.Pi/2 {957					continue958				}959			}960			if norm > i+(thickness/2.0) || norm < (i-thickness/2.0) {961				continue962			}963			proj.Sub(proj, atom)964			projnorm := proj.Norm(2)965			if projnorm <= maxdist {966				selected = append(selected, j)967			}968			if projnorm >= nevercutoff {969				neverselected = append(neverselected, j)970			}971		}972	}973	return selected974}975 976func ext(s string) string {977	s2 := strings.Split(s, ""."")978	return strings.ToLower(s2[len(s2)-1])979}980 981// Attemps to open a using the file extension982// to guess the format among the supported molecule file format. Returns a Molecule983// and an error which will be non nil if the reading fails or if the extension does not984// belong to a supported format.985func MoleculeFileRead(name string) (mol *Molecule, err error) {986	switch ext(name) {987	case ""pdb"":988		mol, err = PDBFileRead(name, true)989	case ""gro"":990		mol, err = GroFileRead(name)991	case ""pdbx"":992		mol, err = PDBxFileRead(name)993	case ""cif"":994		mol, err = PDBxFileRead(name)995	case ""xyz"":996		mol, err = XYZFileRead(name)997	default:998		err = fmt.Errorf(""goChem/MoleculeFileRead: Extension %s not supported"", ext(name))999	}1000	return1001}1002 1003// Creates a map from the old position of a set of atoms, to their new position. It takes the full lenght of the molecule,1004// so the indexes not present in oripos or newpos are assigned their same original position.1005// if an atom A appears in oripos and the atom B appears in the equivalent place in newpos, a map entry A->B will be created.1006// if no entry in oripos contains B, then an additional entry B->A will be created.1007func SwitchMap(oripos, newpos []int, fulllen int) map[int]int {1008	m := make(map[int]int)1009	for i, v := range oripos {1010		m[v] = newpos[i]1011	}1012	//if not given in the list, we also need to ensure that for every a that goes to b, the b element goes to a1013	//of course you might want a go to be, b go to c, c go to a. You need to give those explicitly in that case.1014	for i, v := range newpos {1015		if _, ok := m[v]; ok {1016			continue1017		}1018		m[v] = oripos[i]1019	}1020	//Finally, whatever wasn't mentioned in neither oripos nor newpos, keeps its place1021	for i := 0; i < fulllen; i++ {1022		if _, ok := m[i]; ok {1023			continue1024		}1025		m[i] = i1026	}1027 1028	return m1029}1030 1031// Switches each atom in mol from position i to position switchmap[i] for each i in the switchmap keys.1032// Returns a modified topology. mol itself is not affected (its atoms remain in the same order) but the Indexes and IDs1033// of the atoms, are, since SwitchAtoms doesn't make copies and resets IDs and Indexes to match the atoms's positions1034// in the returned topology.1035func SwitchAtoms(switchmap map[int]int, mol Atomer) *Topology {1036	ret := NewTopology(0, 1)1037	ret.Atoms = make([]*Atom, mol.Len())1038	for i := 0; i < mol.Len(); i++ {1039		at := mol.Atom(i)1040		j, ok := switchmap[i]1041		if !ok {1042			j = i1043		}1044		ret.Atoms[j] = at1045	}1046	for i, v := range ret.Atoms {1047		if v == nil {1048			panic(fmt.Sprintf(""The indexes were such that position %d ended up empty/nil"", i))1049		}1050	}1051	ret.FillIndexes()1052	ret.ResetIDs()1053	return ret1054}1055 1056// Switches each atom in mol from position i to position switchmap[i] for each i in the switchmap keys,1057// returning the modified matrix. The original vec is not affected.1058func SwitchCoords(switchmap map[int]int, vec *v3.Matrix) *v3.Matrix {1059	ret := v3.Zeros(vec.Len())1060	finalpos := make([]int, ret.Len())1061	for i, _ := range finalpos {1062		n, ok := switchmap[i]1063		if ok {1064			finalpos[i] = n1065		} else {1066			finalpos[i] = i1067		}1068 1069	}1070	ret.SetVecs(vec, finalpos)1071	return ret1072}1073","Go"
1074"Computational Biochemistry","rmera/gochem","doc.go",".go","1076","31","/*1075 * doc.go, part of gochem.1076 *1077 * Copyright 2012 Raul Mera <rmera{at}chemDOThelsinkiDOTfi>1078 *1079 * This program is free software; you can redistribute it and/or modify1080 * it under the terms of the GNU Lesser General Public License as1081 * published by the Free Software Foundation; either version 2.1 of the1082 * License, or (at your option) any later version.1083 *1084 * This program is distributed in the hope that it will be useful,1085 * but WITHOUT ANY WARRANTY; without even the implied warranty of1086 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the1087 * GNU General Public License for more details.1088 *1089 * You should have received a copy of the GNU Lesser General1090 * Public License along with this program.  If not, see1091 * <http://www.gnu.org/licenses/>.1092 *1093 */1094 1095/*This is the main package of the goChem library.1096It provides atom and molecule structures, facilities for reading and1097writing some files used in computational chemistry and functions for1098geometric manipulations and shape, among others indicators.1099 1100See www.gochem.org for more information.1101 1102*/1103package chem1104","Go"
1105"Computational Biochemistry","rmera/gochem","chem.go",".go","26233","877","/*1106 * chem.go, part of gochem.1107 *1108 * Copyright 2012 Raul Mera <rmera{at}chemDOThelsinkiDOTfi>1109 *1110 * This program is free software; you can redistribute it and/or modify1111 * it under the terms of the GNU Lesser General Public License as1112 * published by the Free Software Foundation; either version 2.1 of the1113 * License, or (at your option) any later version.1114 *1115 * This program is distributed in the hope that it will be useful,1116 * but WITHOUT ANY WARRANTY; without even the implied warranty of1117 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the1118 * GNU General Public License for more details.1119 *1120 * You should have received a copy of the GNU Lesser General1121 * Public License along with this program.  If not, see1122 * <http://www.gnu.org/licenses/>.1123 *1124 * Gochem is developed at the laboratory for instruction in Swedish, Department of Chemistry,1125 * University of Helsinki, Finland.1126 *1127 */1128 1129package chem1130 1131import (1132	""fmt""1133	""sort""1134 1135	v3 ""github.com/rmera/gochem/v3""1136)1137 1138//import ""strings""1139 1140/* Many funcitons here panic instead of returning errors. This is because they are ""fundamental""1141 * functions. I considered that if something goes wrong here, the program is way-most likely wrong and should1142 * crash. Most panics are related to using the funciton on a nil object or trying to access out-of bounds1143 * fields1144 */1145 1146// Atom contains the information to represent an atom, except for the coordinates, which will be in a separate *v3.Matrix1147// and the b-factors, which are in a separate slice of float64.1148type Atom struct {1149	Name      string  //PDB name of the atom1150	ID        int     //The PDB index of the atom1151	index     int     //The place of the atom in a set. I won't make it accessible to ensure that it does correspond to the ordering.1152	Tag       int     //Just added this for something that someone might want to keep that is not a float.1153	MolName   string  //PDB name of the residue or molecule (3-letter code for residues)1154	MolName1  byte    //the one letter name for residues and nucleotids1155	Char16    byte    //Whatever is in the column 16 (counting from 0) in a PDB file, anything.1156	MolID     int     //PDB index of the corresponding residue or molecule1157	Chain     string  //One-character PDB name for a chain.1158	Mass      float64 //hopefully all these float64 are not too much memory1159	Occupancy float64 //a PDB crystallographic field, often used to store values of interest.1160	Vdw       float64 //radius1161	Charge    float64 //Partial charge on an atom1162	Symbol    string1163	Het       bool    // is the atom an hetatm in the pdb file? (if applicable)1164	Bonds     []*Bond //The bonds connecting the atom to others.1165}1166 1167//Atom methods1168 1169// Copy puts in the receiver a copy of A1170func (N *Atom) Copy(A *Atom) {1171	if A == nil || N == nil {1172		panic(ErrNilAtom)1173	}1174	N.Name = A.Name1175	N.ID = A.ID1176	N.Tag = A.Tag1177	N.MolName = A.MolName1178	N.MolName1 = A.MolName11179	N.MolID = A.MolID1180	N.Chain = A.Chain1181	N.Mass = A.Mass1182	N.Occupancy = A.Occupancy1183	N.Vdw = A.Vdw1184	N.Charge = A.Charge1185	N.Symbol = A.Symbol1186	N.Het = A.Het1187}1188 1189// Index returns the index of the atom1190func (N *Atom) Index() int {1191	return N.index1192}1193 1194// Index returns the index of the atom1195func (N *Atom) SetIndex(i int) {1196	N.index = i1197}1198 1199/*****Topology type***/1200 

Showing the first 1,200 of 44057 lines. Download the file for the rest.