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