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