3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
64 #include "gmx_fatal.h"
69 /* mopac interface routines */
71 F77_FUNC(domldt,DOMLDT)(int *nrqmat, int labels[], char keywords[]);
74 F77_FUNC(domop,DOMOP)(int *nrqmat,double qmcrd[],int *nrmmat,
75 double mmchrg[],double mmcrd[],double qmgrad[],
76 double mmgrad[], double *energy,double qmcharges[]);
80 void init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
82 /* initializes the mopac routines ans sets up the semiempirical
83 * computation by calling moldat(). The inline mopac routines can
84 * only perform gradient operations. If one would like to optimize a
85 * structure or find a transition state at PM3 level, gaussian is
93 if(!qm->bSH){ /* if rerun then grad should not be done! */
94 sprintf(keywords,"PRECISE GEO-OK CHARGE=%d GRAD MMOK ANALYT %s\n",
96 eQMmethod_names[qm->QMmethod]);
99 sprintf(keywords,"PRECISE GEO-OK CHARGE=%d SINGLET GRAD %s C.I.=(%d,%d) root=2 MECI \n",
101 eQMmethod_names[qm->QMmethod],
102 qm->CASorbitals,qm->CASelectrons/2);
103 F77_FUNC(domldt,DOMLDT)(&qm->nrQMatoms,qm->atomicnumberQM,keywords);
104 fprintf(stderr,"keywords are: %s\n",keywords);
109 real call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
110 rvec f[], rvec fshift[])
112 /* do the actual QMMM calculation using directly linked mopac subroutines
114 double /* always double as the MOPAC routines are always compiled in
115 double precission! */
116 *qmcrd=NULL,*qmchrg=NULL,*mmcrd=NULL,*mmchrg=NULL,
117 *qmgrad,*mmgrad=NULL,energy;
122 snew(qmcrd, 3*(qm->nrQMatoms));
123 snew(qmgrad,3*(qm->nrQMatoms));
124 /* copy the data from qr into the arrays that are going to be used
125 * in the fortran routines of MOPAC
127 for(i=0;i<qm->nrQMatoms;i++){
129 qmcrd[3*i+j] = (double)qm->xQM[i][j]*10;
133 /* later we will add the point charges here. There are some
134 * conceptual problems with semi-empirical QM in combination with
135 * point charges that we need to solve first....
137 gmx_fatal(FARGS,"At present only ONIOM is allowed in combination"
138 " with MOPAC QM subroutines\n");
141 /* now compute the energy and the gradients.
144 snew(qmchrg,qm->nrQMatoms);
145 F77_FUNC(domop,DOMOP)(&qm->nrQMatoms,qmcrd,&mm->nrMMatoms,
146 mmchrg,mmcrd,qmgrad,mmgrad,&energy,qmchrg);
147 /* add the gradients to the f[] array, and also to the fshift[].
148 * the mopac gradients are in kCal/angstrom.
150 for(i=0;i<qm->nrQMatoms;i++){
152 f[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
153 fshift[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
156 QMener = (real)CAL2JOULE*energy;
157 /* do we do something with the mulliken charges?? */
166 real call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
167 rvec f[], rvec fshift[])
169 /* do the actual SH QMMM calculation using directly linked mopac
172 double /* always double as the MOPAC routines are always compiled in
173 double precission! */
174 *qmcrd=NULL,*qmchrg=NULL,*mmcrd=NULL,*mmchrg=NULL,
175 *qmgrad,*mmgrad=NULL,energy;
181 snew(qmcrd, 3*(qm->nrQMatoms));
182 snew(qmgrad,3*(qm->nrQMatoms));
183 /* copy the data from qr into the arrays that are going to be used
184 * in the fortran routines of MOPAC
186 for(i=0;i<qm->nrQMatoms;i++){
188 qmcrd[3*i+j] = (double)qm->xQM[i][j]*10;
192 /* later we will add the point charges here. There are some
193 * conceptual problems with semi-empirical QM in combination with
194 * point charges that we need to solve first....
196 gmx_fatal(FARGS,"At present only ONIOM is allowed in combination with MOPAC\n");
199 /* now compute the energy and the gradients.
201 snew(qmchrg,qm->nrQMatoms);
203 F77_FUNC(domop,DOMOP)(&qm->nrQMatoms,qmcrd,&mm->nrMMatoms,
204 mmchrg,mmcrd,qmgrad,mmgrad,&energy,qmchrg);
205 /* add the gradients to the f[] array, and also to the fshift[].
206 * the mopac gradients are in kCal/angstrom.
208 for(i=0;i<qm->nrQMatoms;i++){
210 f[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
211 fshift[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
214 QMener = (real)CAL2JOULE*energy;
219 } /* call_mopac_SH */
223 gmx_qmmm_mopac_empty;