src/gromacs/mdlib/qm_gamess.cpp

   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,2015,2017,2018,2019, 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 #include "gmxpre.h"
  38
  39 #include "qm_gamess.h"
  40
  41 #include "config.h"
  42
  43 #include <stdio.h>
  44 #include <stdlib.h>
  45 #include <string.h>
  46
  47 #include <cmath>
  48
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/gmxlib/network.h"
  51 #include "gromacs/gmxlib/nrnb.h"
  52 #include "gromacs/math/units.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/mdlib/qmmm.h"
  55 #include "gromacs/mdtypes/commrec.h"
  56 #include "gromacs/mdtypes/forcerec.h"
  57 #include "gromacs/mdtypes/md_enums.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 // When not built in a configuration with QMMM support, much of this
  62 // code is unreachable by design. Tell clang not to warn about it.
  63 #pragma GCC diagnostic push
  64 #pragma GCC diagnostic ignored "-Wmissing-noreturn"
  65
  66 /* QMMM sub routines */
  67 /* mopac interface routines */
  68
  69
  70 static void F77_FUNC(inigms, IMIGMS)();
  71
  72 static void F77_FUNC(grads, GRADS)(const int*  nrqmat,
  73                                    real*       qmcrd,
  74                                    const int*  nrmmat,
  75                                    const real* mmchrg,
  76                                    real*       mmcrd,
  77                                    real*       qmgrad,
  78                                    real*       mmgrad,
  79                                    real*       energy);
  80
  81 #if !GMX_QMMM_GAMESS
  82 // Stub definitions to make compilation succeed when not configured
  83 // for GAMESS support. In that case, the module gives a fatal error
  84 // when the initialization function is called, so there is no need to
  85 // issue fatal errors here, because that introduces problems with
  86 // tools suggesting and prohibiting noreturn attributes.
  87
  88 void F77_FUNC(inigms, IMIGMS)(){};
  89 // NOLINTNEXTLINE(readability-named-parameter)
  90 void F77_FUNC(grads, GRADS)(const int*, real*, const int*, const real*, real*, real*, real*, real*){};
  91 #endif
  92
  93 void init_gamess(const t_commrec* cr, t_QMrec* qm, t_MMrec* mm)
  94 {
  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   i, j;
 104     FILE* out;
 105     char  periodic_system[37][3] = { "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ", "O ", "F ",
 106                                     "Ne", "Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar", "K ",
 107                                     "Ca", "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu",
 108                                     "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr" };
 109     if (!GMX_QMMM_GAMESS)
 110     {
 111         gmx_fatal(FARGS,
 112                   "Cannot call GAMESS unless linked against it. Use cmake "
 113                   "-DGMX_QMMM_PROGRAM=GAMESS, and ensure that linking will work correctly.");
 114     }
 115
 116     if (PAR(cr))
 117     {
 118
 119         if (MASTER(cr))
 120         {
 121             out = fopen("FOR009", "w");
 122             /* of these options I am not completely sure....  the overall
 123              * preformance on more than 4 cpu's is rather poor at the moment.
 124              */
 125             fprintf(out, "memory 48000000\nPARALLEL IOMODE SCREENED\n");
 126             fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n", qm->nelectrons,
 127                     qm->multiplicity);
 128             for (i = 0; i < qm->nrQMatoms; i++)
 129             {
 130 #ifdef DOUBLE
 131                 fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n", i / 2., i / 3., i / 4.,
 132                         qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]);
 133 #else
 134                 fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  %2s\n", i / 2., i / 3., i / 4.,
 135                         qm->atomicnumberQM[i] * 1.0, 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", j / 5., j / 6., j / 7., 1.0);
 144 #else
 145                     fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  BQ\n", j / 5., j / 6., j / 7., 2.0);
 146 #endif
 147                 }
 148             }
 149             fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 150                     eQMbasis_names[qm->QMbasis], eQMmethod_names[qm->QMmethod]); /* see enum.h */
 151             fclose(out);
 152         }
 153         gmx_barrier(cr);
 154         F77_FUNC(inigms, IMIGMS)();
 155     }
 156     else /* normal serial run */
 157
 158     {
 159         out = fopen("FOR009", "w");
 160         /* of these options I am not completely sure....  the overall
 161          * preformance on more than 4 cpu's is rather poor at the moment.
 162          */
 163         fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n", qm->nelectrons,
 164                 qm->multiplicity);
 165         for (i = 0; i < qm->nrQMatoms; i++)
 166         {
 167 #ifdef DOUBLE
 168             fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n", i / 2., i / 3., i / 4.,
 169                     qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]);
 170 #else
 171             fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  %2s\n", i / 2., i / 3., i / 4.,
 172                     qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]);
 173 #endif
 174         }
 175         if (mm->nrMMatoms)
 176         {
 177             for (j = i; j < i + 2; j++)
 178             {
 179 #ifdef DOUBLE
 180                 fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n", j / 5., j / 6., j / 7., 1.0);
 181 #else
 182                 fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  BQ\n", j / 5., j / 6., j / 7., 2.0);
 183 #endif
 184             }
 185         }
 186         fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n", eQMbasis_names[qm->QMbasis],
 187                 eQMmethod_names[qm->QMmethod]); /* see enum.h */
 188         F77_FUNC(inigms, IMIGMS)();
 189     }
 190 }
 191
 192 real call_gamess(const t_QMrec* qm, const t_MMrec* mm, rvec f[], rvec fshift[])
 193 {
 194     /* do the actual QMMM calculation using GAMESS-UK. In this
 195      * implementation (3-2001) a system call is made to the GAMESS-UK
 196      * binary. Now we are working to get the electron integral, SCF, and
 197      * gradient routines linked directly
 198      */
 199     int  i, j;
 200     real QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy = 0;
 201
 202     snew(qmcrd, 3 * (qm->nrQMatoms));
 203     snew(mmcrd, 3 * (mm->nrMMatoms));
 204     snew(qmgrad, 3 * (qm->nrQMatoms));
 205     snew(mmgrad, 3 * (mm->nrMMatoms));
 206
 207     /* copy the data from qr into the arrays that are going to be used
 208      * in the fortran routines of gamess
 209      */
 210     for (i = 0; i < qm->nrQMatoms; i++)
 211     {
 212         for (j = 0; j < DIM; j++)
 213         {
 214             qmcrd[DIM * i + j] = 1 / BOHR2NM * qm->xQM[i][j];
 215         }
 216     }
 217     for (i = 0; i < mm->nrMMatoms; i++)
 218     {
 219         for (j = 0; j < DIM; j++)
 220         {
 221             mmcrd[DIM * i + j] = 1 / BOHR2NM * mm->xMM[i][j];
 222         }
 223     }
 224     for (i = 0; i < 3 * qm->nrQMatoms; i += 3)
 225     {
 226         fprintf(stderr, "%8.5f, %8.5f, %8.5f\n", qmcrd[i], qmcrd[i + 1], qmcrd[i + 2]);
 227     }
 228
 229     F77_FUNC(grads, GRADS)
 230     (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mm->MMcharges, mmcrd, qmgrad, mmgrad, &energy);
 231
 232     for (i = 0; i < qm->nrQMatoms; i++)
 233     {
 234         for (j = 0; j < DIM; j++)
 235         {
 236             f[i][j]      = HARTREE_BOHR2MD * qmgrad[3 * i + j];
 237             fshift[i][j] = HARTREE_BOHR2MD * qmgrad[3 * i + j];
 238         }
 239     }
 240     for (i = 0; i < mm->nrMMatoms; i++)
 241     {
 242         for (j = 0; j < DIM; j++)
 243         {
 244             f[i][j]      = HARTREE_BOHR2MD * mmgrad[3 * i + j];
 245             fshift[i][j] = HARTREE_BOHR2MD * mmgrad[3 * i + j];
 246         }
 247     }
 248     /* convert a.u to kJ/mol */
 249     QMener = energy * HARTREE2KJ * AVOGADRO;
 250     return (QMener);
 251 }
 252
 253 #pragma GCC diagnostic pop