src/gromacs/mdlib/qm_gaussian.cpp

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