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
63 #include "gmx_fatal.h"
68 /* mopac interface routines */
70 F77_FUNC(domldt, DOMLDT) (int *nrqmat, int labels[], char keywords[]);
73 F77_FUNC(domop, DOMOP) (int *nrqmat, double qmcrd[], int *nrmmat,
74 double mmchrg[], double mmcrd[], double qmgrad[],
75 double mmgrad[], double *energy, double qmcharges[]);
79 void init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
81 /* initializes the mopac routines ans sets up the semiempirical
82 * computation by calling moldat(). The inline mopac routines can
83 * only perform gradient operations. If one would like to optimize a
84 * structure or find a transition state at PM3 level, gaussian is
92 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]);
100 sprintf(keywords, "PRECISE GEO-OK CHARGE=%d SINGLET GRAD %s C.I.=(%d,%d) root=2 MECI \n",
102 eQMmethod_names[qm->QMmethod],
103 qm->CASorbitals, qm->CASelectrons/2);
105 F77_FUNC(domldt, DOMLDT) (&qm->nrQMatoms, qm->atomicnumberQM, keywords);
106 fprintf(stderr, "keywords are: %s\n", keywords);
111 real call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
112 rvec f[], rvec fshift[])
114 /* do the actual QMMM calculation using directly linked mopac subroutines
116 double /* always double as the MOPAC routines are always compiled in
117 double precission! */
118 *qmcrd = NULL, *qmchrg = NULL, *mmcrd = NULL, *mmchrg = NULL,
119 *qmgrad, *mmgrad = NULL, energy;
124 snew(qmcrd, 3*(qm->nrQMatoms));
125 snew(qmgrad, 3*(qm->nrQMatoms));
126 /* copy the data from qr into the arrays that are going to be used
127 * in the fortran routines of MOPAC
129 for (i = 0; i < qm->nrQMatoms; i++)
131 for (j = 0; j < DIM; j++)
133 qmcrd[3*i+j] = (double)qm->xQM[i][j]*10;
138 /* later we will add the point charges here. There are some
139 * conceptual problems with semi-empirical QM in combination with
140 * point charges that we need to solve first....
142 gmx_fatal(FARGS, "At present only ONIOM is allowed in combination"
143 " with MOPAC QM subroutines\n");
147 /* now compute the energy and the gradients.
150 snew(qmchrg, qm->nrQMatoms);
151 F77_FUNC(domop, DOMOP) (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms,
152 mmchrg, mmcrd, qmgrad, mmgrad, &energy, qmchrg);
153 /* add the gradients to the f[] array, and also to the fshift[].
154 * the mopac gradients are in kCal/angstrom.
156 for (i = 0; i < qm->nrQMatoms; i++)
158 for (j = 0; j < DIM; j++)
160 f[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
161 fshift[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
164 QMener = (real)CAL2JOULE*energy;
165 /* do we do something with the mulliken charges?? */
174 real call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
175 rvec f[], rvec fshift[])
177 /* do the actual SH QMMM calculation using directly linked mopac
180 double /* always double as the MOPAC routines are always compiled in
181 double precission! */
182 *qmcrd = NULL, *qmchrg = NULL, *mmcrd = NULL, *mmchrg = NULL,
183 *qmgrad, *mmgrad = NULL, energy;
189 snew(qmcrd, 3*(qm->nrQMatoms));
190 snew(qmgrad, 3*(qm->nrQMatoms));
191 /* copy the data from qr into the arrays that are going to be used
192 * in the fortran routines of MOPAC
194 for (i = 0; i < qm->nrQMatoms; i++)
196 for (j = 0; j < DIM; j++)
198 qmcrd[3*i+j] = (double)qm->xQM[i][j]*10;
203 /* later we will add the point charges here. There are some
204 * conceptual problems with semi-empirical QM in combination with
205 * point charges that we need to solve first....
207 gmx_fatal(FARGS, "At present only ONIOM is allowed in combination with MOPAC\n");
211 /* now compute the energy and the gradients.
213 snew(qmchrg, qm->nrQMatoms);
215 F77_FUNC(domop, DOMOP) (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms,
216 mmchrg, mmcrd, qmgrad, mmgrad, &energy, qmchrg);
217 /* add the gradients to the f[] array, and also to the fshift[].
218 * the mopac gradients are in kCal/angstrom.
220 for (i = 0; i < qm->nrQMatoms; i++)
222 for (j = 0; j < DIM; j++)
224 f[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
225 fshift[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
228 QMener = (real)CAL2JOULE*energy;
233 } /* call_mopac_SH */
237 gmx_qmmm_mopac_empty;