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