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,2019, 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/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_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
 220 {
 221     int        i;
 222     gmx_bool   bSA;
 223     FILE*      out;
 224     t_QMMMrec* QMMMrec;
 225     QMMMrec = fr->qr;
 226     bSA     = (qm->SAstep > 0);
 227
 228     out = fopen("input.com", "w");
 229     /* write the route */
 230     fprintf(out, "%s", "%scr=input\n");
 231     fprintf(out, "%s", "%rwf=input\n");
 232     fprintf(out, "%s", "%int=input\n");
 233     fprintf(out, "%s", "%d2e=input\n");
 234     /*  if(step)
 235      *   fprintf(out,"%s","%nosave\n");
 236      */
 237     fprintf(out, "%s", "%chk=input\n");
 238     fprintf(out, "%s%d\n", "%mem=", qm->QMmem);
 239     fprintf(out, "%s%3d\n", "%nprocshare=", qm->nQMcpus);
 240
 241     /* use the versions of
 242      * l301 that computes the interaction between MM and QM atoms.
 243      * l510 that can punch the CI coefficients
 244      * l701 that can do gradients on MM atoms
 245      */
 246
 247     /* local version */
 248     fprintf(out, "%s%s%s", "%subst l510 ", qm->devel_dir, "/l510\n");
 249     fprintf(out, "%s%s%s", "%subst l301 ", qm->devel_dir, "/l301\n");
 250     fprintf(out, "%s%s%s", "%subst l701 ", qm->devel_dir, "/l701\n");
 251
 252     fprintf(out, "%s%s%s", "%subst l1003 ", qm->devel_dir, "/l1003\n");
 253     fprintf(out, "%s%s%s", "%subst l9999 ", qm->devel_dir, "/l9999\n");
 254     /* print the nonstandard route
 255      */
 256     fprintf(out, "%s", "#P nonstd\n 1/18=10,20=1,38=1/1;\n");
 257     fprintf(out, "%s", " 2/9=110,15=1,17=6,18=5,40=1/2;\n");
 258     if (mm->nrMMatoms)
 259     {
 260         fprintf(out, " 3/5=%d,6=%d,7=%d,25=1,32=1,43=1,94=-2/1,2,3;\n", qm->SHbasis[0],
 261                 qm->SHbasis[1], qm->SHbasis[2]); /*basisset stuff */
 262     }
 263     else
 264     {
 265         fprintf(out, " 3/5=%d,6=%d,7=%d,25=1,32=1,43=0,94=-2/1,2,3;\n", qm->SHbasis[0],
 266                 qm->SHbasis[1], qm->SHbasis[2]); /*basisset stuff */
 267     }
 268     /* development */
 269     if (step + 1) /* fetch initial guess from check point file */
 270     {             /* hack, to alyays read from chk file!!!!! */
 271         fprintf(out, "%s%d,%s%d%s", " 4/5=1,7=6,17=", qm->CASelectrons, "18=", qm->CASorbitals,
 272                 "/1,5;\n");
 273     }
 274     else /* generate the first checkpoint file */
 275     {
 276         fprintf(out, "%s%d,%s%d%s", " 4/5=0,7=6,17=", qm->CASelectrons, "18=", qm->CASorbitals,
 277                 "/1,5;\n");
 278     }
 279     /* the rest of the input depends on where the system is on the PES
 280      */
 281     if (swap && bSA) /* make a slide to the other surface */
 282     {
 283         if (qm->CASorbitals > 6) /* use direct and no full diag */
 284         {
 285             fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6,97=100/10;\n");
 286         }
 287         else
 288         {
 289             if (qm->cpmcscf)
 290             {
 291                 fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6,97=100/10;\n", qm->accuracy);
 292                 if (mm->nrMMatoms > 0)
 293                 {
 294                     fprintf(out, " 7/7=1,16=-2,30=1/1;\n");
 295                 }
 296                 fprintf(out, " 11/31=1,42=1,45=1/1;\n");
 297                 fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1,97=100/3;\n");
 298                 fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n");
 299             }
 300             else
 301             {
 302                 fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6,97=100/10;\n", qm->accuracy);
 303                 fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 304             }
 305         }
 306     }
 307     else if (bSA) /* do a "state-averaged" CAS calculation */
 308     {
 309         if (qm->CASorbitals > 6) /* no full diag */
 310         {
 311             fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6/10;\n");
 312         }
 313         else
 314         {
 315             if (qm->cpmcscf)
 316             {
 317                 fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6/10;\n", qm->accuracy);
 318                 if (mm->nrMMatoms > 0)
 319                 {
 320                     fprintf(out, " 7/7=1,16=-2,30=1/1;\n");
 321                 }
 322                 fprintf(out, " 11/31=1,42=1,45=1/1;\n");
 323                 fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1/3;\n");
 324                 fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n");
 325             }
 326             else
 327             {
 328                 fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6/10;\n", qm->accuracy);
 329                 fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 330             }
 331         }
 332     }
 333     else if (swap) /* do a "swapped" CAS calculation */
 334     {
 335         if (qm->CASorbitals > 6)
 336         {
 337             fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6,97=100/10;\n");
 338         }
 339         else
 340         {
 341             fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6,97=100/10;\n", qm->accuracy);
 342         }
 343         fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 344     }
 345     else /* do a "normal" CAS calculation */
 346     {
 347         if (qm->CASorbitals > 6)
 348         {
 349             fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6/10;\n");
 350         }
 351         else
 352         {
 353             fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6/10;\n", qm->accuracy);
 354         }
 355         fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
 356     }
 357     fprintf(out, "\ninput-file generated by gromacs\n\n");
 358     fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity);
 359     for (i = 0; i < qm->nrQMatoms; i++)
 360     {
 361         fprintf(out, "%3d %10.7f  %10.7f  %10.7f\n", qm->atomicnumberQM[i],
 362                 qm->xQM[i][XX] / BOHR2NM, qm->xQM[i][YY] / BOHR2NM, qm->xQM[i][ZZ] / BOHR2NM);
 363     }
 364     /* MM point charge data */
 365     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 366     {
 367         fprintf(out, "\n");
 368         for (i = 0; i < mm->nrMMatoms; i++)
 369         {
 370             fprintf(out, "%10.7f  %10.7f  %10.7f %8.4f\n", mm->xMM[i][XX] / BOHR2NM,
 371                     mm->xMM[i][YY] / BOHR2NM, mm->xMM[i][ZZ] / BOHR2NM, mm->MMcharges[i]);
 372         }
 373     }
 374     if (bSA) /* put the SA coefficients at the end of the file */
 375     {
 376         fprintf(out, "\n%10.8f %10.8f\n", qm->SAstep * 0.5 / qm->SAsteps, 1 - qm->SAstep * 0.5 / qm->SAsteps);
 377         fprintf(stderr, "State Averaging level = %d/%d\n", qm->SAstep, qm->SAsteps);
 378     }
 379     fprintf(out, "\n");
 380     fclose(out);
 381 } /* write_gaussian_SH_input */
 382
 383 static void write_gaussian_input(int step, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
 384 {
 385     int        i;
 386     t_QMMMrec* QMMMrec;
 387     FILE*      out;
 388
 389     QMMMrec = fr->qr;
 390     out     = fopen("input.com", "w");
 391     /* write the route */
 392
 393     if (qm->QMmethod >= eQMmethodRHF)
 394     {
 395         fprintf(out, "%s", "%chk=input\n");
 396     }
 397     else
 398     {
 399         fprintf(out, "%s", "%chk=se\n");
 400     }
 401     if (qm->nQMcpus > 1)
 402     {
 403         fprintf(out, "%s%3d\n", "%nprocshare=", qm->nQMcpus);
 404     }
 405     fprintf(out, "%s%d\n", "%mem=", qm->QMmem);
 406     fprintf(out, "%s%s%s", "%subst l701 ", qm->devel_dir, "/l701\n");
 407     fprintf(out, "%s%s%s", "%subst l301 ", qm->devel_dir, "/l301\n");
 408     fprintf(out, "%s%s%s", "%subst l9999 ", qm->devel_dir, "/l9999\n");
 409     if (step)
 410     {
 411         fprintf(out, "%s", "#T ");
 412     }
 413     else
 414     {
 415         fprintf(out, "%s", "#P ");
 416     }
 417     if (qm->QMmethod == eQMmethodB3LYPLAN)
 418     {
 419         fprintf(out, " %s", "B3LYP/GEN Pseudo=Read");
 420     }
 421     else
 422     {
 423         fprintf(out, " %s", eQMmethod_names[qm->QMmethod]);
 424
 425         if (qm->QMmethod >= eQMmethodRHF)
 426         {
 427             if (qm->QMmethod == eQMmethodCASSCF)
 428             {
 429                 /* in case of cas, how many electrons and orbitals do we need?
 430                  */
 431                 fprintf(out, "(%d,%d)", qm->CASelectrons, qm->CASorbitals);
 432             }
 433             fprintf(out, "/%s", eQMbasis_names[qm->QMbasis]);
 434         }
 435     }
 436     if (QMMMrec->QMMMscheme == eQMMMschemenormal && mm->nrMMatoms)
 437     {
 438         fprintf(out, " %s", "Charge ");
 439     }
 440     if (step || qm->QMmethod == eQMmethodCASSCF)
 441     {
 442         /* fetch guess from checkpoint file, always for CASSCF */
 443         fprintf(out, "%s", " guess=read");
 444     }
 445     fprintf(out, "\nNosymm units=bohr\n");
 446
 447     fprintf(out, "FORCE Punch=(Derivatives) ");
 448     fprintf(out, "iop(3/33=1)\n\n");
 449     fprintf(out, "input-file generated by gromacs\n\n");
 450     fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity);
 451     for (i = 0; i < qm->nrQMatoms; i++)
 452     {
 453         fprintf(out, "%3d %10.7f  %10.7f  %10.7f\n", qm->atomicnumberQM[i],
 454                 qm->xQM[i][XX] / BOHR2NM, qm->xQM[i][YY] / BOHR2NM, qm->xQM[i][ZZ] / BOHR2NM);
 455     }
 456
 457     /* Pseudo Potential and ECP are included here if selected (MEthod suffix LAN) */
 458     if (qm->QMmethod == eQMmethodB3LYPLAN)
 459     {
 460         fprintf(out, "\n");
 461         for (i = 0; i < qm->nrQMatoms; i++)
 462         {
 463             if (qm->atomicnumberQM[i] < 21)
 464             {
 465                 fprintf(out, "%d ", i + 1);
 466             }
 467         }
 468         fprintf(out, "\n%s\n****\n", eQMbasis_names[qm->QMbasis]);
 469
 470         for (i = 0; i < qm->nrQMatoms; i++)
 471         {
 472             if (qm->atomicnumberQM[i] > 21)
 473             {
 474                 fprintf(out, "%d ", i + 1);
 475             }
 476         }
 477         fprintf(out, "\n%s\n****\n\n", "lanl2dz");
 478
 479         for (i = 0; i < qm->nrQMatoms; i++)
 480         {
 481             if (qm->atomicnumberQM[i] > 21)
 482             {
 483                 fprintf(out, "%d ", i + 1);
 484             }
 485         }
 486         fprintf(out, "\n%s\n", "lanl2dz");
 487     }
 488
 489
 490     /* MM point charge data */
 491     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 492     {
 493         fprintf(stderr, "nr mm atoms in gaussian.c = %d\n", mm->nrMMatoms);
 494         fprintf(out, "\n");
 495         for (i = 0; i < mm->nrMMatoms; i++)
 496         {
 497             fprintf(out, "%10.7f  %10.7f  %10.7f %8.4f\n", mm->xMM[i][XX] / BOHR2NM,
 498                     mm->xMM[i][YY] / BOHR2NM, mm->xMM[i][ZZ] / BOHR2NM, mm->MMcharges[i]);
 499         }
 500     }
 501     fprintf(out, "\n");
 502
 503
 504     fclose(out);
 505
 506 } /* write_gaussian_input */
 507
 508 static real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], t_QMrec* qm, t_MMrec* mm)
 509 {
 510     int   i;
 511     char  buf[300];
 512     real  QMener;
 513     FILE* in;
 514
 515     in = fopen("fort.7", "r");
 516
 517     /* (There was additional content in the file in case
 518      *    of QM optimizations / transition state search,
 519      *    which was removed.
 520      */
 521     /* the next line is the energy and in the case of CAS, the energy
 522      * difference between the two states.
 523      */
 524     if (nullptr == fgets(buf, 300, in))
 525     {
 526         gmx_fatal(FARGS, "Error reading Gaussian output");
 527     }
 528
 529 #if GMX_DOUBLE
 530     sscanf(buf, "%lf\n", &QMener);
 531 #else
 532     sscanf(buf, "%f\n", &QMener);
 533 #endif
 534     /* next lines contain the gradients of the QM atoms */
 535     for (i = 0; i < qm->nrQMatoms; i++)
 536     {
 537         if (nullptr == fgets(buf, 300, in))
 538         {
 539             gmx_fatal(FARGS, "Error reading Gaussian output");
 540         }
 541 #if GMX_DOUBLE
 542         sscanf(buf, "%lf %lf %lf\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
 543 #else
 544         sscanf(buf, "%f %f %f\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
 545 #endif
 546     }
 547     /* the next lines are the gradients of the MM atoms */
 548     if (qm->QMmethod >= eQMmethodRHF)
 549     {
 550         for (i = 0; i < mm->nrMMatoms; i++)
 551         {
 552             if (nullptr == fgets(buf, 300, in))
 553             {
 554                 gmx_fatal(FARGS, "Error reading Gaussian output");
 555             }
 556 #if GMX_DOUBLE
 557             sscanf(buf, "%lf %lf %lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
 558 #else
 559             sscanf(buf, "%f %f %f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
 560 #endif
 561         }
 562     }
 563     fclose(in);
 564     return (QMener);
 565 }
 566
 567 static real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step, t_QMrec* qm, t_MMrec* mm)
 568 {
 569     int   i;
 570     char  buf[300];
 571     real  QMener, DeltaE;
 572     FILE* in;
 573
 574     in = fopen("fort.7", "r");
 575     /* first line is the energy and in the case of CAS, the energy
 576      * difference between the two states.
 577      */
 578     if (nullptr == fgets(buf, 300, in))
 579     {
 580         gmx_fatal(FARGS, "Error reading Gaussian output");
 581     }
 582
 583 #if GMX_DOUBLE
 584     sscanf(buf, "%lf %lf\n", &QMener, &DeltaE);
 585 #else
 586     sscanf(buf, "%f %f\n", &QMener, &DeltaE);
 587 #endif
 588
 589     /* switch on/off the State Averaging */
 590
 591     if (DeltaE > qm->SAoff)
 592     {
 593         if (qm->SAstep > 0)
 594         {
 595             qm->SAstep--;
 596         }
 597     }
 598     else if (DeltaE < qm->SAon || (qm->SAstep > 0))
 599     {
 600         if (qm->SAstep < qm->SAsteps)
 601         {
 602             qm->SAstep++;
 603         }
 604     }
 605
 606     /* for debugging: */
 607     fprintf(stderr, "Gap = %5f,SA = %3d\n", DeltaE, static_cast<int>(qm->SAstep > 0));
 608     /* next lines contain the gradients of the QM atoms */
 609     for (i = 0; i < qm->nrQMatoms; i++)
 610     {
 611         if (nullptr == fgets(buf, 300, in))
 612         {
 613             gmx_fatal(FARGS, "Error reading Gaussian output");
 614         }
 615
 616 #if GMX_DOUBLE
 617         sscanf(buf, "%lf %lf %lf\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
 618 #else
 619         sscanf(buf, "%f %f %f\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
 620 #endif
 621     }
 622     /* the next lines, are the gradients of the MM atoms */
 623
 624     for (i = 0; i < mm->nrMMatoms; i++)
 625     {
 626         if (nullptr == fgets(buf, 300, in))
 627         {
 628             gmx_fatal(FARGS, "Error reading Gaussian output");
 629         }
 630 #if GMX_DOUBLE
 631         sscanf(buf, "%lf %lf %lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
 632 #else
 633         sscanf(buf, "%f %f %f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
 634 #endif
 635     }
 636
 637     /* the next line contains the two CI eigenvector elements */
 638     if (nullptr == fgets(buf, 300, in))
 639     {
 640         gmx_fatal(FARGS, "Error reading Gaussian output");
 641     }
 642     if (!step)
 643     {
 644         sscanf(buf, "%d", &qm->CIdim);
 645         snew(qm->CIvec1, qm->CIdim);
 646         snew(qm->CIvec1old, qm->CIdim);
 647         snew(qm->CIvec2, qm->CIdim);
 648         snew(qm->CIvec2old, qm->CIdim);
 649     }
 650     else
 651     {
 652         /* before reading in the new current CI vectors, copy the current
 653          * CI vector into the old one.
 654          */
 655         for (i = 0; i < qm->CIdim; i++)
 656         {
 657             qm->CIvec1old[i] = qm->CIvec1[i];
 658             qm->CIvec2old[i] = qm->CIvec2[i];
 659         }
 660     }
 661     /* first vector */
 662     for (i = 0; i < qm->CIdim; i++)
 663     {
 664         if (nullptr == fgets(buf, 300, in))
 665         {
 666             gmx_fatal(FARGS, "Error reading Gaussian output");
 667         }
 668 #if GMX_DOUBLE
 669         sscanf(buf, "%lf\n", &qm->CIvec1[i]);
 670 #else
 671         sscanf(buf, "%f\n", &qm->CIvec1[i]);
 672 #endif
 673     }
 674     /* second vector */
 675     for (i = 0; i < qm->CIdim; i++)
 676     {
 677         if (nullptr == fgets(buf, 300, in))
 678         {
 679             gmx_fatal(FARGS, "Error reading Gaussian output");
 680         }
 681 #if GMX_DOUBLE
 682         sscanf(buf, "%lf\n", &qm->CIvec2[i]);
 683 #else
 684         sscanf(buf, "%f\n", &qm->CIvec2[i]);
 685 #endif
 686     }
 687     fclose(in);
 688     return (QMener);
 689 }
 690
 691 static real inproduct(const real* a, const real* b, int n)
 692 {
 693     int  i;
 694     real dot = 0.0;
 695
 696     /* computes the inner product between two vectors (a.b), both of
 697      * which have length n.
 698      */
 699     for (i = 0; i < n; i++)
 700     {
 701         dot += a[i] * b[i];
 702     }
 703     return (dot);
 704 }
 705
 706 static int hop(int step, t_QMrec* qm)
 707 {
 708     int  swap = 0;
 709     real d11 = 0.0, d12 = 0.0, d21 = 0.0, d22 = 0.0;
 710
 711     /* calculates the inproduct between the current Ci vector and the
 712      * previous CI vector. A diabatic hop will be made if d12 and d21
 713      * are much bigger than d11 and d22. In that case hop returns true,
 714      * otherwise it returns false.
 715      */
 716     if (step) /* only go on if more than one step has been done */
 717     {
 718         d11 = inproduct(qm->CIvec1, qm->CIvec1old, qm->CIdim);
 719         d12 = inproduct(qm->CIvec1, qm->CIvec2old, qm->CIdim);
 720         d21 = inproduct(qm->CIvec2, qm->CIvec1old, qm->CIdim);
 721         d22 = inproduct(qm->CIvec2, qm->CIvec2old, qm->CIdim);
 722     }
 723     fprintf(stderr, "-------------------\n");
 724     fprintf(stderr, "d11 = %13.8f\n", d11);
 725     fprintf(stderr, "d12 = %13.8f\n", d12);
 726     fprintf(stderr, "d21 = %13.8f\n", d21);
 727     fprintf(stderr, "d22 = %13.8f\n", d22);
 728     fprintf(stderr, "-------------------\n");
 729
 730     if ((fabs(d12) > 0.5) && (fabs(d21) > 0.5))
 731     {
 732         swap = 1;
 733     }
 734
 735     return (swap);
 736 }
 737
 738 static void do_gaussian(int step, char* exe)
 739 {
 740     char buf[STRLEN];
 741
 742     /* make the call to the gaussian binary through system()
 743      * The location of the binary will be picked up from the
 744      * environment using getenv().
 745      */
 746     if (step) /* hack to prevent long inputfiles */
 747     {
 748         sprintf(buf, "%s < %s > %s", exe, "input.com", "input.log");
 749     }
 750     else
 751     {
 752         sprintf(buf, "%s < %s > %s", exe, "input.com", "input.log");
 753     }
 754     fprintf(stderr, "Calling '%s'\n", buf);
 755     if (system(buf) != 0)
 756     {
 757         gmx_fatal(FARGS, "Call to '%s' failed\n", buf);
 758     }
 759 }
 760
 761 real call_gaussian(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
 762 {
 763     /* normal gaussian jobs */
 764     static int step = 0;
 765     int        i, j;
 766     real       QMener = 0.0;
 767     rvec *     QMgrad, *MMgrad;
 768     char*      exe;
 769
 770     snew(exe, 30);
 771     sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe);
 772     snew(QMgrad, qm->nrQMatoms);
 773     snew(MMgrad, mm->nrMMatoms);
 774
 775     write_gaussian_input(step, fr, qm, mm);
 776     do_gaussian(step, exe);
 777     QMener = read_gaussian_output(QMgrad, MMgrad, qm, mm);
 778     /* put the QMMM forces in the force array and to the fshift
 779      */
 780     for (i = 0; i < qm->nrQMatoms; i++)
 781     {
 782         for (j = 0; j < DIM; j++)
 783         {
 784             f[i][j]      = HARTREE_BOHR2MD * QMgrad[i][j];
 785             fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
 786         }
 787     }
 788     for (i = 0; i < mm->nrMMatoms; i++)
 789     {
 790         for (j = 0; j < DIM; j++)
 791         {
 792             f[i + qm->nrQMatoms][j]      = HARTREE_BOHR2MD * MMgrad[i][j];
 793             fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
 794         }
 795     }
 796     QMener = QMener * HARTREE2KJ * AVOGADRO;
 797     step++;
 798     free(exe);
 799     return (QMener);
 800
 801 } /* call_gaussian */
 802
 803 real call_gaussian_SH(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
 804 {
 805     /* a gaussian call routine intended for doing diabatic surface
 806      * "sliding". See the manual for the theoretical background of this
 807      * TSH method.
 808      */
 809     static int      step = 0;
 810     int             state, i, j;
 811     real            QMener  = 0.0;
 812     static gmx_bool swapped = FALSE; /* handle for identifying the current PES */
 813     gmx_bool        swap    = FALSE; /* the actual swap */
 814     rvec *          QMgrad, *MMgrad;
 815     char*           buf;
 816     char*           exe;
 817
 818     snew(exe, 30);
 819     sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe);
 820     /* hack to do ground state simulations */
 821     if (!step)
 822     {
 823         snew(buf, 20);
 824         buf = getenv("GMX_QM_GROUND_STATE");
 825         if (buf)
 826         {
 827             sscanf(buf, "%d", &state);
 828         }
 829         else
 830         {
 831             state = 2;
 832         }
 833         if (state == 1)
 834         {
 835             swapped = TRUE;
 836         }
 837     }
 838     /* end of hack */
 839
 840
 841     /* copy the QMMMrec pointer */
 842     snew(QMgrad, qm->nrQMatoms);
 843     snew(MMgrad, mm->nrMMatoms);
 844     /* at step 0 there should be no SA */
 845     /*  if(!step)
 846      * qr->bSA=FALSE;*/
 847     /* temporray set to step + 1, since there is a chk start */
 848     write_gaussian_SH_input(step, swapped, fr, qm, mm);
 849
 850     do_gaussian(step, exe);
 851     QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
 852
 853     /* check for a surface hop. Only possible if we were already state
 854      * averaging.
 855      */
 856     if (qm->SAstep > 0)
 857     {
 858         if (!swapped)
 859         {
 860             swap    = ((step != 0) && (hop(step, qm) != 0));
 861             swapped = swap;
 862         }
 863         else /* already on the other surface, so check if we go back */
 864         {
 865             swap    = ((step != 0) && (hop(step, qm) != 0));
 866             swapped = !swap; /* so swapped shoud be false again */
 867         }
 868         if (swap) /* change surface, so do another call */
 869         {
 870             write_gaussian_SH_input(step, swapped, fr, qm, mm);
 871             do_gaussian(step, exe);
 872             QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
 873         }
 874     }
 875     /* add the QMMM forces to the gmx force array and fshift
 876      */
 877     for (i = 0; i < qm->nrQMatoms; i++)
 878     {
 879         for (j = 0; j < DIM; j++)
 880         {
 881             f[i][j]      = HARTREE_BOHR2MD * QMgrad[i][j];
 882             fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
 883         }
 884     }
 885     for (i = 0; i < mm->nrMMatoms; i++)
 886     {
 887         for (j = 0; j < DIM; j++)
 888         {
 889             f[i + qm->nrQMatoms][j]      = HARTREE_BOHR2MD * MMgrad[i][j];
 890             fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
 891         }
 892     }
 893     QMener = QMener * HARTREE2KJ * AVOGADRO;
 894     fprintf(stderr, "step %5d, SA = %5d, swap = %5d\n", step, static_cast<int>(qm->SAstep > 0),
 895             static_cast<int>(swapped));
 896     step++;
 897     free(exe);
 898     return (QMener);
 899
 900 } /* call_gaussian_SH */
 901
 902 #pragma GCC diagnostic pop