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
39 #ifdef GMX_QMMM_GAMESS
64 #include "gmx_fatal.h"
69 /* QMMM sub routines */
70 /* mopac interface routines */
74 F77_FUNC(inigms, IMIGMS) (void);
77 F77_FUNC(endgms, ENDGMS) (void);
80 F77_FUNC(grads, GRADS) (int *nrqmat, real *qmcrd, int *nrmmat, real *mmchrg,
81 real *mmcrd, real *qmgrad, real *mmgrad, real *energy);
85 void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
87 /* it works hopelessly complicated :-)
88 * first a file is written. Then the standard gamess input/output
89 * routine is called (no system()!) to set up all fortran arrays.
90 * this routine writes a punch file, like in a normal gamess run.
91 * via this punch file the other games routines, needed for gradient
92 * and energy evaluations are called. This setup works fine for
93 * dynamics simulations. 7-6-2002 (London)
100 periodic_system[37][3] = {
101 "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ",
102 "O ", "F ", "Ne", "Na", "Mg", "Al", "Si", "P ",
103 "S ", "Cl", "Ar", "K ", "Ca", "Sc", "Ti", "V ",
104 "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga",
105 "Ge", "As", "Se", "Br", "Kr"
112 out = fopen("FOR009", "w");
113 /* of these options I am not completely sure.... the overall
114 * preformance on more than 4 cpu's is rather poor at the moment.
116 fprintf(out, "memory 48000000\nPARALLEL IOMODE SCREENED\n");
117 fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
118 qm->nelectrons, qm->multiplicity);
119 for (i = 0; i < qm->nrQMatoms; i++)
122 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n",
126 qm->atomicnumberQM[i]*1.0,
127 periodic_system[qm->atomicnumberQM[i]]);
129 fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n",
133 qm->atomicnumberQM[i]*1.0,
134 periodic_system[qm->atomicnumberQM[i]]);
139 for (j = i; j < i+2; j++)
142 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
148 fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n",
158 fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
159 eQMbasis_names[qm->QMbasis],
160 eQMmethod_names[qm->QMmethod]); /* see enum.h */
164 fprintf(out, "END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
165 eQMbasis_names[qm->QMbasis],
166 eQMmethod_names[qm->QMmethod]); /* see enum.h */
171 F77_FUNC(inigms, IMIGMS) ();
173 else /* normal serial run */
176 out = fopen("FOR009", "w");
177 /* of these options I am not completely sure.... the overall
178 * preformance on more than 4 cpu's is rather poor at the moment.
180 fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
181 qm->nelectrons, qm->multiplicity);
182 for (i = 0; i < qm->nrQMatoms; i++)
185 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n",
189 qm->atomicnumberQM[i]*1.0,
190 periodic_system[qm->atomicnumberQM[i]]);
192 fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n",
196 qm->atomicnumberQM[i]*1.0,
197 periodic_system[qm->atomicnumberQM[i]]);
202 for (j = i; j < i+2; j++)
205 fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
211 fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n",
221 fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
222 eQMbasis_names[qm->QMbasis],
223 eQMmethod_names[qm->QMmethod]); /* see enum.h */
227 fprintf(out, "END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
228 eQMbasis_names[qm->QMbasis],
229 eQMmethod_names[qm->QMmethod]); /* see enum.h */
231 F77_FUNC(inigms, IMIGMS) ();
235 real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
236 rvec f[], rvec fshift[])
238 /* do the actual QMMM calculation using GAMESS-UK. In this
239 * implementation (3-2001) a system call is made to the GAMESS-UK
240 * binary. Now we are working to get the electron integral, SCF, and
241 * gradient routines linked directly
246 QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy;
250 /* copy the QMMMrec pointer */
252 snew(qmcrd, 3*(qm->nrQMatoms));
253 snew(mmcrd, 3*(mm->nrMMatoms));
254 snew(qmgrad, 3*(qm->nrQMatoms));
255 snew(mmgrad, 3*(mm->nrMMatoms));
257 /* copy the data from qr into the arrays that are going to be used
258 * in the fortran routines of gamess
260 for (i = 0; i < qm->nrQMatoms; i++)
262 for (j = 0; j < DIM; j++)
264 qmcrd[DIM*i+j] = 1/BOHR2NM*qm->xQM[i][j];
267 for (i = 0; i < mm->nrMMatoms; i++)
269 for (j = 0; j < DIM; j++)
271 mmcrd[DIM*i+j] = 1/BOHR2NM*mm->xMM[i][j];
274 for (i = 0; i < 3*qm->nrQMatoms; i += 3)
276 fprintf(stderr, "%8.5f, %8.5f, %8.5f\n",
282 F77_FUNC(grads, GRADS) (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mm->MMcharges,
283 mmcrd, qmgrad, mmgrad, &energy);
285 for (i = 0; i < qm->nrQMatoms; i++)
287 for (j = 0; j < DIM; j++)
289 f[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
290 fshift[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
293 for (i = 0; i < mm->nrMMatoms; i++)
295 for (j = 0; j < DIM; j++)
297 f[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
298 fshift[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
301 /* convert a.u to kJ/mol */
302 QMener = energy*HARTREE2KJ*AVOGADRO;
308 gmx_qmmm_gamess_empty;