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