src/gromacs/mdlib/qm_gaussian.c

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