src/gromacs/mdlib/qm_gaussian.c

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