src/gromacs/mdlib/qm_orca.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-2008, 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_orca.h"
  40
  41 #include "config.h"
  42
  43 #include <cmath>
  44 #include <cstdio>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include "gromacs/fileio/confio.h"
  49 #include "gromacs/gmxlib/network.h"
  50 #include "gromacs/gmxlib/nrnb.h"
  51 #include "gromacs/math/units.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/mdlib/qmmm.h"
  54 #include "gromacs/mdtypes/forcerec.h"
  55 #include "gromacs/mdtypes/md_enums.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 // When not built in a configuration with QMMM support, much of this
  60 // code is unreachable by design. Tell clang not to warn about it.
  61 #pragma GCC diagnostic push
  62 #pragma GCC diagnostic ignored "-Wmissing-noreturn"
  63
  64 /* ORCA interface routines */
  65
  66 void init_orca(t_QMrec *qm)
  67 {
  68     char *buf;
  69     snew(buf, 200);
  70
  71     if (!GMX_QMMM_ORCA)
  72     {
  73         gmx_fatal(FARGS, "Cannot call ORCA unless linked against it. Use cmake -DGMX_QMMM_PROGRAM=ORCA, and ensure that linking will work correctly.");
  74     }
  75
  76     /* ORCA settings on the system */
  77     buf = getenv("GMX_QM_ORCA_BASENAME");
  78     if (buf)
  79     {
  80         snew(qm->orca_basename, 200);
  81         sscanf(buf, "%s", qm->orca_basename);
  82     }
  83     else
  84     {
  85         gmx_fatal(FARGS, "$GMX_QM_ORCA_BASENAME is not set\n");
  86     }
  87
  88     /* ORCA directory on the system */
  89     snew(buf, 200);
  90     buf = getenv("GMX_ORCA_PATH");
  91
  92     if (buf)
  93     {
  94         snew(qm->orca_dir, 200);
  95         sscanf(buf, "%s", qm->orca_dir);
  96     }
  97     else
  98     {
  99         gmx_fatal(FARGS, "$GMX_ORCA_PATH not set, check manual\n");
 100     }
 101
 102     fprintf(stderr, "Setting ORCA path to: %s...\n", qm->orca_dir);
 103     fprintf(stderr, "ORCA initialised...\n\n");
 104     /* since we append the output to the BASENAME.out file,
 105        we should delete an existent old out-file here. */
 106     sprintf(buf, "%s.out", qm->orca_basename);
 107     remove(buf);
 108 }
 109
 110
 111 static void write_orca_input(const t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
 112 {
 113     int        i;
 114     t_QMMMrec *QMMMrec;
 115     FILE      *out, *pcFile, *addInputFile;
 116     char      *buf, *orcaInput, *addInputFilename, *pcFilename;
 117
 118     QMMMrec = fr->qr;
 119
 120     /* write the first part of the input-file */
 121     snew(orcaInput, 200);
 122     sprintf(orcaInput, "%s.inp", qm->orca_basename);
 123     out = fopen(orcaInput, "w");
 124
 125     snew(addInputFilename, 200);
 126     sprintf(addInputFilename, "%s.ORCAINFO", qm->orca_basename);
 127     addInputFile = fopen(addInputFilename, "r");
 128
 129     fprintf(out, "#input-file generated by GROMACS\n");
 130
 131     fprintf(out, "!EnGrad TightSCF\n");
 132
 133     /* here we include the insertion of the additional orca-input */
 134     snew(buf, 200);
 135     if (addInputFile != nullptr)
 136     {
 137         while (!feof(addInputFile))
 138         {
 139             if (fgets(buf, 200, addInputFile) != nullptr)
 140             {
 141                 fputs(buf, out);
 142             }
 143         }
 144     }
 145     else
 146     {
 147         gmx_fatal(FARGS, "No information on the calculation given in %s\n", addInputFilename);
 148     }
 149
 150     fclose(addInputFile);
 151
 152     /* write charge and multiplicity */
 153     fprintf(out, "*xyz %2d%2d\n", qm->QMcharge, qm->multiplicity);
 154
 155     /* write the QM coordinates */
 156     for (i = 0; i < qm->nrQMatoms; i++)
 157     {
 158         int atomNr;
 159         if (qm->atomicnumberQM[i] == 0)
 160         {
 161             atomNr = 1;
 162         }
 163         else
 164         {
 165             atomNr = qm->atomicnumberQM[i];
 166         }
 167         fprintf(out, "%3d %10.7f  %10.7f  %10.7f\n",
 168                 atomNr,
 169                 qm->xQM[i][XX]/0.1,
 170                 qm->xQM[i][YY]/0.1,
 171                 qm->xQM[i][ZZ]/0.1);
 172     }
 173     fprintf(out, "*\n");
 174
 175     /* write the MM point charge data */
 176     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 177     {
 178         /* name of the point charge file */
 179         snew(pcFilename, 200);
 180         sprintf(pcFilename, "%s.pc", qm->orca_basename);
 181         fprintf(out, "%s%s%s\n", "%pointcharges \"", pcFilename, "\"");
 182         pcFile = fopen(pcFilename, "w");
 183         fprintf(pcFile, "%d\n", mm->nrMMatoms);
 184         for (i = 0; i < mm->nrMMatoms; i++)
 185         {
 186             fprintf(pcFile, "%8.4f %10.7f  %10.7f  %10.7f\n",
 187                     mm->MMcharges[i],
 188                     mm->xMM[i][XX]/0.1,
 189                     mm->xMM[i][YY]/0.1,
 190                     mm->xMM[i][ZZ]/0.1);
 191         }
 192         fprintf(pcFile, "\n");
 193         fclose(pcFile);
 194     }
 195     fprintf(out, "\n");
 196
 197     fclose(out);
 198 }  /* write_orca_input */
 199
 200 static real read_orca_output(rvec QMgrad[], rvec MMgrad[], const t_forcerec *fr,
 201                              t_QMrec *qm, t_MMrec *mm)
 202 {
 203     int
 204         i, j;
 205     char
 206         buf[300], orca_pcgradFilename[300], orca_engradFilename[300];
 207     real
 208         QMener;
 209     FILE
 210        *pcgrad, *engrad;
 211     int k;
 212     t_QMMMrec
 213        *QMMMrec;
 214     QMMMrec = fr->qr;
 215
 216     /* the energy and gradients for the QM part are stored in the engrad file
 217      * and the gradients for the point charges are stored in the pc file.
 218      */
 219     sprintf(orca_engradFilename, "%s.engrad", qm->orca_basename);
 220     engrad = fopen(orca_engradFilename, "r");
 221     /* we read the energy and the gradient for the qm-atoms from the engrad file
 222      */
 223     /* we can skip the first seven lines
 224      */
 225     for (j = 0; j < 7; j++)
 226     {
 227         if (fgets(buf, 300, engrad) == nullptr)
 228         {
 229             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 230         }
 231     }
 232     /* now comes the energy
 233      */
 234     if (fgets(buf, 300, engrad) == nullptr)
 235     {
 236         gmx_fatal(FARGS, "Unexpected end of ORCA output");
 237     }
 238 #if GMX_DOUBLE
 239     sscanf(buf, "%lf\n", &QMener);
 240 #else
 241     sscanf(buf, "%f\n", &QMener);
 242 #endif
 243     /* we can skip the next three lines
 244      */
 245     for (j = 0; j < 3; j++)
 246     {
 247         if (fgets(buf, 300, engrad) == nullptr)
 248         {
 249             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 250         }
 251     }
 252     /* next lines contain the gradients of the QM atoms
 253      * now comes the gradient, one value per line:
 254      * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
 255      */
 256
 257     for (i = 0; i < 3*qm->nrQMatoms; i++)
 258     {
 259         k = i/3;
 260         if (fgets(buf, 300, engrad) == nullptr)
 261         {
 262             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 263         }
 264 #if GMX_DOUBLE
 265         if (i%3 == 0)
 266         {
 267             sscanf(buf, "%lf\n", &QMgrad[k][XX]);
 268         }
 269         else if (i%3 == 1)
 270         {
 271             sscanf(buf, "%lf\n", &QMgrad[k][YY]);
 272         }
 273         else if (i%3 == 2)
 274         {
 275             sscanf(buf, "%lf\n", &QMgrad[k][ZZ]);
 276         }
 277 #else
 278         if (i%3 == 0)
 279         {
 280             sscanf(buf, "%f\n", &QMgrad[k][XX]);
 281         }
 282         else if (i%3 == 1)
 283         {
 284             sscanf(buf, "%f\n", &QMgrad[k][YY]);
 285         }
 286         else if (i%3 == 2)
 287         {
 288             sscanf(buf, "%f\n", &QMgrad[k][ZZ]);
 289         }
 290 #endif
 291     }
 292     fclose(engrad);
 293     /* write the MM point charge data
 294      */
 295     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 296     {
 297         sprintf(orca_pcgradFilename, "%s.pcgrad", qm->orca_basename);
 298         pcgrad = fopen(orca_pcgradFilename, "r");
 299
 300         /* we read the gradient for the mm-atoms from the pcgrad file
 301          */
 302         /* we can skip the first line
 303          */
 304         if (fgets(buf, 300, pcgrad) == nullptr)
 305         {
 306             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 307         }
 308         for (i = 0; i < mm->nrMMatoms; i++)
 309         {
 310             if (fgets(buf, 300, pcgrad) == nullptr)
 311             {
 312                 gmx_fatal(FARGS, "Unexpected end of ORCA output");
 313             }
 314     #if GMX_DOUBLE
 315             sscanf(buf, "%lf%lf%lf\n",
 316                    &MMgrad[i][XX],
 317                    &MMgrad[i][YY],
 318                    &MMgrad[i][ZZ]);
 319     #else
 320             sscanf(buf, "%f%f%f\n",
 321                    &MMgrad[i][XX],
 322                    &MMgrad[i][YY],
 323                    &MMgrad[i][ZZ]);
 324     #endif
 325         }
 326         fclose(pcgrad);
 327     }
 328     return(QMener);
 329 }
 330
 331 static void do_orca(char *orca_dir, char *basename)
 332 {
 333
 334     /* make the call to the orca binary through system()
 335      * The location of the binary is set through the
 336      * environment.
 337      */
 338     char
 339         buf[100];
 340     sprintf(buf, "%s/%s %s.inp >> %s.out",
 341             orca_dir,
 342             "orca",
 343             basename,
 344             basename);
 345     fprintf(stderr, "Calling '%s'\n", buf);
 346     if (system(buf) != 0)
 347     {
 348         gmx_fatal(FARGS, "Call to '%s' failed\n", buf);
 349     }
 350 }
 351
 352 real call_orca(const t_forcerec *fr,
 353                t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 354 {
 355     /* normal orca jobs */
 356     static int
 357         step = 0;
 358     int
 359         i, j;
 360     real
 361         QMener;
 362     rvec
 363        *QMgrad, *MMgrad;
 364     char
 365        *exe;
 366
 367     snew(exe, 30);
 368     sprintf(exe, "%s", "orca");
 369     snew(QMgrad, qm->nrQMatoms);
 370     snew(MMgrad, mm->nrMMatoms);
 371
 372     write_orca_input(fr, qm, mm);
 373     do_orca(qm->orca_dir, qm->orca_basename);
 374     QMener = read_orca_output(QMgrad, MMgrad, fr, qm, mm);
 375     /* put the QMMM forces in the force array and to the fshift
 376      */
 377     for (i = 0; i < qm->nrQMatoms; i++)
 378     {
 379         for (j = 0; j < DIM; j++)
 380         {
 381             f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];
 382             fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
 383         }
 384     }
 385     for (i = 0; i < mm->nrMMatoms; i++)
 386     {
 387         for (j = 0; j < DIM; j++)
 388         {
 389             f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];
 390             fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
 391         }
 392     }
 393     QMener = QMener*HARTREE2KJ*AVOGADRO;
 394     step++;
 395     free(exe);
 396     return(QMener);
 397 } /* call_orca */
 398
 399 /* end of orca sub routines */
 400
 401 #pragma GCC diagnostic pop