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