src/mdlib/qm_gamess.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #ifdef GMX_QMMM_GAMESS
  43
  44 #include <math.h>
  45 #include "sysstuff.h"
  46 #include "typedefs.h"
  47 #include "macros.h"
  48 #include "smalloc.h"
  49 #include "assert.h"
  50 #include "physics.h"
  51 #include "macros.h"
  52 #include "vec.h"
  53 #include "force.h"
  54 #include "invblock.h"
  55 #include "confio.h"
  56 #include "names.h"
  57 #include "network.h"
  58 #include "pbc.h"
  59 #include "ns.h"
  60 #include "nrnb.h"
  61 #include "bondf.h"
  62 #include "mshift.h"
  63 #include "txtdump.h"
  64 #include "copyrite.h"
  65 #include "qmmm.h"
  66 #include <stdio.h>
  67 #include <string.h>
  68 #include "gmx_fatal.h"
  69 #include "typedefs.h"
  70 #include <stdlib.h>
  71
  72
  73 /* QMMM sub routines */
  74 /* mopac interface routines */
  75
  76
  77 void
  78 F77_FUNC(inigms,IMIGMS)(void);
  79
  80 void
  81 F77_FUNC(endgms,ENDGMS)(void);
  82
  83 void
  84 F77_FUNC(grads,GRADS)(int *nrqmat,real *qmcrd,int *nrmmat, real *mmchrg,
  85                       real *mmcrd, real *qmgrad,real *mmgrad, real *energy);
  86
  87
  88
  89 void init_gamess(t_commrec *cr,t_QMrec *qm, t_MMrec *mm){
  90   /* it works hopelessly complicated :-)
  91    * first a file is written. Then the standard gamess input/output
  92    * routine is called (no system()!) to set up all fortran arrays.
  93    * this routine writes a punch file, like in a normal gamess run.
  94    * via this punch file the other games routines, needed for gradient
  95    * and energy evaluations are called. This setup works fine for
  96    * dynamics simulations. 7-6-2002 (London)
  97    */
  98   int
  99     i,j,rank;
 100   FILE
 101     *out;
 102   char
 103     periodic_system[37][3]={"XX","H ","He","Li","Be","B ","C ","N ",
 104                             "O ","F ","Ne","Na","Mg","Al","Si","P ",
 105                             "S ","Cl","Ar","K ","Ca","Sc","Ti","V ",
 106                             "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga",
 107                             "Ge","As","Se","Br","Kr"};
 108
 109   if (PAR(cr)){
 110
 111     if MASTER(cr){
 112       out=fopen("FOR009","w");
 113       /* of these options I am not completely sure....  the overall
 114        * preformance on more than 4 cpu's is rather poor at the moment.
 115        */
 116       fprintf(out,"memory 48000000\nPARALLEL IOMODE SCREENED\n");
 117       fprintf(out,"ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
 118               qm->nelectrons,qm->multiplicity);
 119       for (i=0;i<qm->nrQMatoms;i++){
 120 #ifdef DOUBLE
 121         fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n",
 122                 i/2.,
 123                 i/3.,
 124                 i/4.,
 125                 qm->atomicnumberQM[i]*1.0,
 126                 periodic_system[qm->atomicnumberQM[i]]);
 127 #else
 128         fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  %2s\n",
 129                 i/2.,
 130                 i/3.,
 131                 i/4.,
 132                 qm->atomicnumberQM[i]*1.0,
 133                 periodic_system[qm->atomicnumberQM[i]]);
 134 #endif
 135       }
 136       if(mm->nrMMatoms){
 137         for (j=i;j<i+2;j++){
 138 #ifdef DOUBLE
 139           fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 140                   j/5.,
 141                   j/6.,
 142                   j/7.,
 143                   1.0);
 144 #else
 145           fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 146                   j/5.,
 147                   j/6.,
 148                   j/7.,
 149                   2.0);
 150 #endif
 151         }
 152       }
 153       if(!qm->bTS)
 154         fprintf(out,"END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 155                 eQMbasis_names[qm->QMbasis],
 156                 eQMmethod_names[qm->QMmethod]); /* see enum.h */
 157       else
 158         fprintf(out,"END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
 159                 eQMbasis_names[qm->QMbasis],
 160                 eQMmethod_names[qm->QMmethod]); /* see enum.h */
 161       fclose(out);
 162     }
 163     gmx_barrier(cr);
 164     F77_FUNC(inigms,IMIGMS)();
 165   }
 166   else{ /* normal serial run */
 167
 168     out=fopen("FOR009","w");
 169     /* of these options I am not completely sure....  the overall
 170      * preformance on more than 4 cpu's is rather poor at the moment.
 171      */
 172     fprintf(out,"ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
 173             qm->nelectrons,qm->multiplicity);
 174     for (i=0;i<qm->nrQMatoms;i++){
 175 #ifdef DOUBLE
 176       fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n",
 177               i/2.,
 178               i/3.,
 179               i/4.,
 180               qm->atomicnumberQM[i]*1.0,
 181               periodic_system[qm->atomicnumberQM[i]]);
 182 #else
 183       fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  %2s\n",
 184               i/2.,
 185               i/3.,
 186               i/4.,
 187               qm->atomicnumberQM[i]*1.0,
 188               periodic_system[qm->atomicnumberQM[i]]);
 189 #endif
 190     }
 191     if(mm->nrMMatoms){
 192       for (j=i;j<i+2;j++){
 193 #ifdef DOUBLE
 194         fprintf(out,"%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 195                 j/5.,
 196                 j/6.,
 197                 j/7.,
 198                 1.0);
 199 #else
 200         fprintf(out,"%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 201                 j/5.,
 202                 j/6.,
 203                 j/7.,
 204                 2.0);
 205 #endif
 206       }
 207     }
 208     if(!qm->bTS)
 209       fprintf(out,"END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 210               eQMbasis_names[qm->QMbasis],
 211               eQMmethod_names[qm->QMmethod]); /* see enum.h */
 212     else
 213       fprintf(out,"END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
 214               eQMbasis_names[qm->QMbasis],
 215               eQMmethod_names[qm->QMmethod]); /* see enum.h */
 216     F77_FUNC(inigms,IMIGMS)();
 217   }
 218 }
 219
 220 real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
 221                  rvec f[], rvec fshift[])
 222 {
 223   /* do the actual QMMM calculation using GAMESS-UK. In this
 224    * implementation (3-2001) a system call is made to the GAMESS-UK
 225    * binary. Now we are working to get the electron integral, SCF, and
 226    * gradient routines linked directly
 227    */
 228   int
 229     i,j,rank;
 230   real
 231     QMener=0.0,*qmgrad,*mmgrad,*mmcrd,*qmcrd,energy;
 232   t_QMMMrec
 233     *qr;
 234
 235   /* copy the QMMMrec pointer */
 236   qr = fr->qr;
 237   snew(qmcrd, 3*(qm->nrQMatoms));
 238   snew(mmcrd,3*(mm->nrMMatoms));
 239   snew(qmgrad,3*(qm->nrQMatoms));
 240   snew(mmgrad,3*(mm->nrMMatoms));
 241
 242   /* copy the data from qr into the arrays that are going to be used
 243    * in the fortran routines of gamess
 244    */
 245   for(i=0;i<qm->nrQMatoms;i++){
 246     for (j=0;j<DIM;j++){
 247       qmcrd[DIM*i+j] = 1/BOHR2NM*qm->xQM[i][j];
 248     }
 249   }
 250   for(i=0;i<mm->nrMMatoms;i++){
 251     for (j=0;j<DIM;j++){
 252       mmcrd[DIM*i+j] = 1/BOHR2NM*mm->xMM[i][j];
 253     }
 254   }
 255   for (i=0;i<3*qm->nrQMatoms;i+=3){
 256     fprintf(stderr,"%8.5f, %8.5f, %8.5f\n",
 257             qmcrd[i],
 258             qmcrd[i+1],
 259             qmcrd[i+2]);
 260   }
 261
 262   F77_FUNC(grads,GRADS)(&qm->nrQMatoms,qmcrd,&mm->nrMMatoms,mm->MMcharges,
 263                         mmcrd,qmgrad,mmgrad,&energy);
 264
 265   for(i=0;i<qm->nrQMatoms;i++){
 266     for(j=0;j<DIM;j++){
 267       f[i][j]      = HARTREE_BOHR2MD*qmgrad[3*i+j];
 268       fshift[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
 269     }
 270   }
 271   for(i=0;i<mm->nrMMatoms;i++){
 272     for(j=0;j<DIM;j++){
 273       f[i][j]      = HARTREE_BOHR2MD*mmgrad[3*i+j];
 274       fshift[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
 275     }
 276   }
 277   /* convert a.u to kJ/mol */
 278   QMener=energy*HARTREE2KJ*AVOGADRO;
 279   return(QMener);
 280 }
 281
 282 #else
 283 int
 284 gmx_qmmm_gamess_empty;
 285 #endif
 286