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