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, 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/ns.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
  72     F77_FUNC(inigms, IMIGMS) ();
  73
  74 static void
  75     F77_FUNC(grads, GRADS) (const int *nrqmat, real *qmcrd, const int *nrmmat, const real *mmchrg,
  76                             real *mmcrd, real *qmgrad, real *mmgrad, real *energy);
  77
  78 #if !GMX_QMMM_GAMESS
  79 // Stub definitions to make compilation succeed when not configured
  80 // for GAMESS support. In that case, the module gives a fatal error
  81 // when the initialization function is called, so there is no need to
  82 // issue fatal errors here, because that introduces problems with
  83 // tools suggesting and prohibiting noreturn attributes.
  84
  85 void F77_FUNC(inigms, IMIGMS) ()
  86 {
  87 };
  88 // NOLINTNEXTLINE(readability-named-parameter)
  89 void F77_FUNC(grads, GRADS) (const int *, real *, const int *,
  90                              const real *, real *, real *,
  91                              real *, real *)
  92 {
  93 };
  94 #endif
  95
  96 void init_gamess(const t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
  97 {
  98     /* it works hopelessly complicated :-)
  99      * first a file is written. Then the standard gamess input/output
 100      * routine is called (no system()!) to set up all fortran arrays.
 101      * this routine writes a punch file, like in a normal gamess run.
 102      * via this punch file the other games routines, needed for gradient
 103      * and energy evaluations are called. This setup works fine for
 104      * dynamics simulations. 7-6-2002 (London)
 105      */
 106     int
 107         i, j;
 108     FILE
 109        *out;
 110     char
 111         periodic_system[37][3] = {
 112         "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ",
 113         "O ", "F ", "Ne", "Na", "Mg", "Al", "Si", "P ",
 114         "S ", "Cl", "Ar", "K ", "Ca", "Sc", "Ti", "V ",
 115         "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga",
 116         "Ge", "As", "Se", "Br", "Kr"
 117     };
 118     if (!GMX_QMMM_GAMESS)
 119     {
 120         gmx_fatal(FARGS, "Cannot call GAMESS unless linked against it. Use cmake -DGMX_QMMM_PROGRAM=GAMESS, and ensure that linking will work correctly.");
 121     }
 122
 123     if (PAR(cr))
 124     {
 125
 126         if (MASTER(cr))
 127         {
 128             out = fopen("FOR009", "w");
 129             /* of these options I am not completely sure....  the overall
 130              * preformance on more than 4 cpu's is rather poor at the moment.
 131              */
 132             fprintf(out, "memory 48000000\nPARALLEL IOMODE SCREENED\n");
 133             fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
 134                     qm->nelectrons, qm->multiplicity);
 135             for (i = 0; i < qm->nrQMatoms; i++)
 136             {
 137 #ifdef DOUBLE
 138                 fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n",
 139                         i/2.,
 140                         i/3.,
 141                         i/4.,
 142                         qm->atomicnumberQM[i]*1.0,
 143                         periodic_system[qm->atomicnumberQM[i]]);
 144 #else
 145                 fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  %2s\n",
 146                         i/2.,
 147                         i/3.,
 148                         i/4.,
 149                         qm->atomicnumberQM[i]*1.0,
 150                         periodic_system[qm->atomicnumberQM[i]]);
 151 #endif
 152             }
 153             if (mm->nrMMatoms)
 154             {
 155                 for (j = i; j < i+2; j++)
 156                 {
 157 #ifdef DOUBLE
 158                     fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 159                             j/5.,
 160                             j/6.,
 161                             j/7.,
 162                             1.0);
 163 #else
 164                     fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 165                             j/5.,
 166                             j/6.,
 167                             j/7.,
 168                             2.0);
 169 #endif
 170                 }
 171             }
 172             fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 173                     eQMbasis_names[qm->QMbasis],
 174                     eQMmethod_names[qm->QMmethod]); /* see enum.h */
 175             fclose(out);
 176         }
 177         gmx_barrier(cr);
 178         F77_FUNC(inigms, IMIGMS) ();
 179     }
 180     else /* normal serial run */
 181
 182     {
 183         out = fopen("FOR009", "w");
 184         /* of these options I am not completely sure....  the overall
 185          * preformance on more than 4 cpu's is rather poor at the moment.
 186          */
 187         fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
 188                 qm->nelectrons, qm->multiplicity);
 189         for (i = 0; i < qm->nrQMatoms; i++)
 190         {
 191 #ifdef DOUBLE
 192             fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  %2s\n",
 193                     i/2.,
 194                     i/3.,
 195                     i/4.,
 196                     qm->atomicnumberQM[i]*1.0,
 197                     periodic_system[qm->atomicnumberQM[i]]);
 198 #else
 199             fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  %2s\n",
 200                     i/2.,
 201                     i/3.,
 202                     i/4.,
 203                     qm->atomicnumberQM[i]*1.0,
 204                     periodic_system[qm->atomicnumberQM[i]]);
 205 #endif
 206         }
 207         if (mm->nrMMatoms)
 208         {
 209             for (j = i; j < i+2; j++)
 210             {
 211 #ifdef DOUBLE
 212                 fprintf(out, "%10.7lf  %10.7lf  %10.7lf  %5.3lf  BQ\n",
 213                         j/5.,
 214                         j/6.,
 215                         j/7.,
 216                         1.0);
 217 #else
 218                 fprintf(out, "%10.7f  %10.7f  %10.7f  %5.3f  BQ\n",
 219                         j/5.,
 220                         j/6.,
 221                         j/7.,
 222                         2.0);
 223 #endif
 224             }
 225         }
 226         fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
 227                 eQMbasis_names[qm->QMbasis],
 228                 eQMmethod_names[qm->QMmethod]); /* see enum.h */
 229         F77_FUNC(inigms, IMIGMS) ();
 230     }
 231 }
 232
 233 real call_gamess(const t_QMrec *qm, const t_MMrec *mm,
 234                  rvec f[], rvec fshift[])
 235 {
 236     /* do the actual QMMM calculation using GAMESS-UK. In this
 237      * implementation (3-2001) a system call is made to the GAMESS-UK
 238      * binary. Now we are working to get the electron integral, SCF, and
 239      * gradient routines linked directly
 240      */
 241     int
 242         i, j;
 243     real
 244         QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy = 0;
 245
 246     snew(qmcrd, 3*(qm->nrQMatoms));
 247     snew(mmcrd, 3*(mm->nrMMatoms));
 248     snew(qmgrad, 3*(qm->nrQMatoms));
 249     snew(mmgrad, 3*(mm->nrMMatoms));
 250
 251     /* copy the data from qr into the arrays that are going to be used
 252      * in the fortran routines of gamess
 253      */
 254     for (i = 0; i < qm->nrQMatoms; i++)
 255     {
 256         for (j = 0; j < DIM; j++)
 257         {
 258             qmcrd[DIM*i+j] = 1/BOHR2NM*qm->xQM[i][j];
 259         }
 260     }
 261     for (i = 0; i < mm->nrMMatoms; i++)
 262     {
 263         for (j = 0; j < DIM; j++)
 264         {
 265             mmcrd[DIM*i+j] = 1/BOHR2NM*mm->xMM[i][j];
 266         }
 267     }
 268     for (i = 0; i < 3*qm->nrQMatoms; i += 3)
 269     {
 270         fprintf(stderr, "%8.5f, %8.5f, %8.5f\n",
 271                 qmcrd[i],
 272                 qmcrd[i+1],
 273                 qmcrd[i+2]);
 274     }
 275
 276     F77_FUNC(grads, GRADS) (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mm->MMcharges,
 277                             mmcrd, qmgrad, mmgrad, &energy);
 278
 279     for (i = 0; i < qm->nrQMatoms; i++)
 280     {
 281         for (j = 0; j < DIM; j++)
 282         {
 283             f[i][j]      = HARTREE_BOHR2MD*qmgrad[3*i+j];
 284             fshift[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
 285         }
 286     }
 287     for (i = 0; i < mm->nrMMatoms; i++)
 288     {
 289         for (j = 0; j < DIM; j++)
 290         {
 291             f[i][j]      = HARTREE_BOHR2MD*mmgrad[3*i+j];
 292             fshift[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
 293         }
 294     }
 295     /* convert a.u to kJ/mol */
 296     QMener = energy*HARTREE2KJ*AVOGADRO;
 297     return(QMener);
 298 }
 299
 300 #pragma GCC diagnostic pop