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