src/mdlib/qm_gaussian.c

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