src/gromacs/mdlib/qm_orca.c

   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, 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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "macros.h"
  45 #include "gromacs/utility/smalloc.h"
  46 #include "physics.h"
  47 #include "macros.h"
  48 #include "vec.h"
  49 #include "force.h"
  50 #include "invblock.h"
  51 #include "gromacs/fileio/confio.h"
  52 #include "names.h"
  53 #include "network.h"
  54 #include "pbc.h"
  55 #include "ns.h"
  56 #include "nrnb.h"
  57 #include "bondf.h"
  58 #include "mshift.h"
  59 #include "txtdump.h"
  60 #include "qmmm.h"
  61 #include <stdio.h>
  62 #include <string.h>
  63 #include "gmx_fatal.h"
  64 #include "typedefs.h"
  65 #include <stdlib.h>
  66
  67 /* ORCA interface routines */
  68
  69 void init_orca(t_QMrec *qm)
  70 {
  71     char *buf;
  72     snew(buf, 200);
  73
  74     /* ORCA settings on the system */
  75     buf = getenv("BASENAME");
  76     if (buf)
  77     {
  78         snew(qm->orca_basename, 200);
  79         sscanf(buf, "%s", qm->orca_basename);
  80     }
  81     else
  82     {
  83         gmx_fatal(FARGS, "$BASENAME not set\n");
  84     }
  85
  86     /* ORCA directory on the system */
  87     snew(buf, 200);
  88     buf = getenv("ORCA_PATH");
  89
  90     if (buf)
  91     {
  92         snew(qm->orca_dir, 200);
  93         sscanf(buf, "%s", qm->orca_dir);
  94     }
  95     else
  96     {
  97         gmx_fatal(FARGS, "$ORCA_PATH not set, check manual\n");
  98     }
  99
 100     fprintf(stderr, "Setting ORCA path to: %s...\n", qm->orca_dir);
 101     fprintf(stderr, "ORCA initialised...\n\n");
 102     /* since we append the output to the BASENAME.out file,
 103        we should delete an existent old out-file here. */
 104     sprintf(buf, "%s.out", qm->orca_basename);
 105     remove(buf);
 106 }
 107
 108
 109 void write_orca_input(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
 110 {
 111     int        i;
 112     t_QMMMrec *QMMMrec;
 113     FILE      *out, *pcFile, *addInputFile, *LJCoeff;
 114     char      *buf, *orcaInput, *addInputFilename, *LJCoeffFilename, *pcFilename, *exclInName, *exclOutName;
 115
 116     QMMMrec = fr->qr;
 117
 118     /* write the first part of the input-file */
 119     snew(orcaInput, 200);
 120     sprintf(orcaInput, "%s.inp", qm->orca_basename);
 121     out = fopen(orcaInput, "w");
 122
 123     snew(addInputFilename, 200);
 124     sprintf(addInputFilename, "%s.ORCAINFO", qm->orca_basename);
 125     addInputFile = fopen(addInputFilename, "r");
 126
 127     fprintf(out, "#input-file generated by GROMACS\n");
 128
 129     if (qm->bTS)
 130     {
 131         fprintf(out, "!QMMMOpt TightSCF\n");
 132         fprintf(out, "%s\n", "%geom TS_Search EF end");
 133     }
 134     else if (qm->bOPT)
 135     {
 136         fprintf(out, "!QMMMOpt TightSCF\n");
 137     }
 138     else
 139     {
 140         fprintf(out, "!EnGrad TightSCF\n");
 141     }
 142
 143     /* here we include the insertion of the additional orca-input */
 144     snew(buf, 200);
 145     if (addInputFile != NULL)
 146     {
 147         while (!feof(addInputFile))
 148         {
 149             if (fgets(buf, 200, addInputFile) != NULL)
 150             {
 151                 fputs(buf, out);
 152             }
 153         }
 154     }
 155     else
 156     {
 157         gmx_fatal(FARGS, "No information on the calculation given in %s\n", addInputFilename);
 158     }
 159
 160     fclose(addInputFile);
 161
 162     if (qm->bTS || qm->bOPT)
 163     {
 164         /* freeze the frontier QM atoms and Link atoms. This is
 165          * important only if a full QM subsystem optimization is done
 166          * with a frozen MM environmeent. For dynamics, or gromacs's own
 167          * optimization routines this is not important.
 168          */
 169         /* ORCA reads the exclusions from LJCoeffFilename.Excl,
 170          * so we have to rename the file
 171          */
 172         int didStart = 0;
 173         for (i = 0; i < qm->nrQMatoms; i++)
 174         {
 175             if (qm->frontatoms[i])
 176             {
 177                 if (!didStart)
 178                 {
 179                     fprintf(out, "%s\n", "%geom");
 180                     fprintf(out, "   Constraints \n");
 181                     didStart = 1;
 182                 }
 183                 fprintf(out, "        {C %d C}\n", i); /* counting from 0 */
 184             }
 185         }
 186         if (didStart)
 187         {
 188             fprintf(out, "     end\n   end\n");
 189         }
 190         /* make a file with information on the C6 and C12 coefficients */
 191         if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 192         {
 193             snew(exclInName, 200);
 194             snew(exclOutName, 200);
 195             sprintf(exclInName, "QMMMexcl.dat");
 196             sprintf(exclOutName, "%s.LJ.Excl", qm->orca_basename);
 197             rename(exclInName, exclOutName);
 198             snew(LJCoeffFilename, 200);
 199             sprintf(LJCoeffFilename, "%s.LJ", qm->orca_basename);
 200             fprintf(out, "%s%s%s\n", "%LJCOEFFICIENTS \"", LJCoeffFilename, "\"");
 201             /* make a file with information on the C6 and C12 coefficients */
 202             LJCoeff = fopen(LJCoeffFilename, "w");
 203             fprintf(LJCoeff, "%d\n", qm->nrQMatoms);
 204             for (i = 0; i < qm->nrQMatoms; i++)
 205             {
 206 #ifdef GMX_DOUBLE
 207                 fprintf(LJCoeff, "%10.7lf  %10.7lf\n", qm->c6[i], qm->c12[i]);
 208 #else
 209                 fprintf(LJCoeff, "%10.7f  %10.7f\n", qm->c6[i], qm->c12[i]);
 210 #endif
 211             }
 212             fprintf(LJCoeff, "%d\n", mm->nrMMatoms);
 213             for (i = 0; i < mm->nrMMatoms; i++)
 214             {
 215 #ifdef GMX_DOUBLE
 216                 fprintf(LJCoeff, "%10.7lf  %10.7lf\n", mm->c6[i], mm->c12[i]);
 217 #else
 218                 fprintf(LJCoeff, "%10.7f  %10.7f\n", mm->c6[i], mm->c12[i]);
 219 #endif
 220             }
 221             fclose(LJCoeff);
 222         }
 223     }
 224
 225     /* write charge and multiplicity */
 226     fprintf(out, "*xyz %2d%2d\n", qm->QMcharge, qm->multiplicity);
 227
 228     /* write the QM coordinates */
 229     for (i = 0; i < qm->nrQMatoms; i++)
 230     {
 231         int atomNr;
 232         if (qm->atomicnumberQM[i] == 0)
 233         {
 234             atomNr = 1;
 235         }
 236         else
 237         {
 238             atomNr = qm->atomicnumberQM[i];
 239         }
 240 #ifdef GMX_DOUBLE
 241         fprintf(out, "%3d %10.7lf  %10.7lf  %10.7lf\n",
 242                 atomNr,
 243                 qm->xQM[i][XX]/0.1,
 244                 qm->xQM[i][YY]/0.1,
 245                 qm->xQM[i][ZZ]/0.1);
 246 #else
 247         fprintf(out, "%3d %10.7f  %10.7f  %10.7f\n",
 248                 atomNr,
 249                 qm->xQM[i][XX]/0.1,
 250                 qm->xQM[i][YY]/0.1,
 251                 qm->xQM[i][ZZ]/0.1);
 252 #endif
 253     }
 254     fprintf(out, "*\n");
 255
 256     /* write the MM point charge data */
 257     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 258     {
 259         /* name of the point charge file */
 260         snew(pcFilename, 200);
 261         sprintf(pcFilename, "%s.pc", qm->orca_basename);
 262         fprintf(out, "%s%s%s\n", "%pointcharges \"", pcFilename, "\"");
 263         pcFile = fopen(pcFilename, "w");
 264         fprintf(pcFile, "%d\n", mm->nrMMatoms);
 265         for (i = 0; i < mm->nrMMatoms; i++)
 266         {
 267 #ifdef GMX_DOUBLE
 268             fprintf(pcFile, "%8.4lf %10.7lf  %10.7lf  %10.7lf\n",
 269                     mm->MMcharges[i],
 270                     mm->xMM[i][XX]/0.1,
 271                     mm->xMM[i][YY]/0.1,
 272                     mm->xMM[i][ZZ]/0.1);
 273 #else
 274             fprintf(pcFile, "%8.4f %10.7f  %10.7f  %10.7f\n",
 275                     mm->MMcharges[i],
 276                     mm->xMM[i][XX]/0.1,
 277                     mm->xMM[i][YY]/0.1,
 278                     mm->xMM[i][ZZ]/0.1);
 279 #endif
 280         }
 281         fprintf(pcFile, "\n");
 282         fclose(pcFile);
 283     }
 284     fprintf(out, "\n");
 285
 286     fclose(out);
 287 }  /* write_orca_input */
 288
 289 real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr,
 290                       t_QMrec *qm, t_MMrec *mm)
 291 {
 292     int
 293         i, j, atnum;
 294     char
 295         buf[300], tmp[300], orca_xyzFilename[300], orca_pcgradFilename[300], orca_engradFilename[300];
 296     real
 297         QMener;
 298     FILE
 299        *xyz, *pcgrad, *engrad;
 300     int k;
 301     t_QMMMrec
 302        *QMMMrec;
 303     QMMMrec = fr->qr;
 304     /* in case of an optimization, the coordinates are printed in the
 305      * xyz file, the energy and gradients for the QM part are stored in the engrad file
 306      * and the gradients for the point charges are stored in the pc file.
 307      */
 308
 309     /* we need the new xyz coordinates of the QM atoms only for separate QM-optimization
 310      */
 311
 312     if (qm->bTS || qm->bOPT)
 313     {
 314         sprintf(orca_xyzFilename, "%s.xyz", qm->orca_basename);
 315         xyz = fopen(orca_xyzFilename, "r");
 316         if (fgets(buf, 300, xyz) == NULL)
 317         {
 318             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 319         }
 320         if (fgets(buf, 300, xyz) == NULL)
 321         {
 322             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 323         }
 324         for (i = 0; i < qm->nrQMatoms; i++)
 325         {
 326             if (fgets(buf, 300, xyz) == NULL)
 327             {
 328                 gmx_fatal(FARGS, "Unexpected end of ORCA output");
 329             }
 330 #ifdef GMX_DOUBLE
 331             sscanf(buf, "%s%lf%lf%lf\n",
 332                    tmp,
 333                    &qm->xQM[i][XX],
 334                    &qm->xQM[i][YY],
 335                    &qm->xQM[i][ZZ]);
 336 #else
 337             sscanf(buf, "%d%f%f%f\n",
 338                    &atnum,
 339                    &qm->xQM[i][XX],
 340                    &qm->xQM[i][YY],
 341                    &qm->xQM[i][ZZ]);
 342 #endif
 343             for (j = 0; j < DIM; j++)
 344             {
 345                 qm->xQM[i][j] *= 0.1;
 346             }
 347         }
 348         fclose(xyz);
 349     }
 350     sprintf(orca_engradFilename, "%s.engrad", qm->orca_basename);
 351     engrad = fopen(orca_engradFilename, "r");
 352     /* we read the energy and the gradient for the qm-atoms from the engrad file
 353      */
 354     /* we can skip the first seven lines
 355      */
 356     for (j = 0; j < 7; j++)
 357     {
 358         if (fgets(buf, 300, engrad) == NULL)
 359         {
 360             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 361         }
 362     }
 363     /* now comes the energy
 364      */
 365     if (fgets(buf, 300, engrad) == NULL)
 366     {
 367         gmx_fatal(FARGS, "Unexpected end of ORCA output");
 368     }
 369 #ifdef GMX_DOUBLE
 370     sscanf(buf, "%lf\n", &QMener);
 371 #else
 372     sscanf(buf, "%f\n", &QMener);
 373 #endif
 374     /* we can skip the next three lines
 375      */
 376     for (j = 0; j < 3; j++)
 377     {
 378         if (fgets(buf, 300, engrad) == NULL)
 379         {
 380             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 381         }
 382     }
 383     /* next lines contain the gradients of the QM atoms
 384      * now comes the gradient, one value per line:
 385      * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
 386      */
 387
 388     for (i = 0; i < 3*qm->nrQMatoms; i++)
 389     {
 390         k = i/3;
 391         if (fgets(buf, 300, engrad) == NULL)
 392         {
 393             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 394         }
 395 #ifdef GMX_DOUBLE
 396         if (i%3 == 0)
 397         {
 398             sscanf(buf, "%lf\n", &QMgrad[k][XX]);
 399         }
 400         else if (i%3 == 1)
 401         {
 402             sscanf(buf, "%lf\n", &QMgrad[k][YY]);
 403         }
 404         else if (i%3 == 2)
 405         {
 406             sscanf(buf, "%lf\n", &QMgrad[k][ZZ]);
 407         }
 408 #else
 409         if (i%3 == 0)
 410         {
 411             sscanf(buf, "%f\n", &QMgrad[k][XX]);
 412         }
 413         else if (i%3 == 1)
 414         {
 415             sscanf(buf, "%f\n", &QMgrad[k][YY]);
 416         }
 417         else if (i%3 == 2)
 418         {
 419             sscanf(buf, "%f\n", &QMgrad[k][ZZ]);
 420         }
 421 #endif
 422     }
 423     fclose(engrad);
 424     /* write the MM point charge data
 425      */
 426     if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
 427     {
 428         sprintf(orca_pcgradFilename, "%s.pcgrad", qm->orca_basename);
 429         pcgrad = fopen(orca_pcgradFilename, "r");
 430
 431         /* we read the gradient for the mm-atoms from the pcgrad file
 432          */
 433         /* we can skip the first line
 434          */
 435         if (fgets(buf, 300, pcgrad) == NULL)
 436         {
 437             gmx_fatal(FARGS, "Unexpected end of ORCA output");
 438         }
 439         for (i = 0; i < mm->nrMMatoms; i++)
 440         {
 441             if (fgets(buf, 300, pcgrad) == NULL)
 442             {
 443                 gmx_fatal(FARGS, "Unexpected end of ORCA output");
 444             }
 445     #ifdef GMX_DOUBLE
 446             sscanf(buf, "%lf%lf%lf\n",
 447                    &MMgrad[i][XX],
 448                    &MMgrad[i][YY],
 449                    &MMgrad[i][ZZ]);
 450     #else
 451             sscanf(buf, "%f%f%f\n",
 452                    &MMgrad[i][XX],
 453                    &MMgrad[i][YY],
 454                    &MMgrad[i][ZZ]);
 455     #endif
 456         }
 457         fclose(pcgrad);
 458     }
 459     return(QMener);
 460 }
 461
 462 void do_orca(char *orca_dir, char *basename)
 463 {
 464
 465     /* make the call to the orca binary through system()
 466      * The location of the binary is set through the
 467      * environment.
 468      */
 469     char
 470         buf[100];
 471     sprintf(buf, "%s/%s %s.inp >> %s.out",
 472             orca_dir,
 473             "orca",
 474             basename,
 475             basename);
 476     fprintf(stderr, "Calling '%s'\n", buf);
 477     if (system(buf) != 0)
 478     {
 479         gmx_fatal(FARGS, "Call to '%s' failed\n", buf);
 480     }
 481 }
 482
 483 real call_orca(t_forcerec *fr,
 484                t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 485 {
 486     /* normal orca jobs */
 487     static int
 488         step = 0;
 489     int
 490         i, j;
 491     real
 492         QMener;
 493     rvec
 494        *QMgrad, *MMgrad;
 495     char
 496        *exe;
 497
 498     snew(exe, 30);
 499     sprintf(exe, "%s", "orca");
 500     snew(QMgrad, qm->nrQMatoms);
 501     snew(MMgrad, mm->nrMMatoms);
 502
 503     write_orca_input(fr, qm, mm);
 504     do_orca(qm->orca_dir, qm->orca_basename);
 505     QMener = read_orca_output(QMgrad, MMgrad, fr, qm, mm);
 506     /* put the QMMM forces in the force array and to the fshift
 507      */
 508     for (i = 0; i < qm->nrQMatoms; i++)
 509     {
 510         for (j = 0; j < DIM; j++)
 511         {
 512             f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];
 513             fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
 514         }
 515     }
 516     for (i = 0; i < mm->nrMMatoms; i++)
 517     {
 518         for (j = 0; j < DIM; j++)
 519         {
 520             f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];
 521             fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
 522         }
 523     }
 524     QMener = QMener*HARTREE2KJ*AVOGADRO;
 525     step++;
 526     free(exe);
 527     return(QMener);
 528 } /* call_orca */
 529
 530 /* end of orca sub routines */