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