Remove qm dependence on forcerec in do_force_lowlevel
authorejjordan <e.jjordan12@gmail.com>
Wed, 13 Nov 2019 17:07:52 +0000 (18:07 +0100)
committerMagnus Lundborg <magnus.lundborg@scilifelab.se>
Thu, 14 Nov 2019 16:47:01 +0000 (17:47 +0100)
Also cleaned up some initialization.

Change-Id: Iccc5e80ace52409d734089714c49302c2f0fb994

src/gromacs/mdlib/force.cpp
src/gromacs/mdlib/qm_gaussian.cpp
src/gromacs/mdlib/qm_gaussian.h
src/gromacs/mdlib/qm_orca.cpp
src/gromacs/mdlib/qm_orca.h
src/gromacs/mdlib/qmmm.cpp
src/gromacs/mdlib/qmmm.h

index f659874b92a704244cb4c6b0d92dc9a50f89f176..8b36c27e8c51971d255494c600142101a6e9dcbe 100644 (file)
@@ -126,7 +126,7 @@ void do_force_lowlevel(t_forcerec*                         fr,
     /* do QMMM first if requested */
     if (fr->bQMMM)
     {
-        enerd->term[F_EQM] = calculate_QMMM(cr, &forceOutputs->forceWithShiftForces(), fr);
+        enerd->term[F_EQM] = calculate_QMMM(cr, &forceOutputs->forceWithShiftForces(), fr->qr);
     }
 
     /* Call the short range functions all in one go. */
index b9885bfff340c261b5e4abe841d0b41512292731..60fb57d0f2ddfb73ec7c000fa11a366401ee1a4b 100644 (file)
@@ -52,7 +52,6 @@
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/force.h"
-#include "gromacs/mdlib/forcerec.h"
 #include "gromacs/mdlib/qmmm.h"
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/utility/cstringutil.h"
@@ -216,16 +215,11 @@ void init_gaussian(t_QMrec* qm)
 }
 
 
-static void write_gaussian_SH_input(int step, gmx_bool swap, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
+static void write_gaussian_SH_input(int step, gmx_bool swap, const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
 {
-    int        i;
-    gmx_bool   bSA;
-    FILE*      out;
-    t_QMMMrec* QMMMrec;
-    QMMMrec = fr->qr;
-    bSA     = (qm->SAstep > 0);
-
-    out = fopen("input.com", "w");
+    int   i;
+    bool  bSA = (qm->SAstep > 0);
+    FILE* out = fopen("input.com", "w");
     /* write the route */
     fprintf(out, "%s", "%scr=input\n");
     fprintf(out, "%s", "%rwf=input\n");
@@ -380,14 +374,11 @@ static void write_gaussian_SH_input(int step, gmx_bool swap, const t_forcerec* f
     fclose(out);
 } /* write_gaussian_SH_input */
 
-static void write_gaussian_input(int step, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
+static void write_gaussian_input(int step, const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
 {
-    int        i;
-    t_QMMMrec* QMMMrec;
-    FILE*      out;
+    int i;
 
-    QMMMrec = fr->qr;
-    out     = fopen("input.com", "w");
+    FILE* out = fopen("input.com", "w");
     /* write the route */
 
     if (qm->QMmethod >= eQMmethodRHF)
@@ -758,7 +749,7 @@ static void do_gaussian(int step, char* exe)
     }
 }
 
-real call_gaussian(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
+real call_gaussian(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
 {
     /* normal gaussian jobs */
     static int step = 0;
@@ -772,7 +763,7 @@ real call_gaussian(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rve
     snew(QMgrad, qm->nrQMatoms);
     snew(MMgrad, mm->nrMMatoms);
 
-    write_gaussian_input(step, fr, qm, mm);
+    write_gaussian_input(step, qmmm, qm, mm);
     do_gaussian(step, exe);
     QMener = read_gaussian_output(QMgrad, MMgrad, qm, mm);
     /* put the QMMM forces in the force array and to the fshift
@@ -800,7 +791,7 @@ real call_gaussian(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rve
 
 } /* call_gaussian */
 
-real call_gaussian_SH(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
+real call_gaussian_SH(const t_QMMMrec* qmmm, 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
@@ -845,7 +836,7 @@ real call_gaussian_SH(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[],
     /*  if(!step)
      * qr->bSA=FALSE;*/
     /* temporray set to step + 1, since there is a chk start */
-    write_gaussian_SH_input(step, swapped, fr, qm, mm);
+    write_gaussian_SH_input(step, swapped, qmmm, qm, mm);
 
     do_gaussian(step, exe);
     QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
@@ -867,7 +858,7 @@ real call_gaussian_SH(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[],
         }
         if (swap) /* change surface, so do another call */
         {
-            write_gaussian_SH_input(step, swapped, fr, qm, mm);
+            write_gaussian_SH_input(step, swapped, qmmm, qm, mm);
             do_gaussian(step, exe);
             QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
         }
index 1933fc7893c39d628e4ef264da49f30fc8063a54..9c3088413232c801a9420865f64fd12635312ee3 100644 (file)
@@ -50,23 +50,23 @@ void init_gaussian(t_QMrec* qm);
 /*! \brief
  * Call gaussian to do qm calculation.
  *
- * \param[in] fr Global forcerec.
- * \param[in] qm QM part of forcerec.
- * \param[in] mm mm part of forcerec.
- * \param[in] f  force vector.
+ * \param[in] qmmm   QMMM part forcerec.
+ * \param[in] qm     QM part of forcerec.
+ * \param[in] mm     mm part of forcerec.
+ * \param[in] f      force vector.
  * \param[in] fshift shift of force vector.
  */
-real call_gaussian(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
+real call_gaussian(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
 
 /*! \brief
  * Call gaussian SH(?) to do qm calculation.
  *
- * \param[in] fr Global forcerec.
- * \param[in] qm QM part of forcerec.
- * \param[in] mm mm part of forcerec.
- * \param[in] f  force vector.
+ * \param[in] qmmm   QMMM part forcerec.
+ * \param[in] qm     QM part of forcerec.
+ * \param[in] mm     mm part of forcerec.
+ * \param[in] f      force vector.
  * \param[in] fshift shift of force vector.
  */
-real call_gaussian_SH(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
+real call_gaussian_SH(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
 
 #endif
index 075d8d68340f8e5f296ddb22f31a0fbde91430be..1bd0cca0738566ae0d625a27866189e5d29e3e0e 100644 (file)
@@ -51,7 +51,6 @@
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/qmmm.h"
-#include "gromacs/mdtypes/forcerec.h"
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
@@ -110,19 +109,16 @@ void init_orca(t_QMrec* qm)
 }
 
 
-static void write_orca_input(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
+static void write_orca_input(const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
 {
-    int        i;
-    t_QMMMrec* QMMMrec;
-    FILE *     out, *pcFile, *addInputFile;
-    char *     buf, *orcaInput, *addInputFilename, *pcFilename;
-
-    QMMMrec = fr->qr;
+    int   i;
+    FILE *pcFile, *addInputFile;
+    char *buf, *orcaInput, *addInputFilename, *pcFilename;
 
     /* write the first part of the input-file */
     snew(orcaInput, 200);
     sprintf(orcaInput, "%s.inp", qm->orca_basename);
-    out = fopen(orcaInput, "w");
+    FILE* out = fopen(orcaInput, "w");
 
     snew(addInputFilename, 200);
     sprintf(addInputFilename, "%s.ORCAINFO", qm->orca_basename);
@@ -193,15 +189,13 @@ static void write_orca_input(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
     fclose(out);
 } /* write_orca_input */
 
-static real read_orca_output(rvec QMgrad[], rvec MMgrad[], const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
+static real read_orca_output(rvec QMgrad[], rvec MMgrad[], const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
 {
-    int        i, j;
-    char       buf[300], orca_pcgradFilename[300], orca_engradFilename[300];
-    real       QMener;
-    FILE *     pcgrad, *engrad;
-    int        k;
-    t_QMMMrec* QMMMrec;
-    QMMMrec = fr->qr;
+    int   i, j;
+    char  buf[300], orca_pcgradFilename[300], orca_engradFilename[300];
+    real  QMener;
+    FILE *pcgrad, *engrad;
+    int   k;
 
     /* 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.
@@ -328,7 +322,7 @@ static void do_orca(char* orca_dir, char* basename)
     }
 }
 
-real call_orca(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
+real call_orca(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
 {
     /* normal orca jobs */
     static int step = 0;
@@ -342,9 +336,9 @@ real call_orca(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fs
     snew(QMgrad, qm->nrQMatoms);
     snew(MMgrad, mm->nrMMatoms);
 
-    write_orca_input(fr, qm, mm);
+    write_orca_input(qmmm, qm, mm);
     do_orca(qm->orca_dir, qm->orca_basename);
-    QMener = read_orca_output(QMgrad, MMgrad, fr, qm, mm);
+    QMener = read_orca_output(QMgrad, MMgrad, qmmm, qm, mm);
     /* put the QMMM forces in the force array and to the fshift
      */
     for (i = 0; i < qm->nrQMatoms; i++)
index cf89d941c8fc907e3d15d3331147fc9172253487..c6a3b55917bb8e9038bf5d8838b21440b4cca37d 100644 (file)
@@ -39,6 +39,6 @@
 
 void init_orca(t_QMrec* qm);
 
-real call_orca(const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
+real call_orca(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
 
 #endif
index 2482bf592712f7d86d2c09f694cbc7ddfbea4dab..24b3c75d2adbb5f0268a941f659659b06c3c6b4a 100644 (file)
@@ -95,7 +95,7 @@ static bool struct_comp(const t_j_particle& a, const t_j_particle& b)
 }
 
 static real call_QMroutine(const t_commrec gmx_unused* cr,
-                           const t_forcerec gmx_unused* fr,
+                           const t_QMMMrec gmx_unused* qmmm,
                            t_QMrec gmx_unused* qm,
                            t_MMrec gmx_unused* mm,
                            rvec gmx_unused f[],
@@ -131,7 +131,7 @@ static real call_QMroutine(const t_commrec gmx_unused* cr,
         {
             if (GMX_QMMM_GAUSSIAN)
             {
-                return call_gaussian_SH(fr, qm, mm, f, fshift);
+                return call_gaussian_SH(qmmm, qm, mm, f, fshift);
             }
             else
             {
@@ -146,11 +146,11 @@ static real call_QMroutine(const t_commrec gmx_unused* cr,
             }
             else if (GMX_QMMM_GAUSSIAN)
             {
-                return call_gaussian(fr, qm, mm, f, fshift);
+                return call_gaussian(qmmm, qm, mm, f, fshift);
             }
             else if (GMX_QMMM_ORCA)
             {
-                return call_orca(fr, qm, mm, f, fshift);
+                return call_orca(qmmm, qm, mm, f, fshift);
             }
             else
             {
@@ -827,17 +827,16 @@ void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, co
     }
 } /* update_QMMM_rec */
 
-real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_forcerec* fr)
+real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qr)
 {
     real QMener = 0.0;
     /* a selection for the QM package depending on which is requested
      * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
      * it works through defines.... Not so nice yet
      */
-    t_QMMMrec* qr;
-    t_QMrec *  qm, *qm2;
-    t_MMrec*   mm     = nullptr;
-    rvec *     forces = nullptr, *fshift = nullptr, *forces2 = nullptr,
+    t_QMrec *qm, *qm2;
+    t_MMrec* mm     = nullptr;
+    rvec *   forces = nullptr, *fshift = nullptr, *forces2 = nullptr,
          *fshift2 = nullptr; /* needed for multilayer ONIOM */
     int i, j, k;
 
@@ -848,7 +847,6 @@ real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShi
 
     /* make a local copy the QMMMrec pointer
      */
-    qr = fr->qr;
     mm = qr->mm;
 
     /* now different procedures are carried out for one layer ONION and
@@ -861,7 +859,7 @@ real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShi
         qm = qr->qm[0];
         snew(forces, (qm->nrQMatoms + mm->nrMMatoms));
         snew(fshift, (qm->nrQMatoms + mm->nrMMatoms));
-        QMener = call_QMroutine(cr, fr, qm, mm, forces, fshift);
+        QMener = call_QMroutine(cr, qr, qm, mm, forces, fshift);
         for (i = 0; i < qm->nrQMatoms; i++)
         {
             for (j = 0; j < DIM; j++)
@@ -907,13 +905,13 @@ real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShi
             srenew(fshift, qm->nrQMatoms);
             /* we need to re-initialize the QMroutine every step... */
             init_QMroutine(cr, qm, mm);
-            QMener += call_QMroutine(cr, fr, qm, mm, forces, fshift);
+            QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
 
             /* this layer at the lower level of theory */
             srenew(forces2, qm->nrQMatoms);
             srenew(fshift2, qm->nrQMatoms);
             init_QMroutine(cr, qm2, mm);
-            QMener -= call_QMroutine(cr, fr, qm2, mm, forces2, fshift2);
+            QMener -= call_QMroutine(cr, qr, qm2, mm, forces2, fshift2);
             /* E = E1high-E1low The next layer includes the current layer at
              * the lower level of theory, which provides + E2low
              * this is similar for gradients
@@ -933,7 +931,7 @@ real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShi
         init_QMroutine(cr, qm, mm);
         srenew(forces, qm->nrQMatoms);
         srenew(fshift, qm->nrQMatoms);
-        QMener += call_QMroutine(cr, fr, qm, mm, forces, fshift);
+        QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
         for (i = 0; i < qm->nrQMatoms; i++)
         {
             for (j = 0; j < DIM; j++)
index 8172e9766007e9b1f51e3c48f455e8dbc4a0facd..8de4084926469fa39c73dde1186c99d88bfce919 100644 (file)
@@ -139,7 +139,7 @@ void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, co
  * routine should be called at every step, since it updates the MM
  * elements of the t_QMMMrec struct.
  */
-real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_forcerec* fr);
+real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qmmm);
 
 /* QMMM computes the QM forces. This routine makes either function
  * calls to gmx QM routines (derived from MOPAC7 (semi-emp.) and MPQC