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