src/mdlib/qm_gamess.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_GAMESS
  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 /* QMMM sub routines */
  71 /* mopac interface routines */
  72
  73
  74 #ifndef F77_FUNC
  75 #define F77_FUNC(name,NAME) name ## _
  76 #endif
  77
  78
  79 void
  80 F77_FUNC(inigms,IMIGMS)(void);
  81
  82 void
  83 F77_FUNC(endgms,ENDGMS)(void);
  84
  85 void
  86 F77_FUNC(grads,GRADS)(int *nrqmat,real *qmcrd,int *nrmmat, real *mmchrg,
  87                       real *mmcrd, real *qmgrad,real *mmgrad, real *energy);
  88
  89
  90
  91 void init_gamess(t_commrec *cr,t_QMrec *qm, t_MMrec *mm){
  92   /* it works hopelessly complicated :-)
  93    * first a file is written. Then the standard gamess input/output
  94    * routine is called (no system()!) to set up all fortran arrays.
  95    * this routine writes a punch file, like in a normal gamess run.
  96    * via this punch file the other games routines, needed for gradient
  97    * and energy evaluations are called. This setup works fine for
  98    * dynamics simulations. 7-6-2002 (London)
  99    */
 100   int
 101     i,j,rank;
 102   FILE
 103     *out;
 104   char
 105     periodic_system[37][3]={"XX","H ","He","Li","Be","B ","C ","N ",
 106                             "O ","F ","Ne","Na","Mg","Al","Si","P ",
 107                             "S ","Cl","Ar","K ","Ca","Sc","Ti","V ",
 108                             "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga",
 109                             "Ge","As","Se","Br","Kr"};
 110
 111   if (PAR(cr)){
 112
 113     if MASTER(cr){
 114       out=fopen("FOR009","w");
 115       /* of these options I am not completely sure....  the overall
 116        * preformance on more than 4 cpu's is rather poor at the moment.
 117        */
 118       fprintf(out,"memory 48000000\nPARALLEL IOMODE SCREENED\n");
 119       fprintf(out,"ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
 120               qm->nelectrons,qm->multiplicity);
 121       for (i=0;i<qm->nrQMatoms;i++){
 122 #ifdef DOUBLE
 123         fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n",
 124                 i/2.,
 125                 i/3.,
 126                 i/4.,
 127                 qm->atomicnumberQM[i]*1.0,
 128                 periodic_system[qm->atomicnumberQM[i]]);
 129 #else
 130         fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  %2s\n",
 131                 i/2.,
 132                 i/3.,
 133                 i/4.,
 134                 qm->atomicnumberQM[i]*1.0,
 135                 periodic_system[qm->atomicnumberQM[i]]);
 136 #endif
 137       }
 138       if(mm->nrMMatoms){
 139         for (j=i;j<i+2;j++){
 140 #ifdef DOUBLE
 141           fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 142                   j/5.,
 143                   j/6.,
 144                   j/7.,
 145                   1.0);
 146 #else
 147           fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 148                   j/5.,
 149                   j/6.,
 150                   j/7.,
 151                   2.0);
 152 #endif
 153         }
 154       }
 155       if(!qm->bTS)
 156         fprintf(out,"END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 157                 eQMbasis_names[qm->QMbasis],
 158                 eQMmethod_names[qm->QMmethod]); /* see enum.h */
 159       else
 160         fprintf(out,"END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
 161                 eQMbasis_names[qm->QMbasis],
 162                 eQMmethod_names[qm->QMmethod]); /* see enum.h */
 163       fclose(out);
 164     }
 165     gmx_barrier(cr);
 166     F77_FUNC(inigms,IMIGMS)();
 167   }
 168   else{ /* normal serial run */
 169
 170     out=fopen("FOR009","w");
 171     /* of these options I am not completely sure....  the overall
 172      * preformance on more than 4 cpu's is rather poor at the moment.
 173      */
 174     fprintf(out,"ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
 175             qm->nelectrons,qm->multiplicity);
 176     for (i=0;i<qm->nrQMatoms;i++){
 177 #ifdef DOUBLE
 178       fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n",
 179               i/2.,
 180               i/3.,
 181               i/4.,
 182               qm->atomicnumberQM[i]*1.0,
 183               periodic_system[qm->atomicnumberQM[i]]);
 184 #else
 185       fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  %2s\n",
 186               i/2.,
 187               i/3.,
 188               i/4.,
 189               qm->atomicnumberQM[i]*1.0,
 190               periodic_system[qm->atomicnumberQM[i]]);
 191 #endif
 192     }
 193     if(mm->nrMMatoms){
 194       for (j=i;j<i+2;j++){
 195 #ifdef DOUBLE
 196         fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 197                 j/5.,
 198                 j/6.,
 199                 j/7.,
 200                 1.0);
 201 #else
 202         fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 203                 j/5.,
 204                 j/6.,
 205                 j/7.,
 206                 2.0);
 207 #endif
 208       }
 209     }
 210     if(!qm->bTS)
 211       fprintf(out,"END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 212               eQMbasis_names[qm->QMbasis],
 213               eQMmethod_names[qm->QMmethod]); /* see enum.h */
 214     else
 215       fprintf(out,"END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
 216               eQMbasis_names[qm->QMbasis],
 217               eQMmethod_names[qm->QMmethod]); /* see enum.h */
 218     F77_FUNC(inigms,IMIGMS)();
 219   }
 220 }
 221
 222 real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
 223                  rvec f[], rvec fshift[])
 224 {
 225   /* do the actual QMMM calculation using GAMESS-UK. In this
 226    * implementation (3-2001) a system call is made to the GAMESS-UK
 227    * binary. Now we are working to get the electron integral, SCF, and
 228    * gradient routines linked directly
 229    */
 230   int
 231     i,j,rank;
 232   real
 233     QMener=0.0,*qmgrad,*mmgrad,*mmcrd,*qmcrd,energy;
 234   t_QMMMrec
 235     *qr;
 236
 237   /* copy the QMMMrec pointer */
 238   qr = fr->qr;
 239   snew(qmcrd, 3*(qm->nrQMatoms));
 240   snew(mmcrd,3*(mm->nrMMatoms));
 241   snew(qmgrad,3*(qm->nrQMatoms));
 242   snew(mmgrad,3*(mm->nrMMatoms));
 243
 244   /* copy the data from qr into the arrays that are going to be used
 245    * in the fortran routines of gamess
 246    */
 247   for(i=0;i<qm->nrQMatoms;i++){
 248     for (j=0;j<DIM;j++){
 249       qmcrd[DIM*i+j] = 1/BOHR2NM*qm->xQM[i][j];
 250     }
 251   }
 252   for(i=0;i<mm->nrMMatoms;i++){
 253     for (j=0;j<DIM;j++){
 254       mmcrd[DIM*i+j] = 1/BOHR2NM*mm->xMM[i][j];
 255     }
 256   }
 257   for (i=0;i<3*qm->nrQMatoms;i+=3){
 258     fprintf(stderr,"%8.5f, %8.5f, %8.5f\n",
 259             qmcrd[i],
 260             qmcrd[i+1],
 261             qmcrd[i+2]);
 262   }
 263
 264   F77_FUNC(grads,GRADS)(&qm->nrQMatoms,qmcrd,&mm->nrMMatoms,mm->MMcharges,
 265                         mmcrd,qmgrad,mmgrad,&energy);
 266
 267   for(i=0;i<qm->nrQMatoms;i++){
 268     for(j=0;j<DIM;j++){
 269       f[i][j]      = HARTREE_BOHR2MD*qmgrad[3*i+j];
 270       fshift[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
 271     }
 272   }
 273   for(i=0;i<mm->nrMMatoms;i++){
 274     for(j=0;j<DIM;j++){
 275       f[i][j]      = HARTREE_BOHR2MD*mmgrad[3*i+j];
 276       fshift[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
 277     }
 278   }
 279   /* convert a.u to kJ/mol */
 280   QMener=energy*HARTREE2KJ*AVOGADRO;
 281   return(QMener);
 282 }
 283
 284 #else
 285 int
 286 gmx_qmmm_gamess_empty;
 287 #endif
 288