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 #define F77_FUNC(name,NAME) name ## _
79 F77_FUNC(inigms,IMIGMS)(void);
82 F77_FUNC(endgms,ENDGMS)(void);
85 F77_FUNC(grads,GRADS)(int *nrqmat,real *qmcrd,int *nrmmat, real *mmchrg,
86 real *mmcrd, real *qmgrad,real *mmgrad, real *energy);
90 void init_gamess(t_commrec *cr,t_QMrec *qm, t_MMrec *mm){
91 /* it works hopelessly complicated :-)
92 * first a file is written. Then the standard gamess input/output
93 * routine is called (no system()!) to set up all fortran arrays.
94 * this routine writes a punch file, like in a normal gamess run.
95 * via this punch file the other games routines, needed for gradient
96 * and energy evaluations are called. This setup works fine for
97 * dynamics simulations. 7-6-2002 (London)
104 periodic_system[37][3]={"XX","H ","He","Li","Be","B ","C ","N ",
105 "O ","F ","Ne","Na","Mg","Al","Si","P ",
106 "S ","Cl","Ar","K ","Ca","Sc","Ti","V ",
107 "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga",
108 "Ge","As","Se","Br","Kr"};
113 out=fopen("FOR009","w");
114 /* of these options I am not completely sure.... the overall
115 * preformance on more than 4 cpu's is rather poor at the moment.
117 fprintf(out,"memory 48000000\nPARALLEL IOMODE SCREENED\n");
118 fprintf(out,"ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
119 qm->nelectrons,qm->multiplicity);
120 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]]);
140 fprintf(out,"%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
146 fprintf(out,"%10.7f %10.7f %10.7f %5.3f BQ\n",
155 fprintf(out,"END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
156 eQMbasis_names[qm->QMbasis],
157 eQMmethod_names[qm->QMmethod]); /* see enum.h */
159 fprintf(out,"END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
160 eQMbasis_names[qm->QMbasis],
161 eQMmethod_names[qm->QMmethod]); /* see enum.h */
165 F77_FUNC(inigms,IMIGMS)();
167 else{ /* normal serial run */
169 out=fopen("FOR009","w");
170 /* of these options I am not completely sure.... the overall
171 * preformance on more than 4 cpu's is rather poor at the moment.
173 fprintf(out,"ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
174 qm->nelectrons,qm->multiplicity);
175 for (i=0;i<qm->nrQMatoms;i++){
177 fprintf(out,"%10.7lf %10.7lf %10.7lf %5.3lf %2s\n",
181 qm->atomicnumberQM[i]*1.0,
182 periodic_system[qm->atomicnumberQM[i]]);
184 fprintf(out,"%10.7f %10.7f %10.7f %5.3f %2s\n",
188 qm->atomicnumberQM[i]*1.0,
189 periodic_system[qm->atomicnumberQM[i]]);
195 fprintf(out,"%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
201 fprintf(out,"%10.7f %10.7f %10.7f %5.3f BQ\n",
210 fprintf(out,"END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
211 eQMbasis_names[qm->QMbasis],
212 eQMmethod_names[qm->QMmethod]); /* see enum.h */
214 fprintf(out,"END\nBASIS %s\nRUNTYPE SADDLE\nSCFTYPE %s\n",
215 eQMbasis_names[qm->QMbasis],
216 eQMmethod_names[qm->QMmethod]); /* see enum.h */
217 F77_FUNC(inigms,IMIGMS)();
221 real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
222 rvec f[], rvec fshift[])
224 /* do the actual QMMM calculation using GAMESS-UK. In this
225 * implementation (3-2001) a system call is made to the GAMESS-UK
226 * binary. Now we are working to get the electron integral, SCF, and
227 * gradient routines linked directly
232 QMener=0.0,*qmgrad,*mmgrad,*mmcrd,*qmcrd,energy;
236 /* copy the QMMMrec pointer */
238 snew(qmcrd, 3*(qm->nrQMatoms));
239 snew(mmcrd,3*(mm->nrMMatoms));
240 snew(qmgrad,3*(qm->nrQMatoms));
241 snew(mmgrad,3*(mm->nrMMatoms));
243 /* copy the data from qr into the arrays that are going to be used
244 * in the fortran routines of gamess
246 for(i=0;i<qm->nrQMatoms;i++){
248 qmcrd[DIM*i+j] = 1/BOHR2NM*qm->xQM[i][j];
251 for(i=0;i<mm->nrMMatoms;i++){
253 mmcrd[DIM*i+j] = 1/BOHR2NM*mm->xMM[i][j];
256 for (i=0;i<3*qm->nrQMatoms;i+=3){
257 fprintf(stderr,"%8.5f, %8.5f, %8.5f\n",
263 F77_FUNC(grads,GRADS)(&qm->nrQMatoms,qmcrd,&mm->nrMMatoms,mm->MMcharges,
264 mmcrd,qmgrad,mmgrad,&energy);
266 for(i=0;i<qm->nrQMatoms;i++){
268 f[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
269 fshift[i][j] = HARTREE_BOHR2MD*qmgrad[3*i+j];
272 for(i=0;i<mm->nrMMatoms;i++){
274 f[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
275 fshift[i][j] = HARTREE_BOHR2MD*mmgrad[3*i+j];
278 /* convert a.u to kJ/mol */
279 QMener=energy*HARTREE2KJ*AVOGADRO;
285 gmx_qmmm_gamess_empty;