src/gromacs/applied_forces/qmmm/qmmminputgenerator.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements input generator class for CP2K QMMM
  38  *
  39  * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
  40  * \ingroup module_applied_forces
  41  */
  42
  43 #include "gmxpre.h"
  44
  45 #include "qmmminputgenerator.h"
  46
  47 #include "gromacs/math/units.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/utility/stringutil.h"
  50
  51 namespace gmx
  52 {
  53
  54 QMMMInputGenerator::QMMMInputGenerator(const QMMMParameters& parameters,
  55                                        PbcType               pbcType,
  56                                        const matrix          box,
  57                                        ArrayRef<const real>  q,
  58                                        ArrayRef<const RVec>  x) :
  59     parameters_(parameters),
  60     pbc_(pbcType),
  61     qmBox_{ { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } },
  62     qmCenter_{ 0.0, 0.0, 0.0 },
  63     qmTrans_{ 0.0, 0.0, 0.0 },
  64     q_(q),
  65     x_(x)
  66 {
  67     // Copy box
  68     copy_mat(box, box_);
  69
  70     // Fill qmAtoms_ set
  71     for (const auto& val : parameters_.qmIndices_)
  72     {
  73         qmAtoms_.emplace(val);
  74     }
  75
  76     // Compute QM box
  77     computeQMBox(sc_qmBoxScale, sc_qmBoxMinLength);
  78 }
  79
  80 bool QMMMInputGenerator::isQMAtom(index globalAtomIndex) const
  81 {
  82     return (qmAtoms_.find(globalAtomIndex) != qmAtoms_.end());
  83 }
  84
  85 void QMMMInputGenerator::computeQMBox(real scale, real minNorm)
  86 {
  87     // Init atom numbers
  88     size_t nQm = parameters_.qmIndices_.size();
  89
  90     // If there is only one QM atom, then just copy the box_
  91     if (nQm < 2)
  92     {
  93         copy_mat(box_, qmBox_);
  94         qmCenter_ = x_[parameters_.qmIndices_[0]];
  95         qmTrans_  = RVec(qmBox_[0]) / 2 + RVec(qmBox_[1]) / 2 + RVec(qmBox_[2]) / 2 - qmCenter_;
  96         return;
  97     }
  98
  99     // Initialize pbc
 100     t_pbc pbc;
 101     set_pbc(&pbc, pbc_, box_);
 102
 103     /* To compute qmBox_:
 104      * 1) Compute maximum dimension - maxDist betweeen QM atoms within PBC
 105      * 2) Make projection of the each box_ vector onto the cross prod of other two vectors
 106      * 3) Calculate scales so that norm will be maxDist for each box_ vector
 107      * 4) Apply scales to the box_ to get qmBox_
 108      */
 109     RVec dx(0.0, 0.0, 0.0);
 110     real maxDist = 0.0;
 111
 112     // Search for the maxDist - maximum distance within QM system
 113     for (size_t i = 0; i < nQm - 1; i++)
 114     {
 115         for (size_t j = i + 1; j < nQm; j++)
 116         {
 117             pbc_dx(&pbc, x_[parameters_.qmIndices_[i]], x_[parameters_.qmIndices_[j]], dx);
 118             maxDist = dx.norm() > maxDist ? dx.norm() : maxDist;
 119         }
 120     }
 121
 122     // Apply scale factor: qmBox_ should be *scale times bigger than maxDist
 123     maxDist *= scale;
 124
 125     // Compute QM Box vectors
 126     RVec vec0 = computeQMBoxVec(box_[0], box_[1], box_[2], maxDist, minNorm, norm(box_[0]));
 127     RVec vec1 = computeQMBoxVec(box_[1], box_[0], box_[2], maxDist, minNorm, norm(box_[1]));
 128     RVec vec2 = computeQMBoxVec(box_[2], box_[0], box_[1], maxDist, minNorm, norm(box_[2]));
 129
 130     // Form qmBox_
 131     copy_rvec(vec0, qmBox_[0]);
 132     copy_rvec(vec1, qmBox_[1]);
 133     copy_rvec(vec2, qmBox_[2]);
 134
 135     /* Now we need to also compute translation vector.
 136      * In order to center QM atoms in the computed qmBox_
 137      *
 138      * First compute center of QM system by averaging PBC-aware distance vectors
 139      * with respect to the first QM atom.
 140      */
 141     for (size_t i = 1; i < nQm; i++)
 142     {
 143         // compute pbc-aware distance vector between QM atom 0 and QM atom i
 144         pbc_dx(&pbc, x_[parameters_.qmIndices_[i]], x_[parameters_.qmIndices_[0]], dx);
 145
 146         // add that to the qmCenter_
 147         qmCenter_ += dx;
 148     }
 149
 150     // Average over all nQm atoms and add first atom coordiantes
 151     qmCenter_ = qmCenter_ / nQm + x_[parameters_.qmIndices_[0]];
 152
 153     // Translation vector will be center of qmBox_ - qmCenter_
 154     qmTrans_ = vec0 / 2 + vec1 / 2 + vec2 / 2 - qmCenter_;
 155 }
 156
 157 std::string QMMMInputGenerator::generateGlobalSection()
 158 {
 159     std::string res;
 160
 161     res += "&GLOBAL\n";
 162     res += "  PRINT_LEVEL LOW\n";
 163     res += "  PROJECT GROMACS\n";
 164     res += "  RUN_TYPE ENERGY_FORCE\n";
 165     res += "&END GLOBAL\n";
 166
 167     return res;
 168 }
 169
 170 std::string QMMMInputGenerator::generateDFTSection() const
 171 {
 172     std::string res;
 173
 174     res += "  &DFT\n";
 175
 176     // write charge and multiplicity
 177     res += formatString("    CHARGE %d\n", parameters_.qmCharge_);
 178     res += formatString("    MULTIPLICITY %d\n", parameters_.qmMult_);
 179
 180     // If multiplicity is not 1 then we should use unrestricted Kohn-Sham
 181     if (parameters_.qmMult_ > 1)
 182     {
 183         res += "    UKS\n";
 184     }
 185
 186     // Basis files, Grid setup and SCF parameters
 187     res += "    BASIS_SET_FILE_NAME  BASIS_MOLOPT\n";
 188     res += "    POTENTIAL_FILE_NAME  POTENTIAL\n";
 189     res += "    &MGRID\n";
 190     res += "      NGRIDS 5\n";
 191     res += "      CUTOFF 450\n";
 192     res += "      REL_CUTOFF 50\n";
 193     res += "      COMMENSURATE\n";
 194     res += "    &END MGRID\n";
 195     res += "    &SCF\n";
 196     res += "      SCF_GUESS RESTART\n";
 197     res += "      EPS_SCF 5.0E-8\n";
 198     res += "      MAX_SCF 20\n";
 199     res += "      &OT  T\n";
 200     res += "        MINIMIZER  DIIS\n";
 201     res += "        STEPSIZE   0.15\n";
 202     res += "        PRECONDITIONER FULL_ALL\n";
 203     res += "      &END OT\n";
 204     res += "      &OUTER_SCF  T\n";
 205     res += "        MAX_SCF 20\n";
 206     res += "        EPS_SCF 5.0E-8\n";
 207     res += "      &END OUTER_SCF\n";
 208     res += "    &END SCF\n";
 209
 210     // DFT functional parameters
 211     res += "    &XC\n";
 212     res += "      DENSITY_CUTOFF     1.0E-12\n";
 213     res += "      GRADIENT_CUTOFF    1.0E-12\n";
 214     res += "      TAU_CUTOFF         1.0E-12\n";
 215     res += formatString("      &XC_FUNCTIONAL %s\n", c_qmmmQMMethodNames[parameters_.qmMethod_]);
 216     res += "      &END XC_FUNCTIONAL\n";
 217     res += "    &END XC\n";
 218     res += "    &QS\n";
 219     res += "     METHOD GPW\n";
 220     res += "     EPS_DEFAULT 1.0E-10\n";
 221     res += "     EXTRAPOLATION ASPC\n";
 222     res += "     EXTRAPOLATION_ORDER  4\n";
 223     res += "    &END QS\n";
 224     res += "  &END DFT\n";
 225
 226     return res;
 227 }
 228
 229 std::string QMMMInputGenerator::generateQMMMSection() const
 230 {
 231     std::string res;
 232
 233     // Init some numbers
 234     const size_t nQm = parameters_.qmIndices_.size();
 235
 236     // Count the numbers of individual QM atoms per type
 237     std::vector<int> num_atoms(periodic_system.size(), 0);
 238     for (size_t i = 0; i < nQm; i++)
 239     {
 240         num_atoms[parameters_.atomNumbers_[parameters_.qmIndices_[i]]]++;
 241     }
 242
 243     res += "  &QMMM\n";
 244
 245     // QM cell
 246     res += "    &CELL\n";
 247     res += formatString(
 248             "      A %.3lf %.3lf %.3lf\n", qmBox_[0][0] * 10, qmBox_[0][1] * 10, qmBox_[0][2] * 10);
 249     res += formatString(
 250             "      B %.3lf %.3lf %.3lf\n", qmBox_[1][0] * 10, qmBox_[1][1] * 10, qmBox_[1][2] * 10);
 251     res += formatString(
 252             "      C %.3lf %.3lf %.3lf\n", qmBox_[2][0] * 10, qmBox_[2][1] * 10, qmBox_[2][2] * 10);
 253     res += "      PERIODIC XYZ\n";
 254     res += "    &END CELL\n";
 255
 256     res += "    CENTER EVERY_STEP\n";
 257     res += "    CENTER_GRID TRUE\n";
 258     res += "    &WALLS\n";
 259     res += "      TYPE REFLECTIVE\n";
 260     res += "    &END WALLS\n";
 261
 262     res += "    ECOUPL GAUSS\n";
 263     res += "    USE_GEEP_LIB 12\n";
 264     res += "    &PERIODIC\n";
 265     res += "      GMAX     1.0E+00\n";
 266     res += "      &MULTIPOLE ON\n";
 267     res += "         RCUT     1.0E+01\n";
 268     res += "         EWALD_PRECISION     1.0E-06\n";
 269     res += "      &END\n";
 270     res += "    &END PERIODIC\n";
 271
 272     // Print indicies of QM atoms
 273     // Loop over counter of QM atom types
 274     for (size_t i = 0; i < num_atoms.size(); i++)
 275     {
 276         if (num_atoms[i] > 0)
 277         {
 278             res += formatString("    &QM_KIND %3s\n", periodic_system[i].c_str());
 279             res += "      MM_INDEX";
 280             // Loop over all QM atoms indexes
 281             for (size_t j = 0; j < nQm; j++)
 282             {
 283                 if (parameters_.atomNumbers_[parameters_.qmIndices_[j]] == static_cast<index>(i))
 284                 {
 285                     res += formatString(" %d", static_cast<int>(parameters_.qmIndices_[j] + 1));
 286                 }
 287             }
 288
 289             res += "\n";
 290             res += "    &END QM_KIND\n";
 291         }
 292     }
 293
 294     // Print &LINK groups
 295     // Loop over parameters_.link_
 296     for (size_t i = 0; i < parameters_.link_.size(); i++)
 297     {
 298         res += "    &LINK\n";
 299         res += formatString("      QM_INDEX %d\n", static_cast<int>(parameters_.link_[i].qm) + 1);
 300         res += formatString("      MM_INDEX %d\n", static_cast<int>(parameters_.link_[i].mm) + 1);
 301         res += "    &END LINK\n";
 302     }
 303
 304     res += "  &END QMMM\n";
 305
 306     return res;
 307 }
 308
 309 std::string QMMMInputGenerator::generateMMSection()
 310 {
 311     std::string res;
 312
 313     res += "  &MM\n";
 314     res += "    &FORCEFIELD\n";
 315     res += "      DO_NONBONDED FALSE\n";
 316     res += "    &END FORCEFIELD\n";
 317     res += "    &POISSON\n";
 318     res += "      &EWALD\n";
 319     res += "        EWALD_TYPE NONE\n";
 320     res += "      &END EWALD\n";
 321     res += "    &END POISSON\n";
 322     res += "  &END MM\n";
 323
 324     return res;
 325 }
 326
 327 std::string QMMMInputGenerator::generateSubsysSection() const
 328 {
 329     std::string res;
 330
 331     // Init some numbers
 332     size_t nQm = parameters_.qmIndices_.size();
 333     size_t nMm = parameters_.mmIndices_.size();
 334     size_t nAt = nQm + nMm;
 335
 336     // Count the numbers of individual QM atoms per type
 337     std::vector<int> num_atoms(periodic_system.size(), 0);
 338     for (size_t i = 0; i < nQm; i++)
 339     {
 340         num_atoms[parameters_.atomNumbers_[parameters_.qmIndices_[i]]]++;
 341     }
 342
 343     res += "  &SUBSYS\n";
 344
 345     // Print cell parameters
 346     res += "    &CELL\n";
 347     res += formatString("      A %.3lf %.3lf %.3lf\n", box_[0][0] * 10, box_[0][1] * 10, box_[0][2] * 10);
 348     res += formatString("      B %.3lf %.3lf %.3lf\n", box_[1][0] * 10, box_[1][1] * 10, box_[1][2] * 10);
 349     res += formatString("      C %.3lf %.3lf %.3lf\n", box_[2][0] * 10, box_[2][1] * 10, box_[2][2] * 10);
 350     res += "      PERIODIC XYZ\n";
 351     res += "    &END CELL\n";
 352
 353     // Print topology section
 354     res += "    &TOPOLOGY\n";
 355
 356     // Placeholder for PDB file name that will be replaced with actuall name during mdrun
 357     res += "      COORD_FILE_NAME %s\n";
 358
 359     res += "      COORD_FILE_FORMAT PDB\n";
 360     res += "      CHARGE_EXTENDED TRUE\n";
 361     res += "      CONNECTIVITY OFF\n";
 362     res += "      &GENERATE\n";
 363     res += "         &ISOLATED_ATOMS\n";
 364     res += formatString("            LIST %d..%d\n", 1, static_cast<int>(nAt));
 365     res += "         &END\n";
 366     res += "      &END GENERATE\n";
 367     res += "    &END TOPOLOGY\n";
 368
 369     // Now we will print basises for all types of QM atoms
 370     // Loop over counter of QM atom types
 371     for (size_t i = 0; i < num_atoms.size(); i++)
 372     {
 373         if (num_atoms[i] > 0)
 374         {
 375             res += "    &KIND " + periodic_system[i] + "\n";
 376             res += "      ELEMENT " + periodic_system[i] + "\n";
 377             res += "      BASIS_SET DZVP-MOLOPT-GTH\n";
 378             res += "      POTENTIAL GTH-" + std::string(c_qmmmQMMethodNames[parameters_.qmMethod_]);
 379             res += "\n";
 380             res += "    &END KIND\n";
 381         }
 382     }
 383
 384     // Add element kind X - they are represents virtual sites
 385     res += "    &KIND X\n";
 386     res += "      ELEMENT H\n";
 387     res += "    &END KIND\n";
 388     res += "  &END SUBSYS\n";
 389
 390     return res;
 391 }
 392
 393
 394 std::string QMMMInputGenerator::generateCP2KInput() const
 395 {
 396
 397     std::string inp;
 398
 399     // Begin CP2K input generation
 400     inp += generateGlobalSection();
 401
 402     inp += "&FORCE_EVAL\n";
 403     inp += "  METHOD QMMM\n";
 404
 405     // Add DFT parameters
 406     inp += generateDFTSection();
 407
 408     // Add QMMM parameters
 409     inp += generateQMMMSection();
 410
 411     // Add MM parameters
 412     inp += generateMMSection();
 413
 414     // Add SUBSYS parameters
 415     inp += generateSubsysSection();
 416
 417     inp += "&END FORCE_EVAL\n";
 418
 419     return inp;
 420 }
 421
 422 std::string QMMMInputGenerator::generateCP2KPdb() const
 423 {
 424
 425     std::string pdb;
 426
 427     /* Generate *.pdb formatted lines
 428      * and append to std::string pdb
 429      */
 430     for (size_t i = 0; i < x_.size(); i++)
 431     {
 432         pdb += "ATOM  ";
 433
 434         // Here we need to print i % 100000 because atom counter in *.pdb has only 5 digits
 435         pdb += formatString("%5d ", static_cast<int>(i % 100000));
 436
 437         // Atom name
 438         pdb += formatString(" %3s ", periodic_system[parameters_.atomNumbers_[i]].c_str());
 439
 440         // Lable atom as QM or MM residue
 441         if (isQMAtom(i))
 442         {
 443             pdb += " QM     1     ";
 444         }
 445         else
 446         {
 447             pdb += " MM     2     ";
 448         }
 449
 450         // Coordinates
 451         pdb += formatString("%7.3lf %7.3lf %7.3lf  1.00  0.00         ",
 452                             (x_[i][XX] + qmTrans_[XX]) * 10,
 453                             (x_[i][YY] + qmTrans_[YY]) * 10,
 454                             (x_[i][ZZ] + qmTrans_[ZZ]) * 10);
 455
 456         // Atom symbol
 457         pdb += formatString(" %3s ", periodic_system[parameters_.atomNumbers_[i]].c_str());
 458
 459         // Point charge for MM atoms or 0 for QM atoms
 460         pdb += formatString("%lf\n", q_[i]);
 461     }
 462
 463     return pdb;
 464 }
 465
 466 const RVec& QMMMInputGenerator::qmTrans() const
 467 {
 468     return qmTrans_;
 469 }
 470
 471 const matrix& QMMMInputGenerator::qmBox() const
 472 {
 473     return qmBox_;
 474 }
 475
 476
 477 RVec computeQMBoxVec(const RVec& a, const RVec& b, const RVec& c, real h, real minNorm, real maxNorm)
 478 {
 479     RVec vec0(a);
 480     RVec vec1(b);
 481     RVec vec2(c);
 482     RVec dx(0.0, 0.0, 0.0);
 483     RVec res(0.0, 0.0, 0.0);
 484
 485     // Compute scales sc that will transform a
 486     dx = vec1.cross(vec2);
 487     dx /= dx.norm();
 488
 489     // Transform a
 490     res = h / static_cast<real>(fabs(vec0.dot(dx))) * a;
 491
 492     // If vector is smaller than minL then scale it up
 493     if (res.norm() < minNorm)
 494     {
 495         res *= minNorm / res.norm();
 496     }
 497
 498     // If vector is longer than maxL then scale it down
 499     if (res.norm() > maxNorm)
 500     {
 501         res *= maxNorm / res.norm();
 502     }
 503
 504     return res;
 505 }
 506
 507 } // namespace gmx