2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
41 #ifdef GMX_QMMM_GAMESS
47 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/fileio/confio.h"
65 #include "gmx_fatal.h"
70 /* QMMM sub routines */
71 /* mopac interface routines */
75 F77_FUNC(inigms, IMIGMS) (void);
78 F77_FUNC(endgms, ENDGMS) (void);
81 F77_FUNC(grads, GRADS) (int *nrqmat, real *qmcrd, int *nrmmat, real *mmchrg,
82 real *mmcrd, real *qmgrad, real *mmgrad, real *energy);
86 void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
88 /* it works hopelessly complicated :-)
89 * first a file is written. Then the standard gamess input/output
90 * routine is called (no system()!) to set up all fortran arrays.
91 * this routine writes a punch file, like in a normal gamess run.
92 * via this punch file the other games routines, needed for gradient
93 * and energy evaluations are called. This setup works fine for
94 * dynamics simulations. 7-6-2002 (London)
101 periodic_system[37][3] = {
102 "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ",
103 "O ", "F ", "Ne", "Na", "Mg", "Al", "Si", "P ",
104 "S ", "Cl", "Ar", "K ", "Ca", "Sc", "Ti", "V ",
105 "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga",
106 "Ge", "As", "Se", "Br", "Kr"
114 out = fopen("FOR009", "w");
115 /* of these options I am not completely sure.... the overall
116 * preformance on more than 4 cpu's is rather poor at the moment.
118 fprintf(out, "memory 48000000\nPARALLEL IOMODE SCREENED\n");
119 fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
120 qm->nelectrons, qm->multiplicity);
121 for (i = 0; i < qm->nrQMatoms; i++)
124 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n",
128 qm->atomicnumberQM[i]*1.0,
129 periodic_system[qm->atomicnumberQM[i]]);
131 fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n",
135 qm->atomicnumberQM[i]*1.0,
136 periodic_system[qm->atomicnumberQM[i]]);
141 for (j = i; j < i+2; j++)
144 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
150 fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n",
160 fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
161 eQMbasis_names[qm->QMbasis],
162 eQMmethod_names[qm->QMmethod]); /* see enum.h */
166 fprintf(out, "END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
167 eQMbasis_names[qm->QMbasis],
168 eQMmethod_names[qm->QMmethod]); /* see enum.h */
173 F77_FUNC(inigms, IMIGMS) ();
175 else /* normal serial run */
178 out = fopen("FOR009", "w");
179 /* of these options I am not completely sure.... the overall
180 * preformance on more than 4 cpu's is rather poor at the moment.
182 fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
183 qm->nelectrons, qm->multiplicity);
184 for (i = 0; i < qm->nrQMatoms; i++)
187 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n",
191 qm->atomicnumberQM[i]*1.0,
192 periodic_system[qm->atomicnumberQM[i]]);
194 fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n",
198 qm->atomicnumberQM[i]*1.0,
199 periodic_system[qm->atomicnumberQM[i]]);
204 for (j = i; j < i+2; j++)
207 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
213 fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n",
223 fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
224 eQMbasis_names[qm->QMbasis],
225 eQMmethod_names[qm->QMmethod]); /* see enum.h */
229 fprintf(out, "END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
230 eQMbasis_names[qm->QMbasis],
231 eQMmethod_names[qm->QMmethod]); /* see enum.h */
233 F77_FUNC(inigms, IMIGMS) ();
237 real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
238 rvec f[], rvec fshift[])
240 /* do the actual QMMM calculation using GAMESS-UK. In this
241 * implementation (3-2001) a system call is made to the GAMESS-UK
242 * binary. Now we are working to get the electron integral, SCF, and
243 * gradient routines linked directly
248 QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy;
252 /* copy the QMMMrec pointer */
254 snew(qmcrd, 3*(qm->nrQMatoms));
255 snew(mmcrd, 3*(mm->nrMMatoms));
256 snew(qmgrad, 3*(qm->nrQMatoms));
257 snew(mmgrad, 3*(mm->nrMMatoms));
259 /* copy the data from qr into the arrays that are going to be used
260 * in the fortran routines of gamess
262 for (i = 0; i < qm->nrQMatoms; i++)
264 for (j = 0; j < DIM; j++)
266 qmcrd[DIM*i+j] = 1/BOHR2NM*qm->xQM[i][j];
269 for (i = 0; i < mm->nrMMatoms; i++)
271 for (j = 0; j < DIM; j++)
273 mmcrd[DIM*i+j] = 1/BOHR2NM*mm->xMM[i][j];
276 for (i = 0; i < 3*qm->nrQMatoms; i += 3)
278 fprintf(stderr, "%8.5f, %8.5f, %8.5f\n",
284 F77_FUNC(grads, GRADS) (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mm->MMcharges,
285 mmcrd, qmgrad, mmgrad, &energy);
287 for (i = 0; i < qm->nrQMatoms; i++)
289 for (j = 0; j < DIM; j++)
291 f[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
292 fshift[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
295 for (i = 0; i < mm->nrMMatoms; i++)
297 for (j = 0; j < DIM; j++)
299 f[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
300 fshift[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
303 /* convert a.u to kJ/mol */
304 QMener = energy*HARTREE2KJ*AVOGADRO;
310 gmx_qmmm_gamess_empty;