src/gromacs/mdlib/qm_gaussian.c

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