2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 The GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/gmxlib/network.h"
46 #include "gromacs/linearalgebra/nrjac.h"
47 #include "gromacs/math/do_fit.h"
48 #include "gromacs/math/functions.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdrunutility/multisim.h"
51 #include "gromacs/mdtypes/commrec.h"
52 #include "gromacs/mdtypes/fcdata.h"
53 #include "gromacs/mdtypes/inputrec.h"
54 #include "gromacs/mdtypes/mdatom.h"
55 #include "gromacs/mdtypes/state.h"
56 #include "gromacs/pbcutil/ishift.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/topology/ifunc.h"
59 #include "gromacs/topology/mtop_util.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/utility/arrayref.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/pleasecite.h"
64 #include "gromacs/utility/smalloc.h"
69 // TODO This implementation of ensemble orientation restraints is nasty because
70 // a user can't just do multi-sim with single-sim orientation restraints.
72 void init_orires(FILE* fplog,
73 const gmx_mtop_t* mtop,
76 const gmx_multisim_t* ms,
80 od->nr = gmx_mtop_ftype_count(mtop, F_ORIRES);
83 /* Not doing orientation restraints */
87 const int numFitParams = 5;
88 if (od->nr <= numFitParams)
91 "The system has %d orientation restraints, but at least %d are required, since "
92 "there are %d fitting parameters.",
98 if (ir->bPeriodicMols)
100 /* Since we apply fitting, we need to make molecules whole and this
101 * can not be done when periodic molecules are present.
104 "Orientation restraints can not be applied when periodic molecules are present "
111 "Orientation restraints do not work with MPI parallelization. Choose 1 MPI rank, "
115 GMX_RELEASE_ASSERT(globalState != nullptr, "We need a valid global state in init_orires");
117 od->fc = ir->orires_fc;
124 int* nr_ex = nullptr;
125 int typeMin = INT_MAX;
127 gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
129 while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol))
131 const int numOrires = (*il)[F_ORIRES].size();
132 if (nmol > 1 && numOrires > 0)
135 "Found %d copies of a molecule with orientation restrains while the current "
136 "code only supports a single copy. If you want to ensemble average, run "
137 "multiple copies of the system using the multi-sim feature of mdrun.",
141 for (int i = 0; i < numOrires; i += 3)
143 int type = (*il)[F_ORIRES].iatoms[i];
144 int ex = mtop->ffparams.iparams[type].orires.ex;
147 srenew(nr_ex, ex + 1);
148 for (int j = od->nex; j < ex + 1; j++)
154 GMX_ASSERT(nr_ex, "Check for allocated nr_ex to keep the static analyzer happy");
157 typeMin = std::min(typeMin, type);
158 typeMax = std::max(typeMax, type);
161 /* With domain decomposition we use the type index for indexing in global arrays */
163 typeMax - typeMin + 1 == od->nr,
164 "All orientation restraint parameter entries in the topology should be consecutive");
165 /* Store typeMin so we can index array with the type offset */
166 od->typeMin = typeMin;
168 snew(od->S, od->nex);
169 /* When not doing time averaging, the instaneous and time averaged data
170 * are indentical and the pointers can point to the same memory.
172 snew(od->Dinsl, od->nr);
176 snew(od->Dins, od->nr);
180 od->Dins = od->Dinsl;
183 if (ir->orires_tau == 0)
191 snew(od->Dtav, od->nr);
192 od->edt = std::exp(-ir->delta_t / ir->orires_tau);
193 od->edt_1 = 1.0 - od->edt;
195 /* Extend the state with the orires history */
196 globalState->flags |= enumValueToBitMask(StateEntry::OrireInitF);
197 globalState->hist.orire_initf = 1;
198 globalState->flags |= enumValueToBitMask(StateEntry::OrireDtav);
199 globalState->hist.norire_Dtav = od->nr * 5;
200 snew(globalState->hist.orire_Dtav, globalState->hist.norire_Dtav);
203 snew(od->oinsl, od->nr);
206 snew(od->oins, od->nr);
210 od->oins = od->oinsl;
212 if (ir->orires_tau == 0)
218 snew(od->otav, od->nr);
220 snew(od->tmpEq, od->nex);
223 for (int i = 0; i < mtop->natoms; i++)
225 if (getGroupType(mtop->groups, SimulationAtomGroupType::OrientationRestraintsFit, i) == 0)
230 snew(od->mref, od->nref);
231 snew(od->xref, od->nref);
232 snew(od->xtmp, od->nref);
234 snew(od->eig, od->nex * 12);
236 /* Determine the reference structure on the master node.
237 * Copy it to the other nodes after checking multi compatibility,
238 * so we are sure the subsystems match before copying.
240 auto x = makeArrayRef(globalState->x);
241 rvec com = { 0, 0, 0 };
244 for (const AtomProxy atomP : AtomRange(*mtop))
246 const t_atom& local = atomP.atom();
247 int i = atomP.globalAtomNumber();
248 if (mtop->groups.groupNumbers[SimulationAtomGroupType::OrientationRestraintsFit].empty()
249 || mtop->groups.groupNumbers[SimulationAtomGroupType::OrientationRestraintsFit][i] == 0)
251 /* Not correct for free-energy with changing masses */
252 od->mref[j] = local.m;
253 // Note that only one rank per sim is supported.
256 copy_rvec(x[i], od->xref[j]);
257 for (int d = 0; d < DIM; d++)
259 com[d] += od->mref[j] * x[i][d];
266 svmul(1.0 / mtot, com, com);
269 for (int j = 0; j < od->nref; j++)
271 rvec_dec(od->xref[j], com);
275 fprintf(fplog, "Found %d orientation experiments\n", od->nex);
276 for (int i = 0; i < od->nex; i++)
278 fprintf(fplog, " experiment %d has %d restraints\n", i + 1, nr_ex[i]);
283 fprintf(fplog, " the fit group consists of %d atoms and has total mass %g\n", od->nref, mtot);
287 fprintf(fplog, " the orientation restraints are ensemble averaged over %d systems\n", ms->numSimulations_);
289 check_multi_int(fplog, ms, od->nr, "the number of orientation restraints", FALSE);
290 check_multi_int(fplog, ms, od->nref, "the number of fit atoms for orientation restraining", FALSE);
291 check_multi_int(fplog, ms, ir->nsteps, "nsteps", FALSE);
292 /* Copy the reference coordinates from the master to the other nodes */
293 gmx_sum_sim(DIM * od->nref, od->xref[0], ms);
296 please_cite(fplog, "Hess2003");
299 void diagonalize_orires_tensors(t_oriresdata* od)
301 if (od->M == nullptr)
304 for (int i = 0; i < DIM; i++)
308 snew(od->eig_diag, DIM);
310 for (int i = 0; i < DIM; i++)
316 for (int ex = 0; ex < od->nex; ex++)
318 /* Rotate the S tensor back to the reference frame */
320 mmul(od->R, od->S[ex], TMP);
321 mtmul(TMP, od->R, S);
322 for (int i = 0; i < DIM; i++)
324 for (int j = 0; j < DIM; j++)
326 od->M[i][j] = S[i][j];
331 jacobi(od->M, DIM, od->eig_diag, od->v, &nrot);
334 for (int i = 0; i < DIM; i++)
338 for (int i = 0; i < DIM; i++)
340 for (int j = i + 1; j < DIM; j++)
342 if (gmx::square(od->eig_diag[ord[j]]) > gmx::square(od->eig_diag[ord[i]]))
351 for (int i = 0; i < DIM; i++)
353 od->eig[ex * 12 + i] = od->eig_diag[ord[i]];
355 for (int i = 0; i < DIM; i++)
357 for (int j = 0; j < DIM; j++)
359 od->eig[ex * 12 + 3 + 3 * i + j] = od->v[j][ord[i]];
365 void print_orires_log(FILE* log, t_oriresdata* od)
369 diagonalize_orires_tensors(od);
371 for (int ex = 0; ex < od->nex; ex++)
373 eig = od->eig + ex * 12;
374 fprintf(log, " Orientation experiment %d:\n", ex + 1);
375 fprintf(log, " order parameter: %g\n", eig[0]);
376 for (int i = 0; i < DIM; i++)
379 " eig: %6.3f %6.3f %6.3f %6.3f\n",
380 (eig[0] != 0) ? eig[i] / eig[0] : eig[i],
381 eig[DIM + i * DIM + XX],
382 eig[DIM + i * DIM + YY],
383 eig[DIM + i * DIM + ZZ]);
389 real calc_orires_dev(const gmx_multisim_t* ms,
391 const t_iatom forceatoms[],
392 const t_iparams ip[],
394 ArrayRef<const RVec> xWholeMolecules,
398 const history_t* hist)
401 real edt, edt_1, invn, pfac, r2, invr, corrfac, wsv2, sw, dev;
405 rvec * xref, *xtmp, com, r_unrot, r;
407 const real two_thr = 2.0 / 3.0;
411 /* This means that this is not the master node */
413 "Orientation restraints are only supported on the master rank, use fewer ranks");
416 bTAV = (od->edt != 0);
427 od->exp_min_t_tau = hist->orire_initf * edt;
429 /* Correction factor to correct for the lack of history
432 corrfac = 1.0 / (1.0 - od->exp_min_t_tau);
441 invn = 1.0 / ms->numSimulations_;
451 auto massT = md->massT;
452 auto cORF = md->cORF;
453 for (int i = 0; i < md->nr; i++)
457 copy_rvec(xWholeMolecules[i], xtmp[j]);
459 for (int d = 0; d < DIM; d++)
461 com[d] += mref[j] * xtmp[j][d];
467 svmul(1.0 / mtot, com, com);
468 for (int j = 0; j < nref; j++)
470 rvec_dec(xtmp[j], com);
472 /* Calculate the rotation matrix to rotate x to the reference orientation */
473 calc_fit_R(DIM, nref, mref, xref, xtmp, od->R);
475 for (int fa = 0; fa < nfa; fa += 3)
477 const int type = forceatoms[fa];
478 const int restraintIndex = type - od->typeMin;
481 pbc_dx_aiuc(pbc, x[forceatoms[fa + 1]], x[forceatoms[fa + 2]], r_unrot);
485 rvec_sub(x[forceatoms[fa + 1]], x[forceatoms[fa + 2]], r_unrot);
487 mvmul(od->R, r_unrot, r);
489 invr = gmx::invsqrt(r2);
490 /* Calculate the prefactor for the D tensor, this includes the factor 3! */
491 pfac = ip[type].orires.c * invr * invr * 3;
492 for (int i = 0; i < ip[type].orires.power; i++)
496 rvec5& Dinsl = od->Dinsl[restraintIndex];
497 Dinsl[0] = pfac * (2 * r[0] * r[0] + r[1] * r[1] - r2);
498 Dinsl[1] = pfac * (2 * r[0] * r[1]);
499 Dinsl[2] = pfac * (2 * r[0] * r[2]);
500 Dinsl[3] = pfac * (2 * r[1] * r[1] + r[0] * r[0] - r2);
501 Dinsl[4] = pfac * (2 * r[1] * r[2]);
505 for (int i = 0; i < 5; i++)
507 od->Dins[restraintIndex][i] = Dinsl[i] * invn;
514 gmx_sum_sim(5 * od->nr, od->Dins[0], ms);
517 /* Calculate the order tensor S for each experiment via optimization */
518 for (int ex = 0; ex < od->nex; ex++)
520 for (int i = 0; i < 5; i++)
522 matEq[ex].rhs[i] = 0;
523 for (int j = 0; j <= i; j++)
525 matEq[ex].mat[i][j] = 0;
530 for (int fa = 0; fa < nfa; fa += 3)
532 const int type = forceatoms[fa];
533 const int restraintIndex = type - od->typeMin;
534 rvec5& Dtav = od->Dtav[restraintIndex];
537 /* Here we update Dtav in t_fcdata using the data in history_t.
538 * Thus the results stay correct when this routine
539 * is called multiple times.
541 for (int i = 0; i < 5; i++)
543 Dtav[i] = edt * hist->orire_Dtav[restraintIndex * 5 + i]
544 + edt_1 * od->Dins[restraintIndex][i];
548 int ex = ip[type].orires.ex;
549 real weight = ip[type].orires.kfac;
550 /* Calculate the vector rhs and half the matrix T for the 5 equations */
551 for (int i = 0; i < 5; i++)
553 matEq[ex].rhs[i] += Dtav[i] * ip[type].orires.obs * weight;
554 for (int j = 0; j <= i; j++)
556 matEq[ex].mat[i][j] += Dtav[i] * Dtav[j] * weight;
560 /* Now we have all the data we can calculate S */
561 for (int ex = 0; ex < od->nex; ex++)
563 OriresMatEq& eq = matEq[ex];
564 /* Correct corrfac and copy one half of T to the other half */
565 for (int i = 0; i < 5; i++)
567 eq.rhs[i] *= corrfac;
568 eq.mat[i][i] *= gmx::square(corrfac);
569 for (int j = 0; j < i; j++)
571 eq.mat[i][j] *= gmx::square(corrfac);
572 eq.mat[j][i] = eq.mat[i][j];
575 m_inv_gen(&eq.mat[0][0], 5, &eq.mat[0][0]);
576 /* Calculate the orientation tensor S for this experiment */
577 matrix& S = od->S[ex];
583 for (int i = 0; i < 5; i++)
585 S[0][0] += 1.5 * eq.mat[0][i] * eq.rhs[i];
586 S[0][1] += 1.5 * eq.mat[1][i] * eq.rhs[i];
587 S[0][2] += 1.5 * eq.mat[2][i] * eq.rhs[i];
588 S[1][1] += 1.5 * eq.mat[3][i] * eq.rhs[i];
589 S[1][2] += 1.5 * eq.mat[4][i] * eq.rhs[i];
594 S[2][2] = -S[0][0] - S[1][1];
597 const matrix* S = od->S;
602 for (int fa = 0; fa < nfa; fa += 3)
604 const int type = forceatoms[fa];
605 const int restraintIndex = type - od->typeMin;
606 const int ex = ip[type].orires.ex;
608 const rvec5& Dtav = od->Dtav[restraintIndex];
609 od->otav[restraintIndex] =
611 * (S[ex][0][0] * Dtav[0] + S[ex][0][1] * Dtav[1] + S[ex][0][2] * Dtav[2]
612 + S[ex][1][1] * Dtav[3] + S[ex][1][2] * Dtav[4]);
615 const rvec5& Dins = od->Dins[restraintIndex];
616 od->oins[restraintIndex] =
618 * (S[ex][0][0] * Dins[0] + S[ex][0][1] * Dins[1] + S[ex][0][2] * Dins[2]
619 + S[ex][1][1] * Dins[3] + S[ex][1][2] * Dins[4]);
623 /* When ensemble averaging is used recalculate the local orientation
624 * for output to the energy file.
626 const rvec5& Dinsl = od->Dinsl[restraintIndex];
627 od->oinsl[restraintIndex] =
629 * (S[ex][0][0] * Dinsl[0] + S[ex][0][1] * Dinsl[1] + S[ex][0][2] * Dinsl[2]
630 + S[ex][1][1] * Dinsl[3] + S[ex][1][2] * Dinsl[4]);
633 dev = od->otav[restraintIndex] - ip[type].orires.obs;
635 wsv2 += ip[type].orires.kfac * gmx::square(dev);
636 sw += ip[type].orires.kfac;
638 od->rmsdev = std::sqrt(wsv2 / sw);
640 /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
641 for (int ex = 0; ex < od->nex; ex++)
644 tmmul(od->R, od->S[ex], RS);
645 mmul(RS, od->R, od->S[ex]);
650 /* Approx. 120*nfa/3 flops */
654 const t_iatom forceatoms[],
655 const t_iparams ip[],
660 real gmx_unused lambda,
661 real gmx_unused* dvdlambda,
662 gmx::ArrayRef<const real> /*charge*/,
663 t_fcdata gmx_unused* fcd,
664 t_disresdata gmx_unused* disresdata,
665 t_oriresdata* oriresdata,
666 int gmx_unused* global_atom_index)
668 int ex, power, ki = CENTRAL;
669 real r2, invr, invr2, fc, smooth_fc, dev, devins, pfac;
676 if (oriresdata->fc != 0)
678 bTAV = (oriresdata->edt != 0);
680 smooth_fc = oriresdata->fc;
683 /* Smoothly switch on the restraining when time averaging is used */
684 smooth_fc *= (1.0 - oriresdata->exp_min_t_tau);
687 for (int fa = 0; fa < nfa; fa += 3)
689 const int type = forceatoms[fa];
690 const int ai = forceatoms[fa + 1];
691 const int aj = forceatoms[fa + 2];
692 const int restraintIndex = type - oriresdata->typeMin;
695 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], r);
699 rvec_sub(x[ai], x[aj], r);
702 invr = gmx::invsqrt(r2);
704 ex = ip[type].orires.ex;
705 power = ip[type].orires.power;
706 fc = smooth_fc * ip[type].orires.kfac;
707 dev = oriresdata->otav[restraintIndex] - ip[type].orires.obs;
710 * there is no real potential when time averaging is applied
712 vtot += 0.5 * fc * gmx::square(dev);
716 /* Calculate the force as the sqrt of tav times instantaneous */
717 devins = oriresdata->oins[restraintIndex] - ip[type].orires.obs;
718 if (dev * devins <= 0)
724 dev = std::sqrt(dev * devins);
732 pfac = fc * ip[type].orires.c * invr2;
733 for (int i = 0; i < power; i++)
737 mvmul(oriresdata->S[ex], r, Sr);
738 for (int i = 0; i < DIM; i++)
740 fij[i] = -pfac * dev * (4 * Sr[i] - 2 * (2 + power) * invr2 * iprod(Sr, r) * r[i]);
743 for (int i = 0; i < DIM; i++)
749 fshift[ki][i] += fij[i];
750 fshift[CENTRAL][i] -= fij[i];
758 /* Approx. 80*nfa/3 flops */
761 void update_orires_history(const t_oriresdata& od, history_t* hist)
765 /* Copy the new time averages that have been calculated
766 * in calc_orires_dev.
768 hist->orire_initf = od.exp_min_t_tau;
769 for (int pair = 0; pair < od.nr; pair++)
771 for (int i = 0; i < 5; i++)
773 hist->orire_Dtav[pair * 5 + i] = od.Dtav[pair][i];