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