Project

General

Profile

qm_orca.c

Grzegorz Wieczorek, 04/01/2016 03:42 AM

 
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], *filename, *filename_ab, *err_msg1 = "Unexpected end of ORCA output";
257
    real QMener;
258
    FILE *ifile;
259
    int k;
260
    t_QMMMrec *QMMMrec;
261

    
262
    QMMMrec = fr->qr;
263

    
264
    orca_b_len = strlen(qm->orca_basename);
265
    snew(filename, orca_b_len+10);
266
    filename_ab = filename + orca_b_len; /* pointer in the filename just after the basename */
267
    memcpy(filename, qm->orca_basename, orca_b_len);
268

    
269
    /* in case of an optimization, the coordinates are printed in the
270
     * xyz file, the energy and gradients for the QM part are stored in the engrad file
271
     * and the gradients for the point charges are stored in the pc file.
272
     */
273

    
274
    /* we need the new xyz coordinates of the QM atoms only for separate QM-optimization */
275

    
276
    if (qm->bTS || qm->bOPT) {
277
        sprintf(filename_ab, ".xyz");
278
        ifile = fopen(filename, "r");
279
        for(i=0;i<2;i++) {
280
            if (fgets(buf, 300, ifile) == NULL) {
281
                gmx_fatal(FARGS, err_msg1);
282
            }
283
        }
284
        for (i = 0; i < qm->nrQMatoms; i++) {
285
            if(fscanf(ifile,
286
#ifdef GMX_DOUBLE
287
                      "%s%lf%lf%lf",
288
#else
289
                      "%s%f%f%f",
290
#endif
291
                       buf,&qm->xQM[i][XX],&qm->xQM[i][YY],&qm->xQM[i][ZZ])!=4) {
292
                break;
293
            }
294
            for (j = 0; j < DIM; j++) {
295
                qm->xQM[i][j] *= 0.1;
296
            }
297
        }
298
        fclose(ifile);
299
        if(i<qm->nrQMatoms) {
300
            gmx_fatal(FARGS, err_msg1);
301
        }
302
    }
303
    sprintf(filename_ab, ".engrad");
304
    ifile = fopen(filename, "r");
305
    /* we read the energy and the gradient for the qm-atoms from the engrad file */
306
    /* we can skip the first seven lines */
307
    for (j = 0; j < 7; j++) {
308
        if (fgets(buf, 300, ifile) == NULL) {
309
            gmx_fatal(FARGS, err_msg1);
310
        }
311
    }
312
    /* now comes the energy */
313
    if(fscanf(ifile,
314
#ifdef GMX_DOUBLE
315
              "%lf",
316
#else
317
              "%f",
318
#endif
319
              &QMener)!=1) {
320
        gmx_fatal(FARGS, err_msg1);
321
    }
322
    /* we can skip the next three lines (plus EOL from the previous one)*/
323
    for (j = 0; j < 4; j++) {
324
        if (fgets(buf, 300, ifile) == NULL) {
325
            gmx_fatal(FARGS, err_msg1);
326
        }
327
    }
328
    /* next lines contain the gradients of the QM atoms
329
     * now comes the gradient, one value per line:
330
     * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
331
     */
332

    
333
    for (i = 0; i < qm->nrQMatoms; i++) {
334
        if(fscanf(ifile,
335
#ifdef GMX_DOUBLE
336
                  "%lf%lf%lf",
337
#else
338
                  "%f%f%f",
339
#endif
340
                  &QMgrad[i][XX],&QMgrad[i][YY],&QMgrad[i][ZZ])!=3) {
341
            break;
342
        }
343
    }
344
    fclose(ifile);
345
    if(i<qm->nrQMatoms) {
346
        gmx_fatal(FARGS, err_msg1);
347
    }
348
    /* write the MM point charge data */
349
    if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) {
350
        sprintf(filename_ab, ".pcgrad");
351
        ifile = fopen(filename, "r");
352

    
353
        /* we read the gradient for the mm-atoms from the pcgrad file */
354
        /* we can skip the first line */
355
        if (fgets(buf, 300, ifile) == NULL) {
356
            gmx_fatal(FARGS, err_msg1);
357
        }
358
        for (i = 0; i < mm->nrMMatoms; i++) {
359
            if(fscanf(ifile,
360
#ifdef GMX_DOUBLE
361
                      "%lf%lf%lf",
362
#else
363
                      "%f%f%f",
364
#endif
365
                      &MMgrad[i][XX],&MMgrad[i][YY],&MMgrad[i][ZZ])!=3) {
366
                break;
367
            }
368
        }
369
        fclose(ifile);
370
        if(i<qm->nrQMatoms) {
371
            gmx_fatal(FARGS, err_msg1);
372
        }
373
    }
374
    sfree(filename);
375
    return QMener;
376
}
377

    
378
void do_orca(char *orca_dir, char *basename)
379
{
380

    
381
    /* make the call to the orca binary through system()
382
     * The location of the binary is set through the
383
     * environment.
384
     */
385
    char *buf;
386
    int len = strlen(orca_dir) + strlen(basename)*2 + 20;
387
    snew(buf,len);
388
    sprintf(buf,"%s/%s %s.inp >> %s.out",orca_dir,"orca",basename,basename);
389
    fprintf(stderr,"Calling '%s'\n",buf);
390
    if (system(buf)!= 0) { /* ahhh... this pain again...*/
391
        gmx_fatal(FARGS,"Call to '%s' failed\n",buf);
392
    }
393
    sfree(buf);
394
}
395

    
396
real call_orca(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
397
{
398
    /* normal orca jobs */
399
    static int step = 0;
400
    int i, j;
401
    real QMener;
402
    rvec *QMgrad, *MMgrad;
403

    
404
    snew(QMgrad, qm->nrQMatoms);
405
    snew(MMgrad, mm->nrMMatoms);
406

    
407
    write_orca_input(fr, qm, mm);
408
    do_orca(qm->orca_dir, qm->orca_basename);
409
    QMener = read_orca_output(QMgrad, MMgrad, fr, qm, mm);
410
    /* put the QMMM forces in the force array and to the fshift
411
     */
412
    for (i = 0; i < qm->nrQMatoms; i++) {
413
        for (j = 0; j < DIM; j++) {
414
            f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];
415
            fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
416
        }
417
    }
418
    for (i = 0; i < mm->nrMMatoms; i++) {
419
        for (j = 0; j < DIM; j++) {
420
            f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];
421
            fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
422
        }
423
    }
424
    QMener = QMener*HARTREE2KJ*AVOGADRO;
425
    step++;
426
    return(QMener);
427
} /* call_orca */
428

    
429
/* end of orca sub routines */