github.com/rmera/gochem@v0.7.1/chem.go (about) 1 /* 2 * chem.go, part of gochem. 3 * 4 * Copyright 2012 Raul Mera <rmera{at}chemDOThelsinkiDOTfi> 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU Lesser General Public License as 8 * published by the Free Software Foundation; either version 2.1 of the 9 * License, or (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU Lesser General 17 * Public License along with this program. If not, see 18 * <http://www.gnu.org/licenses/>. 19 * 20 * Gochem is developed at the laboratory for instruction in Swedish, Department of Chemistry, 21 * University of Helsinki, Finland. 22 * 23 */ 24 25 package chem 26 27 import ( 28 "fmt" 29 "sort" 30 31 v3 "github.com/rmera/gochem/v3" 32 ) 33 34 //import "strings" 35 36 /* Many funcitons here panic instead of returning errors. This is because they are "fundamental" 37 * functions. I considered that if something goes wrong here, the program is way-most likely wrong and should 38 * crash. Most panics are related to using the funciton on a nil object or trying to access out-of bounds 39 * fields 40 */ 41 42 // Atom contains the information to represent an atom, except for the coordinates, which will be in a separate *v3.Matrix 43 // and the b-factors, which are in a separate slice of float64. 44 type Atom struct { 45 Name string //PDB name of the atom 46 ID int //The PDB index of the atom 47 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. 48 Tag int //Just added this for something that someone might want to keep that is not a float. 49 MolName string //PDB name of the residue or molecule (3-letter code for residues) 50 MolName1 byte //the one letter name for residues and nucleotids 51 Char16 byte //Whatever is in the column 16 (counting from 0) in a PDB file, anything. 52 MolID int //PDB index of the corresponding residue or molecule 53 Chain string //One-character PDB name for a chain. 54 Mass float64 //hopefully all these float64 are not too much memory 55 Occupancy float64 //a PDB crystallographic field, often used to store values of interest. 56 Vdw float64 //radius 57 Charge float64 //Partial charge on an atom 58 Symbol string 59 Het bool // is the atom an hetatm in the pdb file? (if applicable) 60 Bonds []*Bond //The bonds connecting the atom to others. 61 } 62 63 //Atom methods 64 65 // Copy returns a copy of the Atom object. 66 // puts the copy into the 67 func (N *Atom) Copy(A *Atom) { 68 if A == nil || N == nil { 69 panic(ErrNilAtom) 70 } 71 N.Name = A.Name 72 N.ID = A.ID 73 N.Tag = A.Tag 74 N.MolName = A.MolName 75 N.MolName1 = A.MolName1 76 N.MolID = A.MolID 77 N.Chain = A.Chain 78 N.Mass = A.Mass 79 N.Occupancy = A.Occupancy 80 N.Vdw = A.Vdw 81 N.Charge = A.Charge 82 N.Symbol = A.Symbol 83 N.Het = A.Het 84 } 85 86 // Index returns the index of the atom 87 func (N *Atom) Index() int { 88 return N.index 89 } 90 91 // Index returns the index of the atom 92 func (N *Atom) SetIndex(i int) { 93 N.index = i 94 } 95 96 /*****Topology type***/ 97 98 // Topology contains information about a molecule which is not expected to change in time (i.e. everything except for coordinates and b-factors) 99 type Topology struct { 100 Atoms []*Atom 101 charge int 102 multi int 103 } 104 105 // NewTopology returns topology with ats atoms, 106 // charge charge and multi multiplicity. 107 // It doesnt check for consitency across slices, correct charge 108 // or unpaired electrons. 109 func NewTopology(charge, multi int, ats ...[]*Atom) *Topology { 110 top := new(Topology) 111 if len(ats) == 0 || ats[0] == nil { 112 top.Atoms = make([]*Atom, 0, 0) //return nil, fmt.Errorf("Supplied a nil Topology") 113 } else { 114 top.Atoms = ats[0] 115 } 116 top.charge = charge 117 top.multi = multi 118 return top 119 } 120 121 /*Topology methods*/ 122 123 // Charge returns the total charge of the topology 124 func (T *Topology) Charge() int { 125 return T.charge 126 } 127 128 // Multi returns the multiplicity in the topology 129 func (T *Topology) Multi() int { 130 return T.multi 131 } 132 133 // SetCharge sets the total charge of the topology to i 134 func (T *Topology) SetCharge(i int) { 135 T.charge = i 136 } 137 138 // SetMulti sets the multiplicity in the topology to i 139 func (T *Topology) SetMulti(i int) { 140 T.multi = i 141 } 142 143 // FillMasses tries to get fill the masses for atom that don't have one 144 // by getting it from the symbol. Only a few common elements are supported 145 func (T *Topology) FillMasses() { 146 for _, val := range T.Atoms { 147 if val.Symbol != "" && val.Mass == 0 { 148 val.Mass = symbolMass[val.Symbol] //Not error checking 149 } 150 } 151 } 152 153 // FillsIndexes sets the Index value of each atom to that cooresponding to its 154 // place in the molecule. 155 func (T *Topology) FillIndexes() { 156 for key, val := range T.Atoms { 157 val.index = key 158 } 159 160 } 161 162 // FillVdw tries to get fill the van der Waals radii for the atoms in the molecule 163 // from a symbol->radii map. Only a few common elements are supported 164 func (T *Topology) FillVdw() { 165 for _, val := range T.Atoms { 166 if val.Symbol != "" && val.Vdw == 0 { 167 val.Vdw = symbolVdwrad[val.Symbol] //Not error checking 168 //I'm not super sure about this. 169 } 170 } 171 } 172 173 // ResetIDs sets the current order of atoms as ID and the order of molecules as 174 // MolID for all atoms 175 func (T *Topology) ResetIDs() { 176 currid := 1 177 currid2 := 1 178 for key, val := range T.Atoms { 179 T.Atoms[key].ID = key + 1 180 if currid == val.MolID { 181 continue 182 } 183 if currid == val.MolID-1 { //We hit a new molecule 184 currid2++ 185 currid++ 186 continue 187 } 188 //change of residue after fixing one that didnt match position 189 if currid2 != val.MolID { 190 currid2 = T.Atoms[key].MolID 191 T.Atoms[key].MolID = currid + 1 192 currid = currid + 1 193 continue 194 } 195 //A residue's index doesnt match its position 196 T.Atoms[key].MolID = currid 197 198 } 199 } 200 201 // CopyAtoms copies the atoms form A into the receiver topology. This is a deep copy, so the receiver must have 202 // at least as many atoms as A. 203 func (T *Topology) CopyAtoms(A Atomer) { 204 //T := new(Topology) 205 if len(T.Atoms) < A.Len() { 206 T.Atoms = make([]*Atom, A.Len()) 207 for i, _ := range T.Atoms { 208 T.Atoms[i] = new(Atom) 209 } 210 } 211 for key := 0; key < A.Len(); key++ { 212 T.Atoms[key].Copy(A.Atom(key)) 213 } 214 } 215 216 // Atom returns the Atom corresponding to the index i 217 // of the Atom slice in the Topology. Panics if 218 // out of range. 219 func (T *Topology) Atom(i int) *Atom { 220 if i >= T.Len() { 221 panic(ErrAtomOutOfRange) 222 } 223 return T.Atoms[i] 224 } 225 226 // SetAtom sets the (i+1)th Atom of the topology to aT. 227 // Panics if out of range 228 func (T *Topology) SetAtom(i int, at *Atom) { 229 if i >= T.Len() { 230 panic(ErrAtomOutOfRange) 231 } 232 T.Atoms[i] = at 233 } 234 235 // AppendAtom appends an atom at the end of the reference 236 func (T *Topology) AppendAtom(at *Atom) { 237 /*newmol, ok := T.CopyAtoms().(*Topology) 238 if !ok { 239 panic("cant happen") 240 } 241 newmol.Atoms = append(newmol.Atoms, at)*/ 242 T.Atoms = append(T.Atoms, at) 243 } 244 245 // SelectAtoms puts the subset of atoms in T that have 246 // indexes in atomlist into the receiver. Panics if problem. 247 func (R *Topology) SomeAtoms(T Atomer, atomlist []int) { 248 var ret []*Atom 249 lenatoms := T.Len() 250 for k, j := range atomlist { 251 if j > lenatoms-1 { 252 err := fmt.Sprintf("goChem: Atom requested (Number: %d, value: %d) out of range", k, j) 253 panic(PanicMsg(err)) 254 } 255 ret = append(ret, T.Atom(j)) 256 } 257 R.Atoms = ret 258 } 259 260 // SelectAtomsSafe puts the atoms of T 261 // with indexes in atomlist into the receiver. Returns error if problem. 262 func (R *Topology) SomeAtomsSafe(T Atomer, atomlist []int) error { 263 f := func() { R.SomeAtoms(T, atomlist) } 264 return gnMaybe(gnPanicker(f)) 265 } 266 267 // DelAtom Deletes atom i by reslicing. 268 // This means that the copy still uses as much memory as the original T. 269 func (T *Topology) DelAtom(i int) { 270 if i >= T.Len() { 271 panic(ErrAtomOutOfRange) 272 } 273 if i == T.Len()-1 { 274 T.Atoms = T.Atoms[:i] 275 } else { 276 T.Atoms = append(T.Atoms[:i], T.Atoms[i+1:]...) 277 } 278 } 279 280 // Len returns the number of atoms in the topology. 281 func (T *Topology) Len() int { 282 //if T.Atoms is nil, return len(T.Atoms) will panic, so I will let that happen for now. 283 // if T.Atoms == nil { 284 // panic(ErrNilAtoms) 285 // } 286 return len(T.Atoms) 287 } 288 289 // Masses returns a slice of float64 with the masses of the atoms in the topology, or nil and an error if they have not been calculated 290 func (T *Topology) Masses() ([]float64, error) { 291 mass := make([]float64, T.Len()) 292 for i := 0; i < T.Len(); i++ { 293 thisatom := T.Atom(i) 294 if thisatom.Mass == 0 { 295 return nil, CError{fmt.Sprintf("goChem: Not all the masses have been obtained: %d %v", i, thisatom), []string{"Topology.Masses"}} 296 } 297 mass[i] = thisatom.Mass 298 } 299 return mass, nil 300 } 301 302 // AssignBonds assigns bonds to a molecule based on a simple distance 303 // criterium, similar to that described in DOI:10.1186/1758-2946-3-33 304 func (T *Topology) AssignBonds(coord *v3.Matrix) error { 305 // might get slow for 306 //large systems. It's really not thought 307 //for proteins or macromolecules. 308 //For this reason, this method is not called automatically when building a new topology. 309 //Well, that and that it requires a Matrix object, which would mean changing the 310 //signature of the NewTopology function. 311 var t1, t2 *v3.Matrix 312 var at1, at2 *Atom 313 T.FillIndexes() 314 t3 := v3.Zeros(1) 315 bonds := make([]*Bond, 0, 10) 316 tot := T.Len() 317 var nextIndex int 318 for i := 0; i < tot; i++ { 319 t1 = coord.VecView(i) 320 at1 = T.Atoms[i] 321 cov1 := symbolCovrad[at1.Symbol] 322 if cov1 == 0 { 323 err := new(CError) 324 err.msg = fmt.Sprintf("Couldn't find the covalent radii for %s %d", at1.Symbol, i) 325 err.Decorate("AssignBonds") 326 return err 327 } 328 for j := i + 1; j < tot; j++ { 329 t2 = coord.VecView(j) 330 at2 = T.Atoms[j] 331 cov2 := symbolCovrad[at2.Symbol] 332 if cov2 == 0 { 333 err := new(CError) 334 err.msg = fmt.Sprintf("Couldn't find the covalent radii for %s %d", at2.Symbol, j) 335 err.Decorate("AssignBonds") 336 return err 337 } 338 339 t3.Sub(t2, t1) 340 d := t3.Norm(2) 341 if d < cov1+cov2+bondtol && d > tooclose { 342 b := &Bond{Index: nextIndex, Dist: d, At1: at1, At2: at2} 343 at1.Bonds = append(at1.Bonds, b) 344 at2.Bonds = append(at2.Bonds, b) 345 bonds = append(bonds, b) //just to easily keep track of them. 346 nextIndex++ 347 } 348 349 } 350 } 351 352 //Now we check that no atom has too many bonds. 353 for i := 0; i < tot; i++ { 354 at := T.Atoms[i] 355 max := symbolMaxBonds[at.Symbol] 356 if max == 0 { //means there is not a specified number of bonds for this atom. 357 continue 358 } 359 sort.Slice(at.Bonds, func(i, j int) bool { return at.Bonds[i].Dist < at.Bonds[j].Dist }) 360 //I am hoping this will remove bonds until len(at.Bonds) is not 361 //greater than max. 362 for i := len(at.Bonds); i > max; i = len(at.Bonds) { 363 err := at.Bonds[len(at.Bonds)-1].Remove() //we remove the longest bond 364 if err != nil { 365 return errDecorate(err, "AssignBonds") 366 } 367 } 368 369 } 370 371 return nil 372 } 373 374 /**Type Molecule**/ 375 376 // Molecule contains all the info for a molecule in many states. The info that is expected to change between states, 377 // Coordinates and b-factors are stored separately from other atomic info. 378 type Molecule struct { 379 *Topology 380 Coords []*v3.Matrix 381 Bfactors [][]float64 382 XYZFileData []string //This can be anything. The main rationale for including it is that XYZ files have a "comment" 383 //line after the first one. This line is sometimes used to write the energy of the structure. 384 //So here the line can be kept for each XYZ frame, and parse later 385 current int 386 } 387 388 // NewMolecule makes a molecule with ats atoms, coords coordinates, bfactors b-factors 389 // charge charge and unpaired unpaired electrons, and returns it. It doesnt check for 390 // consitency across slices or correct charge or unpaired electrons. 391 func NewMolecule(coords []*v3.Matrix, ats Atomer, bfactors [][]float64) (*Molecule, error) { 392 if ats == nil { 393 return nil, CError{"Supplied a nil Reference", []string{"NewMolCule"}} 394 } 395 if coords == nil { 396 return nil, CError{"Supplied a nil Coord slice", []string{"NewMolCule"}} 397 } 398 // if bfactors == nil { 399 // return nil, fmt.Errorf("Supplied a nil Bfactors slice") 400 // } 401 mol := new(Molecule) 402 atcopy := func() { 403 mol.Topology = NewTopology(9999, -1, make([]*Atom, 0, ats.Len())) //I use 9999 for charge and -1 or multi to indicate that they are not truly set. So far NewTopology never actually returns any error so it's safe to ignore them. 404 for i := 0; i < ats.Len(); i++ { 405 mol.Atoms = append(mol.Atoms, ats.Atom(i)) 406 } 407 } 408 switch ats := ats.(type) { //for speed 409 case *Topology: 410 mol.Topology = ats 411 case AtomMultiCharger: 412 atcopy() 413 mol.SetMulti(ats.Multi()) 414 mol.SetCharge(ats.Charge()) 415 default: 416 atcopy() 417 } 418 mol.Coords = coords 419 mol.Bfactors = bfactors 420 return mol, nil 421 } 422 423 //The molecule methods: 424 425 // DelCoord deletes the coodinate i from every frame of the molecule. 426 func (M *Molecule) DelCoord(i int) error { 427 //note: Maybe this shouldn't be exported. Unexporting it could be a reasonable API change. 428 r, _ := M.Coords[0].Dims() 429 var err error 430 for j := 0; j < len(M.Coords); j++ { 431 tmp := v3.Zeros(r - 1) 432 tmp.DelVec(M.Coords[j], i) 433 M.Coords[j] = tmp 434 if err != nil { 435 return err 436 } 437 } 438 return nil 439 } 440 441 // Del Deletes atom i and its coordinates from the molecule. 442 func (M *Molecule) Del(i int) error { 443 if i >= M.Len() { 444 panic(ErrAtomOutOfRange) 445 } 446 M.DelAtom(i) 447 err := M.DelCoord(i) 448 return err 449 } 450 451 // Copy puts in the receiver a copy of the molecule A including coordinates 452 func (M *Molecule) Copy(A *Molecule) { 453 if err := A.Corrupted(); err != nil { 454 panic(err.Error()) 455 } 456 r, _ := A.Coords[0].Dims() 457 //mol := new(Molecule) 458 M.Topology = new(Topology) 459 for i := 0; i < A.Len(); i++ { 460 at := new(Atom) 461 at.Copy(A.Atom(i)) 462 M.Topology.Atoms = append(M.Topology.Atoms, at) 463 464 } 465 // M.CopyAtoms(A) 466 M.Coords = make([]*v3.Matrix, 0, len(A.Coords)) 467 M.Bfactors = make([][]float64, 0, len(A.Bfactors)) 468 for key, val := range A.Coords { 469 tmp := v3.Zeros(r) 470 tmp.Copy(val) 471 M.Coords = append(M.Coords, tmp) 472 tmp2 := copyB(A.Bfactors[key]) 473 M.Bfactors = append(M.Bfactors, tmp2) 474 } 475 if err := M.Corrupted(); err != nil { 476 panic(PanicMsg(fmt.Sprintf("goChem: Molecule creation error: %s", err.Error()))) 477 } 478 } 479 480 func copyB(b []float64) []float64 { 481 r := make([]float64, len(b), len(b)) 482 for k, v := range b { 483 r[k] = v 484 } 485 return r 486 } 487 488 // AddFrame akes a matrix of coordinates and appends them at the end of the Coords. 489 // It checks that the number of coordinates matches the number of atoms. 490 func (M *Molecule) AddFrame(newframe *v3.Matrix) { 491 if newframe == nil { 492 panic(ErrNilFrame) 493 } 494 r, c := newframe.Dims() 495 if c != 3 { 496 panic(ErrNotXx3Matrix) 497 } 498 if M.Len() != r { 499 panic(PanicMsg(fmt.Sprintf("goChem: Wrong number of coordinates (%d)", newframe.NVecs()))) 500 } 501 if M.Coords == nil { 502 M.Coords = make([]*v3.Matrix, 1, 1) 503 } 504 M.Coords = append(M.Coords, newframe) 505 } 506 507 // AddManyFrames adds the array of matrices newfames to the molecule. It checks that 508 // the number of coordinates matches the number of atoms. 509 func (M *Molecule) AddManyFrames(newframes []*v3.Matrix) { 510 if newframes == nil { 511 panic(ErrNilFrame) 512 } 513 if M.Coords == nil { 514 M.Coords = make([]*v3.Matrix, 1, len(newframes)) 515 } 516 for key, val := range newframes { 517 f := func() { M.AddFrame(val) } 518 err := gnMaybe(gnPanicker(f)) 519 if err != nil { 520 panic(PanicMsg(fmt.Sprintf("goChem: %s in frame %d", err.Error(), key))) 521 } 522 } 523 } 524 525 // Coord returns the coords for the atom atom in the frame frame. 526 // panics if frame or coords are out of range. 527 func (M *Molecule) Coord(atom, frame int) *v3.Matrix { 528 if frame >= len(M.Coords) { 529 panic(PanicMsg(fmt.Sprintf("goChem: Frame requested (%d) out of range", frame))) 530 } 531 r, _ := M.Coords[frame].Dims() 532 if atom >= r { 533 panic(PanicMsg(fmt.Sprintf("goChem: Requested coordinate (%d) out of bounds (%d)", atom, M.Coords[frame].NVecs()))) 534 } 535 ret := v3.Zeros(1) 536 empt := M.Coords[frame].VecView(atom) 537 ret.Copy(empt) 538 return ret 539 } 540 541 // Current returns the number of the next read frame 542 func (M *Molecule) Current() int { 543 if M == nil { 544 return -1 545 } 546 return M.current 547 } 548 549 // SetCurrent sets the value of the frame nex to be read 550 // to i. 551 func (M *Molecule) SetCurrent(i int) { 552 if i < 0 || i >= len(M.Coords) { 553 panic(PanicMsg(fmt.Sprintf("goChem: Invalid new value for current frame: %d Current frames: %d", i, len(M.Coords)))) 554 } 555 M.current = i 556 } 557 558 /* 559 //SetCoords replaces the coordinates of atoms in the positions given by atomlist with the gnOnes in NewVecs (in order) 560 //If atomlist contains a single element, it replaces as many coordinates as given in NewVecs, starting 561 //at the element in atomlist. In the latter case, the function checks that there are enough coordinates to 562 //replace and returns an error if not. 563 func (M *Molecule) SetCoords(NewVecs *CoordMA, atomlist []int, frame int) { 564 if frame >= len(M.Coords) { 565 panic(fmt.Sprintf("Frame (%d) out of range!", frame)) 566 } 567 //If supplies a list with one number, the NewVecs will replace the old coords 568 //Starting that number. We do check that you don't put more coords than spaces we have. 569 if len(atomlist) == 1 { 570 if NewVecs.Rows() > M.Coords[frame].Rows()-atomlist[0]-1 { 571 panic(fmt.Sprintf("Cant replace starting from position %d: Not enough atoms in molecule", atomlist[0])) 572 } 573 M.Coords[frame].SetMatrix(atomlist[0], 0, NewVecs) 574 return 575 } 576 //If the list has more than one atom 577 lenatoms := M.Len() 578 for k, j := range atomlist { 579 if j > lenatoms-1 { 580 panic(fmt.Sprintf("Requested position number: %d (%d) out of range", k, j)) 581 } 582 M.Coords[frame].SetMatrix(j, 0, NewVecs.GetRowVector(k)) 583 } 584 } 585 586 */ 587 588 // Corrupted checks whether the molecule is corrupted, i.e. the 589 // coordinates don't match the number of atoms. It also checks 590 // That the coordinate matrices have 3 columns. 591 func (M *Molecule) Corrupted() error { 592 var err error 593 if M.Bfactors == nil { 594 M.Bfactors = make([][]float64, 0, len(M.Coords)) 595 M.Bfactors = append(M.Bfactors, make([]float64, M.Len())) 596 } 597 lastbfac := len(M.Bfactors) - 1 598 for i := range M.Coords { 599 r, c := M.Coords[i].Dims() 600 if M.Len() != r || c != 3 { 601 err = CError{fmt.Sprintf("Inconsistent coordinates/atoms in frame %d: Atoms %d, coords: %d", i, M.Len(), M.Coords[i].NVecs()), []string{"Molecule.Corrupted"}} 602 break 603 } 604 //Since bfactors are not as important as coordinates, we will just fill with 605 //zeroes anything that is lacking or incomplete instead of returning an error. 606 607 if lastbfac < i { 608 bfacs := make([]float64, M.Len()) 609 M.Bfactors = append(M.Bfactors, bfacs) 610 } 611 bfr := len(M.Bfactors[i]) 612 if bfr < M.Len() { 613 M.Bfactors[i] = make([]float64, M.Len(), 1) 614 } 615 } 616 return err 617 } 618 619 // NFrames returns the number of frames in the molecule 620 func (M *Molecule) NFrames() int { 621 return len(M.Coords) 622 } 623 624 //Implementaiton of the sort.Interface 625 626 // Swap function, as demanded by sort.Interface. It swaps atoms, coordinates 627 // (all frames) and bfactors of the molecule. 628 func (M *Molecule) Swap(i, j int) { 629 M.Atoms[i], M.Atoms[j] = M.Atoms[j], M.Atoms[i] 630 for k := 0; k < len(M.Coords); k++ { 631 M.Coords[k].SwapVecs(i, j) 632 t1 := M.Bfactors[k][i] 633 t2 := M.Bfactors[k][j] 634 M.Bfactors[k][i] = t2 635 M.Bfactors[k][j] = t1 636 } 637 } 638 639 // Less returns true if the value in the Bfactors for 640 // atom i are less than that for atom j, and false otherwise. 641 func (M *Molecule) Less(i, j int) bool { 642 return M.Bfactors[0][i] < M.Bfactors[0][j] 643 } 644 645 //Len is implemented in Topology 646 647 //End sort.Interface 648 649 /****************************************** 650 //The following implement the Traj interface 651 **********************************************/ 652 653 // Checks that the molecule exists and has some existent 654 // Coordinates, in which case returns true. 655 // It returns false otherwise. 656 // The coordinates could still be empty 657 func (M *Molecule) Readable() bool { 658 if M != nil || M.Coords != nil || M.current < len(M.Coords) { 659 660 return true 661 } 662 return false 663 } 664 665 // Next puts the next frame into V and returns an error or nil 666 // The box argument is never used. 667 func (M *Molecule) Next(V *v3.Matrix, box ...[]float64) error { 668 if M.current >= len(M.Coords) { 669 return newlastFrameError("", len(M.Coords)-1) 670 } 671 // fmt.Println("CURR", M.current, len(M.Coords), V.NVecs(), M.Coords[M.current].NVecs()) //////////////// 672 M.current++ 673 if V == nil { 674 return nil 675 } 676 V.Copy(M.Coords[M.current-1]) 677 return nil 678 } 679 680 // InitRead initializes molecule to be read as a traj (not tested!) 681 func (M *Molecule) InitRead() error { 682 if M == nil || len(M.Coords) == 0 { 683 return CError{"Bad molecule", []string{"InitRead"}} 684 } 685 M.current = 0 686 return nil 687 } 688 689 // NextConc takes a slice of bools and reads as many frames as elements the list has 690 // form the trajectory. The frames are discarted if the corresponding elemetn of the slice 691 // is false. The function returns a slice of channels through each of each of which 692 // a *matrix.DenseMatrix will be transmited 693 func (M *Molecule) NextConc(frames []*v3.Matrix) ([]chan *v3.Matrix, error) { 694 toreturn := make([]chan *v3.Matrix, 0, len(frames)) 695 used := false 696 697 for _, val := range frames { 698 if val == nil { 699 M.current++ 700 toreturn = append(toreturn, nil) 701 continue 702 } 703 if M.current >= len(M.Coords) { 704 lastframe := newlastFrameError("", len(M.Coords)-1) 705 if used == false { 706 return nil, lastframe 707 } else { 708 return toreturn, lastframe 709 } 710 } 711 used = true 712 toreturn = append(toreturn, make(chan *v3.Matrix)) 713 go func(a *v3.Matrix, pipe chan *v3.Matrix) { 714 pipe <- a 715 }(M.Coords[M.current], toreturn[len(toreturn)-1]) 716 M.current++ 717 } 718 return toreturn, nil 719 } 720 721 // Close just sets the "current" counter to 0. 722 // If you are using it as a trajectory, you can always just discard the molecule 723 // and let the CG take care of it, as there is nothing on disk linked to it.. 724 func (M *Molecule) Close() { 725 M.current = 0 726 } 727 728 /**End Traj interface implementation***********/ 729 730 //End Molecule methods 731 732 //Traj Error 733 734 type lastFrameError struct { 735 fileName string 736 frame int 737 deco []string 738 } 739 740 // Error returns an error message string. 741 func (E *lastFrameError) Error() string { 742 return "EOF" //: Last frame in mol-based trajectory from file %10s reached at frame %10d", E.fileName, E.frame) 743 } 744 745 // Format returns the format used by the trajectory that returned the error. 746 func (E *lastFrameError) Format() string { 747 return "mol" 748 } 749 750 // Frame returns the frame at which the error was detected. 751 func (E *lastFrameError) Frame() int { 752 return E.frame 753 } 754 755 func (E *lastFrameError) Critical() bool { 756 return false 757 } 758 759 // FileName returns the name of the file from where the trajectory that gave the error is read. 760 func (E *lastFrameError) FileName() string { 761 return E.fileName 762 } 763 764 // NormalLastFrameTermination does nothing, it is there so we can have an interface unifying all 765 // "normal termination" errors so they can be filtered out by type switch. 766 func (E *lastFrameError) NormalLastFrameTermination() { 767 } 768 769 // Decorate will add the dec string to the decoration slice of strings of the error, 770 // and return the resulting slice. 771 func (E *lastFrameError) Decorate(dec string) []string { 772 if dec == "" { 773 return E.deco 774 } 775 E.deco = append(E.deco, dec) 776 return E.deco 777 } 778 779 func newlastFrameError(filename string, frame int) *lastFrameError { 780 e := new(lastFrameError) 781 e.fileName = filename 782 e.frame = frame 783 return e 784 } 785 786 //End Traj Error 787 788 //The general concrete error type for the package 789 790 // CError (Concrete Error) is the concrete error type 791 // for the chem package, that implements chem.Error 792 type CError struct { 793 msg string 794 deco []string 795 } 796 797 //func (err CError) NonCritical() bool { return err.noncritical } 798 799 func (err CError) Error() string { return err.msg } 800 801 // Decorate will add the dec string to the decoration slice of strings of the error, 802 // and return the resulting slice. 803 func (err CError) Decorate(dec string) []string { 804 if dec == "" { 805 return err.deco 806 } 807 err.deco = append(err.deco, dec) 808 return err.deco 809 } 810 811 // errDecorate will decorate a chem.Error, or use the message of the error 812 // and the name of the called to produce a new chem.Error (if the original error does not 813 // implement chem.Error 814 func errDecorate(err error, caller string) Error { 815 if err == nil { 816 return nil 817 } 818 err2, ok := err.(Error) //I know that is the type returned byt initRead 819 if ok { 820 err2.Decorate(caller) 821 return err2 822 } 823 err3 := CError{err.Error(), []string{caller}} 824 return err3 825 } 826 827 // PanicMsg is the type used for all the panics raised in the chem package 828 // so they can be easily recovered if needed. goChem only raises panics 829 // for programming errors: Attempts to to out of a matrice's bounds, 830 // dimension mismatches in matrices, etc. 831 type PanicMsg string 832 833 // Error returns a string with an error message 834 func (v PanicMsg) Error() string { return string(v) } 835 836 const ( 837 ErrNilData = PanicMsg("goChem: Nil data given ") 838 ErrInconsistentData = PanicMsg("goChem: Inconsistent data length ") 839 ErrNilMatrix = PanicMsg("goChem: Attempted to access nil v3.Matrix or gonum/mat64.Dense") 840 ErrNilAtoms = PanicMsg("goChem: Topology has a nil []*Atom slice") 841 ErrNilAtom = PanicMsg("goChem: Attempted to copy from or to a nil Atom") 842 ErrAtomOutOfRange = PanicMsg("goChem: Requested/Attempted setting Atom out of range") 843 ErrNilFrame = PanicMsg("goChem: Attempted to acces nil frame") 844 ErrNotXx3Matrix = PanicMsg("goChem: A v3.Matrix should have 3 columns") 845 ErrCliffordRotation = PanicMsg("goChem-Clifford: Target and Result matrices must have the same dimensions. They cannot reference the same matrix") //the only panic that the Clifford functions throw. 846 )