Project

General

Profile

qm_orca.c

Grzegorz Wieczorek, 04/01/2016 08:55 PM

 
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
    {
70
        len = strlen(buf);
71
        snew(qm->orca_basename, len);
72
        memcpy(qm->orca_basename, buf, len);
73
    /* since we append the output to the BASENAME.out file,
74
       we should delete an existent old out-file here. */
75
        snew(buf, len+4);
76
        sprintf(buf, "%s.out", qm->orca_basename);
77
        remove(buf);
78
        sfree(buf);
79
    }
80
    else
81
    {
82
        gmx_fatal(FARGS, "$GMX_QM_ORCA_BASENAME is not set\n");
83
    }
84

    
85
    /* ORCA directory on the system */
86
    buf = getenv("GMX_ORCA_PATH");
87
    if (buf)
88
    {
89
        len = strlen(buf);
90
        snew(qm->orca_dir, len);
91
        memcpy(qm->orca_dir, buf, len);
92
    }
93
    else
94
    {
95
        gmx_fatal(FARGS, "$GMX_ORCA_PATH not set, check manual\n");
96
    }
97

    
98
    fprintf(stderr, "Setting ORCA path to: %s...\n", qm->orca_dir);
99
    fprintf(stderr, "ORCA initialised...\n\n");
100
}
101

    
102
void write_orca_input(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
103
{
104
    int        i, orca_b_len;
105
    t_QMMMrec *QMMMrec;
106
    FILE      *out, *iofile;
107
    char      *buf, *filename, *filename_ab, *exclInName = "QMMMexcl.dat";
108

    
109
    QMMMrec = fr->qr;
110

    
111
    orca_b_len = strlen(qm->orca_basename);
112
    snew(filename, orca_b_len+10);
113
    filename_ab = filename + orca_b_len; /* pointer in the filename just after the basename */
114
    memcpy(filename, qm->orca_basename, orca_b_len);
115

    
116
    /* write the first part of the input-file */
117
    sprintf(filename_ab, ".inp");
118
    out = fopen(filename, "w");
119
    if (out == NULL)
120
    {
121
        gmx_fatal(FARGS, "Cannot open ORCA input (%s) for writing\n", filename);
122
    }
123

    
124
    sprintf(filename_ab, ".ORCAINFO");
125
    iofile = fopen(filename, "r");
126
    if (iofile == NULL)
127
    {
128
        gmx_fatal(FARGS, "No information on the calculation given in %s\n", filename);
129
    }
130

    
131
    fprintf(out, "#input-file generated by GROMACS\n");
132

    
133
    if (qm->bTS)
134
    {
135
        fprintf(out, "!QMMMOpt TightSCF\n");
136
        fprintf(out, "%s\n", "%geom TS_Search EF end");
137
    }
138
    else if (qm->bOPT)
139
    {
140
        fprintf(out, "!QMMMOpt TightSCF\n");
141
    }
142
    else
143
    {
144
        fprintf(out, "!EnGrad TightSCF\n");
145
    }
146

    
147
    /* here we include the insertion of the additional orca-input */
148
    snew(buf, 200);
149
    while (!feof(iofile))
150
    {
151
        if (fgets(buf, 200, iofile) != NULL)
152
        {
153
            fputs(buf, out);
154
        }
155
    }
156
    sfree(buf);
157

    
158
    fclose(iofile);
159

    
160
    if (qm->bTS || qm->bOPT)
161
    {
162
        /* freeze the frontier QM atoms and Link atoms. This is
163
         * important only if a full QM subsystem optimization is done
164
         * with a frozen MM environmeent. For dynamics, or gromacs's own
165
         * optimization routines this is not important.
166
         */
167
        /* ORCA reads the exclusions from LJCoeffFilename.Excl,
168
         * so we have to rename the file
169
         */
170
        int didStart = 0;
171
        for (i = 0; i < qm->nrQMatoms; i++)
172
        {
173
            if (qm->frontatoms[i])
174
            {
175
                if (!didStart)
176
                {
177
                    fprintf(out, "%s\n", "%geom");
178
                    fprintf(out, "   Constraints \n");
179
                    didStart = 1;
180
                }
181
                fprintf(out, "        {C %d C}\n", i); /* counting from 0 */
182
            }
183
        }
184
        if (didStart)
185
        {
186
            fprintf(out, "     end\n   end\n");
187
        }
188
        /* make a file with information on the C6 and C12 coefficients */
189
        if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
190
        {
191
            sprintf(filename_ab, ".LJ.Excl");
192
            rename(exclInName, filename);
193
            sprintf(filename_ab, ".LJ");
194
            fprintf(out, "%s%s%s\n", "%LJCOEFFICIENTS \"", filename, "\"");
195
            /* make a file with information on the C6 and C12 coefficients */
196
            iofile = fopen(filename, "w");
197
            fprintf(iofile, "%d\n", qm->nrQMatoms);
198
            for (i = 0; i < qm->nrQMatoms; i++)
199
            {
200
                fprintf(iofile,
201
#ifdef GMX_DOUBLE
202
                        "%10.7lf  %10.7lf\n",
203
#else
204
                        "%10.7f  %10.7f\n",
205
#endif
206
                        qm->c6[i], qm->c12[i]);
207
            }
208
            fprintf(iofile, "%d\n", mm->nrMMatoms);
209
            for (i = 0; i < mm->nrMMatoms; i++)
210
            {
211
                fprintf(iofile,
212
#ifdef GMX_DOUBLE
213
                        "%10.7lf  %10.7lf\n",
214
#else
215
                        "%10.7f  %10.7f\n",
216
#endif
217
                        mm->c6[i], mm->c12[i]);
218
            }
219
            fclose(iofile);
220
        }
221
    }
222

    
223
    /* write charge and multiplicity */
224
    fprintf(out, "*xyz %2d%2d\n", qm->QMcharge, qm->multiplicity);
225

    
226
    /* write the QM coordinates */
227
    for (i = 0; i < qm->nrQMatoms; i++)
228
    {
229
        int atomNr;
230
        if (qm->atomicnumberQM[i] == 0)
231
        {
232
            atomNr = 1;
233
        }
234
        else
235
        {
236
            atomNr = qm->atomicnumberQM[i];
237
        }
238
        fprintf(out,
239
#ifdef GMX_DOUBLE
240
                 "%3d %10.7lf  %10.7lf  %10.7lf\n",
241
#else
242
                 "%3d %10.7f  %10.7f  %10.7f\n",
243
#endif
244
                atomNr,
245
                qm->xQM[i][XX]/0.1,
246
                qm->xQM[i][YY]/0.1,
247
                qm->xQM[i][ZZ]/0.1);
248
    }
249
    fprintf(out, "*\n");
250

    
251
    /* write the MM point charge data */
252
    if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
253
    {
254
        /* name of the point charge file */
255
        sprintf(filename_ab, ".pc");
256
        fprintf(out, "%s%s%s\n", "%pointcharges \"", filename, "\"");
257
        iofile = fopen(filename, "w");
258
        fprintf(iofile, "%d\n", mm->nrMMatoms);
259
        for (i = 0; i < mm->nrMMatoms; i++)
260
        {
261
            fprintf(iofile,
262
#ifdef GMX_DOUBLE
263
                    "%8.4lf %10.7lf  %10.7lf  %10.7lf\n",
264
#else
265
                    "%8.4f %10.7f  %10.7f  %10.7f\n",
266
#endif
267
                    mm->MMcharges[i],
268
                    mm->xMM[i][XX]/0.1,
269
                    mm->xMM[i][YY]/0.1,
270
                    mm->xMM[i][ZZ]/0.1);
271
        }
272
        fprintf(iofile, "\n");
273
        fclose(iofile);
274
    }
275
    fprintf(out, "\n");
276

    
277
    sfree(filename);
278
    fclose(out);
279
}  /* write_orca_input */
280

    
281
real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr,
282
                      t_QMrec *qm, t_MMrec *mm)
283
{
284
    int i, j, orca_b_len;
285
    char buf[300], *filename, *filename_ab, *err_msg1 = "Unexpected end of ORCA output";
286
    real QMener;
287
    FILE *ifile;
288
    int k;
289
    t_QMMMrec *QMMMrec;
290

    
291
    QMMMrec = fr->qr;
292

    
293
    orca_b_len = strlen(qm->orca_basename);
294
    snew(filename, orca_b_len+10);
295
    filename_ab = filename + orca_b_len; /* pointer in the filename just after the basename */
296
    memcpy(filename, qm->orca_basename, orca_b_len);
297

    
298
    /* in case of an optimization, the coordinates are printed in the
299
     * xyz file, the energy and gradients for the QM part are stored in the engrad file
300
     * and the gradients for the point charges are stored in the pc file.
301
     */
302

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

    
305
    if (qm->bTS || qm->bOPT)
306
    {
307
        sprintf(filename_ab, ".xyz");
308
        ifile = fopen(filename, "r");
309
        for(i=0;i<2;i++)
310
        {
311
            if (fgets(buf, 300, ifile) == NULL)
312
            {
313
                gmx_fatal(FARGS, err_msg1);
314
            }
315
        }
316
        for (i = 0; i < qm->nrQMatoms; i++)
317
        {
318
            if(fscanf(ifile,
319
#ifdef GMX_DOUBLE
320
                      "%s%lf%lf%lf",
321
#else
322
                      "%s%f%f%f",
323
#endif
324
                       buf,&qm->xQM[i][XX],&qm->xQM[i][YY],&qm->xQM[i][ZZ])!=4)
325
            {
326
                break;
327
            }
328
            for (j = 0; j < DIM; j++)
329
            {
330
                qm->xQM[i][j] *= 0.1;
331
            }
332
        }
333
        fclose(ifile);
334
        if(i<qm->nrQMatoms)
335
        {
336
            gmx_fatal(FARGS, err_msg1);
337
        }
338
    }
339
    sprintf(filename_ab, ".engrad");
340
    ifile = fopen(filename, "r");
341
    /* we read the energy and the gradient for the qm-atoms from the engrad file */
342
    /* we can skip the first seven lines */
343
    for (j = 0; j < 7; j++)
344
    {
345
        if (fgets(buf, 300, ifile) == NULL)
346
        {
347
            gmx_fatal(FARGS, err_msg1);
348
        }
349
    }
350
    /* now comes the energy */
351
    if(fscanf(ifile,
352
#ifdef GMX_DOUBLE
353
              "%lf",
354
#else
355
              "%f",
356
#endif
357
              &QMener)!=1)
358
    {
359
        gmx_fatal(FARGS, err_msg1);
360
    }
361
    /* we can skip the next three lines (plus EOL from the previous one)*/
362
    for (j = 0; j < 4; j++)
363
    {
364
        if (fgets(buf, 300, ifile) == NULL)
365
        {
366
            gmx_fatal(FARGS, err_msg1);
367
        }
368
    }
369
    /* next lines contain the gradients of the QM atoms
370
     * now comes the gradient, one value per line:
371
     * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
372
     */
373

    
374
    for (i = 0; i < qm->nrQMatoms; i++)
375
    {
376
        if(fscanf(ifile,
377
#ifdef GMX_DOUBLE
378
                  "%lf%lf%lf",
379
#else
380
                  "%f%f%f",
381
#endif
382
                  &QMgrad[i][XX],&QMgrad[i][YY],&QMgrad[i][ZZ])!=3)
383
        {
384
            break;
385
        }
386
    }
387
    fclose(ifile);
388
    if(i<qm->nrQMatoms)
389
    {
390
        gmx_fatal(FARGS, err_msg1);
391
    }
392
    /* write the MM point charge data */
393
    if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
394
    {
395
        sprintf(filename_ab, ".pcgrad");
396
        ifile = fopen(filename, "r");
397

    
398
        /* we read the gradient for the mm-atoms from the pcgrad file */
399
        /* we can skip the first line */
400
        if (fgets(buf, 300, ifile) == NULL)
401
        {
402
            gmx_fatal(FARGS, err_msg1);
403
        }
404
        for (i = 0; i < mm->nrMMatoms; i++)
405
        {
406
            if(fscanf(ifile,
407
#ifdef GMX_DOUBLE
408
                      "%lf%lf%lf",
409
#else
410
                      "%f%f%f",
411
#endif
412
                      &MMgrad[i][XX],&MMgrad[i][YY],&MMgrad[i][ZZ])!=3)
413
            {
414
                break;
415
            }
416
        }
417
        fclose(ifile);
418
        if(i<qm->nrQMatoms)
419
        {
420
            gmx_fatal(FARGS, err_msg1);
421
        }
422
    }
423
    sfree(filename);
424
    return QMener;
425
}
426

    
427
void do_orca(char *orca_dir, char *basename)
428
{
429

    
430
    /* make the call to the orca binary through system()
431
     * The location of the binary is set through the
432
     * environment.
433
     */
434
    char *buf;
435
    int len = strlen(orca_dir) + strlen(basename)*2 + 20;
436
    snew(buf,len);
437
    sprintf(buf,"%s/%s %s.inp >> %s.out",orca_dir,"orca",basename,basename);
438
    fprintf(stderr,"Calling '%s'\n",buf);
439
    if (system(buf)!= 0)
440
    {
441
        gmx_fatal(FARGS,"Call to '%s' failed\n",buf);
442
    }
443
    sfree(buf);
444
}
445

    
446
real call_orca(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
447
{
448
    /* normal orca jobs */
449
    static int step = 0;
450
    int i, j;
451
    real QMener;
452
    rvec *QMgrad, *MMgrad;
453

    
454
    snew(QMgrad, qm->nrQMatoms);
455
    snew(MMgrad, mm->nrMMatoms);
456

    
457
    write_orca_input(fr, qm, mm);
458
    do_orca(qm->orca_dir, qm->orca_basename);
459
    QMener = read_orca_output(QMgrad, MMgrad, fr, qm, mm);
460
    /* put the QMMM forces in the force array and to the fshift
461
     */
462
    for (i = 0; i < qm->nrQMatoms; i++)
463
    {
464
        for (j = 0; j < DIM; j++)
465
        {
466
            f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];
467
            fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
468
        }
469
    }
470
    for (i = 0; i < mm->nrMMatoms; i++)
471
    {
472
        for (j = 0; j < DIM; j++)
473
        {
474
            f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];
475
            fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
476
        }
477
    }
478
    QMener = QMener*HARTREE2KJ*AVOGADRO;
479
    step++;
480
    return(QMener);
481
} /* call_orca */
482

    
483
/* end of orca sub routines */