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 "qmmm.h"
  61 #include <stdio.h>
  62 #include <string.h>
  63 #include "gmx_fatal.h"
  64 #include "typedefs.h"
  65 #include <stdlib.h>
  66
  67
  68 /* QMMM sub routines */
  69 /* mopac interface routines */
  70
  71
  72 void
  73 F77_FUNC(inigms, IMIGMS) (void);
  74
  75 void
  76 F77_FUNC(endgms, ENDGMS) (void);
  77
  78 void
  79 F77_FUNC(grads, GRADS) (int *nrqmat, real *qmcrd, int *nrmmat, real *mmchrg,
  80                         real *mmcrd, real *qmgrad, real *mmgrad, real *energy);
  81
  82
  83
  84 void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
  85 {
  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] = {
 100         "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ",
 101         "O ", "F ", "Ne", "Na", "Mg", "Al", "Si", "P ",
 102         "S ", "Cl", "Ar", "K ", "Ca", "Sc", "Ti", "V ",
 103         "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga",
 104         "Ge", "As", "Se", "Br", "Kr"
 105     };
 106
 107     if (PAR(cr))
 108     {
 109
 110         if (MASTER(cr))
 111         {
 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             {
 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             {
 139                 for (j = i; j < i+2; j++)
 140                 {
 141 #ifdef DOUBLE
 142                     fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 143                             j/5.,
 144                             j/6.,
 145                             j/7.,
 146                             1.0);
 147 #else
 148                     fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 149                             j/5.,
 150                             j/6.,
 151                             j/7.,
 152                             2.0);
 153 #endif
 154                 }
 155             }
 156             if (!qm->bTS)
 157             {
 158                 fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 159                         eQMbasis_names[qm->QMbasis],
 160                         eQMmethod_names[qm->QMmethod]); /* see enum.h */
 161             }
 162             else
 163             {
 164                 fprintf(out, "END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
 165                         eQMbasis_names[qm->QMbasis],
 166                         eQMmethod_names[qm->QMmethod]); /* see enum.h */
 167             }
 168             fclose(out);
 169         }
 170         gmx_barrier(cr);
 171         F77_FUNC(inigms, IMIGMS) ();
 172     }
 173     else /* normal serial run */
 174
 175     {
 176         out = fopen("FOR009", "w");
 177         /* of these options I am not completely sure....  the overall
 178          * preformance on more than 4 cpu's is rather poor at the moment.
 179          */
 180         fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
 181                 qm->nelectrons, qm->multiplicity);
 182         for (i = 0; i < qm->nrQMatoms; i++)
 183         {
 184 #ifdef DOUBLE
 185             fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n",
 186                     i/2.,
 187                     i/3.,
 188                     i/4.,
 189                     qm->atomicnumberQM[i]*1.0,
 190                     periodic_system[qm->atomicnumberQM[i]]);
 191 #else
 192             fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  %2s\n",
 193                     i/2.,
 194                     i/3.,
 195                     i/4.,
 196                     qm->atomicnumberQM[i]*1.0,
 197                     periodic_system[qm->atomicnumberQM[i]]);
 198 #endif
 199         }
 200         if (mm->nrMMatoms)
 201         {
 202             for (j = i; j < i+2; j++)
 203             {
 204 #ifdef DOUBLE
 205                 fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 206                         j/5.,
 207                         j/6.,
 208                         j/7.,
 209                         1.0);
 210 #else
 211                 fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 212                         j/5.,
 213                         j/6.,
 214                         j/7.,
 215                         2.0);
 216 #endif
 217             }
 218         }
 219         if (!qm->bTS)
 220         {
 221             fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 222                     eQMbasis_names[qm->QMbasis],
 223                     eQMmethod_names[qm->QMmethod]); /* see enum.h */
 224         }
 225         else
 226         {
 227             fprintf(out, "END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
 228                     eQMbasis_names[qm->QMbasis],
 229                     eQMmethod_names[qm->QMmethod]); /* see enum.h */
 230         }
 231         F77_FUNC(inigms, IMIGMS) ();
 232     }
 233 }
 234
 235 real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
 236                  rvec f[], rvec fshift[])
 237 {
 238     /* do the actual QMMM calculation using GAMESS-UK. In this
 239      * implementation (3-2001) a system call is made to the GAMESS-UK
 240      * binary. Now we are working to get the electron integral, SCF, and
 241      * gradient routines linked directly
 242      */
 243     int
 244         i, j, rank;
 245     real
 246         QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy;
 247     t_QMMMrec
 248        *qr;
 249
 250     /* copy the QMMMrec pointer */
 251     qr = fr->qr;
 252     snew(qmcrd, 3*(qm->nrQMatoms));
 253     snew(mmcrd, 3*(mm->nrMMatoms));
 254     snew(qmgrad, 3*(qm->nrQMatoms));
 255     snew(mmgrad, 3*(mm->nrMMatoms));
 256
 257     /* copy the data from qr into the arrays that are going to be used
 258      * in the fortran routines of gamess
 259      */
 260     for (i = 0; i < qm->nrQMatoms; i++)
 261     {
 262         for (j = 0; j < DIM; j++)
 263         {
 264             qmcrd[DIM*i+j] = 1/BOHR2NM*qm->xQM[i][j];
 265         }
 266     }
 267     for (i = 0; i < mm->nrMMatoms; i++)
 268     {
 269         for (j = 0; j < DIM; j++)
 270         {
 271             mmcrd[DIM*i+j] = 1/BOHR2NM*mm->xMM[i][j];
 272         }
 273     }
 274     for (i = 0; i < 3*qm->nrQMatoms; i += 3)
 275     {
 276         fprintf(stderr, "%8.5f, %8.5f, %8.5f\n",
 277                 qmcrd[i],
 278                 qmcrd[i+1],
 279                 qmcrd[i+2]);
 280     }
 281
 282     F77_FUNC(grads, GRADS) (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mm->MMcharges,
 283                             mmcrd, qmgrad, mmgrad, &energy);
 284
 285     for (i = 0; i < qm->nrQMatoms; i++)
 286     {
 287         for (j = 0; j < DIM; j++)
 288         {
 289             f[i][j]      = HARTREE_BOHR2MD*qmgrad[3*i+j];
 290             fshift[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
 291         }
 292     }
 293     for (i = 0; i < mm->nrMMatoms; i++)
 294     {
 295         for (j = 0; j < DIM; j++)
 296         {
 297             f[i][j]      = HARTREE_BOHR2MD*mmgrad[3*i+j];
 298             fshift[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
 299         }
 300     }
 301     /* convert a.u to kJ/mol */
 302     QMener = energy*HARTREE2KJ*AVOGADRO;
 303     return(QMener);
 304 }
 305
 306 #else
 307 int
 308     gmx_qmmm_gamess_empty;
 309 #endif