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