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