src/gromacs/mdlib/qm_gaussian.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "qm_gaussian.h"
  40
  41 #include "config.h"
  42
  43 #include <stdio.h>
  44 #include <stdlib.h>
  45 #include <string.h>
  46
  47 #include <cmath>
  48
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/gmxlib/network.h"
  51 #include "gromacs/gmxlib/nrnb.h"
  52 #include "gromacs/math/units.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/mdlib/force.h"
  55 #include "gromacs/mdlib/forcerec.h"
  56 #include "gromacs/mdlib/ns.h"
  57 #include "gromacs/mdlib/qmmm.h"
  58 #include "gromacs/mdtypes/md_enums.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63 // When not built in a configuration with QMMM support, much of this
  64 // code is unreachable by design. Tell clang not to warn about it.
  65 #pragma GCC diagnostic push
  66 #pragma GCC diagnostic ignored "-Wmissing-noreturn"
  67
  68 /* TODO: this should be made thread-safe */
  69
  70 /* Gaussian interface routines */
  71
  72 void init_gaussian(t_QMrec *qm)
  73 {
  74     ivec
  75         basissets[eQMbasisNR] = {{0, 3, 0},
  76                                  {0, 3, 0}, /*added for double sto-3g entry in names.c*/
  77                                  {5, 0, 0},
  78                                  {5, 0, 1},
  79                                  {5, 0, 11},
  80                                  {5, 6, 0},
  81                                  {1, 6, 0},
  82                                  {1, 6, 1},
  83                                  {1, 6, 11},
  84                                  {4, 6, 0}};
  85     char
  86        *buf = nullptr;
  87     int
  88         i;
  89
  90     if (!GMX_QMMM_GAUSSIAN)
  91     {
  92         gmx_fatal(FARGS, "Cannot call GAUSSIAN unless linked against it. Use cmake -DGMX_QMMM_PROGRAM=GAUSSIAN, and ensure that linking will work correctly.");
  93     }
  94
  95     /* using the ivec above to convert the basis read form the mdp file
  96      * in a human readable format into some numbers for the gaussian
  97      * route. This is necessary as we are using non standard routes to
  98      * do SH.
  99      */
 100
 101     /* per layer we make a new subdir for integral file, checkpoint
 102      * files and such. These dirs are stored in the QMrec for
 103      * convenience
 104      */
 105
 106
 107     if (!qm->nQMcpus) /* this we do only once per layer
 108                        * as we call g01 externally
 109                        */
 110
 111     {
 112         for (i = 0; i < DIM; i++)
 113         {
 114             qm->SHbasis[i] = basissets[qm->QMbasis][i];
 115         }
 116
 117         /* init gradually switching on of the SA */
 118         qm->SAstep = 0;
 119         /* we read the number of cpus and environment from the environment
 120          * if set.
 121          */
 122         buf = getenv("GMX_QM_GAUSSIAN_NCPUS");
 123         if (buf)
 124         {
 125             sscanf(buf, "%d", &qm->nQMcpus);
 126         }
 127         else
 128         {
 129             qm->nQMcpus = 1;
 130         }
 131         fprintf(stderr, "number of CPUs for gaussian = %d\n", qm->nQMcpus);
 132         buf = getenv("GMX_QM_GAUSSIAN_MEMORY");
 133         if (buf)
 134         {
 135             sscanf(buf, "%d", &qm->QMmem);
 136         }
 137         else
 138         {
 139             qm->QMmem = 50000000;
 140         }
 141         fprintf(stderr, "memory for gaussian = %d\n", qm->QMmem);
 142         buf = getenv("GMX_QM_ACCURACY");
 143         if (buf)
 144         {
 145             sscanf(buf, "%d", &qm->accuracy);
 146         }
 147         else
 148         {
 149             qm->accuracy = 8;
 150         }
 151         fprintf(stderr, "accuracy in l510 = %d\n", qm->accuracy);
 152
 153         buf = getenv("GMX_QM_CPMCSCF");
 154         if (buf)
 155         {
 156             sscanf(buf, "%d", &i);
 157             qm->cpmcscf = (i != 0);
 158         }
 159         else
 160         {
 161             qm->cpmcscf = FALSE;
 162         }
 163         if (qm->cpmcscf)
 164         {
 165             fprintf(stderr, "using cp-mcscf in l1003\n");
 166         }
 167         else
 168         {
 169             fprintf(stderr, "NOT using cp-mcscf in l1003\n");
 170         }
 171         buf = getenv("GMX_QM_SA_STEP");
 172         if (buf)
 173         {
 174             sscanf(buf, "%d", &qm->SAstep);
 175         }
 176         else
 177         {
 178             /* init gradually switching on of the SA */
 179             qm->SAstep = 0;
 180         }
 181         /* we read the number of cpus and environment from the environment
 182          * if set.
 183          */
 184         fprintf(stderr, "Level of SA at start = %d\n", qm->SAstep);
 185         /* gaussian settings on the system */
 186         buf = getenv("GMX_QM_GAUSS_DIR");
 187
 188         if (buf)
 189         {
 190             qm->gauss_dir = gmx_strdup(buf);
 191         }
 192         else
 193         {
 194             gmx_fatal(FARGS, "no $GMX_QM_GAUSS_DIR, check gaussian manual\n");
 195         }
 196
 197         buf = getenv("GMX_QM_GAUSS_EXE");
 198         if (buf)
 199         {
 200             qm->gauss_exe = gmx_strdup(buf);
 201         }
 202         else
 203         {
 204             gmx_fatal(FARGS, "no $GMX_QM_GAUSS_EXE set, check gaussian manual\n");
 205         }
 206         buf = getenv("GMX_QM_MODIFIED_LINKS_DIR");
 207         if (buf)
 208         {
 209             qm->devel_dir = gmx_strdup (buf);
 210         }
 211         else
 212         {
 213             gmx_fatal(FARGS, "no $GMX_QM_MODIFIED_LINKS_DIR, this is were the modified links reside.\n");
 214         }
 215
 216         /*  if(fr->bRF){*/
 217         /* reactionfield, file is needed using gaussian */
 218         /*    rffile=fopen("rf.dat","w");*/
 219         /*   fprintf(rffile,"%f %f\n",fr->epsilon_r,fr->rcoulomb/BOHR2NM);*/
 220         /* fclose(rffile);*/
 221         /*  }*/
 222     }
 223     fprintf(stderr, "gaussian initialised...\n");
 224 }
 225
 226
 227 static void write_gaussian_SH_input(int step, gmx_bool swap,
 228                                     const t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
 229 {
 230     int
 231         i;
 232     gmx_bool
 233         bSA;
 234     FILE
 235        *out;
 236     t_QMMMrec
 237        *QMMMrec;
 238     QMMMrec = fr->qr;
 239     bSA     = (qm->SAstep > 0);
 240
 241     out = fopen("input.com", "w");
 242     /* write the route */
 243     fprintf(out, "%s", "%scr=input\n");
 244     fprintf(out, "%s", "%rwf=input\n");
 245     fprintf(out, "%s", "%int=input\n");
 246     fprintf(out, "%s", "%d2e=input\n");
 247 /*  if(step)
 248  *   fprintf(out,"%s","%nosave\n");
 249  */
 250     fprintf(out, "%s", "%chk=input\n");
 251     fprintf(out, "%s%d\n", "%mem=", qm->QMmem);
 252     fprintf(out, "%s%3d\n", "%nprocshare=", qm->nQMcpus);
 253
 254     /* use the versions of
 255      * l301 that computes the interaction between MM and QM atoms.
 256      * l510 that can punch the CI coefficients
 257      * l701 that can do gradients on MM atoms
 258      */
 259
 260     /* local version */
 261     fprintf(out, "%s%s%s",
 262             "%subst l510 ",
 263             qm->devel_dir,
 264             "/l510\n");
 265     fprintf(out, "%s%s%s",
 266             "%subst l301 ",
 267             qm->devel_dir,
 268             "/l301\n");
 269     fprintf(out, "%s%s%s",
 270             "%subst l701 ",
 271             qm->devel_dir,
 272             "/l701\n");
 273
 274     fprintf(out, "%s%s%s",
 275             "%subst l1003 ",
 276             qm->devel_dir,
 277             "/l1003\n");
 278     fprintf(out, "%s%s%s",
 279             "%subst l9999 ",
 280             qm->devel_dir,
 281             "/l9999\n");
 282     /* print the nonstandard route
 283      */
 284     fprintf(out, "%s",
 285             "#P nonstd\n 1/18=10,20=1,38=1/1;\n");
 286     fprintf(out, "%s",
 287             " 2/9=110,15=1,17=6,18=5,40=1/2;\n");
 288     if (mm->nrMMatoms)
 289     {
 290         fprintf(out,
 291                 " 3/5=%d,6=%d,7=%d,25=1,32=1,43=1,94=-2/1,2,3;\n",
 292                 qm->SHbasis[0],
 293                 qm->SHbasis[1],
 294                 qm->SHbasis[2]); /*basisset stuff */
 295     }
 296     else
 297     {
 298         fprintf(out,
 299                 " 3/5=%d,6=%d,7=%d,25=1,32=1,43=0,94=-2/1,2,3;\n",
 300                 qm->SHbasis[0],
 301                 qm->SHbasis[1],
 302                 qm->SHbasis[2]); /*basisset stuff */
 303     }
 304     /* development */
 305     if (step+1) /* fetch initial guess from check point file */
 306     {           /* hack, to alyays read from chk file!!!!! */
 307         fprintf(out, "%s%d,%s%d%s", " 4/5=1,7=6,17=",
 308                 qm->CASelectrons,
 309                 "18=", qm->CASorbitals, "/1,5;\n");
 310     }
 311     else /* generate the first checkpoint file */
 312     {
 313         fprintf(out, "%s%d,%s%d%s", " 4/5=0,7=6,17=",
 314                 qm->CASelectrons,
 315                 "18=", qm->CASorbitals, "/1,5;\n");
 316     }
 317     /* the rest of the input depends on where the system is on the PES
 318      */
 319     if (swap && bSA)             /* make a slide to the other surface */
 320     {
 321         if (qm->CASorbitals > 6) /* use direct and no full diag */
 322         {
 323             fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6,97=100/10;\n");
 324         }
 325         else
 326         {
 327             if (qm->cpmcscf)
 328             {
 329                 fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6,97=100/10;\n",
 330                         qm->accuracy);
 331                 if (mm->nrMMatoms > 0)
 332                 {
 333                     fprintf(out, " 7/7=1,16=-2,30=1/1;\n");
 334                 }
 335                 fprintf(out, " 11/31=1,42=1,45=1/1;\n");
 336                 fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1,97=100/3;\n");
 337                 fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n");
 338             }
 339             else
 340             {
 341                 fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6,97=100/10;\n",
 342                         qm->accuracy);
 343                 fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 344             }
 345         }
 346     }
 347     else if (bSA)                /* do a "state-averaged" CAS calculation */
 348     {
 349         if (qm->CASorbitals > 6) /* no full diag */
 350         {
 351             fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6/10;\n");
 352         }
 353         else
 354         {
 355             if (qm->cpmcscf)
 356             {
 357                 fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6/10;\n",
 358                         qm->accuracy);
 359                 if (mm->nrMMatoms > 0)
 360                 {
 361                     fprintf(out, " 7/7=1,16=-2,30=1/1;\n");
 362                 }
 363                 fprintf(out, " 11/31=1,42=1,45=1/1;\n");
 364                 fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1/3;\n");
 365                 fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n");
 366             }
 367             else
 368             {
 369                 fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6/10;\n",
 370                         qm->accuracy);
 371                 fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 372             }
 373         }
 374     }
 375     else if (swap) /* do a "swapped" CAS calculation */
 376     {
 377         if (qm->CASorbitals > 6)
 378         {
 379             fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6,97=100/10;\n");
 380         }
 381         else
 382         {
 383             fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6,97=100/10;\n",
 384                     qm->accuracy);
 385         }
 386         fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 387     }
 388     else /* do a "normal" CAS calculation */
 389     {
 390         if (qm->CASorbitals > 6)
 391         {
 392             fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6/10;\n");
 393         }
 394         else
 395         {
 396             fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6/10;\n",
 397                     qm->accuracy);
 398         }
 399         fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 400     }
 401     fprintf(out, "\ninput-file generated by gromacs\n\n");
 402     fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity);
 403     for (i = 0; i < qm->nrQMatoms; i++)
 404     {
 405         fprintf(out, "%3d %10.7f  %10.7f  %10.7f\n",
 406                 qm->atomicnumberQM[i],
 407                 qm->xQM[i][XX]/BOHR2NM,
 408                 qm->xQM[i][YY]/BOHR2NM,
 409                 qm->xQM[i][ZZ]/BOHR2NM);
 410     }
 411     /* MM point charge data */
 412     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 413     {
 414         fprintf(out, "\n");
 415         for (i = 0; i < mm->nrMMatoms; i++)
 416         {
 417             fprintf(out, "%10.7f  %10.7f  %10.7f %8.4f\n",
 418                     mm->xMM[i][XX]/BOHR2NM,
 419                     mm->xMM[i][YY]/BOHR2NM,
 420                     mm->xMM[i][ZZ]/BOHR2NM,
 421                     mm->MMcharges[i]);
 422         }
 423     }
 424     if (bSA) /* put the SA coefficients at the end of the file */
 425     {
 426         fprintf(out, "\n%10.8f %10.8f\n",
 427                 qm->SAstep*0.5/qm->SAsteps,
 428                 1-qm->SAstep*0.5/qm->SAsteps);
 429         fprintf(stderr, "State Averaging level = %d/%d\n", qm->SAstep, qm->SAsteps);
 430     }
 431     fprintf(out, "\n");
 432     fclose(out);
 433 }  /* write_gaussian_SH_input */
 434
 435 static void write_gaussian_input(int step, const t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
 436 {
 437     int
 438         i;
 439     t_QMMMrec
 440        *QMMMrec;
 441     FILE
 442        *out;
 443
 444     QMMMrec = fr->qr;
 445     out     = fopen("input.com", "w");
 446     /* write the route */
 447
 448     if (qm->QMmethod >= eQMmethodRHF)
 449     {
 450         fprintf(out, "%s",
 451                 "%chk=input\n");
 452     }
 453     else
 454     {
 455         fprintf(out, "%s",
 456                 "%chk=se\n");
 457     }
 458     if (qm->nQMcpus > 1)
 459     {
 460         fprintf(out, "%s%3d\n",
 461                 "%nprocshare=", qm->nQMcpus);
 462     }
 463     fprintf(out, "%s%d\n",
 464             "%mem=", qm->QMmem);
 465     fprintf(out, "%s%s%s",
 466             "%subst l701 ", qm->devel_dir, "/l701\n");
 467     fprintf(out, "%s%s%s",
 468             "%subst l301 ", qm->devel_dir, "/l301\n");
 469     fprintf(out, "%s%s%s",
 470             "%subst l9999 ", qm->devel_dir, "/l9999\n");
 471     if (step)
 472     {
 473         fprintf(out, "%s",
 474                 "#T ");
 475     }
 476     else
 477     {
 478         fprintf(out, "%s",
 479                 "#P ");
 480     }
 481     if (qm->QMmethod == eQMmethodB3LYPLAN)
 482     {
 483         fprintf(out, " %s",
 484                 "B3LYP/GEN Pseudo=Read");
 485     }
 486     else
 487     {
 488         fprintf(out, " %s",
 489                 eQMmethod_names[qm->QMmethod]);
 490
 491         if (qm->QMmethod >= eQMmethodRHF)
 492         {
 493             if (qm->QMmethod == eQMmethodCASSCF)
 494             {
 495                 /* in case of cas, how many electrons and orbitals do we need?
 496                  */
 497                 fprintf(out, "(%d,%d)",
 498                         qm->CASelectrons, qm->CASorbitals);
 499             }
 500             fprintf(out, "/%s",
 501                     eQMbasis_names[qm->QMbasis]);
 502         }
 503     }
 504     if (QMMMrec->QMMMscheme == eQMMMschemenormal && mm->nrMMatoms)
 505     {
 506         fprintf(out, " %s",
 507                 "Charge ");
 508     }
 509     if (step || qm->QMmethod == eQMmethodCASSCF)
 510     {
 511         /* fetch guess from checkpoint file, always for CASSCF */
 512         fprintf(out, "%s", " guess=read");
 513     }
 514     fprintf(out, "\nNosymm units=bohr\n");
 515
 516     fprintf(out, "FORCE Punch=(Derivatives) ");
 517     fprintf(out, "iop(3/33=1)\n\n");
 518     fprintf(out, "input-file generated by gromacs\n\n");
 519     fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity);
 520     for (i = 0; i < qm->nrQMatoms; i++)
 521     {
 522         fprintf(out, "%3d %10.7f  %10.7f  %10.7f\n",
 523                 qm->atomicnumberQM[i],
 524                 qm->xQM[i][XX]/BOHR2NM,
 525                 qm->xQM[i][YY]/BOHR2NM,
 526                 qm->xQM[i][ZZ]/BOHR2NM);
 527     }
 528
 529     /* Pseudo Potential and ECP are included here if selected (MEthod suffix LAN) */
 530     if (qm->QMmethod == eQMmethodB3LYPLAN)
 531     {
 532         fprintf(out, "\n");
 533         for (i = 0; i < qm->nrQMatoms; i++)
 534         {
 535             if (qm->atomicnumberQM[i] < 21)
 536             {
 537                 fprintf(out, "%d ", i+1);
 538             }
 539         }
 540         fprintf(out, "\n%s\n****\n", eQMbasis_names[qm->QMbasis]);
 541
 542         for (i = 0; i < qm->nrQMatoms; i++)
 543         {
 544             if (qm->atomicnumberQM[i] > 21)
 545             {
 546                 fprintf(out, "%d ", i+1);
 547             }
 548         }
 549         fprintf(out, "\n%s\n****\n\n", "lanl2dz");
 550
 551         for (i = 0; i < qm->nrQMatoms; i++)
 552         {
 553             if (qm->atomicnumberQM[i] > 21)
 554             {
 555                 fprintf(out, "%d ", i+1);
 556             }
 557         }
 558         fprintf(out, "\n%s\n", "lanl2dz");
 559     }
 560
 561
 562
 563     /* MM point charge data */
 564     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 565     {
 566         fprintf(stderr, "nr mm atoms in gaussian.c = %d\n", mm->nrMMatoms);
 567         fprintf(out, "\n");
 568         for (i = 0; i < mm->nrMMatoms; i++)
 569         {
 570             fprintf(out, "%10.7f  %10.7f  %10.7f %8.4f\n",
 571                     mm->xMM[i][XX]/BOHR2NM,
 572                     mm->xMM[i][YY]/BOHR2NM,
 573                     mm->xMM[i][ZZ]/BOHR2NM,
 574                     mm->MMcharges[i]);
 575         }
 576     }
 577     fprintf(out, "\n");
 578
 579
 580     fclose(out);
 581
 582 }  /* write_gaussian_input */
 583
 584 static real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], t_QMrec *qm, t_MMrec *mm)
 585 {
 586     int
 587         i;
 588     char
 589         buf[300];
 590     real
 591         QMener;
 592     FILE
 593        *in;
 594
 595     in = fopen("fort.7", "r");
 596
 597     /* (There was additional content in the file in case
 598      *    of QM optimizations / transition state search,
 599      *    which was removed.
 600      */
 601     /* the next line is the energy and in the case of CAS, the energy
 602      * difference between the two states.
 603      */
 604     if (nullptr == fgets(buf, 300, in))
 605     {
 606         gmx_fatal(FARGS, "Error reading Gaussian output");
 607     }
 608
 609 #if GMX_DOUBLE
 610     sscanf(buf, "%lf\n", &QMener);
 611 #else
 612     sscanf(buf, "%f\n", &QMener);
 613 #endif
 614     /* next lines contain the gradients of the QM atoms */
 615     for (i = 0; i < qm->nrQMatoms; i++)
 616     {
 617         if (nullptr == fgets(buf, 300, in))
 618         {
 619             gmx_fatal(FARGS, "Error reading Gaussian output");
 620         }
 621 #if GMX_DOUBLE
 622         sscanf(buf, "%lf %lf %lf\n",
 623                &QMgrad[i][XX],
 624                &QMgrad[i][YY],
 625                &QMgrad[i][ZZ]);
 626 #else
 627         sscanf(buf, "%f %f %f\n",
 628                &QMgrad[i][XX],
 629                &QMgrad[i][YY],
 630                &QMgrad[i][ZZ]);
 631 #endif
 632     }
 633     /* the next lines are the gradients of the MM atoms */
 634     if (qm->QMmethod >= eQMmethodRHF)
 635     {
 636         for (i = 0; i < mm->nrMMatoms; i++)
 637         {
 638             if (nullptr == fgets(buf, 300, in))
 639             {
 640                 gmx_fatal(FARGS, "Error reading Gaussian output");
 641             }
 642 #if GMX_DOUBLE
 643             sscanf(buf, "%lf %lf %lf\n",
 644                    &MMgrad[i][XX],
 645                    &MMgrad[i][YY],
 646                    &MMgrad[i][ZZ]);
 647 #else
 648             sscanf(buf, "%f %f %f\n",
 649                    &MMgrad[i][XX],
 650                    &MMgrad[i][YY],
 651                    &MMgrad[i][ZZ]);
 652 #endif
 653         }
 654     }
 655     fclose(in);
 656     return(QMener);
 657 }
 658
 659 static real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step, t_QMrec *qm, t_MMrec *mm)
 660 {
 661     int
 662         i;
 663     char
 664         buf[300];
 665     real
 666         QMener, DeltaE;
 667     FILE
 668        *in;
 669
 670     in = fopen("fort.7", "r");
 671     /* first line is the energy and in the case of CAS, the energy
 672      * difference between the two states.
 673      */
 674     if (nullptr == fgets(buf, 300, in))
 675     {
 676         gmx_fatal(FARGS, "Error reading Gaussian output");
 677     }
 678
 679 #if GMX_DOUBLE
 680     sscanf(buf, "%lf %lf\n", &QMener, &DeltaE);
 681 #else
 682     sscanf(buf, "%f %f\n",  &QMener, &DeltaE);
 683 #endif
 684
 685     /* switch on/off the State Averaging */
 686
 687     if (DeltaE > qm->SAoff)
 688     {
 689         if (qm->SAstep > 0)
 690         {
 691             qm->SAstep--;
 692         }
 693     }
 694     else if (DeltaE < qm->SAon || (qm->SAstep > 0))
 695     {
 696         if (qm->SAstep < qm->SAsteps)
 697         {
 698             qm->SAstep++;
 699         }
 700     }
 701
 702     /* for debugging: */
 703     fprintf(stderr, "Gap = %5f,SA = %3d\n", DeltaE, static_cast<int>(qm->SAstep > 0));
 704     /* next lines contain the gradients of the QM atoms */
 705     for (i = 0; i < qm->nrQMatoms; i++)
 706     {
 707         if (nullptr == fgets(buf, 300, in))
 708         {
 709             gmx_fatal(FARGS, "Error reading Gaussian output");
 710         }
 711
 712 #if GMX_DOUBLE
 713         sscanf(buf, "%lf %lf %lf\n",
 714                &QMgrad[i][XX],
 715                &QMgrad[i][YY],
 716                &QMgrad[i][ZZ]);
 717 #else
 718         sscanf(buf, "%f %f %f\n",
 719                &QMgrad[i][XX],
 720                &QMgrad[i][YY],
 721                &QMgrad[i][ZZ]);
 722 #endif
 723     }
 724     /* the next lines, are the gradients of the MM atoms */
 725
 726     for (i = 0; i < mm->nrMMatoms; i++)
 727     {
 728         if (nullptr == fgets(buf, 300, in))
 729         {
 730             gmx_fatal(FARGS, "Error reading Gaussian output");
 731         }
 732 #if GMX_DOUBLE
 733         sscanf(buf, "%lf %lf %lf\n",
 734                &MMgrad[i][XX],
 735                &MMgrad[i][YY],
 736                &MMgrad[i][ZZ]);
 737 #else
 738         sscanf(buf, "%f %f %f\n",
 739                &MMgrad[i][XX],
 740                &MMgrad[i][YY],
 741                &MMgrad[i][ZZ]);
 742 #endif
 743     }
 744
 745     /* the next line contains the two CI eigenvector elements */
 746     if (nullptr == fgets(buf, 300, in))
 747     {
 748         gmx_fatal(FARGS, "Error reading Gaussian output");
 749     }
 750     if (!step)
 751     {
 752         sscanf(buf, "%d", &qm->CIdim);
 753         snew(qm->CIvec1, qm->CIdim);
 754         snew(qm->CIvec1old, qm->CIdim);
 755         snew(qm->CIvec2, qm->CIdim);
 756         snew(qm->CIvec2old, qm->CIdim);
 757     }
 758     else
 759     {
 760         /* before reading in the new current CI vectors, copy the current
 761          * CI vector into the old one.
 762          */
 763         for (i = 0; i < qm->CIdim; i++)
 764         {
 765             qm->CIvec1old[i] = qm->CIvec1[i];
 766             qm->CIvec2old[i] = qm->CIvec2[i];
 767         }
 768     }
 769     /* first vector */
 770     for (i = 0; i < qm->CIdim; i++)
 771     {
 772         if (nullptr == fgets(buf, 300, in))
 773         {
 774             gmx_fatal(FARGS, "Error reading Gaussian output");
 775         }
 776 #if GMX_DOUBLE
 777         sscanf(buf, "%lf\n", &qm->CIvec1[i]);
 778 #else
 779         sscanf(buf, "%f\n", &qm->CIvec1[i]);
 780 #endif
 781     }
 782     /* second vector */
 783     for (i = 0; i < qm->CIdim; i++)
 784     {
 785         if (nullptr == fgets(buf, 300, in))
 786         {
 787             gmx_fatal(FARGS, "Error reading Gaussian output");
 788         }
 789 #if GMX_DOUBLE
 790         sscanf(buf, "%lf\n", &qm->CIvec2[i]);
 791 #else
 792         sscanf(buf, "%f\n", &qm->CIvec2[i]);
 793 #endif
 794     }
 795     fclose(in);
 796     return(QMener);
 797 }
 798
 799 static real inproduct(const real *a, const real *b, int n)
 800 {
 801     int
 802         i;
 803     real
 804         dot = 0.0;
 805
 806     /* computes the inner product between two vectors (a.b), both of
 807      * which have length n.
 808      */
 809     for (i = 0; i < n; i++)
 810     {
 811         dot += a[i]*b[i];
 812     }
 813     return(dot);
 814 }
 815
 816 static int hop(int step, t_QMrec *qm)
 817 {
 818     int
 819         swap = 0;
 820     real
 821         d11 = 0.0, d12 = 0.0, d21 = 0.0, d22 = 0.0;
 822
 823     /* calculates the inproduct between the current Ci vector and the
 824      * previous CI vector. A diabatic hop will be made if d12 and d21
 825      * are much bigger than d11 and d22. In that case hop returns true,
 826      * otherwise it returns false.
 827      */
 828     if (step) /* only go on if more than one step has been done */
 829     {
 830         d11 = inproduct(qm->CIvec1, qm->CIvec1old, qm->CIdim);
 831         d12 = inproduct(qm->CIvec1, qm->CIvec2old, qm->CIdim);
 832         d21 = inproduct(qm->CIvec2, qm->CIvec1old, qm->CIdim);
 833         d22 = inproduct(qm->CIvec2, qm->CIvec2old, qm->CIdim);
 834     }
 835     fprintf(stderr, "-------------------\n");
 836     fprintf(stderr, "d11 = %13.8f\n", d11);
 837     fprintf(stderr, "d12 = %13.8f\n", d12);
 838     fprintf(stderr, "d21 = %13.8f\n", d21);
 839     fprintf(stderr, "d22 = %13.8f\n", d22);
 840     fprintf(stderr, "-------------------\n");
 841
 842     if ((fabs(d12) > 0.5) && (fabs(d21) > 0.5))
 843     {
 844         swap = 1;
 845     }
 846
 847     return(swap);
 848 }
 849
 850 static void do_gaussian(int step, char *exe)
 851 {
 852     char
 853         buf[STRLEN];
 854
 855     /* make the call to the gaussian binary through system()
 856      * The location of the binary will be picked up from the
 857      * environment using getenv().
 858      */
 859     if (step) /* hack to prevent long inputfiles */
 860     {
 861         sprintf(buf, "%s < %s > %s",
 862                 exe,
 863                 "input.com",
 864                 "input.log");
 865     }
 866     else
 867     {
 868         sprintf(buf, "%s < %s > %s",
 869                 exe,
 870                 "input.com",
 871                 "input.log");
 872     }
 873     fprintf(stderr, "Calling '%s'\n", buf);
 874     if (system(buf) != 0)
 875     {
 876         gmx_fatal(FARGS, "Call to '%s' failed\n", buf);
 877     }
 878 }
 879
 880 real call_gaussian(const t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 881 {
 882     /* normal gaussian jobs */
 883     static int
 884         step = 0;
 885     int
 886         i, j;
 887     real
 888         QMener = 0.0;
 889     rvec
 890        *QMgrad, *MMgrad;
 891     char
 892        *exe;
 893
 894     snew(exe, 30);
 895     sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe);
 896     snew(QMgrad, qm->nrQMatoms);
 897     snew(MMgrad, mm->nrMMatoms);
 898
 899     write_gaussian_input(step, fr, qm, mm);
 900     do_gaussian(step, exe);
 901     QMener = read_gaussian_output(QMgrad, MMgrad, qm, mm);
 902     /* put the QMMM forces in the force array and to the fshift
 903      */
 904     for (i = 0; i < qm->nrQMatoms; i++)
 905     {
 906         for (j = 0; j < DIM; j++)
 907         {
 908             f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];
 909             fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
 910         }
 911     }
 912     for (i = 0; i < mm->nrMMatoms; i++)
 913     {
 914         for (j = 0; j < DIM; j++)
 915         {
 916             f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];
 917             fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
 918         }
 919     }
 920     QMener = QMener*HARTREE2KJ*AVOGADRO;
 921     step++;
 922     free(exe);
 923     return(QMener);
 924
 925 } /* call_gaussian */
 926
 927 real call_gaussian_SH(const t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 928 {
 929     /* a gaussian call routine intended for doing diabatic surface
 930      * "sliding". See the manual for the theoretical background of this
 931      * TSH method.
 932      */
 933     static int
 934         step = 0;
 935     int
 936         state, i, j;
 937     real
 938         QMener = 0.0;
 939     static  gmx_bool
 940         swapped = FALSE; /* handle for identifying the current PES */
 941     gmx_bool
 942         swap = FALSE;    /* the actual swap */
 943     rvec
 944        *QMgrad, *MMgrad;
 945     char
 946        *buf;
 947     char
 948        *exe;
 949
 950     snew(exe, 30);
 951     sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe);
 952     /* hack to do ground state simulations */
 953     if (!step)
 954     {
 955         snew(buf, 20);
 956         buf = getenv("GMX_QM_GROUND_STATE");
 957         if (buf)
 958         {
 959             sscanf(buf, "%d", &state);
 960         }
 961         else
 962         {
 963             state = 2;
 964         }
 965         if (state == 1)
 966         {
 967             swapped = TRUE;
 968         }
 969     }
 970     /* end of hack */
 971
 972
 973     /* copy the QMMMrec pointer */
 974     snew(QMgrad, qm->nrQMatoms);
 975     snew(MMgrad, mm->nrMMatoms);
 976     /* at step 0 there should be no SA */
 977     /*  if(!step)
 978      * qr->bSA=FALSE;*/
 979     /* temporray set to step + 1, since there is a chk start */
 980     write_gaussian_SH_input(step, swapped, fr, qm, mm);
 981
 982     do_gaussian(step, exe);
 983     QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
 984
 985     /* check for a surface hop. Only possible if we were already state
 986      * averaging.
 987      */
 988     if (qm->SAstep > 0)
 989     {
 990         if (!swapped)
 991         {
 992             swap    = ((step != 0) && (hop(step, qm) != 0));
 993             swapped = swap;
 994         }
 995         else /* already on the other surface, so check if we go back */
 996         {
 997             swap    = ((step != 0) && (hop(step, qm) != 0));
 998             swapped = !swap; /* so swapped shoud be false again */
 999         }
1000         if (swap)            /* change surface, so do another call */
1001         {
1002             write_gaussian_SH_input(step, swapped, fr, qm, mm);
1003             do_gaussian(step, exe);
1004             QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
1005         }
1006     }
1007     /* add the QMMM forces to the gmx force array and fshift
1008      */
1009     for (i = 0; i < qm->nrQMatoms; i++)
1010     {
1011         for (j = 0; j < DIM; j++)
1012         {
1013             f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];
1014             fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
1015         }
1016     }
1017     for (i = 0; i < mm->nrMMatoms; i++)
1018     {
1019         for (j = 0; j < DIM; j++)
1020         {
1021             f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];
1022             fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
1023         }
1024     }
1025     QMener = QMener*HARTREE2KJ*AVOGADRO;
1026     fprintf(stderr, "step %5d, SA = %5d, swap = %5d\n",
1027             step, static_cast<int>(qm->SAstep > 0), static_cast<int>(swapped));
1028     step++;
1029     free(exe);
1030     return(QMener);
1031
1032 } /* call_gaussian_SH */
1033
1034 #pragma GCC diagnostic pop