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