GROMACS

/*

 * This file is part of the GROMACS molecular simulation package.

 *

 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.

 * Copyright (c) 2001-2008, The GROMACS development team.

 * Copyright (c) 2013,2014, by the GROMACS development team, led by

 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,

 * and including many others, as listed in the AUTHORS file in the

 * top-level source directory and at http://www.gromacs.org.

 *

 * GROMACS is free software; you can redistribute it and/or

 * modify it under the terms of the GNU Lesser General Public License

 * as published by the Free Software Foundation; either version 2.1

 * of the License, or (at your option) any later version.

 *

 * GROMACS is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * Lesser General Public License for more details.

 *

 * You should have received a copy of the GNU Lesser General Public

 * License along with GROMACS; if not, see

 * http://www.gnu.org/licenses, or write to the Free Software Foundation,

 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.

 *

 * If you want to redistribute modifications to GROMACS, please

 * consider that scientific software is very special. Version

 * control is crucial - bugs must be traceable. We will be happy to

 * consider code for inclusion in the official distribution, but

 * derived work must not be called official GROMACS. Details are found

 * in the README & COPYING files - if they are missing, get the

 * official version at http://www.gromacs.org.

 *

 * To help us fund GROMACS development, we humbly ask that you cite

 * the research papers on the package. Check out http://www.gromacs.org.

 */

#include "gmxpre.h"

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "gromacs/fileio/confio.h"

#include "gromacs/legacyheaders/force.h"

#include "gromacs/legacyheaders/macros.h"

#include "gromacs/legacyheaders/names.h"

#include "gromacs/legacyheaders/network.h"

#include "gromacs/legacyheaders/nrnb.h"

#include "gromacs/legacyheaders/ns.h"

#include "gromacs/legacyheaders/qmmm.h"

#include "gromacs/legacyheaders/txtdump.h"

#include "gromacs/legacyheaders/typedefs.h"

#include "gromacs/math/units.h"

#include "gromacs/math/vec.h"

#include "gromacs/utility/fatalerror.h"

#include "gromacs/utility/smalloc.h"

/* ORCA interface routines */

void init_orca(t_QMrec *qm)

{

    char *buf;

    int len;

    /* ORCA settings on the system */

    buf = getenv("GMX_QM_ORCA_BASENAME");

    if (buf) {

        len = strlen(buf);

        snew(qm->orca_basename, len);

        memcpy(qm->orca_basename, buf, len);

    /* since we append the output to the BASENAME.out file,

       we should delete an existent old out-file here. */

        snew(buf, len+4);

        sprintf(buf, "%s.out", qm->orca_basename);

        remove(buf);

        sfree(buf);

    } else {

        gmx_fatal(FARGS, "$GMX_QM_ORCA_BASENAME is not set\n");

    }

    /* ORCA directory on the system */

    buf = getenv("GMX_ORCA_PATH");

    if (buf) {

        len = strlen(buf);

        snew(qm->orca_dir, len);

        memcpy(qm->orca_dir, buf, len);

    } else {

        gmx_fatal(FARGS, "$GMX_ORCA_PATH not set, check manual\n");

    }

    fprintf(stderr, "Setting ORCA path to: %s...\n", qm->orca_dir);

    fprintf(stderr, "ORCA initialised...\n\n");

}

void write_orca_input(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)

{

    int        i, orca_b_len;

    t_QMMMrec *QMMMrec;

    FILE      *out, *iofile;

    char      *buf, *filename, *filename_ab, *exclInName = "QMMMexcl.dat";

    QMMMrec = fr->qr;

    orca_b_len = strlen(qm->orca_basename);

    snew(filename, orca_b_len+10);

    filename_ab = filename + orca_b_len; /* pointer in the filename just after the basename */

    memcpy(filename, qm->orca_basename, orca_b_len);

    /* write the first part of the input-file */

    sprintf(filename_ab, ".inp");

    out = fopen(filename, "w");

    if (out == NULL) {

        gmx_fatal(FARGS, "Cannot open ORCA input (%s) for writing\n", filename);

    }

    sprintf(filename_ab, ".ORCAINFO");

    iofile = fopen(filename, "r");

    if (iofile == NULL) {

        gmx_fatal(FARGS, "No information on the calculation given in %s\n", filename);

    }

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

    if (qm->bTS) {

        fprintf(out, "!QMMMOpt TightSCF\n");

        fprintf(out, "%s\n", "%geom TS_Search EF end");

    } else if (qm->bOPT) {

        fprintf(out, "!QMMMOpt TightSCF\n");

    } else {

        fprintf(out, "!EnGrad TightSCF\n");

    }

    /* here we include the insertion of the additional orca-input */

    snew(buf, 200);

    while (!feof(iofile)) {

        if (fgets(buf, 200, iofile) != NULL) {

            fputs(buf, out);

        }

    }

    sfree(buf);

    fclose(iofile);

    if (qm->bTS || qm->bOPT) {

        /* freeze the frontier QM atoms and Link atoms. This is

         * important only if a full QM subsystem optimization is done

         * with a frozen MM environmeent. For dynamics, or gromacs's own

         * optimization routines this is not important.

         */

        /* ORCA reads the exclusions from LJCoeffFilename.Excl,

         * so we have to rename the file

         */

        int didStart = 0;

        for (i = 0; i < qm->nrQMatoms; i++) {

            if (qm->frontatoms[i]) {

                if (!didStart) {

                    fprintf(out, "%s\n", "%geom");

                    fprintf(out, "   Constraints \n");

                    didStart = 1;

                }

                fprintf(out, "        {C %d C}\n", i); /* counting from 0 */

            }

        }

        if (didStart) {

            fprintf(out, "     end\n   end\n");

        }

        /* make a file with information on the C6 and C12 coefficients */

        if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) {

            sprintf(filename_ab, ".LJ.Excl");

            rename(exclInName, filename);

            sprintf(filename_ab, ".LJ");

            fprintf(out, "%s%s%s\n", "%LJCOEFFICIENTS \"", filename, "\"");

            /* make a file with information on the C6 and C12 coefficients */

            iofile = fopen(filename, "w");

            fprintf(iofile, "%d\n", qm->nrQMatoms);

            for (i = 0; i < qm->nrQMatoms; i++) {

                fprintf(iofile,

#ifdef GMX_DOUBLE

                        "%10.7lf  %10.7lf\n",

#else

                        "%10.7f  %10.7f\n",

#endif

                        qm->c6[i], qm->c12[i]);

            }

            fprintf(iofile, "%d\n", mm->nrMMatoms);

            for (i = 0; i < mm->nrMMatoms; i++) {

                fprintf(iofile,

#ifdef GMX_DOUBLE

                        "%10.7lf  %10.7lf\n",

#else

                        "%10.7f  %10.7f\n",

#endif

                        mm->c6[i], mm->c12[i]);

            }

            fclose(iofile);

        }

    }

    /* write charge and multiplicity */

    fprintf(out, "*xyz %2d%2d\n", qm->QMcharge, qm->multiplicity);

    /* write the QM coordinates */

    for (i = 0; i < qm->nrQMatoms; i++) {

        int atomNr;

        if (qm->atomicnumberQM[i] == 0) {

            atomNr = 1;

        } else {

            atomNr = qm->atomicnumberQM[i];

        }

        fprintf(out,

#ifdef GMX_DOUBLE

                 "%3d %10.7lf  %10.7lf  %10.7lf\n",

#else

                 "%3d %10.7f  %10.7f  %10.7f\n",

#endif

                atomNr,

                qm->xQM[i][XX]/0.1,

                qm->xQM[i][YY]/0.1,

                qm->xQM[i][ZZ]/0.1);

    }

    fprintf(out, "*\n");

    /* write the MM point charge data */

    if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) {

        /* name of the point charge file */

        sprintf(filename_ab, ".pc");

        fprintf(out, "%s%s%s\n", "%pointcharges \"", filename, "\"");

        iofile = fopen(filename, "w");

        fprintf(iofile, "%d\n", mm->nrMMatoms);

        for (i = 0; i < mm->nrMMatoms; i++) {

            fprintf(iofile,

#ifdef GMX_DOUBLE

                    "%8.4lf %10.7lf  %10.7lf  %10.7lf\n",

#else

                    "%8.4f %10.7f  %10.7f  %10.7f\n",

#endif

                    mm->MMcharges[i],

                    mm->xMM[i][XX]/0.1,

                    mm->xMM[i][YY]/0.1,

                    mm->xMM[i][ZZ]/0.1);

        }

        fprintf(iofile, "\n");

        fclose(iofile);

    }

    fprintf(out, "\n");

    sfree(filename);

    fclose(out);

}  /* write_orca_input */

real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr,

                      t_QMrec *qm, t_MMrec *mm)

{

    int i, j, orca_b_len;

    char buf[300], tmp[300];

    char *filename, *filename_ab, *err_msg1 = "Unexpected end of ORCA output";

    real QMener;

    FILE *ifile;

    int k;

    t_QMMMrec *QMMMrec;

    QMMMrec = fr->qr;

    orca_b_len = strlen(qm->orca_basename);

    snew(filename, orca_b_len+10);

    filename_ab = filename + orca_b_len; /* pointer in the filename just after the basename */

    memcpy(filename, qm->orca_basename, orca_b_len);

    /* in case of an optimization, the coordinates are printed in the

     * xyz file, the energy and gradients for the QM part are stored in the engrad file

     * and the gradients for the point charges are stored in the pc file.

     */

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

    if (qm->bTS || qm->bOPT) {

        sprintf(filename_ab, ".xyz");

        ifile = fopen(filename, "r");

        for(i=0;i<2;i++) {

            if (fgets(buf, 300, ifile) == NULL) {

                gmx_fatal(FARGS, err_msg1);

            }

        }

        for (i = 0; i < qm->nrQMatoms; i++) {

            if(fscanf(ifile,

#ifdef GMX_DOUBLE

                      "%s%lf%lf%lf",

#else

                      "%s%f%f%f",

#endif

                       tmp,&qm->xQM[i][XX],&qm->xQM[i][YY],&qm->xQM[i][ZZ])!=4) {

                break;

            }

            for (j = 0; j < DIM; j++) {

                qm->xQM[i][j] *= 0.1;

            }

        }

        fclose(ifile);

        if(i<qm->nrQMatoms) {

            gmx_fatal(FARGS, err_msg1);

        }

    }

    sprintf(filename_ab, ".engrad");

    ifile = fopen(filename, "r");

    /* we read the energy and the gradient for the qm-atoms from the engrad file */

    /* we can skip the first seven lines */

    for (j = 0; j < 7; j++) {

        if (fgets(buf, 300, ifile) == NULL) {

            gmx_fatal(FARGS, err_msg1);

        }

    }

    /* now comes the energy */

    if(fscanf(ifile,

#ifdef GMX_DOUBLE

              "%lf",

#else

              "%f",

#endif

              &QMener)!=1) {

        gmx_fatal(FARGS, err_msg1);

    }

    /* we can skip the next three lines (plus EOL from the previous one)*/

    for (j = 0; j < 4; j++) {

        if (fgets(buf, 300, ifile) == NULL) {

            gmx_fatal(FARGS, err_msg1);

        }

    }

    /* next lines contain the gradients of the QM atoms

     * now comes the gradient, one value per line:

     * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...

     */

    for (i = 0; i < qm->nrQMatoms; i++) {

        if(fscanf(ifile,

#ifdef GMX_DOUBLE

                  "%lf%lf%lf",

#else

                  "%f%f%f",

#endif

                  &QMgrad[i][XX],&QMgrad[i][YY],&QMgrad[i][ZZ])!=3) {

            break;

        }

    }

    fclose(ifile);

    if(i<qm->nrQMatoms) {

        gmx_fatal(FARGS, err_msg1);

    }

    /* write the MM point charge data */

    if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms) {

        sprintf(filename_ab, ".pcgrad");

        ifile = fopen(filename, "r");

        /* we read the gradient for the mm-atoms from the pcgrad file */

        /* we can skip the first line */

        if (fgets(buf, 300, ifile) == NULL) {

            gmx_fatal(FARGS, err_msg1);

        }

        for (i = 0; i < mm->nrMMatoms; i++) {

            if(fscanf(ifile,

#ifdef GMX_DOUBLE

                      "%lf%lf%lf",

#else

                      "%f%f%f",

#endif

                      &MMgrad[i][XX],&MMgrad[i][YY],&MMgrad[i][ZZ])!=3) {

                break;

            }

        }

        fclose(ifile);

        if(i<qm->nrQMatoms) {

            gmx_fatal(FARGS, err_msg1);

        }

    }

    sfree(filename);

    return QMener;

}

void do_orca(char *orca_dir, char *basename)

{

    /* make the call to the orca binary through system()

     * The location of the binary is set through the

     * environment.

     */

    char *buf;

    int len = strlen(orca_dir) + strlen(basename)*2 + 20;

    snew(buf,len);

    sprintf(buf,"%s/%s %s.inp >> %s.out",orca_dir,"orca",basename,basename);

    fprintf(stderr,"Calling '%s'\n",buf);

    if (system(buf)!= 0) { /* ahhh... this pain again...*/

        gmx_fatal(FARGS,"Call to '%s' failed\n",buf);

    }

    sfree(buf);

}

real call_orca(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])

{

    /* normal orca jobs */

    static int step = 0;

    int i, j;

    real QMener;

    rvec *QMgrad, *MMgrad;

    snew(QMgrad, qm->nrQMatoms);

    snew(MMgrad, mm->nrMMatoms);

    write_orca_input(fr, qm, mm);

    do_orca(qm->orca_dir, qm->orca_basename);

    QMener = read_orca_output(QMgrad, MMgrad, fr, qm, mm);

    /* put the QMMM forces in the force array and to the fshift

     */

    for (i = 0; i < qm->nrQMatoms; i++) {

        for (j = 0; j < DIM; j++) {

            f[i][j]      = HARTREE_BOHR2MD*QMgrad[i][j];

            fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];

        }

    }

    for (i = 0; i < mm->nrMMatoms; i++) {

        for (j = 0; j < DIM; j++) {

            f[i+qm->nrQMatoms][j]      = HARTREE_BOHR2MD*MMgrad[i][j];

            fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];

        }

    }

    QMener = QMener*HARTREE2KJ*AVOGADRO;

    step++;

    return(QMener);

} /* call_orca */

/* end of orca sub routines */