From: Erik Lindahl Date: Tue, 7 Jul 2015 20:50:41 +0000 (+0200) Subject: Move QMMM source to C++ X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?p=alexxy%2Fgromacs.git;a=commitdiff_plain;h=f7e76b5e296e23aa8de0103140811473b5515013 Move QMMM source to C++ Tested by enabling the four different QM/MM interfaces, but not linking to the actual programs since those calls are still unaltered. Change-Id: I86163a0db56003687e6226a78169141af98940b7 --- diff --git a/src/gromacs/mdlib/qm_gamess.c b/src/gromacs/mdlib/qm_gamess.cpp similarity index 95% rename from src/gromacs/mdlib/qm_gamess.c rename to src/gromacs/mdlib/qm_gamess.cpp index d913f7f489..b7c729b29d 100644 --- a/src/gromacs/mdlib/qm_gamess.c +++ b/src/gromacs/mdlib/qm_gamess.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, 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. @@ -55,6 +55,7 @@ #include "gromacs/legacyheaders/qmmm.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/typedefs.h" +#include "gromacs/legacyheaders/types/commrec.h" #include "gromacs/math/units.h" #include "gromacs/math/vec.h" #include "gromacs/utility/fatalerror.h" @@ -66,14 +67,14 @@ void -F77_FUNC(inigms, IMIGMS) (void); + F77_FUNC(inigms, IMIGMS) (void); void -F77_FUNC(endgms, ENDGMS) (void); + F77_FUNC(endgms, ENDGMS) (void); void -F77_FUNC(grads, GRADS) (int *nrqmat, real *qmcrd, int *nrmmat, real *mmchrg, - real *mmcrd, real *qmgrad, real *mmgrad, real *energy); + F77_FUNC(grads, GRADS) (int *nrqmat, real *qmcrd, int *nrmmat, real *mmchrg, + real *mmcrd, real *qmgrad, real *mmgrad, real *energy); @@ -88,7 +89,7 @@ void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm) * dynamics simulations. 7-6-2002 (London) */ int - i, j, rank; + i, j; FILE *out; char @@ -228,7 +229,7 @@ void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm) } } -real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, +real call_gamess(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]) { /* do the actual QMMM calculation using GAMESS-UK. In this @@ -237,7 +238,7 @@ real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, * gradient routines linked directly */ int - i, j, rank; + i, j; real QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy; t_QMMMrec diff --git a/src/gromacs/mdlib/qm_gaussian.c b/src/gromacs/mdlib/qm_gaussian.cpp similarity index 97% rename from src/gromacs/mdlib/qm_gaussian.c rename to src/gromacs/mdlib/qm_gaussian.cpp index d90232c7e8..6f519c26b3 100644 --- a/src/gromacs/mdlib/qm_gaussian.c +++ b/src/gromacs/mdlib/qm_gaussian.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, 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. @@ -57,20 +57,19 @@ #include "gromacs/legacyheaders/typedefs.h" #include "gromacs/math/units.h" #include "gromacs/math/vec.h" +#include "gromacs/utility/cstringutil.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" - /* TODO: this should be made thread-safe */ /* Gaussian interface routines */ -void init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm) +void init_gaussian(t_QMrec *qm) { - FILE - *rffile = NULL, *out = NULL; + FILE *out = NULL; ivec - basissets[eQMbasisNR] = {{0, 3, 0}, + basissets[eQMbasisNR] = {{0, 3, 0}, {0, 3, 0}, /*added for double sto-3g entry in names.c*/ {5, 0, 0}, {5, 0, 1}, @@ -81,9 +80,9 @@ void init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm) {1, 6, 11}, {4, 6, 0}}; char - *buf = NULL; + *buf = NULL; int - i; + i; /* using the ivec above to convert the basis read form the mdp file * in a human readable format into some numbers for the gaussian @@ -695,8 +694,7 @@ void write_gaussian_input(int step, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm) } /* write_gaussian_input */ -real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], int step, - t_QMrec *qm, t_MMrec *mm) +real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], t_QMrec *qm, t_MMrec *mm) { int i, j, atnum; @@ -801,8 +799,7 @@ real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], int step, return(QMener); } -real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step, - gmx_bool swapped, t_QMrec *qm, t_MMrec *mm) +real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step, t_QMrec *qm, t_MMrec *mm) { int i; @@ -1023,8 +1020,7 @@ void do_gaussian(int step, char *exe) } } -real call_gaussian(t_commrec *cr, t_forcerec *fr, - t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]) +real call_gaussian(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]) { /* normal gaussian jobs */ static int @@ -1045,7 +1041,7 @@ real call_gaussian(t_commrec *cr, t_forcerec *fr, write_gaussian_input(step, fr, qm, mm); do_gaussian(step, exe); - QMener = read_gaussian_output(QMgrad, MMgrad, step, qm, mm); + QMener = read_gaussian_output(QMgrad, MMgrad, qm, mm); /* put the QMMM forces in the force array and to the fshift */ for (i = 0; i < qm->nrQMatoms; i++) @@ -1071,8 +1067,7 @@ real call_gaussian(t_commrec *cr, t_forcerec *fr, } /* call_gaussian */ -real call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, - rvec f[], rvec fshift[]) +real call_gaussian_SH(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]) { /* a gaussian call routine intended for doing diabatic surface * "sliding". See the manual for the theoretical background of this @@ -1128,7 +1123,7 @@ real call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, write_gaussian_SH_input(step, swapped, fr, qm, mm); do_gaussian(step, exe); - QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, swapped, qm, mm); + QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm); /* check for a surface hop. Only possible if we were already state * averaging. @@ -1149,7 +1144,7 @@ real call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, { write_gaussian_SH_input(step, swapped, fr, qm, mm); do_gaussian(step, exe); - QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, swapped, qm, mm); + QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm); } } /* add the QMMM forces to the gmx force array and fshift diff --git a/src/gromacs/mdlib/qm_mopac.c b/src/gromacs/mdlib/qm_mopac.cpp similarity index 91% rename from src/gromacs/mdlib/qm_mopac.c rename to src/gromacs/mdlib/qm_mopac.cpp index 6eb689cd3e..7496965bf0 100644 --- a/src/gromacs/mdlib/qm_mopac.c +++ b/src/gromacs/mdlib/qm_mopac.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, 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. @@ -63,16 +63,16 @@ /* mopac interface routines */ void -F77_FUNC(domldt, DOMLDT) (int *nrqmat, int labels[], char keywords[]); + F77_FUNC(domldt, DOMLDT) (int *nrqmat, int labels[], char keywords[]); void -F77_FUNC(domop, DOMOP) (int *nrqmat, double qmcrd[], int *nrmmat, - double mmchrg[], double mmcrd[], double qmgrad[], - double mmgrad[], double *energy, double qmcharges[]); + F77_FUNC(domop, DOMOP) (int *nrqmat, double qmcrd[], int *nrmmat, + double mmchrg[], double mmcrd[], double qmgrad[], + double mmgrad[], double *energy, double qmcharges[]); -void init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm) +void init_mopac(t_QMrec *qm) { /* initializes the mopac routines ans sets up the semiempirical * computation by calling moldat(). The inline mopac routines can @@ -104,8 +104,7 @@ void init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm) } /* init_mopac */ -real call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, - rvec f[], rvec fshift[]) +real call_mopac(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]) { /* do the actual QMMM calculation using directly linked mopac subroutines */ @@ -167,8 +166,7 @@ real call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, return (QMener); } -real call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, - rvec f[], rvec fshift[]) +real call_mopac_SH(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]) { /* do the actual SH QMMM calculation using directly linked mopac subroutines */ diff --git a/src/gromacs/mdlib/qm_orca.c b/src/gromacs/mdlib/qm_orca.cpp similarity index 98% rename from src/gromacs/mdlib/qm_orca.c rename to src/gromacs/mdlib/qm_orca.cpp index 11ddfd325b..2f121806fd 100644 --- a/src/gromacs/mdlib/qm_orca.c +++ b/src/gromacs/mdlib/qm_orca.cpp @@ -3,7 +3,7 @@ * * 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 + * Copyright (c) 2013,2014,2015, 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. @@ -284,7 +284,7 @@ real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr, int i, j, atnum; char - buf[300], tmp[300], orca_xyzFilename[300], orca_pcgradFilename[300], orca_engradFilename[300]; + buf[300], orca_xyzFilename[300], orca_pcgradFilename[300], orca_engradFilename[300]; real QMener; FILE @@ -320,8 +320,8 @@ real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr, gmx_fatal(FARGS, "Unexpected end of ORCA output"); } #ifdef GMX_DOUBLE - sscanf(buf, "%s%lf%lf%lf\n", - tmp, + sscanf(buf, "%d%lf%lf%lf\n", + &atnum, &qm->xQM[i][XX], &qm->xQM[i][YY], &qm->xQM[i][ZZ]); diff --git a/src/gromacs/mdlib/qmmm.c b/src/gromacs/mdlib/qmmm.cpp similarity index 95% rename from src/gromacs/mdlib/qmmm.c rename to src/gromacs/mdlib/qmmm.cpp index bb83405a55..cf7f7b3440 100644 --- a/src/gromacs/mdlib/qmmm.c +++ b/src/gromacs/mdlib/qmmm.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, 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. @@ -45,6 +45,10 @@ #include #include +#include + +#include + #include "gromacs/fileio/confio.h" #include "gromacs/legacyheaders/force.h" #include "gromacs/legacyheaders/macros.h" @@ -73,36 +77,32 @@ void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm); real -call_gamess(t_commrec *cr, t_forcerec *fr, +call_gamess(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]); #elif defined GMX_QMMM_MOPAC /* MOPAC interface */ void -init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm); +init_mopac(t_QMrec *qm); real -call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, - t_MMrec *mm, rvec f[], rvec fshift[]); +call_mopac(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]); real -call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, - t_MMrec *mm, rvec f[], rvec fshift[]); +call_mopac_SH(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]); #elif defined GMX_QMMM_GAUSSIAN /* GAUSSIAN interface */ void -init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm); +init_gaussian(t_QMrec *qm); real -call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, - t_MMrec *mm, rvec f[], rvec fshift[]); +call_gaussian_SH(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]); real -call_gaussian(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, - t_MMrec *mm, rvec f[], rvec fshift[]); +call_gaussian(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]); #elif defined GMX_QMMM_ORCA /* ORCA interface */ @@ -165,11 +165,11 @@ real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec #ifdef GMX_QMMM_MOPAC if (qm->bSH) { - QMener = call_mopac_SH(cr, fr, qm, mm, f, fshift); + QMener = call_mopac_SH(qm, mm, f, fshift); } else { - QMener = call_mopac(cr, fr, qm, mm, f, fshift); + QMener = call_mopac(qm, mm, f, fshift); } #else gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac."); @@ -181,7 +181,7 @@ real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec if (qm->bSH && qm->QMmethod == eQMmethodCASSCF) { #ifdef GMX_QMMM_GAUSSIAN - QMener = call_gaussian_SH(cr, fr, qm, mm, f, fshift); + QMener = call_gaussian_SH(fr, qm, mm, f, fshift); #else gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian."); #endif @@ -189,9 +189,9 @@ real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec else { #ifdef GMX_QMMM_GAMESS - QMener = call_gamess(cr, fr, qm, mm, f, fshift); + QMener = call_gamess(fr, qm, mm, f, fshift); #elif defined GMX_QMMM_GAUSSIAN - QMener = call_gaussian(cr, fr, qm, mm, f, fshift); + QMener = call_gaussian(fr, qm, mm, f, fshift); #elif defined GMX_QMMM_ORCA QMener = call_orca(fr, qm, mm, f, fshift); #else @@ -210,7 +210,7 @@ void init_QMroutine(t_commrec gmx_unused *cr, t_QMrec gmx_unused *qm, t_MMrec gm { #ifdef GMX_QMMM_MOPAC /* do a semi-empiprical calculation */ - init_mopac(cr, qm, mm); + init_mopac(qm); #else gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac."); #endif @@ -221,7 +221,7 @@ void init_QMroutine(t_commrec gmx_unused *cr, t_QMrec gmx_unused *qm, t_MMrec gm #ifdef GMX_QMMM_GAMESS init_gamess(cr, qm, mm); #elif defined GMX_QMMM_GAUSSIAN - init_gaussian(cr, qm, mm); + init_gaussian(qm); #elif defined GMX_QMMM_ORCA init_orca(qm); #else @@ -262,7 +262,7 @@ static void punch_QMMM_excl(t_QMrec *qm, t_MMrec *mm, t_blocka *excls) FILE *out = NULL; int - i, j, k, nrexcl = 0, *excluded = NULL, max = 0; + i, j, k, nrexcl = 0, *excluded = NULL, max_excl = 0; out = fopen("QMMMexcl.dat", "w"); @@ -280,10 +280,10 @@ static void punch_QMMM_excl(t_QMrec *qm, t_MMrec *mm, t_blocka *excls) { if (mm->indexMM[k] == excls->a[j]) /* the excluded MM atom */ { - if (nrexcl >= max) + if (nrexcl >= max_excl) { - max += 1000; - srenew(excluded, max); + max_excl += 1000; + srenew(excluded, max_excl); } excluded[nrexcl++] = k; continue; @@ -485,8 +485,8 @@ void init_QMMMrec(t_commrec *cr, t_ilist *ilist_mol; gmx_mtop_atomlookup_t alook; - c6au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 6)); - c12au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 12)); + c6au = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 6)); + c12au = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 12)); /* issue a fatal if the user wants to run with more than one node */ if (PAR(cr)) { @@ -753,7 +753,7 @@ void init_QMMMrec(t_commrec *cr, { #ifdef GMX_QMMM_MOPAC /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */ - init_mopac(cr, qr->qm[0], qr->mm); + init_mopac(qr->qm[0]); #else gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac."); #endif @@ -764,7 +764,7 @@ void init_QMMMrec(t_commrec *cr, #ifdef GMX_QMMM_GAMESS init_gamess(cr, qr->qm[0], qr->mm); #elif defined GMX_QMMM_GAUSSIAN - init_gaussian(cr, qr->qm[0], qr->mm); + init_gaussian(qr->qm[0]); #elif defined GMX_QMMM_ORCA init_orca(qr->qm[0]); #else @@ -797,8 +797,6 @@ void update_QMMMrec(t_commrec *cr, QMMMlist; rvec dx, crd; - int - *MMatoms; t_QMrec *qm; t_MMrec @@ -810,8 +808,8 @@ void update_QMMMrec(t_commrec *cr, real c12au, c6au; - c6au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 6)); - c12au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 12)); + c6au = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 6)); + c12au = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 12)); /* every cpu has this array. On every processor we fill this array * with 1's and 0's. 1's indicate the atoms is a QM atom on the @@ -873,7 +871,7 @@ void update_QMMMrec(t_commrec *cr, crd[0] = IS2X(QMMMlist.shift[i]) + IS2X(qm_i_particles[i].shift); crd[1] = IS2Y(QMMMlist.shift[i]) + IS2Y(qm_i_particles[i].shift); crd[2] = IS2Z(QMMMlist.shift[i]) + IS2Z(qm_i_particles[i].shift); - is = XYZ2IS(crd[0], crd[1], crd[2]); + is = static_cast(XYZ2IS(crd[0], crd[1], crd[2])); for (j = QMMMlist.jindex[i]; j < QMMMlist.jindex[i+1]; j++) @@ -897,9 +895,13 @@ void update_QMMMrec(t_commrec *cr, qsort(qm_i_particles, QMMMlist.nri, (size_t)sizeof(qm_i_particles[0]), struct_comp); - qsort(mm_j_particles, mm_nr, - (size_t)sizeof(mm_j_particles[0]), - struct_comp); + /* The mm_j_particles argument to qsort is not allowed to be NULL */ + if (mm_nr > 0) + { + qsort(mm_j_particles, mm_nr, + (size_t)sizeof(mm_j_particles[0]), + struct_comp); + } /* remove multiples in the QM shift array, since in init_QMMM() we * went through the atom numbers from 0 to md.nr, the order sorted * here matches the one of QMindex already.