github.com/rmera/gochem@v0.7.1/gochem_test.go (about) 1 /* 2 * gochem_test.go 3 * 4 * Copyright 2013 <rmera@Holmes> 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation; either version 2 of the License, or 9 * (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 General Public License 17 * along with this program; if not, write to the Free Software 18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 19 * MA 02110-1301, USA. 20 * 21 * 22 */ 23 24 package chem 25 26 import ( 27 "fmt" 28 "os" 29 "runtime" 30 "testing" 31 32 v3 "github.com/rmera/gochem/v3" 33 "gonum.org/v1/gonum/mat" 34 ) 35 36 //import "runtime" 37 38 func TestSymbolError(Te *testing.T) { 39 _, err := PDBFileRead("test/SymbolErrorTest.pdb") 40 if err != nil { 41 Te.Error(err) 42 } 43 } 44 45 func TestGROIO(Te *testing.T) { 46 mol, err := GroFileRead("test/test.gro") 47 if err != nil { 48 Te.Error(err) 49 } 50 fmt.Println("NO WEI", mol.Len(), len(mol.Coords)) 51 fmt.Println(mol.Atom(3), mol.Coords[0].VecView(3)) 52 err = PDBFileWrite("test/testgro.pdb", mol.Coords[0], mol, nil) 53 if err != nil { 54 Te.Error(err) 55 } 56 err = GroFileWrite("test/testgro.gro", mol.Coords, mol) 57 if err != nil { 58 Te.Error(err) 59 } 60 61 } 62 63 //TestMultiXYZ tests that multi-XYZ files are opened and read correctly. 64 func TestXYZIO(Te *testing.T) { 65 mol, err := XYZFileRead("test/sample.xyz") 66 if err != nil { 67 fmt.Println("There was an error!", err.Error()) 68 Te.Error(err) 69 } 70 fmt.Println("XYZ read!") 71 XYZFileWrite("test/sampleFirst.xyz", mol.Coords[0], mol) 72 } 73 74 func TestPDBIO(Te *testing.T) { 75 mol, err := PDBFileRead("test/2c9v.pdb", true) 76 if err != nil { 77 Te.Error(err) 78 } 79 fmt.Println("NO WEI") 80 err = PDBFileWrite("test/2c9vIO.pdb", mol.Coords[0], mol, mol.Bfactors[0]) 81 if err != nil { 82 Te.Error(err) 83 } 84 //for the 3 residue I should get -131.99, 152.49. 85 } 86 87 //TestChangeAxis reads the PDB 2c9v.pdb from the test directory, collects 88 //The CA and CB of residue D124 of the chain A, and uses Clifford algebra to rotate the 89 //whole molecule such as the vector defined by these 2 atoms is 90 //aligned with the Z axis. The new molecule is written 91 //as 2c9v_aligned.pdb to the test folder. 92 func TestChangeAxis(Te *testing.T) { 93 //runtime.GOMAXPROCS(2) /////////////////////////// 94 mol, err := PDBFileRead("test/2c9v.pdb", true) 95 if err != nil { 96 Te.Error(err) 97 } 98 PDBFileWrite("test/2c9v-Readtest.pdb", mol.Coords[0], mol, nil) 99 //The selection thing 100 orient_atoms := [2]int{0, 0} 101 for index := 0; index < mol.Len(); index++ { 102 atom := mol.Atom(index) 103 if atom.Chain == "A" && atom.MolID == 124 { 104 if atom.Name == "CA" { 105 orient_atoms[0] = index 106 } else if atom.Name == "CB" { 107 orient_atoms[1] = index 108 } 109 } 110 } 111 //Get the axis of rotation 112 //ov1:=mol.Coord(orient_atoms[0], 0) 113 ov2 := mol.Coord(orient_atoms[1], 0) 114 //now we center the thing in the beta carbon of D124 115 mol.Coords[0].SubVec(mol.Coords[0], ov2) 116 PDBFileWrite("test/2c9v-translated.pdb", mol.Coords[0], mol, nil) 117 //Now the rotation 118 ov1 := mol.Coord(orient_atoms[0], 0) //make sure we have the correct versions 119 ov2 = mol.Coord(orient_atoms[1], 0) //same 120 orient := v3.Zeros(ov2.NVecs()) 121 orient.Sub(ov2, ov1) 122 // PDBWrite(mol,"test/2c9v-124centered.pdb") 123 Z, _ := v3.NewMatrix([]float64{0, 0, 1}) 124 axis := cross(orient, Z) 125 angle := Angle(orient, Z) 126 oldcoords := v3.Zeros(mol.Coords[0].NVecs()) 127 oldcoords.Copy(mol.Coords[0]) 128 mol.Coords[0] = Rotate(oldcoords, mol.Coords[0], axis, angle) 129 if err != nil { 130 Te.Error(err) 131 } 132 PDBFileWrite("test/2c9v-aligned.pdb", mol.Coords[0], mol, nil) 133 fmt.Println("bench1") 134 } 135 136 func TestMolidNameChain2Index(Te *testing.T) { 137 //runtime.GOMAXPROCS(2) /////////////////////////// 138 mol, err := PDBFileRead("test/2c9v.pdb", true) 139 if err != nil { 140 Te.Error(err) 141 } 142 index, err := MolIDNameChain2Index(mol, 46, "ND1", "A") 143 if err != nil { 144 Te.Error(err) 145 } 146 fmt.Println("this should print the index and the Atom object for H46, chain A: ", index, mol.Atom(index)) 147 148 } 149 150 //TestOldChangeAxis reads the PDB 2c9v.pdb from the test directory, collects 151 //The CA and CB of residue D124 of the chain A, and rotates the 152 //whole molecule such as the vector defined by these 2 atoms is 153 //aligned with the Z axis. The new molecule is written 154 //as 2c9v_aligned.pdb to the test folder. 155 func TestOldChangeAxis(Te *testing.T) { 156 viej, _ := os.Open("test/2c9v.pdb") 157 mol, err := PDBRead(viej, true) 158 if err != nil { 159 Te.Error(err) 160 } 161 orient_atoms := [2]int{0, 0} 162 for index := 0; index < mol.Len(); index++ { 163 atom := mol.Atom(index) 164 if atom.Chain == "A" && atom.MolID == 124 { 165 if atom.Name == "CA" { 166 orient_atoms[0] = index 167 } else if atom.Name == "CB" { 168 orient_atoms[1] = index 169 } 170 } 171 } 172 ov1 := mol.Coord(orient_atoms[0], 0) 173 ov2 := mol.Coord(orient_atoms[1], 0) 174 //now we center the thing in the beta carbon of D124 175 mol.Coords[0].SubVec(mol.Coords[0], ov2) 176 //Now the rotation 177 ov1 = mol.Coord(orient_atoms[0], 0) //make sure we have the correct versions 178 ov2 = mol.Coord(orient_atoms[1], 0) //same 179 orient := v3.Zeros(ov2.NVecs()) 180 orient.Sub(ov2, ov1) 181 rotation := RotatorToNewZ(orient) 182 cr, cc := mol.Coords[0].Dims() 183 fmt.Println("rotation: ", rotation, cr, cc) //////////////////////////////////////////////////////// 184 mol.Coords[0].Mul(mol.Coords[0], rotation) 185 // fmt.Println(orient_atoms[1], mol.Atom(orient_atoms[1]),mol.Atom(orient_atoms[0])) 186 if err != nil { 187 Te.Error(err) 188 } 189 PDBFileWrite("test/2c9v-old-aligned.pdb", mol.Coords[0], mol, nil) 190 fmt.Println("bench2") 191 } 192 193 //Aligns the main plane of a molecule with the XY-plane. 194 //Here XYZRead and XYZWrite are tested 195 func TestPutInXYPlane(Te *testing.T) { 196 myxyz, _ := os.Open("test/sample_plane.xyz") 197 mol, err := XYZRead(myxyz) 198 if err != nil { 199 Te.Error(err) 200 } 201 indexes := []int{0, 1, 2, 3, 23, 22, 21, 20, 25, 44, 39, 40, 41, 42, 61, 60, 59, 58, 63, 5} 202 some := v3.Zeros(len(indexes)) 203 some.SomeVecs(mol.Coords[0], indexes) 204 //for most rotation things it is good to have the molecule centered on its mean. 205 mol.Coords[0], _, _ = MassCenter(mol.Coords[0], some, nil) 206 //The test molecule is not completely planar so we use a subset of atoms that are contained in a plane 207 //These are the atoms given in the indexes slice. 208 some.SomeVecs(mol.Coords[0], indexes) 209 //The strategy is: Take the normal to the plane of the molecule (os molecular subset), and rotate it until it matches the Z-axis 210 //This will mean that the plane of the molecule will now match the XY-plane. 211 best, err := BestPlane(some, nil) 212 if err != nil { 213 err2 := err.(Error) 214 fmt.Println(err2.Decorate("")) 215 Te.Errorf(err.Error()) 216 // panic(err.Error()) 217 } 218 z, _ := v3.NewMatrix([]float64{0, 0, 1}) 219 zero, _ := v3.NewMatrix([]float64{0, 0, 0}) 220 fmt.Println("Best Plane", best, z) 221 axis := v3.Zeros(1) 222 axis.Cross(best, z) 223 fmt.Println("axis", axis) 224 //The main part of the program, where the rotation actually happens. Note that we rotate the whole 225 //molecule, not just the planar subset, this is only used to calculate the rotation angle. 226 // fmt.Println("DATA", mol.Coords[0], zero, axis, Angle(best, z)) 227 mol.Coords[0], err = RotateAbout(mol.Coords[0], zero, axis, Angle(best, z)) 228 if err != nil { 229 Te.Error(err) 230 } 231 // fmt.Println("after!", mol.Coords[0], err) 232 //Now we write the rotated result. 233 outxyz, _ := os.Create("test/Rotated.xyz") //This is the XYZWrite written file 234 XYZWrite(outxyz, mol.Coords[0], mol) 235 } 236 237 func TestDelete(Te *testing.T) { 238 mol, err := XYZFileRead("test/ethanol.xyz") 239 if err != nil { 240 Te.Error(err) 241 } 242 fmt.Println("Calling with 8") 243 mol.Del(8) 244 XYZFileWrite("test/ethanolDel8.xyz", mol.Coords[0], mol) 245 mol2, err := XYZFileRead("test/ethanol.xyz") 246 if err != nil { 247 Te.Error(err) 248 } 249 fmt.Println("Calling with 4") 250 mol2.Del(4) 251 XYZFileWrite("test/ethanolDel4.xyz", mol2.Coords[0], mol2) 252 253 } 254 255 func TestWater(Te *testing.T) { 256 // runtime.GOMAXPROCS(2) /////////////////////////// 257 mol, err := XYZFileRead("test/sample.xyz") 258 if err != nil { 259 Te.Error(err) 260 } 261 for i := 0; i < 6; i++ { 262 s := new(Atom) 263 if i == 0 || i == 3 { 264 s.Symbol = "O" 265 } else { 266 s.Symbol = "H" 267 } 268 mol.AppendAtom(s) 269 } 270 mol.SetCharge(1) 271 mol.SetMulti(1) 272 c2 := v3.Zeros(mol.Len()) 273 v := v3.Zeros(6) 274 l, _ := mol.Coords[0].Dims() 275 fmt.Println(l, mol.Len()) 276 c2.Stack(mol.Coords[0], v) 277 mol.Coords[0] = c2 278 c := mol.Coords[0].VecView(43) 279 h1 := mol.Coords[0].VecView(42) 280 coords := v3.Zeros(mol.Len()) 281 coords.Copy(mol.Coords[0]) 282 w1 := MakeWater(c, h1, 2, Deg2Rad*30, true) 283 w2 := MakeWater(c, h1, 2, Deg2Rad*-30, false) 284 tmp := v3.Zeros(6) 285 tmp.Stack(w1, w2) 286 fmt.Println("tmp water", w1, w2, tmp, c, h1) 287 coords.SetMatrix(mol.Len()-6, 0, tmp) 288 XYZFileWrite("test/WithWater.xyz", coords, mol) 289 fmt.Println("Done TestWater") 290 } 291 292 func TesstFixPDB(Te *testing.T) { 293 mol, err := PDBFileRead("test/2c9vbroken.pdb", true) 294 if err != nil { 295 Te.Error(err) 296 } 297 FixNumbering(mol) 298 PDBFileWrite("test/2c9vfixed.pdb", mol.Coords[0], mol, nil) 299 fmt.Println("DoneTestFixPDB") 300 } 301 302 //will fail if reduce is not installed! 303 func TTTestReduce(Te *testing.T) { //silenced 304 fmt.Println("Start TestReduce") 305 mol, err := PDBFileRead("test/2c9v.pdb", true) 306 if err != nil { 307 Te.Error(err) 308 } 309 logger, err := os.Create("test/reducereport.log") 310 if err != nil { 311 Te.Error(err) 312 } 313 defer logger.Close() 314 mol2, err := Reduce(mol, mol.Coords[0], 2, logger, "") 315 if err != nil { 316 Te.Error(err) 317 } 318 PDBFileWrite("test/2c9vHReduce.pdb", mol2.Coords[0], mol2, nil) 319 fmt.Println("END TestReduce") 320 } 321 322 func TTestShape(Te *testing.T) { 323 myhandle, _ := os.Open("test/2c9v.pdb") 324 mol1, err := PDBRead(myhandle, true) //true means that we try to read the symbol from the PDB file. 325 masses, err := mol1.Masses() 326 if err != nil { 327 Te.Error(err) 328 } 329 moment, err := MomentTensor(mol1.Coords[0], masses) 330 if err != nil { 331 Te.Error(err) 332 } 333 rhos, err := Rhos(moment) 334 if err != nil { 335 Te.Error(err) 336 } 337 linear, circular, err := RhoShapeIndexes(rhos) 338 if err != nil { 339 Te.Error(err) 340 } 341 fmt.Println("liner,circular distortion:", linear, circular) 342 lin2, circ2, err := EasyShape(mol1.Coords[0], -1, mol1) 343 if err != nil { 344 Te.Error(err) 345 } 346 fmt.Println("Easier way linear,circular:", lin2, circ2) 347 mol2, _ := XYZFileRead("test/sample_plane.xyz") 348 lin3, circ3, err := EasyShape(mol2.Coords[0], -1, mol2) 349 if err != nil { 350 Te.Error(err) 351 } 352 fmt.Println("sample_plane.xyz shape indicators; linear,circular:", lin3, circ3) 353 //now the shapetests batterty! 354 mol2, _ = XYZFileRead("test/shapetests/porphyrin.xyz") 355 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 356 if err != nil { 357 Te.Error(err) 358 } 359 fmt.Println("porphyrin.xyz shape indicators; linear,circular:", lin3, circ3) 360 mol2, _ = XYZFileRead("test/shapetests/2-mesoporphyrin.xyz") 361 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 362 if err != nil { 363 Te.Error(err) 364 } 365 fmt.Println("2-mesoporphyrin.xyz shape indicators; linear,circular:", lin3, circ3) 366 mol2, _ = XYZFileRead("test/shapetests/4-mesoporphyrin.xyz") 367 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 368 if err != nil { 369 Te.Error(err) 370 } 371 fmt.Println("4-mesoporphyrin.xyz shape indicators; linear,circular:", lin3, circ3) 372 mol2, _ = XYZFileRead("test/shapetests/heptane.xyz") 373 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 374 if err != nil { 375 Te.Error(err) 376 } 377 fmt.Println("heptane.xyz shape indicators; linear,circular:", lin3, circ3) 378 mol2, _ = XYZFileRead("test/shapetests/decane.xyz") 379 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 380 if err != nil { 381 Te.Error(err) 382 } 383 fmt.Println("decane.xyz shape indicators; linear,circular:", lin3, circ3) 384 mol2, _ = XYZFileRead("test/shapetests/phenantrene.xyz") 385 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 386 if err != nil { 387 Te.Error(err) 388 } 389 fmt.Println("phenantrene.xyz shape indicators; linear,circular:", lin3, circ3) 390 mol2, _ = XYZFileRead("test/shapetests/methylphenantrene.xyz") 391 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 392 if err != nil { 393 Te.Error(err) 394 } 395 fmt.Println("methylphenantrene shape indicators; linear,circular:", lin3, circ3) 396 mol2, _ = XYZFileRead("test/shapetests/tbutylphenantrene.xyz") 397 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 398 if err != nil { 399 Te.Error(err) 400 } 401 fmt.Println("tbutylphenantrene shape indicators; linear,circular:", lin3, circ3) 402 mol2, _ = XYZFileRead("test/shapetests/fullerene20.xyz") 403 lin3, circ3, err = EasyShape(mol2.Coords[0], -1, mol2) 404 if err != nil { 405 Te.Error(err) 406 } 407 fmt.Println("fullerene20.xyz shape indicators; linear,circular:", lin3, circ3) 408 mol2, _ = XYZFileRead("test/shapetests/fullerene60.xyz") 409 lin3, circ3, err = EasyShape(mol2.Coords[0], 0.0001, mol2) //maybe it's too symmetrical for the default epsilon? 410 if err != nil { 411 Te.Error(err) 412 } 413 fmt.Println("fullerene60.xyz shape indicators; linear,circular:", lin3, circ3) 414 415 } 416 417 //Here PDBRead and PDBWrite are tested 418 func TTestSuper(Te *testing.T) { 419 backbone := []string{"CA", "C", "N"} //The PDB name of the atoms in the backbone. 420 myhandle, _ := os.Open("test/2c9v.pdb") 421 mol1, err := PDBRead(myhandle, true) //true means that we try to read the symbol from the PDB file. 422 mol2, err2 := PDBFileRead("test/1uxm.pdb", true) 423 if err != nil || err2 != nil { 424 panic("Unable to open input files!") 425 } 426 mols := []*Molecule{mol1, mol2} 427 superlist := make([][]int, 2, 2) 428 //We collect the atoms that are part of the backbone. 429 for molnumber, mol := range mols { 430 for atomindex, atom := range mol.Atoms { 431 if isInString(backbone, atom.Name) && atom.Chain == "A" { 432 superlist[molnumber] = append(superlist[molnumber], atomindex) 433 // fmt.Println(atom) 434 } 435 } 436 } 437 fmt.Println("superlists!!", len(superlist[0]), len(superlist[1])) 438 mol1.Coords[0], err = Super(mol1.Coords[0], mol2.Coords[0], superlist[0], superlist[1]) 439 rmsd1, _ := rMSD(mol1.Coords[0], mol2.Coords[0], superlist[0], superlist[1]) 440 rmsd2, _ := RMSD(mol1.Coords[0], mol2.Coords[0], superlist[0], superlist[1]) 441 fmt.Println("RMSDs for proteins!", rmsd2, rmsd1) 442 fmt.Println("Atoms superimposed:", len(superlist[0])) 443 if err != nil { 444 panic(err.Error()) 445 } 446 newname := "test/2c9v_super.pdb" //This is the PDBWrite written file 447 pdbout, _ := os.Create(newname) 448 PDBWrite(pdbout, mol1.Coords[0], mol1, nil) 449 //Now for a full molecule 450 ptest, _ := XYZFileRead("test/Rotated.xyz") 451 ptempla, _ := XYZFileRead("test/sample_plane.xyz") 452 newp, err := Super(ptest.Coords[0], ptempla.Coords[0], nil, nil) 453 rmsd2, _ = RMSD(newp, ptempla.Coords[0]) 454 rmsd3, _ := RMSD(newp, ptempla.Coords[0], nil, nil) 455 rmsd1, _ = rMSD(newp, ptempla.Coords[0], nil, nil) 456 fmt.Println("RMSD mol (should be 0):", rmsd1, rmsd2, rmsd3) 457 if err != nil { 458 panic(err.Error()) 459 } 460 XYZFileWrite("test/SuperPlane.xyz", newp, ptest) 461 462 } 463 464 func TestRotateBz(Te *testing.T) { 465 runtime.GOMAXPROCS(2) 466 fmt.Println("Here we go!") 467 mol, err := XYZFileRead("test/BZ.xyz") 468 if err != nil { 469 panic(err.Error()) 470 } 471 carbonIn := []int{} 472 bzIn := []int{} 473 for i := 0; i < mol.Len(); i++ { 474 at := mol.Atom(i) 475 if at.Symbol == "C" { 476 bzIn = append(bzIn, i) 477 carbonIn = append(carbonIn, i) 478 } else if at.Symbol == "H" { 479 bzIn = append(bzIn, i) 480 } 481 } 482 coordsI := mol.Coords[0] 483 carbons := v3.Zeros(6) 484 bz := v3.Zeros(12) 485 carbons.SomeVecs(coordsI, carbonIn) 486 coords := v3.Zeros(mol.Len()) 487 coords, _, _ = MassCenter(coordsI, carbons, nil) 488 bz.SomeVecs(coords, bzIn) 489 carbons.SomeVecs(coords, carbonIn) 490 planevec, err := BestPlane(carbons, nil) 491 if err != nil { 492 if e, ok := err.(Error); ok { 493 fmt.Println("DEcoration:", e.Decorate("")) 494 } 495 Te.Errorf(err.Error()) 496 } 497 basename := "BZ" 498 newcoords := v3.Zeros(mol.Len()) 499 origin := v3.Zeros(1) 500 bzcopy := v3.Zeros(12) 501 bzcopy2 := v3.Zeros(12) //testing 502 rot := v3.Zeros(12) 503 rot3 := v3.Zeros(12) 504 for _, angle := range []float64{0, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 180} { 505 bzcopy.Copy(bz) 506 bzcopy2.Copy(bz) //testing 507 rot = Rotate(bzcopy, rot, planevec, Deg2Rad*angle) 508 rot3 = RotateSer(bzcopy, rot, planevec, Deg2Rad*angle) 509 rot2, _ := EulerRotateAbout(bzcopy2, origin, planevec, Deg2Rad*angle) //should be the same as the previous 510 if !mat.EqualApprox(rot, rot2, 0.01) { 511 Te.Errorf("Rotors Rotate and EulerRotate not equal for angle %3.2f", angle) 512 } else if !mat.EqualApprox(rot2, rot3, 0.01) { 513 Te.Errorf("Rotors RotateSer and EulerRotate not equal for angle %3.2f", angle) 514 515 } else { 516 fmt.Println("Rotors EQUAL for angle", angle) 517 518 } 519 fmt.Println("rot", rot, "rot2", rot2) 520 newcoords.Copy(coords) 521 newcoords.SetVecs(rot, bzIn) 522 //test 523 // tempcoords.Stack(planevec,origin) 524 // testxyz.Stack(newcoords,tempcoords) 525 //end 526 XYZFileWrite(fmt.Sprintf("test/%s-%3.1f.xyz", basename, angle), newcoords, mol) 527 528 } 529 // fmt.Println(mol, planevec) 530 } 531 532 func TestProjectionAndAntiProjection(Te *testing.T) { 533 A := v3.Zeros(1) 534 A.Set(0, 0, 2.0) 535 B, _ := v3.NewMatrix([]float64{1, 1, 0}) 536 C := AntiProjection(A, B) 537 D := Projection(B, A) 538 fmt.Println("Projection of B on A (D)", D) 539 fmt.Println("Anti-projection of A on B (C):", C) 540 fmt.Println("Norm of C: ", C.Norm(2), " Norm of A,B: ", A.Norm(2), B.Norm(2), "Norm of D:", D.Norm(2)) 541 } 542 543 func TestBondsBz(Te *testing.T) { 544 runtime.GOMAXPROCS(2) 545 fmt.Println("Bonds a la carga!") 546 mol, err := XYZFileRead("test/BZBonds.xyz") 547 if err != nil { 548 panic(err.Error()) 549 } 550 var C *Atom 551 552 for i := 0; i < mol.Len(); i++ { 553 C = mol.Atoms[i] 554 if C.Symbol == "C" { 555 break 556 } 557 } 558 mol.AssignBonds(mol.Coords[0]) 559 for i, v := range mol.Atoms { 560 fmt.Printf("Atom %d index %d %s has %d bonds. It's bonded to atoms:\n", i, v.index, v.Symbol, len(v.Bonds)) 561 for _, w := range v.Bonds { 562 a := w.Cross(v) 563 fmt.Printf("Atom: %d %s\n", a.index, a.Symbol) 564 } 565 fmt.Printf("\n") 566 567 } 568 569 path := BondedPaths(C, C.index, &BondedOptions{OnlyShortest: true}) 570 if path == nil { 571 Te.Errorf("No path found!") 572 } 573 fmt.Printf("Path %v has %d nodes\n", path, len(path)) 574 575 } 576 577 func TestBondsPorf(Te *testing.T) { 578 runtime.GOMAXPROCS(2) 579 fmt.Println("A por las tienoporfi!") 580 mol, err := XYZFileRead("test/sample_plane.xyz") 581 if err != nil { 582 panic(err.Error()) 583 } 584 var S *Atom 585 586 for i := 0; i < mol.Len(); i++ { 587 S = mol.Atoms[i] 588 if S.Symbol == "S" { 589 break 590 } 591 } 592 mol.AssignBonds(mol.Coords[0]) 593 594 spath := BondedPaths(S, S.index) 595 if spath == nil { 596 Te.Errorf("No short path found!") 597 } 598 fmt.Printf("Short path %v has %d nodes\n", spath[0], len(spath[0])) 599 paths := BondedPaths(S, S.index) 600 fmt.Printf("Long path %v has %d nodes\n", paths[len(paths)-1], len(paths[len(paths)-1])) 601 602 }