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