github.com/rmera/gochem@v0.7.1/qm/xtb.go (about) 1 /* 2 * xtb.go, part of gochem. 3 * 4 * 5 * Copyright 2016 Raul Mera <rmera{at}chemDOThelsinkiDOTfi> 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU Lesser General Public License as 9 * published by the Free Software Foundation; either version 2.1 of the 10 * License, or (at your option) any later version. 11 * 12 * This program is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General 18 * Public License along with this program. If not, see 19 * <http://www.gnu.org/licenses/>. 20 * 21 * 22 * 23 */ 24 //In order to use this part of the library you need the xtb program, which must be obtained from Prof. Stefan Grimme's group. 25 //Please cite the the xtb references if you used the program. 26 27 package qm 28 29 import ( 30 // "bufio" 31 "bufio" 32 "fmt" 33 "math" 34 "os" 35 "os/exec" 36 "runtime" 37 "strconv" 38 "strings" 39 40 chem "github.com/rmera/gochem" 41 v3 "github.com/rmera/gochem/v3" 42 ) 43 44 // XTBHandle represents an xtb calculation 45 type XTBHandle struct { 46 //Note that the default methods and basis vary with each program, and even 47 //for a given program they are NOT considered part of the API, so they can always change. 48 //This is unavoidable, as methods change with time 49 command string 50 inputname string 51 nCPU int 52 options []string 53 gfnff bool 54 relconstraints bool 55 force float64 56 wrkdir string 57 inputfile string 58 } 59 60 // NewXTBHandle initializes and returns an xtb handle 61 // with values set to their defaults. Defaults might change 62 // as new methods appear, so they are not part of the API. 63 func NewXTBHandle() *XTBHandle { 64 run := new(XTBHandle) 65 run.SetDefaults() 66 return run 67 } 68 69 //XTBHandle methods 70 71 // SetnCPU sets the number of CPU to be used 72 func (O *XTBHandle) SetnCPU(cpu int) { 73 O.nCPU = cpu 74 } 75 76 // Command returns the path and name for the xtb excecutable 77 func (O *XTBHandle) Command() string { 78 return O.command 79 } 80 81 // SetName sets the name for the calculations 82 // which is defines the input and output file names 83 func (O *XTBHandle) SetName(name string) { 84 O.inputname = name 85 } 86 87 // SetCommand sets the path and name for the xtb excecutable 88 func (O *XTBHandle) SetCommand(name string) { 89 O.command = name 90 } 91 92 // SetWorkDir sets the name of the working directory for the calculations 93 func (O *XTBHandle) SetWorkDir(d string) { 94 O.wrkdir = d 95 } 96 97 // RelConstraints sets the use of relative contraints 98 // instead of absolute position restraints 99 // with the force force constant. If force is 100 // less than 0, the default value is employed. 101 func (O *XTBHandle) RelConstraints(force float64) { 102 if force > 0 { 103 //goChem units (Kcal/A^2) to xtb units (Eh/Bohr^2) 104 O.force = force * (chem.Kcal2H / (chem.A2Bohr * chem.A2Bohr)) 105 } 106 O.relconstraints = true 107 108 } 109 110 // SetDefaults sets calculations parameters to their defaults. 111 // Defaults might change 112 // as new methods appear, so they are not part of the API. 113 func (O *XTBHandle) SetDefaults() { 114 O.command = os.ExpandEnv("xtb") 115 // if O.command == "/xtb" { //if XTBHOME was not defined 116 // O.command = "./xtb" 117 // } 118 cpu := runtime.NumCPU() / 2 119 O.nCPU = cpu 120 121 } 122 123 func (O *XTBHandle) seticonstraints(Q *Calc, xcontrol []string) []string { 124 //Here xtb expects distances in A and angles in deg, so no conversion needed. 125 var g2x = map[byte]string{ 126 'B': "distance: ", 127 'A': "angle: ", 128 'D': "dihedral: ", 129 } 130 131 for _, v := range Q.IConstraints { 132 constra := g2x[v.Class] 133 134 for _, w := range v.CAtoms { 135 constra += strconv.Itoa(w+1) + ", " //1-based indexes for xtb 136 } 137 strings.TrimRight(constra, ",") 138 if v.UseVal { 139 constra += fmt.Sprintf(" %4.2f\n", v.Val) 140 } else { 141 constra += " auto\n" 142 } 143 xcontrol = append(xcontrol, constra) 144 145 } 146 return xcontrol 147 } 148 149 // BuildInput builds an input for XTB. Right now it's very limited, only singlets are allowed and 150 // only unconstrained optimizations and single-points. 151 func (O *XTBHandle) BuildInput(coords *v3.Matrix, atoms chem.AtomMultiCharger, Q *Calc) error { 152 //Now lets write the thing 153 if O.wrkdir != "" { 154 O.wrkdir += "/" 155 } 156 w := O.wrkdir 157 if O.inputname == "" { 158 O.inputname = "gochem" 159 } 160 //Only error so far 161 if atoms == nil || coords == nil { 162 return Error{ErrMissingCharges, "XTB", O.inputname, "", []string{"BuildInput"}, true} 163 } 164 err := chem.XYZFileWrite(w+O.inputname+".xyz", coords, atoms) 165 if err != nil { 166 return Error{ErrCantInput, "XTB", O.inputname, "", []string{"BuildInput"}, true} 167 } 168 // mem := "" 169 if Q.Memory != 0 { 170 //Here we can adjust memory if needed 171 } 172 173 xcontroltxt := make([]string, 0, 10) 174 O.options = make([]string, 0, 6) 175 O.options = append(O.options, O.command) 176 if Q.Method == "gfnff" { 177 O.gfnff = true 178 } 179 O.options = append(O.options, O.inputname+".xyz") 180 O.options = append(O.options, fmt.Sprintf("-c %d", atoms.Charge())) 181 O.options = append(O.options, fmt.Sprintf("-u %d", (atoms.Multi()-1))) 182 if O.nCPU > 1 { 183 O.options = append(O.options, fmt.Sprintf("-P %d", O.nCPU)) 184 } 185 //Added new things to select a method in xtb 186 if !isInString([]string{"gfn1", "gfn2", "gfn0", "gfnff"}, Q.Method) { 187 O.options = append(O.options, "--gfn 2") //default method 188 } else if Q.Method != "gfnff" { 189 m := strings.ReplaceAll(Q.Method, "gfn", "") //so m should be "0", "1" or "2" 190 O.options = append(O.options, "--gfn "+m) //default method 191 } 192 193 if Q.Dielectric > 0 && Q.Method != "gfn0" { //as of the current version, gfn0 doesn't support implicit solvation 194 solvent, ok := dielectric2Solvent[int(Q.Dielectric)] 195 if ok { 196 O.options = append(O.options, "--alpb "+solvent) 197 } 198 } 199 //O.options = append(O.options, "-gfn") 200 fixed := "" 201 fixtoken := "$fix\n" 202 if O.relconstraints { 203 force := "" 204 if O.force > 0 { 205 force = fmt.Sprintf("force constant = %4.2f\n", O.force) 206 } 207 fixtoken = "$constrain\n" + force 208 } 209 if Q.CConstraints != nil || Q.IConstraints != nil { 210 xcontroltxt = append(xcontroltxt, fixtoken) 211 if Q.CConstraints != nil { 212 fixed = "atoms: " 213 for _, v := range Q.CConstraints { 214 fixed = fixed + strconv.Itoa(v+1) + ", " //1-based indexes 215 } 216 strings.TrimRight(fixed, ",") 217 fixed = fixed + "\n" 218 xcontroltxt = append(xcontroltxt, fixed) 219 } 220 if Q.IConstraints != nil { 221 if !O.relconstraints { 222 xcontroltxt = append(xcontroltxt, ("$end\n$constrain\n")) 223 } 224 xcontroltxt = O.seticonstraints(Q, xcontroltxt) 225 226 } 227 xcontroltxt = append(xcontroltxt, "$end\n") 228 229 } 230 jc := jobChoose{} 231 jc.opti = func() { 232 add := "-o normal" 233 if Q.OptTightness == 2 { 234 add = "-o tight" 235 } else if Q.OptTightness > 2 { 236 add = "-o verytight" 237 } 238 O.options = append(O.options, add) 239 } 240 jc.forces = func() { 241 O.options = append(O.options, "--ohess") 242 } 243 244 jc.md = func() { 245 O.options = append(O.options, "--md") 246 //There are specific settings needed with gfnff, mainly, a shorter timestep 247 //The restart=false option doesn't have any effect, but it's added so it's easier later to use sed or whatever to change it to true, and restart 248 //a calculation. 249 if Q.Method == "gfnff" { 250 xcontroltxt = append(xcontroltxt, fmt.Sprintf("$md\n temp=%5.3f\n time=%d\n velo=false\n nvt=true\n step=2.0\n hmass=4.0\n shake=0\n restart=false\n$end", Q.MDTemp, Q.MDTime)) 251 } else { 252 xcontroltxt = append(xcontroltxt, fmt.Sprintf("$md\n temp=%5.3f\n time=%d\n velo=false\n nvt=true\n restart=false\n$end", Q.MDTemp, Q.MDTime)) 253 } 254 255 } 256 // O.options = append(O.options, "--input xcontrol") 257 O.options = append(O.options, Q.Others) 258 Q.Job.Do(jc) 259 if len(xcontroltxt) == 0 { 260 return nil //no need to write a control file 261 } 262 O.inputfile = O.inputname + ".inp" //if not input file was written 263 //this will just be an empty string. 264 xcontrol, err := os.Create(w + O.inputfile) 265 if err != nil { 266 return err 267 } 268 for _, v := range xcontroltxt { 269 xcontrol.WriteString(v) 270 271 } 272 xcontrol.Close() 273 return nil 274 } 275 276 // Run runs the command given by the string O.command 277 // it waits or not for the result depending on wait. 278 // Not waiting for results works 279 // only for unix-compatible systems, as it uses bash and nohup. 280 func (O *XTBHandle) Run(wait bool) (err error) { 281 var com string 282 extraoptions := "" 283 if len(O.options) >= 3 { 284 extraoptions = strings.Join(O.options[2:], " ") 285 } 286 inputfile := "" 287 if O.inputfile != "" { 288 inputfile = fmt.Sprintf("--input %s", O.inputfile) 289 } 290 291 if O.gfnff { 292 com = fmt.Sprintf(" --gfnff %s.xyz %s %s > %s.out 2>&1", O.inputname, inputfile, extraoptions, O.inputname) 293 } else { 294 295 com = fmt.Sprintf(" %s.xyz %s %s > %s.out 2>&1", O.inputname, inputfile, extraoptions, O.inputname) 296 } 297 if wait { 298 //It would be nice to have this logging as an option. 299 //log.Printf(O.command + com) //this is stderr, I suppose 300 command := exec.Command("sh", "-c", O.command+com) 301 command.Dir = O.wrkdir 302 err = command.Run() 303 304 } else { 305 command := exec.Command("sh", "-c", "nohup "+O.command+com) 306 command.Dir = O.wrkdir 307 err = command.Start() 308 } 309 if err != nil { 310 err = Error{ErrNotRunning, XTB, O.inputname, err.Error(), []string{"exec.Start", "Run"}, true} 311 } 312 if err != nil { 313 return err 314 } 315 os.Remove("xtbrestart") 316 return nil 317 } 318 319 // OptimizedGeometry returns the latest geometry from an XTB optimization. It doesn't actually need the chem.Atomer 320 // but requires it so XTBHandle fits with the QM interface. 321 func (O *XTBHandle) OptimizedGeometry(atoms chem.Atomer) (*v3.Matrix, error) { 322 inp := O.wrkdir + O.inputname 323 if !O.normalTermination() { 324 return nil, Error{ErrNoGeometry, XTB, inp, "Calculation didn't end normally", []string{"OptimizedGeometry"}, true} 325 } 326 mol, err := chem.XYZFileRead(O.wrkdir + "xtbopt.xyz") //Trying to run several calculations in parallel in the same directory will fail as the output has always the same name. 327 if err != nil { 328 return nil, Error{ErrNoGeometry, XTB, inp, "", []string{"OptimizedGeometry"}, true} 329 } 330 return mol.Coords[0], nil 331 } 332 333 // This checks that an xtb calculation has terminated normally 334 // I know this duplicates code, I wrote this one first and then the other one. 335 func (O *XTBHandle) normalTermination() bool { 336 inp := O.wrkdir + O.inputname 337 if searchBackwards("normal termination of x", fmt.Sprintf("%s.out", inp)) != "" || searchBackwards("abnormal termination of x", fmt.Sprintf("%s.out", inp)) == "" { 338 return true 339 } 340 // fmt.Println(fmt.Sprintf("%s.out", O.inputname), searchBackwards("normal termination of x",fmt.Sprintf("%s.out", O.inputname))) //////////////////// 341 return false 342 } 343 344 // search a file backwards, i.e., starting from the end, for a string. Returns the line that contains the string, or an empty string. 345 // I really really should have commented this one. 346 func searchBackwards(str, filename string) string { 347 var ini int64 = 0 348 var end int64 = 0 349 var first bool 350 first = true 351 buf := make([]byte, 1) 352 f, err := os.Open(filename) 353 if err != nil { 354 return "" 355 } 356 defer f.Close() 357 var i int64 = 1 358 for ; ; i++ { 359 if _, err := f.Seek(-1*i, 2); err != nil { 360 return "" 361 } 362 if _, err := f.Read(buf); err != nil { 363 return "" 364 } 365 if buf[0] == byte('\n') && first == false { 366 first = true 367 } else if buf[0] == byte('\n') && end == 0 { 368 end = i 369 } else if buf[0] == byte('\n') && ini == 0 { 370 i-- 371 ini = i 372 f.Seek(-1*(ini), 2) 373 bufF := make([]byte, ini-end) 374 f.Read(bufF) 375 if strings.Contains(string(bufF), str) { 376 return string(bufF) 377 } 378 // first=false 379 end = 0 380 ini = 0 381 } 382 383 } 384 } 385 386 // Energy returns the energy of a previous XTB calculations, in kcal/mol. 387 // Returns error if problem, and also if the energy returned that is product of an 388 // abnormally-terminated ORCA calculation. (in this case error is "Probable problem 389 // in calculation") 390 func (O *XTBHandle) Energy() (float64, error) { 391 inp := O.wrkdir + O.inputname 392 var err error 393 var energy float64 394 energyline := searchBackwards("TOTAL ENERGY", fmt.Sprintf("%s.out", inp)) 395 if energyline == "" { 396 return 0, Error{ErrNoEnergy, XTB, inp, fmt.Sprintf("%s.out", inp), []string{"searchBackwards", "Energy"}, true} 397 } 398 split := strings.Fields(energyline) 399 if len(split) < 5 { 400 return 0, Error{ErrNoEnergy, XTB, inp, err.Error(), []string{"Energy"}, true} 401 402 } 403 energy, err = strconv.ParseFloat(split[3], 64) 404 if err != nil { 405 return 0, Error{ErrNoEnergy, XTB, inp, err.Error(), []string{"strconv.ParseFloat", "Energy"}, true} 406 } 407 408 return energy * chem.H2Kcal, err //dummy thin 409 } 410 411 // LargestImaginary returns the absolute value of the wave number (in 1/cm) for the largest imaginary mode in the vibspectrum file 412 // produced by a forces calculation with xtb. Returns an error and -1 if unable to check. 413 func (O *XTBHandle) LargestImaginary() (float64, error) { 414 largestimag := 0.0 415 vibf, err := os.Open(O.wrkdir + "vibspectrum") 416 if err != nil { 417 //fmt.Println("Unable to open file!!") 418 return -1, Error{ErrCantValue, XTB, "vibspectrum", err.Error(), []string{"os.Open", "LargestImaginary"}, true} 419 } 420 vib := bufio.NewReader(vibf) 421 for i := 0; i < 3; i++ { 422 _, err := vib.ReadString('\n') //The text in "data" could be anything, including just "\n" 423 if err != nil { 424 return -1, Error{ErrCantValue, XTB, "vibspectrum", err.Error(), []string{"ReadString", "LargestImaginary"}, true} 425 426 } 427 } 428 for { 429 line, err := vib.ReadString('\n') 430 if err != nil { //inefficient, (errs[1] can be checked once before), but clearer. 431 if strings.Contains(err.Error(), "EOF") { 432 err = nil //it's not an actual error 433 break 434 } else { 435 return -1, Error{ErrCantValue, XTB, "vibspectrum", err.Error(), []string{"ReadString", "LargestImaginary"}, true} 436 437 } 438 } 439 440 fields := strings.Fields(line) 441 if len(fields) < 5 { 442 return -1, Error{ErrCantValue, XTB, "vibspectrum", "Can't parse vibspectrum", []string{"ReadString", "LargestImaginary"}, true} 443 } 444 wave, err := strconv.ParseFloat(fields[len(fields)-4], 64) 445 if err != nil { 446 return -1, Error{ErrCantValue, XTB, "vibspectrum", "Can't parse vibspectrum", []string{"strconv.ParseFloat", "LargestImaginary"}, true} 447 } 448 if wave > 0.0 { 449 return largestimag, nil //no more imaginary frequencies so we just return the largest so far. 450 } else if math.Abs(wave) > largestimag { 451 largestimag = math.Abs(wave) 452 } 453 } 454 return largestimag, nil 455 } 456 457 // FixImaginary prepares and runs a calculation on a geometry, produced by xtb on a previous Hessian calculation, which 458 // is distorted along the main imaginary mode found, if any. It such mode was not found, and thus the geometry was not 459 // produced by xtb, FixImaginary returns an error. 460 func (O *XTBHandle) FixImaginary(wait bool) error { 461 var com string 462 var err error 463 if _, err := os.Stat("xtbhess.coord"); os.IsNotExist(err) { 464 return fmt.Errorf("xtbhess.coord doesn't exist. There is likely no significant imaginary mode") 465 } 466 if O.gfnff { 467 com = fmt.Sprintf(" --gfnff xtbhess.coord --input %s.inp %s > %s.out 2>&1", O.inputname, strings.Join(O.options[2:], " "), O.inputname) 468 } else { 469 470 com = fmt.Sprintf(" xtbhess.coord --input %s.inp %s > %s.out 2>&1", O.inputname, strings.Join(O.options[2:], " "), O.inputname) 471 } 472 if wait == true { 473 //log.Printf(com) //this is stderr, I suppose 474 command := exec.Command("sh", "-c", O.command+com) 475 command.Dir = O.wrkdir 476 err = command.Run() 477 478 } else { 479 command := exec.Command("sh", "-c", "nohup "+O.command+com) 480 command.Dir = O.wrkdir 481 err = command.Start() 482 } 483 if err != nil { 484 err = Error{ErrNotRunning, XTB, O.inputname, err.Error(), []string{"exec.Start", "Run"}, true} 485 } 486 if err != nil { 487 return err 488 } 489 os.Remove("xtbrestart") 490 return nil 491 492 } 493 494 // FreeEnergy returns the Gibbs free energy of a previous XTB calculations. 495 // A frequencies/solvation calculation is needed for this to work. FreeEnergy does _not_ check that the structure was at a minimum. You can check that with 496 // the LargestIm 497 func (O *XTBHandle) FreeEnergy() (float64, error) { 498 var err error 499 var energy float64 500 inp := O.wrkdir + O.inputname 501 energyline := searchBackwards("total free energy", fmt.Sprintf("%s.out", inp)) 502 if energyline == "" { 503 return 0, Error{ErrNoFreeEnergy, XTB, inp, fmt.Sprintf("%s.out", inp), []string{"searchBackwards", "FreeEnergy"}, true} 504 } 505 split := strings.Fields(energyline) 506 if len(split) < 4 { 507 return 0, Error{ErrNoFreeEnergy, XTB, inp, err.Error(), []string{"Energy"}, true} 508 509 } 510 energy, err = strconv.ParseFloat(split[4], 64) 511 if err != nil { 512 return 0, Error{ErrNoFreeEnergy, XTB, inp, err.Error(), []string{"strconv.ParseFloat", "Energy"}, true} 513 } 514 515 return energy * chem.H2Kcal, err //err should be nil at this point. 516 } 517 518 // MDAverageEnergy gets the average potential and kinetic energy along a trajectory. 519 func (O *XTBHandle) MDAverageEnergy(start, skip int) (float64, float64, error) { 520 inp := O.wrkdir + O.inputname 521 var potential, kinetic float64 522 if !O.normalTermination() { 523 return 0, 0, Error{ErrNoEnergy, XTB, inp, "Calculation didn't end normally", []string{"MDAverageEnergy"}, true} 524 } 525 outname := fmt.Sprintf("%s.out", inp) 526 outfile, err := os.Open(outname) 527 if err != nil { 528 return 0, 0, Error{ErrNoEnergy, XTB, inp, "Couldn't open output file", []string{"MDAverageEnergy"}, true} 529 } 530 out := bufio.NewReader(outfile) 531 reading := false 532 cont := 0 533 read := 0 534 for { 535 line, err := out.ReadString('\n') 536 // fmt.Println("LINE", line) ///////// 537 if err != nil && strings.Contains(err.Error(), "EOF") { 538 break 539 } else if err != nil { 540 return 0, 0, Error{ErrNoEnergy, XTB, inp, "Error while iterating through output file", []string{"MDAverageEnergy"}, true} 541 } 542 if strings.Contains(line, "time (ps) <Epot> Ekin <T> T Etot") { 543 reading = true 544 continue 545 } 546 if !reading { 547 continue 548 } 549 fields := strings.Fields(line) 550 if len(fields) != 7 { 551 continue 552 } 553 cont++ 554 if (cont-1)%skip != 0 || (cont-1) < start { 555 continue 556 } 557 K, err := strconv.ParseFloat(fields[3], 64) 558 if err != nil { 559 return 0, 0, Error{ErrNoEnergy, XTB, inp, fmt.Sprintf("Error while retrieving %d th kinetic energy", cont), []string{"MDAverageEnergy"}, true} 560 561 } 562 V, err := strconv.ParseFloat(fields[3], 64) 563 if err != nil { 564 return 0, 0, Error{ErrNoEnergy, XTB, inp, fmt.Sprintf("Error while retrieving %d th potential energy", cont), []string{"MDAverageEnergy"}, true} 565 566 } 567 fmt.Println("potential", V) ////////// 568 kinetic += K 569 potential += V 570 read++ 571 } 572 N := float64(read) 573 if math.IsNaN(potential/N) || math.IsNaN(kinetic/N) { //note that we still return whatever we got here, in addition to the error. The user can decide. 574 return potential / N, kinetic / N, Error{ErrProbableProblem, XTB, inp, "At least one of the energies is NaN", []string{"MDAverageEnergy"}, true} 575 } 576 return potential / N, kinetic / N, nil 577 } 578 579 var dielectric2Solvent = map[int]string{ 580 80: "h2o", 581 5: "chcl3", 582 9: "ch2cl2", 583 10: "octanol", 584 21: "acetone", 585 37: "acetonitrile", 586 33: "methanol", 587 2: "toluene", 588 1: "hexadecane", //not quite 1 589 7: "thf", 590 47: "dmso", 591 38: "dmf", 592 1000: "woctanol", //This is a hackish way to have both dry and wet octanol. I gave wet octanol an fake epsilon that won't be used by anything else. 593 //really, what I should do is to add to the API a way to specify either epsilon or solvent name. FIX 594 } 595 596 //old code 597 598 // out, err := os.Create(fmt.Sprintf("%s.out", O.inputname)) 599 // if err != nil { 600 // return Error{ErrNotRunning, XTB, O.inputname, "", []string{"Run"}, true} 601 // } 602 // ferr, err := os.Create(fmt.Sprintf("%s.err", O.inputname)) 603 // 604 // if err != nil { 605 // return Error{ErrNotRunning, XTB, O.inputname, "", []string{"Run"}, true} 606 // } 607 // defer out.Close() 608 // defer ferr.Close() 609 // fullCommand:=strings.Join(O.options," ") 610 // fmt.Println(fullCommand) //("Command", O.command, O.options) //////////////////////// 611 // command := exec.Command(fullCommand) //, O.options...) 612 // command.Stdout = out 613 // command.Stderr = ferr 614 // err = command.Run() 615 // fmt.Println(O.command+fmt.Sprintf(" %s.xyz %s > %s.out &", O.inputname, strings.Join(O.options[2:]," "), O.inputname)) ////////////////////////