src/mdlib/qm_orca.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 4.5
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2008, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Groningen Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "sysstuff.h"
  41 #include "typedefs.h"
  42 #include "macros.h"
  43 #include "smalloc.h"
  44 #include "assert.h"
  45 #include "physics.h"
  46 #include "macros.h"
  47 #include "vec.h"
  48 #include "force.h"
  49 #include "invblock.h"
  50 #include "confio.h"
  51 #include "names.h"
  52 #include "network.h"
  53 #include "pbc.h"
  54 #include "ns.h"
  55 #include "nrnb.h"
  56 #include "bondf.h"
  57 #include "mshift.h"
  58 #include "txtdump.h"
  59 #include "copyrite.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_commrec *cr, t_QMrec *qm, t_MMrec *mm)
  70 {
  71  char
  72    *buf;
  73  snew(buf,200);
  74  /* ORCA settings on the system */
  75  buf = getenv("BASENAME");
  76  if (buf){
  77      snew(qm->orca_basename,200);
  78      sscanf(buf,"%s",qm->orca_basename);
  79  }
  80  else
  81      gmx_fatal(FARGS,"no $BASENAME\n");
  82
  83  /* ORCA directory on the system */
  84  snew(buf,200);
  85  buf = getenv("ORCA_PATH");
  86  fprintf(stderr,"%s",buf);
  87
  88  if (buf){
  89    snew(qm->orca_dir,200);
  90    sscanf(buf,"%s",qm->orca_dir);
  91  }
  92  else
  93    gmx_fatal(FARGS,"no $ORCA_PATH, check manual\n");
  94
  95  fprintf(stderr,"%s...\n",qm->orca_dir);
  96  fprintf(stderr,"orca initialised...\n");
  97  /* since we append the output to the BASENAME.out file,
  98  we should delete an existent old out-file here. */
  99  sprintf(buf,"%s.out",qm->orca_basename);
 100  remove(buf);
 101 }
 102
 103
 104 void write_orca_input(int step ,t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
 105 {
 106  int
 107    i;
 108  t_QMMMrec
 109    *QMMMrec;
 110  FILE
 111    *out, *pcFile, *addInputFile, *LJCoeff;
 112  char
 113    *buf,*orcaInput,*addInputFilename,*LJCoeffFilename,
 114    *pcFilename,*exclInName,*exclOutName;
 115  QMMMrec = fr->qr;
 116  /* write the first part of the input-file */
 117  snew(orcaInput,200);
 118  sprintf(orcaInput,"%s.inp",qm->orca_basename);
 119  out = fopen(orcaInput,"w");
 120  snew(addInputFilename,200);
 121  sprintf(addInputFilename,"%s.ORCAINFO",qm->orca_basename);
 122  addInputFile = fopen(addInputFilename,"r");
 123  fprintf(out, "#input-file generated by gromacs\n");
 124  if(qm->bTS){
 125      fprintf(out,"!QMMMOpt TightSCF\n");
 126      fprintf(out,"%s\n","%geom TS_Search EF end");
 127  }
 128  else if (qm->bOPT){
 129      fprintf(out,"!QMMMOpt TightSCF\n");
 130  }
 131  else{
 132      fprintf(out,"!EnGrad TightSCF\n");
 133  }
 134  /* here we include the insertion of the additional orca-input */
 135  snew(buf,200);
 136  if (addInputFile!=NULL) {
 137      while (!feof(addInputFile)) {
 138          if (fgets(buf, 200, addInputFile) != NULL)
 139              fputs(buf, out);
 140      }
 141  }
 142  else {
 143      fprintf(stderr,"No information on the calculation given in <%s>\n",addInputFilename);
 144      gmx_call("qm_orca.c");
 145  }
 146  fclose(addInputFile);
 147  if(qm->bTS||qm->bOPT){
 148      /* freeze the frontier QM atoms and Link atoms. This is
 149       * important only if a full QM subsystem optimization is done
 150       * with a frozen MM environmeent. For dynamics, or gromacs's own
 151       * optimization routines this is not important.
 152       */
 153      /* ORCA reads the exclusions from LJCoeffFilename.Excl,
 154       *so we have to rename the file
 155       */
 156      int didStart = 0;
 157      for(i=0;i<qm->nrQMatoms;i++){
 158           if(qm->frontatoms[i]){
 159              if (!didStart){
 160                  fprintf(out,"%s\n","%geom");
 161                  fprintf(out,"   Constraints \n");
 162                  didStart = 1;
 163              }
 164               fprintf(out,"        {C %d C}\n",i); /* counting from 0 */
 165           }
 166      }
 167      if (didStart) fprintf(out,"     end\n   end\n");
 168      /* make a file with information on the C6 and C12 coefficients */
 169      if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
 170          snew(exclInName,200);
 171          snew(exclOutName,200);
 172          sprintf(exclInName,"QMMMexcl.dat");
 173          sprintf(exclOutName,"%s.LJ.Excl",qm->orca_basename);
 174          rename(exclInName,exclOutName);
 175          snew(LJCoeffFilename,200);
 176          sprintf(LJCoeffFilename,"%s.LJ",qm->orca_basename);
 177          fprintf(out,"%s%s%s\n","%LJCOEFFICIENTS \"",LJCoeffFilename,"\"");
 178          /* make a file with information on the C6 and C12 coefficients */
 179          LJCoeff = fopen(LJCoeffFilename,"w");
 180          fprintf(LJCoeff,"%d\n",qm->nrQMatoms);
 181          for (i=0;i<qm->nrQMatoms;i++){
 182 #ifdef GMX_DOUBLE
 183              fprintf(LJCoeff,"%10.7lf  %10.7lf\n",qm->c6[i],qm->c12[i]);
 184 #else
 185              fprintf(LJCoeff,"%10.7f  %10.7f\n",qm->c6[i],qm->c12[i]);
 186 #endif
 187          }
 188          fprintf(LJCoeff,"%d\n",mm->nrMMatoms);
 189          for (i=0;i<mm->nrMMatoms;i++){
 190 #ifdef GMX_DOUBLE
 191              fprintf(LJCoeff,"%10.7lf  %10.7lf\n",mm->c6[i],mm->c12[i]);
 192 #else
 193              fprintf(LJCoeff,"%10.7f  %10.7f\n",mm->c6[i],mm->c12[i]);
 194 #endif
 195          }
 196          fclose(LJCoeff);
 197      }
 198  }
 199  /* write charge and multiplicity
 200   */
 201  fprintf(out,"*xyz %2d%2d\n",qm->QMcharge,qm->multiplicity);
 202  /* write the QM coordinates
 203   */
 204  for (i=0;i<qm->nrQMatoms;i++){
 205      int atomNr;
 206      if (qm->atomicnumberQM[i]==0)
 207          atomNr = 1;
 208      else
 209          atomNr = qm->atomicnumberQM[i];
 210 #ifdef GMX_DOUBLE
 211      fprintf(out,"%3d %10.7lf  %10.7lf  %10.7lf\n",
 212      atomNr,
 213      qm->xQM[i][XX]/0.1,
 214      qm->xQM[i][YY]/0.1,
 215      qm->xQM[i][ZZ]/0.1);
 216 #else
 217      fprintf(out,"%3d %10.7f  %10.7f  %10.7f\n",
 218      atomNr,
 219      qm->xQM[i][XX]/0.1,
 220      qm->xQM[i][YY]/0.1,
 221      qm->xQM[i][ZZ]/0.1);
 222 #endif
 223  }
 224  fprintf(out,"*\n");
 225  /* write the MM point charge data
 226   */
 227  if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
 228      /* name of the point charge file */
 229      snew(pcFilename,200);
 230      sprintf(pcFilename,"%s.pc",qm->orca_basename);
 231      fprintf(out,"%s%s%s\n","%pointcharges \"",pcFilename,"\"");
 232      pcFile = fopen(pcFilename,"w");
 233      fprintf(pcFile,"%d\n",mm->nrMMatoms);
 234      for(i=0;i<mm->nrMMatoms;i++){
 235 #ifdef GMX_DOUBLE
 236          fprintf(pcFile,"%8.4lf %10.7lf  %10.7lf  %10.7lf\n",
 237                         mm->MMcharges[i],
 238                         mm->xMM[i][XX]/0.1,
 239                         mm->xMM[i][YY]/0.1,
 240                         mm->xMM[i][ZZ]/0.1);
 241 #else
 242          fprintf(pcFile,"%8.4f %10.7f  %10.7f  %10.7f\n",
 243                         mm->MMcharges[i],
 244                         mm->xMM[i][XX]/0.1,
 245                         mm->xMM[i][YY]/0.1,
 246                         mm->xMM[i][ZZ]/0.1);
 247 #endif
 248      }
 249      fprintf(pcFile,"\n");
 250      fclose(pcFile);
 251  }
 252  fprintf(out,"\n");
 253
 254  fclose(out);
 255 }  /* write_orca_input */
 256
 257 real read_orca_output(rvec QMgrad[],rvec MMgrad[],int step,t_forcerec *fr,
 258                           t_QMrec *qm, t_MMrec *mm)
 259 {
 260  int
 261    i,j,atnum;
 262  char
 263    buf[300], tmp[300], orca_xyzFilename[300], orca_pcgradFilename[300], orca_engradFilename[300];
 264  real
 265    QMener;
 266  FILE
 267    *xyz, *pcgrad, *engrad;
 268  int k;
 269  t_QMMMrec
 270    *QMMMrec;
 271  QMMMrec = fr->qr;
 272  /* in case of an optimization, the coordinates are printed in the
 273   * xyz file, the energy and gradients for the QM part are stored in the engrad file
 274   * and the gradients for the point charges are stored in the pc file.
 275   */
 276
 277  /* we need the new xyz coordinates of the QM atoms only for separate QM-optimization
 278   */
 279
 280  if(qm->bTS||qm->bOPT){
 281      sprintf(orca_xyzFilename,"%s.xyz",qm->orca_basename);
 282      xyz=fopen(orca_xyzFilename,"r");
 283      if (fgets(buf,300,xyz) == NULL)
 284          gmx_fatal(FARGS, "Unexpected end of ORCA output");
 285      if (fgets(buf,300,xyz) == NULL)
 286          gmx_fatal(FARGS, "Unexpected end of ORCA output");
 287      for(i=0;i<qm->nrQMatoms;i++){
 288          if (fgets(buf,300,xyz) == NULL)
 289              gmx_fatal(FARGS, "Unexpected end of ORCA output");
 290 #ifdef GMX_DOUBLE
 291          sscanf(buf,"%s%lf%lf%lf\n",
 292                     tmp,
 293                     &qm->xQM[i][XX],
 294                     &qm->xQM[i][YY],
 295                     &qm->xQM[i][ZZ]);
 296 #else
 297          sscanf(buf,"%d%f%f%f\n",
 298                     &atnum,
 299                     &qm->xQM[i][XX],
 300                     &qm->xQM[i][YY],
 301                     &qm->xQM[i][ZZ]);
 302 #endif
 303          for(j=0;j<DIM;j++){
 304              qm->xQM[i][j]*=0.1;
 305          }
 306      }
 307      fclose(xyz);
 308  }
 309  sprintf(orca_engradFilename,"%s.engrad",qm->orca_basename);
 310  engrad=fopen(orca_engradFilename,"r");
 311  /* we read the energy and the gradient for the qm-atoms from the engrad file
 312   */
 313  /* we can skip the first seven lines
 314   */
 315  for (j=0;j<7;j++){
 316      if (fgets(buf,300,engrad) == NULL)
 317          gmx_fatal(FARGS, "Unexpected end of ORCA output");
 318  }
 319  /* now comes the energy
 320   */
 321  if (fgets(buf,300,engrad) == NULL)
 322      gmx_fatal(FARGS, "Unexpected end of ORCA output");
 323 #ifdef GMX_DOUBLE
 324  sscanf(buf,"%lf\n",&QMener);
 325 #else
 326  sscanf(buf,"%f\n", &QMener);
 327 #endif
 328  /* we can skip the next three lines
 329   */
 330  for (j=0;j<3;j++){
 331      if (fgets(buf,300,engrad) == NULL)
 332          gmx_fatal(FARGS, "Unexpected end of ORCA output");
 333  }
 334  /* next lines contain the gradients of the QM atoms
 335   * now comes the gradient, one value per line:
 336   * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
 337   */
 338
 339  for(i=0;i<3*qm->nrQMatoms;i++){
 340      k = i/3;
 341      if (fgets(buf,300,engrad) == NULL)
 342          gmx_fatal(FARGS, "Unexpected end of ORCA output");
 343 #ifdef GMX_DOUBLE
 344      if (i%3==0)
 345          sscanf(buf,"%lf\n", &QMgrad[k][XX]);
 346      else if (i%3==1)
 347          sscanf(buf,"%lf\n", &QMgrad[k][YY]);
 348      else if (i%3==2)
 349          sscanf(buf,"%lf\n", &QMgrad[k][ZZ]);
 350 #else
 351      if (i%3==0)
 352          sscanf(buf,"%f\n", &QMgrad[k][XX]);
 353      else if (i%3==1)
 354          sscanf(buf,"%f\n", &QMgrad[k][YY]);
 355      else if (i%3==2)
 356          sscanf(buf,"%f\n", &QMgrad[k][ZZ]);
 357 #endif
 358  }
 359  fclose(engrad);
 360  /* write the MM point charge data
 361   */
 362  if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
 363      sprintf(orca_pcgradFilename,"%s.pcgrad",qm->orca_basename);
 364      pcgrad=fopen(orca_pcgradFilename,"r");
 365
 366      /* we read the gradient for the mm-atoms from the pcgrad file
 367       */
 368      /* we can skip the first line
 369       */
 370      if (fgets(buf,300,pcgrad) == NULL)
 371          gmx_fatal(FARGS, "Unexpected end of ORCA output");
 372      for(i=0;i<mm->nrMMatoms;i++){
 373          if (fgets(buf,300,pcgrad) == NULL)
 374              gmx_fatal(FARGS, "Unexpected end of ORCA output");
 375     #ifdef GMX_DOUBLE
 376          sscanf(buf,"%lf%lf%lf\n",
 377                     &MMgrad[i][XX],
 378                     &MMgrad[i][YY],
 379                     &MMgrad[i][ZZ]);
 380     #else
 381          sscanf(buf,"%f%f%f\n",
 382                     &MMgrad[i][XX],
 383                     &MMgrad[i][YY],
 384                     &MMgrad[i][ZZ]);
 385     #endif
 386      }
 387      fclose(pcgrad);
 388  }
 389  return(QMener);
 390 }
 391
 392 void do_orca(int step,char *exe, char *orca_dir, char *basename)
 393 {
 394
 395  /* make the call to the orca binary through system()
 396   * The location of the binary is set through the
 397   * environment.
 398   */
 399  char
 400    buf[100];
 401  sprintf(buf,"%s/%s %s.inp >> %s.out",
 402              orca_dir,
 403              "orca",
 404              basename,
 405              basename);
 406  fprintf(stderr,"Calling '%s'\n",buf);
 407  if ( system(buf) != 0 )
 408      gmx_fatal(FARGS,"Call to '%s' failed\n",buf);
 409 }
 410
 411 real call_orca(t_commrec *cr,  t_forcerec *fr,
 412                    t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 413 {
 414  /* normal orca jobs */
 415  static int
 416    step=0;
 417  int
 418    i,j;
 419  real
 420    QMener=0.0;
 421  rvec
 422    *QMgrad,*MMgrad;
 423  char
 424    *exe;
 425
 426  snew(exe,30);
 427  sprintf(exe,"%s","orca");
 428  snew(QMgrad,qm->nrQMatoms);
 429  snew(MMgrad,mm->nrMMatoms);
 430
 431  write_orca_input(step,fr,qm,mm);
 432  do_orca(step,exe,qm->orca_dir,qm->orca_basename);
 433  QMener = read_orca_output(QMgrad,MMgrad,step,fr,qm,mm);
 434  /* put the QMMM forces in the force array and to the fshift
 435   */
 436  for(i=0;i<qm->nrQMatoms;i++){
 437      for(j=0;j<DIM;j++){
 438          f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];
 439          fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
 440      }
 441  }
 442  for(i=0;i<mm->nrMMatoms;i++){
 443      for(j=0;j<DIM;j++){
 444          f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];
 445          fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
 446      }
 447  }
 448  QMener = QMener*HARTREE2KJ*AVOGADRO;
 449  step++;
 450  free(exe);
 451  return(QMener);
 452 } /* call_orca */
 453
 454 /* end of orca sub routines */