src/gromacs/mdlib/qm_gaussian.c

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