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