src/mdlib/qm_mopac.c

   1 /*
   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 3.2.0
  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-2004, 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  * GROwing Monsters And Cloning Shrimps
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #ifdef GMX_QMMM_MOPAC
  40
  41 #include <math.h>
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "macros.h"
  45 #include "smalloc.h"
  46 #include "assert.h"
  47 #include "physics.h"
  48 #include "macros.h"
  49 #include "vec.h"
  50 #include "force.h"
  51 #include "invblock.h"
  52 #include "confio.h"
  53 #include "names.h"
  54 #include "network.h"
  55 #include "pbc.h"
  56 #include "ns.h"
  57 #include "nrnb.h"
  58 #include "bondf.h"
  59 #include "mshift.h"
  60 #include "txtdump.h"
  61 #include "copyrite.h"
  62 #include "qmmm.h"
  63 #include <stdio.h>
  64 #include <string.h>
  65 #include "gmx_fatal.h"
  66 #include "typedefs.h"
  67 #include <stdlib.h>
  68
  69
  70 /* mopac interface routines */
  71
  72 #ifndef F77_FUNC
  73 #define F77_FUNC(name,NAME) name ## _
  74 #endif
  75
  76
  77 void
  78 F77_FUNC(domldt,DOMLDT)(int *nrqmat, int labels[], char keywords[]);
  79
  80 void
  81 F77_FUNC(domop,DOMOP)(int *nrqmat,double qmcrd[],int *nrmmat,
  82                       double mmchrg[],double mmcrd[],double qmgrad[],
  83                       double mmgrad[], double *energy,double qmcharges[]);
  84
  85
  86
  87 void init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
  88 {
  89   /* initializes the mopac routines ans sets up the semiempirical
  90    * computation by calling moldat(). The inline mopac routines can
  91    * only perform gradient operations. If one would like to optimize a
  92    * structure or find a transition state at PM3 level, gaussian is
  93    * used instead.
  94    */
  95   char
  96     *keywords;
  97
  98   snew(keywords,240);
  99
 100   if(!qm->bSH){    /* if rerun then grad should not be done! */
 101     sprintf(keywords,"PRECISE GEO-OK CHARGE=%d GRAD MMOK ANALYT %s\n",
 102             qm->QMcharge,
 103             eQMmethod_names[qm->QMmethod]);
 104   }
 105   else
 106     sprintf(keywords,"PRECISE GEO-OK CHARGE=%d SINGLET GRAD %s C.I.=(%d,%d) root=2 MECI \n",
 107             qm->QMcharge,
 108             eQMmethod_names[qm->QMmethod],
 109             qm->CASorbitals,qm->CASelectrons/2);
 110   F77_FUNC(domldt,DOMLDT)(&qm->nrQMatoms,qm->atomicnumberQM,keywords);
 111   fprintf(stderr,"keywords are: %s\n",keywords);
 112   free(keywords);
 113
 114 } /* init_mopac */
 115
 116 real call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
 117                 rvec f[], rvec fshift[])
 118 {
 119   /* do the actual QMMM calculation using directly linked mopac subroutines
 120    */
 121   double /* always double as the MOPAC routines are always compiled in
 122             double precission! */
 123     *qmcrd=NULL,*qmchrg=NULL,*mmcrd=NULL,*mmchrg=NULL,
 124     *qmgrad,*mmgrad=NULL,energy;
 125   int
 126     i,j;
 127   real
 128     QMener=0.0;
 129   snew(qmcrd, 3*(qm->nrQMatoms));
 130   snew(qmgrad,3*(qm->nrQMatoms));
 131   /* copy the data from qr into the arrays that are going to be used
 132    * in the fortran routines of MOPAC
 133    */
 134   for(i=0;i<qm->nrQMatoms;i++){
 135     for (j=0;j<DIM;j++){
 136       qmcrd[3*i+j] = (double)qm->xQM[i][j]*10;
 137     }
 138   }
 139   if(mm->nrMMatoms){
 140     /* later we will add the point charges here. There are some
 141      * conceptual problems with semi-empirical QM in combination with
 142      * point charges that we need to solve first....
 143      */
 144     gmx_fatal(FARGS,"At present only ONIOM is allowed in combination"
 145                 " with MOPAC QM subroutines\n");
 146   }
 147   else {
 148     /* now compute the energy and the gradients.
 149      */
 150
 151     snew(qmchrg,qm->nrQMatoms);
 152     F77_FUNC(domop,DOMOP)(&qm->nrQMatoms,qmcrd,&mm->nrMMatoms,
 153            mmchrg,mmcrd,qmgrad,mmgrad,&energy,qmchrg);
 154     /* add the gradients to the f[] array, and also to the fshift[].
 155      * the mopac gradients are in kCal/angstrom.
 156      */
 157     for(i=0;i<qm->nrQMatoms;i++){
 158       for(j=0;j<DIM;j++){
 159         f[i][j]       = (real)10*CAL2JOULE*qmgrad[3*i+j];
 160         fshift[i][j]  = (real)10*CAL2JOULE*qmgrad[3*i+j];
 161       }
 162     }
 163     QMener = (real)CAL2JOULE*energy;
 164     /* do we do something with the mulliken charges?? */
 165
 166     free(qmchrg);
 167 }
 168   free(qmgrad);
 169   free(qmcrd);
 170   return (QMener);
 171 }
 172
 173 real call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
 174                    rvec f[], rvec fshift[])
 175 {
 176   /* do the actual SH QMMM calculation using directly linked mopac
 177    subroutines */
 178
 179   double /* always double as the MOPAC routines are always compiled in
 180             double precission! */
 181     *qmcrd=NULL,*qmchrg=NULL,*mmcrd=NULL,*mmchrg=NULL,
 182     *qmgrad,*mmgrad=NULL,energy;
 183   int
 184     i,j;
 185   real
 186     QMener=0.0;
 187
 188   snew(qmcrd, 3*(qm->nrQMatoms));
 189   snew(qmgrad,3*(qm->nrQMatoms));
 190   /* copy the data from qr into the arrays that are going to be used
 191    * in the fortran routines of MOPAC
 192    */
 193   for(i=0;i<qm->nrQMatoms;i++){
 194     for (j=0;j<DIM;j++){
 195       qmcrd[3*i+j] = (double)qm->xQM[i][j]*10;
 196     }
 197   }
 198   if(mm->nrMMatoms){
 199     /* later we will add the point charges here. There are some
 200      * conceptual problems with semi-empirical QM in combination with
 201      * point charges that we need to solve first....
 202      */
 203     gmx_fatal(FARGS,"At present only ONIOM is allowed in combination with MOPAC\n");
 204   }
 205   else {
 206     /* now compute the energy and the gradients.
 207      */
 208     snew(qmchrg,qm->nrQMatoms);
 209
 210     F77_FUNC(domop,DOMOP)(&qm->nrQMatoms,qmcrd,&mm->nrMMatoms,
 211            mmchrg,mmcrd,qmgrad,mmgrad,&energy,qmchrg);
 212     /* add the gradients to the f[] array, and also to the fshift[].
 213      * the mopac gradients are in kCal/angstrom.
 214      */
 215     for(i=0;i<qm->nrQMatoms;i++){
 216       for(j=0;j<DIM;j++){
 217         f[i][j]      = (real)10*CAL2JOULE*qmgrad[3*i+j];
 218         fshift[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
 219       }
 220     }
 221     QMener = (real)CAL2JOULE*energy;
 222   }
 223   free(qmgrad);
 224   free(qmcrd);
 225   return (QMener);
 226 } /* call_mopac_SH */
 227
 228 #else
 229 int
 230 gmx_qmmm_mopac_empty;
 231 #endif