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