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,2015,2017,2018, 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.
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/gmxlib/network.h"
51 #include "gromacs/gmxlib/nrnb.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdlib/ns.h"
55 #include "gromacs/mdlib/qmmm.h"
56 #include "gromacs/mdtypes/commrec.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/smalloc.h"
62 /* QMMM sub routines */
63 /* mopac interface routines */
67 F77_FUNC(inigms
, IMIGMS
) (void);
70 F77_FUNC(endgms
, ENDGMS
) (void);
73 F77_FUNC(grads
, GRADS
) (int *nrqmat
, real
*qmcrd
, int *nrmmat
, real
*mmchrg
,
74 real
*mmcrd
, real
*qmgrad
, real
*mmgrad
, real
*energy
);
78 void init_gamess(t_commrec
*cr
, t_QMrec
*qm
, t_MMrec
*mm
)
80 /* it works hopelessly complicated :-)
81 * first a file is written. Then the standard gamess input/output
82 * routine is called (no system()!) to set up all fortran arrays.
83 * this routine writes a punch file, like in a normal gamess run.
84 * via this punch file the other games routines, needed for gradient
85 * and energy evaluations are called. This setup works fine for
86 * dynamics simulations. 7-6-2002 (London)
93 periodic_system
[37][3] = {
94 "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ",
95 "O ", "F ", "Ne", "Na", "Mg", "Al", "Si", "P ",
96 "S ", "Cl", "Ar", "K ", "Ca", "Sc", "Ti", "V ",
97 "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga",
98 "Ge", "As", "Se", "Br", "Kr"
106 out
= fopen("FOR009", "w");
107 /* of these options I am not completely sure.... the overall
108 * preformance on more than 4 cpu's is rather poor at the moment.
110 fprintf(out
, "memory 48000000\nPARALLEL IOMODE SCREENED\n");
111 fprintf(out
, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
112 qm
->nelectrons
, qm
->multiplicity
);
113 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
116 fprintf(out
, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n",
120 qm
->atomicnumberQM
[i
]*1.0,
121 periodic_system
[qm
->atomicnumberQM
[i
]]);
123 fprintf(out
, "%10.7f %10.7f %10.7f %5.3f %2s\n",
127 qm
->atomicnumberQM
[i
]*1.0,
128 periodic_system
[qm
->atomicnumberQM
[i
]]);
133 for (j
= i
; j
< i
+2; j
++)
136 fprintf(out
, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
142 fprintf(out
, "%10.7f %10.7f %10.7f %5.3f BQ\n",
150 fprintf(out
, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
151 eQMbasis_names
[qm
->QMbasis
],
152 eQMmethod_names
[qm
->QMmethod
]); /* see enum.h */
156 F77_FUNC(inigms
, IMIGMS
) ();
158 else /* normal serial run */
161 out
= fopen("FOR009", "w");
162 /* of these options I am not completely sure.... the overall
163 * preformance on more than 4 cpu's is rather poor at the moment.
165 fprintf(out
, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n",
166 qm
->nelectrons
, qm
->multiplicity
);
167 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
170 fprintf(out
, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n",
174 qm
->atomicnumberQM
[i
]*1.0,
175 periodic_system
[qm
->atomicnumberQM
[i
]]);
177 fprintf(out
, "%10.7f %10.7f %10.7f %5.3f %2s\n",
181 qm
->atomicnumberQM
[i
]*1.0,
182 periodic_system
[qm
->atomicnumberQM
[i
]]);
187 for (j
= i
; j
< i
+2; j
++)
190 fprintf(out
, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n",
196 fprintf(out
, "%10.7f %10.7f %10.7f %5.3f BQ\n",
204 fprintf(out
, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
205 eQMbasis_names
[qm
->QMbasis
],
206 eQMmethod_names
[qm
->QMmethod
]); /* see enum.h */
207 F77_FUNC(inigms
, IMIGMS
) ();
211 real
call_gamess(t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
,
212 rvec f
[], rvec fshift
[])
214 /* do the actual QMMM calculation using GAMESS-UK. In this
215 * implementation (3-2001) a system call is made to the GAMESS-UK
216 * binary. Now we are working to get the electron integral, SCF, and
217 * gradient routines linked directly
222 QMener
= 0.0, *qmgrad
, *mmgrad
, *mmcrd
, *qmcrd
, energy
;
226 /* copy the QMMMrec pointer */
228 snew(qmcrd
, 3*(qm
->nrQMatoms
));
229 snew(mmcrd
, 3*(mm
->nrMMatoms
));
230 snew(qmgrad
, 3*(qm
->nrQMatoms
));
231 snew(mmgrad
, 3*(mm
->nrMMatoms
));
233 /* copy the data from qr into the arrays that are going to be used
234 * in the fortran routines of gamess
236 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
238 for (j
= 0; j
< DIM
; j
++)
240 qmcrd
[DIM
*i
+j
] = 1/BOHR2NM
*qm
->xQM
[i
][j
];
243 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
245 for (j
= 0; j
< DIM
; j
++)
247 mmcrd
[DIM
*i
+j
] = 1/BOHR2NM
*mm
->xMM
[i
][j
];
250 for (i
= 0; i
< 3*qm
->nrQMatoms
; i
+= 3)
252 fprintf(stderr
, "%8.5f, %8.5f, %8.5f\n",
258 F77_FUNC(grads
, GRADS
) (&qm
->nrQMatoms
, qmcrd
, &mm
->nrMMatoms
, mm
->MMcharges
,
259 mmcrd
, qmgrad
, mmgrad
, &energy
);
261 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
263 for (j
= 0; j
< DIM
; j
++)
265 f
[i
][j
] = HARTREE_BOHR2MD
*qmgrad
[3*i
+j
];
266 fshift
[i
][j
] = HARTREE_BOHR2MD
*qmgrad
[3*i
+j
];
269 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
271 for (j
= 0; j
< DIM
; j
++)
273 f
[i
][j
] = HARTREE_BOHR2MD
*mmgrad
[3*i
+j
];
274 fshift
[i
][j
] = HARTREE_BOHR2MD
*mmgrad
[3*i
+j
];
277 /* convert a.u to kJ/mol */
278 QMener
= energy
*HARTREE2KJ
*AVOGADRO
;
284 gmx_qmmm_gamess_empty
;