qm_orca.c
1 |
/*
|
---|---|
2 |
* This file is part of the GROMACS molecular simulation package.
|
3 |
*
|
4 |
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
|
5 |
* Copyright (c) 2001-2008, 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.
|
10 |
*
|
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.
|
15 |
*
|
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.
|
20 |
*
|
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.
|
25 |
*
|
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.
|
33 |
*
|
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.
|
36 |
*/
|
37 |
#include "gmxpre.h" |
38 |
|
39 |
#include <math.h> |
40 |
#include <stdio.h> |
41 |
#include <stdlib.h> |
42 |
#include <string.h> |
43 |
|
44 |
#include "gromacs/fileio/confio.h" |
45 |
#include "gromacs/legacyheaders/force.h" |
46 |
#include "gromacs/legacyheaders/macros.h" |
47 |
#include "gromacs/legacyheaders/names.h" |
48 |
#include "gromacs/legacyheaders/network.h" |
49 |
#include "gromacs/legacyheaders/nrnb.h" |
50 |
#include "gromacs/legacyheaders/ns.h" |
51 |
#include "gromacs/legacyheaders/qmmm.h" |
52 |
#include "gromacs/legacyheaders/txtdump.h" |
53 |
#include "gromacs/legacyheaders/typedefs.h" |
54 |
#include "gromacs/math/units.h" |
55 |
#include "gromacs/math/vec.h" |
56 |
#include "gromacs/utility/fatalerror.h" |
57 |
#include "gromacs/utility/smalloc.h" |
58 |
|
59 |
/* ORCA interface routines */
|
60 |
|
61 |
void init_orca(t_QMrec *qm)
|
62 |
{ |
63 |
char *buf;
|
64 |
int len;
|
65 |
|
66 |
/* ORCA settings on the system */
|
67 |
buf = getenv("GMX_QM_ORCA_BASENAME");
|
68 |
if (buf) {
|
69 |
len = strlen(buf); |
70 |
snew(qm->orca_basename, len); |
71 |
memcpy(qm->orca_basename, buf, len); |
72 |
/* since we append the output to the BASENAME.out file,
|
73 |
we should delete an existent old out-file here. */
|
74 |
snew(buf, len+4);
|
75 |
sprintf(buf, "%s.out", qm->orca_basename);
|
76 |
remove(buf); |
77 |
sfree(buf); |
78 |
} else {
|
79 |
gmx_fatal(FARGS, "$GMX_QM_ORCA_BASENAME is not set\n");
|
80 |
} |
81 |
|
82 |
/* ORCA directory on the system */
|
83 |
buf = getenv("GMX_ORCA_PATH");
|
84 |
if (buf) {
|
85 |
len = strlen(buf); |
86 |
snew(qm->orca_dir, len); |
87 |
memcpy(qm->orca_dir, buf, len); |
88 |
} else {
|
89 |
gmx_fatal(FARGS, "$GMX_ORCA_PATH not set, check manual\n");
|
90 |
} |
91 |
|
92 |
fprintf(stderr, "Setting ORCA path to: %s...\n", qm->orca_dir);
|
93 |
fprintf(stderr, "ORCA initialised...\n\n");
|
94 |
} |
95 |
|
96 |
void write_orca_input(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
|
97 |
{ |
98 |
int i, orca_b_len;
|
99 |
t_QMMMrec *QMMMrec; |
100 |
FILE *out, *iofile; |
101 |
char *buf, *filename, *filename_ab, *exclInName = "QMMMexcl.dat"; |
102 |
|
103 |
QMMMrec = fr->qr; |
104 |
|
105 |
orca_b_len = strlen(qm->orca_basename); |
106 |
snew(filename, orca_b_len+10);
|
107 |
filename_ab = filename + orca_b_len; /* pointer in the filename just after the basename */
|
108 |
memcpy(filename, qm->orca_basename, orca_b_len); |
109 |
|
110 |
/* write the first part of the input-file */
|
111 |
sprintf(filename_ab, ".inp");
|
112 |
out = fopen(filename, "w");
|
113 |
if (out == NULL) { |
114 |
gmx_fatal(FARGS, "Cannot open ORCA input (%s) for writing\n", filename);
|
115 |
} |
116 |
|
117 |
sprintf(filename_ab, ".ORCAINFO");
|
118 |
iofile = fopen(filename, "r");
|
119 |
if (iofile == NULL) { |
120 |
gmx_fatal(FARGS, "No information on the calculation given in %s\n", filename);
|
121 |
} |
122 |
|
123 |
fprintf(out, "#input-file generated by GROMACS\n");
|
124 |
|
125 |
if (qm->bTS) {
|
126 |
fprintf(out, "!QMMMOpt TightSCF\n");
|
127 |
fprintf(out, "%s\n", "%geom TS_Search EF end"); |
128 |
} else if (qm->bOPT) { |
129 |
fprintf(out, "!QMMMOpt TightSCF\n");
|
130 |
} else {
|
131 |
fprintf(out, "!EnGrad TightSCF\n");
|
132 |
} |
133 |
|
134 |
/* here we include the insertion of the additional orca-input */
|
135 |
snew(buf, 200);
|
136 |
while (!feof(iofile)) {
|
137 |
if (fgets(buf, 200, iofile) != NULL) { |
138 |
fputs(buf, out); |
139 |
} |
140 |
} |
141 |
sfree(buf); |
142 |
|
143 |
fclose(iofile); |
144 |
|
145 |
if (qm->bTS || qm->bOPT) {
|
146 |
/* freeze the frontier QM atoms and Link atoms. This is
|
147 |
* important only if a full QM subsystem optimization is done
|
148 |
* with a frozen MM environmeent. For dynamics, or gromacs's own
|
149 |
* optimization routines this is not important.
|
150 |
*/
|
151 |
/* ORCA reads the exclusions from LJCoeffFilename.Excl,
|
152 |
* so we have to rename the file
|
153 |
*/
|
154 |
int didStart = 0; |
155 |
for (i = 0; i < qm->nrQMatoms; i++) { |
156 |
if (qm->frontatoms[i]) {
|
157 |
if (!didStart) {
|
158 |
fprintf(out, "%s\n", "%geom"); |
159 |
fprintf(out, " Constraints \n");
|
160 |
didStart = 1;
|
161 |
} |
162 |
fprintf(out, " {C %d C}\n", i); /* counting from 0 */ |
163 |
} |
164 |
} |
165 |
if (didStart) {
|
166 |
fprintf(out, " end\n end\n");
|
167 |
} |
168 |
/* make a file with information on the C6 and C12 coefficients */
|
169 |
if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) {
|
170 |
sprintf(filename_ab, ".LJ.Excl");
|
171 |
rename(exclInName, filename); |
172 |
sprintf(filename_ab, ".LJ");
|
173 |
fprintf(out, "%s%s%s\n", "%LJCOEFFICIENTS \"", filename, "\""); |
174 |
/* make a file with information on the C6 and C12 coefficients */
|
175 |
iofile = fopen(filename, "w");
|
176 |
fprintf(iofile, "%d\n", qm->nrQMatoms);
|
177 |
for (i = 0; i < qm->nrQMatoms; i++) { |
178 |
fprintf(iofile, |
179 |
#ifdef GMX_DOUBLE
|
180 |
"%10.7lf %10.7lf\n",
|
181 |
#else
|
182 |
"%10.7f %10.7f\n",
|
183 |
#endif
|
184 |
qm->c6[i], qm->c12[i]); |
185 |
} |
186 |
fprintf(iofile, "%d\n", mm->nrMMatoms);
|
187 |
for (i = 0; i < mm->nrMMatoms; i++) { |
188 |
fprintf(iofile, |
189 |
#ifdef GMX_DOUBLE
|
190 |
"%10.7lf %10.7lf\n",
|
191 |
#else
|
192 |
"%10.7f %10.7f\n",
|
193 |
#endif
|
194 |
mm->c6[i], mm->c12[i]); |
195 |
} |
196 |
fclose(iofile); |
197 |
} |
198 |
} |
199 |
|
200 |
/* write charge and multiplicity */
|
201 |
fprintf(out, "*xyz %2d%2d\n", qm->QMcharge, qm->multiplicity);
|
202 |
|
203 |
/* write the QM coordinates */
|
204 |
for (i = 0; i < qm->nrQMatoms; i++) { |
205 |
int atomNr;
|
206 |
if (qm->atomicnumberQM[i] == 0) { |
207 |
atomNr = 1;
|
208 |
} else {
|
209 |
atomNr = qm->atomicnumberQM[i]; |
210 |
} |
211 |
fprintf(out, |
212 |
#ifdef GMX_DOUBLE
|
213 |
"%3d %10.7lf %10.7lf %10.7lf\n",
|
214 |
#else
|
215 |
"%3d %10.7f %10.7f %10.7f\n",
|
216 |
#endif
|
217 |
atomNr, |
218 |
qm->xQM[i][XX]/0.1, |
219 |
qm->xQM[i][YY]/0.1, |
220 |
qm->xQM[i][ZZ]/0.1); |
221 |
} |
222 |
fprintf(out, "*\n");
|
223 |
|
224 |
/* write the MM point charge data */
|
225 |
if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) {
|
226 |
/* name of the point charge file */
|
227 |
sprintf(filename_ab, ".pc");
|
228 |
fprintf(out, "%s%s%s\n", "%pointcharges \"", filename, "\""); |
229 |
iofile = fopen(filename, "w");
|
230 |
fprintf(iofile, "%d\n", mm->nrMMatoms);
|
231 |
for (i = 0; i < mm->nrMMatoms; i++) { |
232 |
fprintf(iofile, |
233 |
#ifdef GMX_DOUBLE
|
234 |
"%8.4lf %10.7lf %10.7lf %10.7lf\n",
|
235 |
#else
|
236 |
"%8.4f %10.7f %10.7f %10.7f\n",
|
237 |
#endif
|
238 |
mm->MMcharges[i], |
239 |
mm->xMM[i][XX]/0.1, |
240 |
mm->xMM[i][YY]/0.1, |
241 |
mm->xMM[i][ZZ]/0.1); |
242 |
} |
243 |
fprintf(iofile, "\n");
|
244 |
fclose(iofile); |
245 |
} |
246 |
fprintf(out, "\n");
|
247 |
|
248 |
sfree(filename); |
249 |
fclose(out); |
250 |
} /* write_orca_input */
|
251 |
|
252 |
real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr, |
253 |
t_QMrec *qm, t_MMrec *mm) |
254 |
{ |
255 |
int i, j, orca_b_len;
|
256 |
char buf[300], tmp[300]; |
257 |
char *filename, *filename_ab, *err_msg1 = "Unexpected end of ORCA output"; |
258 |
real QMener; |
259 |
FILE *ifile; |
260 |
int k;
|
261 |
t_QMMMrec *QMMMrec; |
262 |
|
263 |
QMMMrec = fr->qr; |
264 |
|
265 |
orca_b_len = strlen(qm->orca_basename); |
266 |
snew(filename, orca_b_len+10);
|
267 |
filename_ab = filename + orca_b_len; /* pointer in the filename just after the basename */
|
268 |
memcpy(filename, qm->orca_basename, orca_b_len); |
269 |
|
270 |
/* in case of an optimization, the coordinates are printed in the
|
271 |
* xyz file, the energy and gradients for the QM part are stored in the engrad file
|
272 |
* and the gradients for the point charges are stored in the pc file.
|
273 |
*/
|
274 |
|
275 |
/* we need the new xyz coordinates of the QM atoms only for separate QM-optimization */
|
276 |
|
277 |
if (qm->bTS || qm->bOPT) {
|
278 |
sprintf(filename_ab, ".xyz");
|
279 |
ifile = fopen(filename, "r");
|
280 |
for(i=0;i<2;i++) { |
281 |
if (fgets(buf, 300, ifile) == NULL) { |
282 |
gmx_fatal(FARGS, err_msg1); |
283 |
} |
284 |
} |
285 |
for (i = 0; i < qm->nrQMatoms; i++) { |
286 |
if(fscanf(ifile,
|
287 |
#ifdef GMX_DOUBLE
|
288 |
"%s%lf%lf%lf",
|
289 |
#else
|
290 |
"%s%f%f%f",
|
291 |
#endif
|
292 |
tmp,&qm->xQM[i][XX],&qm->xQM[i][YY],&qm->xQM[i][ZZ])!=4) {
|
293 |
break;
|
294 |
} |
295 |
for (j = 0; j < DIM; j++) { |
296 |
qm->xQM[i][j] *= 0.1; |
297 |
} |
298 |
} |
299 |
fclose(ifile); |
300 |
if(i<qm->nrQMatoms) {
|
301 |
gmx_fatal(FARGS, err_msg1); |
302 |
} |
303 |
} |
304 |
sprintf(filename_ab, ".engrad");
|
305 |
ifile = fopen(filename, "r");
|
306 |
/* we read the energy and the gradient for the qm-atoms from the engrad file */
|
307 |
/* we can skip the first seven lines */
|
308 |
for (j = 0; j < 7; j++) { |
309 |
if (fgets(buf, 300, ifile) == NULL) { |
310 |
gmx_fatal(FARGS, err_msg1); |
311 |
} |
312 |
} |
313 |
/* now comes the energy */
|
314 |
if(fscanf(ifile,
|
315 |
#ifdef GMX_DOUBLE
|
316 |
"%lf",
|
317 |
#else
|
318 |
"%f",
|
319 |
#endif
|
320 |
&QMener)!=1) {
|
321 |
gmx_fatal(FARGS, err_msg1); |
322 |
} |
323 |
/* we can skip the next three lines (plus EOL from the previous one)*/
|
324 |
for (j = 0; j < 4; j++) { |
325 |
if (fgets(buf, 300, ifile) == NULL) { |
326 |
gmx_fatal(FARGS, err_msg1); |
327 |
} |
328 |
} |
329 |
/* next lines contain the gradients of the QM atoms
|
330 |
* now comes the gradient, one value per line:
|
331 |
* (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
|
332 |
*/
|
333 |
|
334 |
for (i = 0; i < qm->nrQMatoms; i++) { |
335 |
if(fscanf(ifile,
|
336 |
#ifdef GMX_DOUBLE
|
337 |
"%lf%lf%lf",
|
338 |
#else
|
339 |
"%f%f%f",
|
340 |
#endif
|
341 |
&QMgrad[i][XX],&QMgrad[i][YY],&QMgrad[i][ZZ])!=3) {
|
342 |
break;
|
343 |
} |
344 |
} |
345 |
fclose(ifile); |
346 |
if(i<qm->nrQMatoms) {
|
347 |
gmx_fatal(FARGS, err_msg1); |
348 |
} |
349 |
/* write the MM point charge data */
|
350 |
if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) {
|
351 |
sprintf(filename_ab, ".pcgrad");
|
352 |
ifile = fopen(filename, "r");
|
353 |
|
354 |
/* we read the gradient for the mm-atoms from the pcgrad file */
|
355 |
/* we can skip the first line */
|
356 |
if (fgets(buf, 300, ifile) == NULL) { |
357 |
gmx_fatal(FARGS, err_msg1); |
358 |
} |
359 |
for (i = 0; i < mm->nrMMatoms; i++) { |
360 |
if(fscanf(ifile,
|
361 |
#ifdef GMX_DOUBLE
|
362 |
"%lf%lf%lf",
|
363 |
#else
|
364 |
"%f%f%f",
|
365 |
#endif
|
366 |
&MMgrad[i][XX],&MMgrad[i][YY],&MMgrad[i][ZZ])!=3) {
|
367 |
break;
|
368 |
} |
369 |
} |
370 |
fclose(ifile); |
371 |
if(i<qm->nrQMatoms) {
|
372 |
gmx_fatal(FARGS, err_msg1); |
373 |
} |
374 |
} |
375 |
sfree(filename); |
376 |
return QMener;
|
377 |
} |
378 |
|
379 |
void do_orca(char *orca_dir, char *basename) |
380 |
{ |
381 |
|
382 |
/* make the call to the orca binary through system()
|
383 |
* The location of the binary is set through the
|
384 |
* environment.
|
385 |
*/
|
386 |
char *buf;
|
387 |
int len = strlen(orca_dir) + strlen(basename)*2 + 20; |
388 |
snew(buf,len); |
389 |
sprintf(buf,"%s/%s %s.inp >> %s.out",orca_dir,"orca",basename,basename); |
390 |
fprintf(stderr,"Calling '%s'\n",buf);
|
391 |
if (system(buf)!= 0) { /* ahhh... this pain again...*/ |
392 |
gmx_fatal(FARGS,"Call to '%s' failed\n",buf);
|
393 |
} |
394 |
sfree(buf); |
395 |
} |
396 |
|
397 |
real call_orca(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]) |
398 |
{ |
399 |
/* normal orca jobs */
|
400 |
static int step = 0; |
401 |
int i, j;
|
402 |
real QMener; |
403 |
rvec *QMgrad, *MMgrad; |
404 |
|
405 |
snew(QMgrad, qm->nrQMatoms); |
406 |
snew(MMgrad, mm->nrMMatoms); |
407 |
|
408 |
write_orca_input(fr, qm, mm); |
409 |
do_orca(qm->orca_dir, qm->orca_basename); |
410 |
QMener = read_orca_output(QMgrad, MMgrad, fr, qm, mm); |
411 |
/* put the QMMM forces in the force array and to the fshift
|
412 |
*/
|
413 |
for (i = 0; i < qm->nrQMatoms; i++) { |
414 |
for (j = 0; j < DIM; j++) { |
415 |
f[i][j] = HARTREE_BOHR2MD*QMgrad[i][j]; |
416 |
fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j]; |
417 |
} |
418 |
} |
419 |
for (i = 0; i < mm->nrMMatoms; i++) { |
420 |
for (j = 0; j < DIM; j++) { |
421 |
f[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j]; |
422 |
fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j]; |
423 |
} |
424 |
} |
425 |
QMener = QMener*HARTREE2KJ*AVOGADRO; |
426 |
step++; |
427 |
return(QMener);
|
428 |
} /* call_orca */
|
429 |
|
430 |
/* end of orca sub routines */
|