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