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